SEM and Multilevel Models with IRW Data • irw

Not every measurement model wants a response matrix. This article fits a categorical confirmatory factor analysis with lavaan, which needs the wide form, and then a Rasch model as a generalized linear mixed model with lme4, which consumes IRW long format directly with no reshape at all.

library(irw)
library(lavaan)
library(lme4)

The data

df <- irw_fetch("4thgrade_math_sirt")

That call needs Redivis credentials, which this documentation build does not have, so the article runs on a simulated table of the same shape:

df <- irw_simdata(n_id = 500, n_item = 8, model = "2PL", seed = 3)
head(df, 3)
#>   id item resp
#> 1  1    1    1
#> 2  2    1    1
#> 3  3    1    0

Categorical CFA with lavaan

lavaan needs one column per item, so pivot first:

wide <- irw_long2resp(df)
resp <- as.data.frame(wide[, setdiff(names(wide), "id"), drop = FALSE])

Because the item columns are already named item_1 … item_8, the model syntax can be built from names(resp) rather than typed out:

items <- names(resp)
model <- paste0("F =~ ", paste(items, collapse = " + "))
model
#> [1] "F =~ item_1 + item_2 + item_3 + item_4 + item_5 + item_6 + item_7 + item_8"

Binary items must be declared ordered, which switches lavaan to a diagonally weighted least squares estimator with a probit link. A CFA fitted this way is equivalent to a two-parameter normal-ogive IRT model:

fit <- cfa(model, data = resp, ordered = items)

Declaring the items ordered switches the estimator to DWLS with a scaled and shifted test statistic, so the .scaled fit measures are the ones to read; the unscaled pvalue is reported as NA under this test:

fitMeasures(fit, c("chisq.scaled", "df.scaled", "pvalue.scaled",
                   "cfi.scaled", "tli.scaled", "rmsea.scaled"))
#>  chisq.scaled     df.scaled pvalue.scaled    cfi.scaled    tli.scaled 
#>        15.659        20.000         0.738         1.000         1.043 
#>  rmsea.scaled 
#>         0.000

The fit is close to perfect, which is what simulated unidimensional data should produce. Real IRW tables rarely look this tidy.

head(parameterEstimates(fit, standardized = TRUE)[, c("lhs", "op", "rhs", "est", "std.all")], 8)
#>   lhs op    rhs   est std.all
#> 1   F =~ item_1 1.000   0.290
#> 2   F =~ item_2 0.539   0.156
#> 3   F =~ item_3 2.147   0.623
#> 4   F =~ item_4 1.401   0.406
#> 5   F =~ item_5 1.496   0.434
#> 6   F =~ item_6 1.763   0.511
#> 7   F =~ item_7 1.545   0.448
#> 8   F =~ item_8 1.426   0.414

The standardized loadings are a monotone transformation of IRT discriminations, so items that load highly here are the same items mirt would give a high a parameter.

Rasch as a mixed model with lme4

This is the one workflow where IRW’s native shape is exactly what the package wants. A Rasch model is a logistic mixed model with a random person intercept and fixed item effects, so the long data frame goes straight in:

df$id <- factor(df$id)
df$item <- factor(df$item)

fit_glmm <- glmer(resp ~ 0 + item + (1 | id), data = df, family = binomial)

The fixed effects are item easiness, so negate them to read them as the item difficulties an IRT package would report:

difficulty <- -fixef(fit_glmm)
head(round(difficulty, 3))
#>  item1  item2  item3  item4  item5  item6 
#> -0.206 -0.702 -0.071 -1.141 -0.334 -0.935

The person variance is the squared latent trait standard deviation:

VarCorr(fit_glmm)
#>  Groups Name        Std.Dev.
#>  id     (Intercept) 0.71708

Why bother

Writing the Rasch model this way buys you everything lme4 can do. Person covariates enter as fixed effects, giving an explanatory IRT model in one line, and crossed random effects handle raters without a separate package. IRW tables that carry a rater column suit the latter directly:

# for a table with a rater column, e.g. irw_fetch("swmd_mokken")
glmer(resp ~ 0 + item + (1 | id) + (1 | rater), data = df, family = binomial)

See vignette("dif-linking") for adding a person covariate, which is where irw_covariates() comes in.