Using eirm for estimating explanatory IRT models

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:

data("VerbAgg")
head(VerbAgg)

##   Anger Gender        item    resp id btype  situ mode r2
## 1    20      M S1WantCurse      no  1 curse other want  N
## 2    11      M S1WantCurse      no  2 curse other want  N
## 3    17      F S1WantCurse perhaps  3 curse other want  Y
## 4    21      F S1WantCurse perhaps  4 curse other want  Y
## 5    17      F S1WantCurse perhaps  5 curse other want  Y
## 6    21      F S1WantCurse     yes  6 curse other want  Y

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.
• itemis 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:

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:

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:

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.

mod2a <- eirm(formula = "r2 ~ -1 + situ + btype + mode + (1|id)", data = VerbAgg)

print(mod2a)

## EIRM formula: "r2 ~ -1 + situ + btype + mode + (1|id)" 
## 
## Number of persons: 316 
## 
## Number of observations: 7584 
## 
## Number of predictors: 5 
## 
## Parameter Estimates:
## 
##              Easiness       S.E.   z-value       p-value
## situother   1.7437027 0.10145090  17.18765  3.285796e-66
## situself    0.7158127 0.09775376   7.32261  2.431943e-13
## btypescold -1.0551642 0.06803581 -15.50895  3.017583e-54
## btypeshout -2.0421978 0.07292877 -28.00263 1.509082e-172
## modedo     -0.6715770 0.05621055 -11.94753  6.688613e-33
## 
## Note: The estimated parameters above represent 'easiness'.
## Use difficulty = TRUE to get difficulty parameters.

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).

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).

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).

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:

VerbAgg2 <- polyreformat(data=VerbAgg, id.var = "id", long.format = FALSE, 
                         var.name = "item", val.name = "resp")

head(VerbAgg2)

##   Anger Gender        item    resp id btype  situ mode r2 polycategory
## 1    20      M S1WantCurse      no  1 curse other want  N  cat_perhaps
## 2    11      M S1WantCurse      no  2 curse other want  N  cat_perhaps
## 3    17      F S1WantCurse perhaps  3 curse other want  Y  cat_perhaps
## 4    21      F S1WantCurse perhaps  4 curse other want  Y  cat_perhaps
## 5    17      F S1WantCurse perhaps  5 curse other want  Y  cat_perhaps
## 6    21      F S1WantCurse     yes  6 curse other want  Y  cat_perhaps
##   polyresponse                polyitem
## 1            0 S1WantCurse.cat_perhaps
## 2            0 S1WantCurse.cat_perhaps
## 3            1 S1WantCurse.cat_perhaps
## 4            1 S1WantCurse.cat_perhaps
## 5            1 S1WantCurse.cat_perhaps
## 6           NA S1WantCurse.cat_perhaps

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:

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:

mod4 <- eirm(formula = "polyresponse ~ -1 + btype + situ + mode + 
             polycategory + polycategory:btype + (1|id)", data = VerbAgg2)