EIRM Tutorial - LLTM
In this tutorial, we’ll run a Rasch model and a LLTM based on the same data as De Boeck and Wilson, Chapter 2.
We’ll use the eRm package for modeling. lme4 is loaded just to get access to a dataset.
library(tidyverse)
library(eRm)
library(lme4) # just to get a dataset
data("VerbAgg")
summary(VerbAgg)## Anger Gender item resp id
## Min. :11 F:5832 S1WantCurse: 316 no :3973 1 : 24
## 1st Qu.:17 M:1752 S1WantScold: 316 perhaps:2081 2 : 24
## Median :19 S1WantShout: 316 yes :1530 3 : 24
## Mean :20 S2WantCurse: 316 4 : 24
## 3rd Qu.:23 S2WantScold: 316 5 : 24
## Max. :39 S2WantShout: 316 6 : 24
## (Other) :5688 (Other):7440
## btype situ mode r2
## curse:2528 other:3792 want:3792 N:3973
## scold:2528 self :3792 do :3792 Y:3611
## shout:2528
##
##
##
## We’ll need to create a dichotomous version of the response variable to follow along with Ch. 2. We’ll also make the data wide, which is how eRm operates.
#dichotomize for illustrative purposes
VerbAgg <- VerbAgg %>% mutate(respn = case_when(resp == "no" ~ 0,
resp == "perhaps" ~ 1,
resp == "yes" ~ 1))
#make data wide
resp_wide <- VerbAgg %>% select(id, item, respn) %>%
pivot_wider(names_from = item, values_from = respn)Basic Rasch Model
Now we can run a Rasch model. This should all be review at this point!
#basic rasch model
RM <- RM(resp_wide[2:25], se = T, sum0 = T)
summary(RM)##
## Results of RM estimation:
##
## Call: RM(X = resp_wide[2:25], se = T, sum0 = T)
##
## Conditional log-likelihood: -3049.923
## Number of iterations: 31
## Number of parameters: 23
##
## Item (Category) Difficulty Parameters (eta): with 0.95 CI:
## Estimate Std. Error lower CI upper CI
## S1WantScold -0.731 0.131 -0.987 -0.475
## S1WantShout -0.249 0.128 -0.501 0.003
## S2WantCurse -1.909 0.153 -2.210 -1.608
## S2WantScold -0.873 0.132 -1.132 -0.614
## S2WantShout -0.181 0.128 -0.433 0.070
## S3WantCurse -0.696 0.130 -0.951 -0.440
## S3WantScold 0.514 0.132 0.254 0.773
## S3WantShout 1.358 0.149 1.065 1.650
## S4wantCurse -1.245 0.137 -1.514 -0.976
## S4WantScold 0.178 0.129 -0.076 0.432
## S4WantShout 0.871 0.138 0.601 1.141
## S1DoCurse -1.383 0.140 -1.658 -1.109
## S1DoScold -0.557 0.129 -0.810 -0.303
## S1DoShout 0.698 0.135 0.434 0.963
## S2DoCurse -1.037 0.134 -1.300 -0.774
## S2DoScold -0.113 0.128 -0.365 0.138
## S2DoShout 1.312 0.148 1.022 1.602
## S3DoCurse 0.040 0.129 -0.212 0.293
## S3DoScold 1.335 0.149 1.044 1.626
## S3DoShout 2.871 0.222 2.436 3.306
## S4DoCurse -0.873 0.132 -1.132 -0.614
## S4DoScold 0.213 0.130 -0.042 0.467
## S4DoShout 1.840 0.165 1.516 2.165
##
## Item Easiness Parameters (beta) with 0.95 CI:
## Estimate Std. Error lower CI upper CI
## beta S1WantCurse 1.383 0.140 1.109 1.658
## beta S1WantScold 0.731 0.131 0.475 0.987
## beta S1WantShout 0.249 0.128 -0.003 0.501
## beta S2WantCurse 1.909 0.153 1.608 2.210
## beta S2WantScold 0.873 0.132 0.614 1.132
## beta S2WantShout 0.181 0.128 -0.070 0.433
## beta S3WantCurse 0.696 0.130 0.440 0.951
## beta S3WantScold -0.514 0.132 -0.773 -0.254
## beta S3WantShout -1.358 0.149 -1.650 -1.065
## beta S4wantCurse 1.245 0.137 0.976 1.514
## beta S4WantScold -0.178 0.129 -0.432 0.076
## beta S4WantShout -0.871 0.138 -1.141 -0.601
## beta S1DoCurse 1.383 0.140 1.109 1.658
## beta S1DoScold 0.557 0.129 0.303 0.810
## beta S1DoShout -0.698 0.135 -0.963 -0.434
## beta S2DoCurse 1.037 0.134 0.774 1.300
## beta S2DoScold 0.113 0.128 -0.138 0.365
## beta S2DoShout -1.312 0.148 -1.602 -1.022
## beta S3DoCurse -0.040 0.129 -0.293 0.212
## beta S3DoScold -1.335 0.149 -1.626 -1.044
## beta S3DoShout -2.871 0.222 -3.306 -2.436
## beta S4DoCurse 0.873 0.132 0.614 1.132
## beta S4DoScold -0.213 0.130 -0.467 0.042
## beta S4DoShout -1.840 0.165 -2.165 -1.516#person parameters
RM.pps <- person.parameter(RM)
#separation reliability
summary(SepRel(RM.pps))## Separation Reliability: 0.8821
##
## Observed Variance: 2.4882 (Squared Standard Deviation)
## Mean Square Measurement Error: 0.2933 (Model Error Variance)summary(gofIRT(RM.pps))##
## Goodness-of-Fit Tests
## value df p-value
## Collapsed Deviance 588.043 552 0.140
## Hosmer-Lemeshow 7.820 8 0.451
## Rost Deviance 4455.435 16777192 1.000
## Casewise Deviance 7072.594 7322 0.981
##
## R-Squared Measures
## Pearson R2: 0.366
## Sum-of-Squares R2: 0.366
## McFadden R2: 0.441
##
## Classifier Results - Confusion Matrix (relative frequencies)
## observed
## predicted 0 1
## 0 0.412 0.120
## 1 0.114 0.354
##
## Accuracy: 0.766
## Sensitivity: 0.747
## Specificity: 0.784
## Area under ROC: 0.849
## Gini coefficient: 0.699Linear Logisitic Test Model (LLTM)
Now for the fun part - we’ll estimate the difficulties associated with characteristics of the items and then compare to our doubly-descriptive Rasch model.
First, we’ll need to create a matrix of item characteristics. eRm needs a matrix of item responses (in this case, the same 0/1 data as used in our basic Rasch model) and a numeric matrix of item characteristics.
Let’s grab our item info first:
q_prep <- select(VerbAgg, item, btype, situ, mode) %>%
distinct()What we’ve done here is select just the columns with information about the items, and then saved only the distinct rows (we don’t need all 7,000+ rows of repeated info in the long data).
Now we’ll recode the factors to numeric values, following what De Boeck and Wilson described. We’ll save this as a matrix, which you can think of as a simpler type of data frame. eRm likes matrices, because matrices can be mathematically operated directly on in lots of useful ways.
q_prep <- q_prep %>% mutate(mode = case_when(mode == "do" ~ 1,
mode == "want" ~ 0),
situ = case_when(situ == "other" ~ 1,
situ == "self" ~ 0),
btype_v_shout = case_when(btype == "curse" ~ 1/2,
btype == "scold" ~ 1/2,
btype == "shout" ~ -1),
btype_v_scold = case_when(btype == "curse" ~ 1/2,
btype == "shout" ~ 1/2,
btype == "scold" ~ -1))
qmat <- select(q_prep, mode, situ, btype_v_shout, btype_v_scold) %>%
as.matrix()What we’ve done here is dummy-coded mode and situation, and then done sum coding for behavior type, which has three levels.
Now we can estimate our LLTM:
#Step 2: run LLTM
LLTM <- LLTM(resp_wide[2:25], qmat)
summary(LLTM)##
## Results of LLTM estimation:
##
## Call: LLTM(X = resp_wide[2:25], W = qmat)
##
## Conditional log-likelihood: -3130.414
## Number of iterations: 7
## Number of parameters: 4
##
## Basic Parameters eta with 0.95 CI:
## Estimate Std. Error lower CI upper CI
## mode -0.671 0.057 -0.783 -0.559
## situ 1.027 0.058 0.913 1.141
## btype_v_shout 1.359 0.050 1.261 1.457
## btype_v_scold 0.701 0.046 0.611 0.792
##
## Item Easiness Parameters (beta) with 0.95 CI:
## Estimate Std. Error lower CI upper CI
## beta S1WantCurse 2.057 0.076 1.908 2.207
## beta S1WantScold 1.005 0.070 0.868 1.142
## beta S1WantShout 0.018 0.067 -0.113 0.150
## beta S2WantCurse 2.057 0.076 1.908 2.207
## beta S2WantScold 1.005 0.070 0.868 1.142
## beta S2WantShout 0.018 0.067 -0.113 0.150
## beta S3WantCurse 1.030 0.042 0.948 1.113
## beta S3WantScold -0.022 0.039 -0.099 0.055
## beta S3WantShout -1.009 0.043 -1.092 -0.925
## beta S4wantCurse 1.030 0.042 0.948 1.113
## beta S4WantScold -0.022 0.039 -0.099 0.055
## beta S4WantShout -1.009 0.043 -1.092 -0.925
## beta S1DoCurse 1.386 0.090 1.209 1.563
## beta S1DoScold 0.334 0.088 0.162 0.506
## beta S1DoShout -0.653 0.088 -0.825 -0.480
## beta S2DoCurse 1.386 0.090 1.209 1.563
## beta S2DoScold 0.334 0.088 0.162 0.506
## beta S2DoShout -0.653 0.088 -0.825 -0.480
## beta S3DoCurse 0.359 0.068 0.227 0.492
## beta S3DoScold -0.693 0.069 -0.829 -0.557
## beta S3DoShout -1.680 0.074 -1.825 -1.534
## beta S4DoCurse 0.359 0.068 0.227 0.492
## beta S4DoScold -0.693 0.069 -0.829 -0.557
## beta S4DoShout -1.680 0.074 -1.825 -1.534## check the effects of item features
LLTM$etapar*-1## mode situ btype_v_shout btype_v_scold
## 0.6712050 -1.0269967 -1.3592355 -0.7014138Out output is identical to what De Boeck and Wilson reported:
situation (other/self) is -1.03 (easier to blame others than self)
mode (want/do) is 0.67 (going from wanting to doing is harder)
btype has complicated effects - see p. 65 in the reading for details
In the output above, item difficulties (well, easiness) are calculated by summing up the effects of item characteristics. We can compare these values to the Rasch model estimates.
#Step 3: Compared to basic RM
##correlation and r2
cor(RM$betapar*-1, LLTM$betapar*-1)## [1] 0.9418711cor(RM$betapar*-1, LLTM$betapar*-1)^2## [1] 0.8871211So uh, yeah, this is a pretty solid accounting of what makes these aggression items easy/difficult to endorse.
A graph:
ggplot(data = NULL, aes(LLTM$betapar*-1, y = RM$betapar*-1))+
geom_abline(slope = 1, intercept = 0)+
geom_smooth(method = "lm")+
geom_point()+
labs(x = "LLTM Item Difficulty", y = "Rasch Item Difficulty")+
theme_bw()## `geom_smooth()` using formula 'y ~ x'