---
title: "Using eirm for estimating explanatory IRT models"
output: rmarkdown::html_vignette
description: >
  How does eirm can be used for estimating dichotomous and polytomous 
  explanatory IRT models in R?
vignette: >
  %\VignetteIndexEntry{Using eirm for estimating explanatory IRT models}
  %\VignetteEngine{knitr::rmarkdown}
  \usepackage[utf8]{inputenc}
---

```{r setup, echo = TRUE, eval = TRUE, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, eval = TRUE, message = FALSE, warning = FALSE)
```

### Example 1: Rasch model

The Rasch model (i.e., a fully descriptive model) can be estimated using `eirm`. The following example shows how to estimate Rasch item parameters for the verbal aggression data set (see `?VerbAgg` for further details). A preview of the `VerbAgg` data set is shown below:

```{r, echo=FALSE}
library("eirm")
```

```{r}
data("VerbAgg")
head(VerbAgg)
```

To estimate the Rasch model, a regression-like formula must be defined: `formula = "r2 ~ -1 + item + (1|id)"`. In the formula, 

* `r2` is the variable for dichotomous item responses
* `-1` removes the intercept from the model and yields parameter estimates for all items in the data set. With `1` (instead of `-1`), an intercept representing the parameter of the first item and relative parameters for the remaining items (i.e., distance from the parameter of the first item) would be estimated.
* `item`is the variable representing item IDs in the data set
* `(1|id)` refers to the random effects for persons represented by the `id` column in the data set.

The output for the Rasch model is shown below:

```{r, echo = TRUE, eval = FALSE}
mod1 <- eirm(formula = "r2 ~ -1 + item + (1|id)", data = VerbAgg)
print(mod1)
```

By default, the `eirm` function returns the **easiness** parameters because the function uses a regression model parameterization where positive parameters indicate positive association with the dependent variable. In order to print the difficulty parameters (instead of easiness), `print(mod1, difficulty = TRUE)` must be used:

```{r, echo = TRUE, eval = FALSE}
print(mod1, difficulty = TRUE)
```

The `mod1` object is essentially a `glmerMod`-class object from the `lme4` package (Bates, Maechler, Bolker, & Walker (2015)). All `glmerMod` results for the estimated model can seen with `mod1$model`. For example, estimated random effects for persons (i.e., theta estimates) can be obtained using:

```{r, echo = TRUE, eval = FALSE}
theta <- ranef(mod1$model)$id
head(theta)
```

***
### Example 2: EIRM for dichotomous responses 

The following example shows how to use item-related and person-related explanatory variables to explain dichotomous responses in the verbal aggression data set. 

```{r, echo = TRUE, eval = TRUE}
mod2a <- eirm(formula = "r2 ~ -1 + situ + btype + mode + (1|id)", data = VerbAgg)

print(mod2a)
```

It is possible to visualize the parameters using an item-person map using `plot(mod2a)`, which returns the following plot. Note that this plot is a modified version of the `plotPImap` function from the `eRm` package). 

```{r, echo = TRUE, eval = FALSE}
plot(mod2a)
```

Aesthetic elements such as axis labels and plot title can be added to the plot. For example, the following code updates the x-axis label and the main plot title (see `?plot.eirm` for further details).

```{r, echo = TRUE, eval = FALSE}
plot(mod2a, difficulty = TRUE, main = "Verbal Aggression Example", latdim = "Verbal Aggression")
```

This plot wpuld show the difficulty parameters (instead of easiness), change the main title above the plot, and change the x-axis -- the name for the latent trait being measured. 

Also, it is possible to compare nested explanatory models with each other. The following example shows the estimation of a more compact version of `mod2a` with one less variable and compares the two models (i.e., `mod2a` vs. `mod2b`).

```{r, echo = TRUE, eval = FALSE}
mod2b <- eirm(formula = "r2 ~ -1 + situ + btype + (1|id)", data = VerbAgg)

anova(mod2a$model, mod2b$model)
```

***
### Example 3: EIRM for polytomous responses 

It is also possible to use the `eirm` function with polytomous item responses as well. Because the function only accepts dichotomous responses (i.e., binomial distribution), polytomous data must be reformatted first. To reformat the data, the `polyreformat` function can be used. The following example demonstrates how polytomous responses (no, perhaps, and yes) in the verbal aggression data set can be reformatted into a dichotomous form:

```{r, echo = TRUE, eval = TRUE}
VerbAgg2 <- polyreformat(data=VerbAgg, id.var = "id", long.format = FALSE, 
                         var.name = "item", val.name = "resp")

head(VerbAgg2)
```

In the reformatted data, `polyresponse` is the new dependent variable (i.e., pseudo-dichotomous version of the original response variable `resp`) and `polycategory` represents the response categories. Based on the reformatted data, each item has two rows (number of response categories - 1) based on the following rules (Stanke & Bulut (2019)) for further details on this parameterization):

* If `polycategory` = "cat_perhaps" and `resp` = "no", then `polyresponse` = 0
* If `polycategory` = "cat_perhaps" and `resp` = "perhaps", then `polyresponse` = 1
* If `polycategory` = "cat_perhaps" and `resp` = "yes", then `polyresponse` = NA

and

* If `polycategory` = "cat_yes" and `resp` = "no", then `polyresponse` = NA
* If `polycategory` = "cat_yes" and `resp` = "perhaps", then `polyresponse` = 0
* If `polycategory` = "cat_yes" and `resp` = "yes", then `polyresponse` = 1


**NOTE:** Although `polyreformat` is capable of reshaping wide-format data into long-format and reformat the long data for the analysis with `eirm`, a safer option is to transform the data from wide to long format before using `polyreformat`. The `melt` function from the `reshape2` package can easily transform wide data to long data. 

Several polytomous models can be estimated using the reformatted data:

**Model 1:** It explains only the first threshold (i.e., threshold from no to perhaps) based on explanatory variables:

```{r, echo = TRUE, eval = FALSE}
mod3 <- eirm(formula = "polyresponse ~ -1 + situ + btype + mode + (1|id)", 
             data = VerbAgg2)
```

**Model 2:** It explains the first threshold (i.e., threshold from no to perhaps) and second threshold (perhaps to yes) based on explanatory variables:

```{r, echo = TRUE, eval = FALSE}
mod4 <- eirm(formula = "polyresponse ~ -1 + btype + situ + mode + 
             polycategory + polycategory:btype + (1|id)", data = VerbAgg2)
```