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 <- (rasch_fit$residuals-mean(rasch_fit$residuals))/sd(rasch_fit$residuals)

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

Now some plots! First, the eigenvalues. We’ll zoom in on the first 5 contrasts, like Winsteps does. You’re looking for a sharp drop and elbow.

plot(pca$Vaccounted[1,][1:5])

Nothing much; mostly looks like a gradual decline. No eigenvalues above 2.0, either.

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

plot(pca$Vaccounted[2,][1:5])

Proportion of Unexpectedly large Residuals

First, make the residual matrix long.

st.res.long <- as_tibble(st.res) %>%
  rownames_to_column() %>%
  pivot_longer(Q1:Q23, names_to = "item", values_to = "std.resid")

Then, tally how many residuals are greater than |2.0| and |3.0|, divide by total number of observations.

2.0:

st.res.long %>% group_by(abs(std.resid) > 2.0) %>%
  count() %>%
  mutate(proportion = n/4232)
## # A tibble: 2 × 3
## # Groups:   abs(std.resid) > 2 [2]
##   `abs(std.resid) > 2`     n proportion
##   <lgl>                <int>      <dbl>
## 1 FALSE                 4214    0.996  
## 2 TRUE                    18    0.00425

3.0:

st.res.long %>% group_by(abs(std.resid) > 3.0) %>%
  count() %>%
  mutate(proportion = n/4232)
## # A tibble: 1 × 3
## # Groups:   abs(std.resid) > 3 [1]
##   `abs(std.resid) > 3`     n proportion
##   <lgl>                <int>      <dbl>
## 1 FALSE                 4232          1