---
##title: "IRT functions 1"
author: "Module created by Larry Nelson (l.nelson@curtin.edu.au)"
date: "RStudio code of 7 December 2021"
#output: word_document
output: html_document
#output: html_notebook
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, comment=NA)

library("psych")
library("ltm")
##library("TAM")
##library("WrightMap")

#get R to set the working directory to the folder containing this RMD file.
#(this works ok if this Rmd file has been opened from within the folder it's saved in)
this.dir <- function(directory)
setwd(file.path(getwd(),directory))
```
***
#### Module 'UWO-1': dichotomous IRT analyses

This report has been created by an [R program](https://en.wikipedia.org/wiki/R_(programming_language)) written to read item scores data from a csv file, such as the Omega-IScores.csv file often used with Lertap 5, and then to use functions in selected R packages in order to get a variety of IRT statistics and graphs.

Large sections of the R code used to generate the report are from a paper by Smyth and Johnson, University of Western Ontario, UWO. A copy of their paper, 'Item Response Theory for Dichotomous Items', might be [available here](http://www.uwo.ca/fhs/tc/labs/11.DichotomousIRT.pdf). Their code was integrated into an Rmd code module assembled by Larry Nelson, Curtin University, Western Australia (l.nelson@curtin.edu.au), last updated 28 June 2018.

'Omega-IScores.csv' files are created by Lertap. A reference to the respective Lertap help page is [here](http://www.lertap5.com/HTMLHelp/HTML/index.html?omega1.htm). It is **NOT** necessary to use Omega-IScores.csv files -- replace the 'scored <- read.csv("Omega-IScores.csv")' below with whatever csv file is appropriate.

The csv file used in this analysis was obtained from this folder:    

```{r echo=FALSE}

getwd() # Okay, to confirm: R's working directory has been set to:
 
