To get setup, first we’ll load the following packages and read in our data.

library(tidyverse)
library(TAM)
library(psych)
library(WrightMap)

mcr <- read_csv("multiple choice reading.csv")

We’re using data from Ch. 3 in the McNamara et al. (2019) textbook. This data is from a multiple-choice reading test. The data feature response choices which haven’t been scored. To score these, i.e. recode responses to 0/1, we’ll use a function from the psych package. First we’ll need to save the answer key somewhere handy!

key <- c(3,2,1,1,4,1,4,1,4,4,1,2,3,1,3,2,3,1,2,2,2,3,2)

This function, if you follow the defaults, also runs basic CTT reliability item analyses. But we don’t really want that- just the 0/1 scores are fine. The following code saves our score data in a usable format.

mcrs <- as_tibble(score.multiple.choice(key, mcr[-1], missing = F, score = F))

Rasch Model

We’ll be using the TAM package to run our Rasch analysis. The function tam estimates a Rasch model using Marginal Maximum Likelihood (MML) estimation; this differs slightly from Winsteps.

rasch <- tam(resp = mcrs)

TAM makes it pretty easy to get person reliability:

rasch$EAP.rel # this is person reliability
## [1] 0.7763395

Gathering Item Stats

Now we’re going to start pulling together item stats from various parts of the TAM model object.

rasch_items <- tam.fit(rasch)$itemfit
## Item fit calculation based on 100 simulations
## |**********|
## |----------|
rasch_items$beta <- rasch$xsi$xsi
rasch_items$se <- rasch$xsi$se.xsi

Local dependence

Now we’ll gather some stats on model-data fit. We can use some of it for dimensionality analysis, too…

rasch_fit <- tam.modelfit(rasch)
summary(rasch_fit)

In more detail…

rasch_fit$Q3.matr

Linacre says Q3 of < |.40| nothing much to worry about; in some situations people use a value of |.20| to screen for potential problems with local dependence (Chen, 1997).

rasch_fit$stat.itempair

As we can see here, only a handful of items exceed Q3 > |.20|. This is where you might want to take a look at the individual items themselves and see whether there might be a plausible reason for dependence. In this data, Q18 and Q19 have one of the highest Q3 values; it may be that the content of the two items is rather closely related.

PCA of Residuals

We’ll use the pca() function from psych.

st.res <- scale(rasch_fit$residuals)

pca <- pca(st.res, nfactors = ncol(st.res), rotate = "none")

Now some plots! First, the eigenvalues. You’re looking for a sharp drop and elbow.

plot(pca$Vaccounted[1,])

Nothing much.

We can also plot the variance accounted for by each contrast (it’s basically an identical plot…).

plot(pca$Vaccounted[2,])