```
The files in this folder are listed below:

```{r}
dir()
```

The first step gets R to read the csv file into an internal working dataset called "scored". Then the first few records of the dataset are displayed as a check. **You should confirm that these records are as expected!** The numbers in the first column have been inserted by R. *The second column should contain scores for the first item being processed.*

Note: lines that begin with ## are ignored by the processor.

```{r}
scored <- read.csv("Omega-IScores.csv")
##scored <- read.csv("D1_scored.csv") ## from the FIMS study
##scored <- read.csv("EIRT-binaries.csv") ## from EIRT's main dataset
##scored <- LSAT  ## a standard dataset from the Psych package
head (scored) 
```
***
We'll try fitting the data to the Rasch model, with item discrimination fixed ('constrained') to a value of 1 -- this will be referred to in this report as **'Model1'**.

A reference for the Rasch model is [here](https://en.wikipedia.org/wiki/Rasch_model). A basic introduction to IRT, with a focus on item characteristic curves is [this one](https://support.sas.com/resources/papers/proceedings14/SAS364-2014.pdf) from SAS.

```{r }
fit1 <- rasch(scored, constraint=cbind(length(scored)+1,1))
```

We can use the 'summary' function to get fit statistics for an IRT model. Results appear below under 'Model Summary'. In this case, the results apply to the constrained Rasch model (fit1).

```{r}
summary(fit1)
```

The 'coef' function can be used to get a table of item difficulty, item discrimination (fixed at a value of 1 in this case), and the probablility that an 'average' individual (z-score of 0) will answer an item correctly.

```{r}
coef(fit1,prob=TRUE)
```
A goodness-of-fit measure, 'GoF', is next. A significant p-value (such as 0.050 or less) is taken to suggest that the model is **not** a good fit to the data. BE PATIENT with this step as a FEW MINUTES might be required for the computation cycle to conclude.

(Note: lines with ## at the front are not processed by R, they're ignored. Computing GoF can be computationally demanding, possibly taking more than a minute or two to finish -- for this reason GoF will sometimes be deactivated by the use of ## at the start of some of the lines below.)

```{r}
##GoF.rasch(fit1,B=199)
##GoF.rasch(fit1)
```
***
Following the suggestions in the Smyth and Johnson paper referenced above, we can now see what might happen if the constraint on the discrimination parameter in the Rasch model is relaxed (above it was constrained to a value of 1; now it'll be free to assume any value -- the value appears as 'Dscrmn' below). This will be referred to here as **'Model2'**.

```{r}
fit2 <- rasch(scored)
summary(fit2)
```
***
Is Model2 a better fit to the Rasch model? The 'anova' function in the ltm package may be used to test the hypothesis that there's no significant difference between models. If the p-value below is less than 0.050, the hypothesis is rejected, and the second model might be regarded as providing a better fit.

A reference for these fit statistics is [this Wikipedia page](https://en.wikipedia.org/wiki/Akaike_information_criterion).

```{r}
anova(fit1,fit2)
```
***
Okay, now let's try a third model by allowing items to have different discrimination values, we'll refer to it here as **'Model3'**. This model is called the **'2PL'** model in the literature, the 'two-parameter logistic model'. The table of coefficients (below) has the difficulty estimate ('Dffclt') for each item, followed by the discrimination estimate ('Dscrmn') for each item.

```{r}
fit3 <-ltm(scored ~ z1)
summary(fit3)
```
***
Is Model3, 2PL, better than Model2, the unconstrained Rasch model? That is, does Model3 appear to fit the data better than Model2? If the p-value below is significant, less than 0.050, then yes, it would seem so.

```{r}
anova(fit2, fit3)
```
***
The Rasch and 2PL models, as used above, might be seen as incomplete as they do not reflect the possible effects of guessing. The **'3PL'** IRT model is used to estimate item difficultly, item discrimination, and an item guessing parameter.

We'll get 3PL results next, and, since this is the fourth 'model' used in this report, we'll refer to here it as **'Model4'** (the first model was the constained Rasch model, the second was the unconstrained Rasch model, the third was the 2PL model, and now we'll have 3PL, 'Model4'). In the table below, estimates of item guessing parameters are labelled as 'Gussng'.

```{r}
fit4 <-tpm(scored, type="latent.trait", max.guessing=1)
summary(fit4)
```
***
Let's now see if the data fit the '3PL' IRT model better than the '2PL' IRT model. If the p-value below is significant (less than 0.050), we will have evidence to suggest that it does.

```{r}
anova(fit3, fit4)
```
***
How about ability estimates? It's possible to get them, to be sure, but note: we have used ## markers below to turn off the estimates as the output can be quite lenghty. If you want theta estimates, remove the ## in front of factor.scores(fit3) .

```{r}
##factor.scores(fit3)
```

The factor.scores(fit1) line below can be used (if wanted) to get ability estimates based on the constrained Rasch model:

```{r}
##factor.scores(fit1)
```
***
Graphics are always great to have, a definite aid in the interpretation of IRT results. Here are some **'ICCs'**, item characteristic curves, for **Model1**, the Rasch model with discrimination constrained to have a value of 1.

```{r}
plot(fit1,
legend = TRUE,
     cx = "bottomright",
     lwd = 3,
     cex.main = 1.5,
     cex.lab = 1.3,
cex = 0.8)
```

Next, ICCs for the **2PL** model:

```{r}
plot(fit3,
legend = TRUE,
     cx = "bottomright",
     lwd = 3,
     cex.main = 1.5,
     cex.lab = 1.3,
cex = 0.8)
```

Next: ICCs for the **3PL** model:

```{r}
plot(fit4,
legend = TRUE,
     cx = "bottomright",
     lwd = 3,
     cex.main = 1.5,
     cex.lab = 1.3,
cex = 0.8)
```

**Item information curves** are shown below. They're used to indicate where each item best measures the latent trait. The curves are based on results from the **2PL** model, 'fit3' in this case.

```{r}
plot(fit3,
legend = TRUE,
     cx = "topright",
     cex = 0.8,
     type = "IIC",
     annot = FALSE,
     lwd = 3,
     cex.main = 1.5,
     cex.lab = 1.3)
```

Last but not least we come to a graph of the **'test information function'**. It's used as an indicator of where the items, taken as a whole (that is, as a 'test'), provide the most precise information on the latent trait ('ability'). In this example, 'fit3' has been employed, the 2PL IRT model.

An excellent and classic reference to test information and item information functions (or 'curves') is [this one](http://echo.edres.org:8080/irt/baker/chapter6.pdf) by Frank Baker.

```{r}
plot(fit3, type="IIC",
     items = 0,
     lwd = 3,
     cex.main = 1.5,
     cex.lab = 1.3)
```

***
***
##### THE END!
