Dichotomous Rasch Model
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)
score.multiple.choice(key, mcr[-1], missing = F)## Call: score.multiple.choice(key = key, data = mcr[-1], missing = F)
##
## (Unstandardized) Alpha:
## [1] 0.78
##
## Average item correlation:
## [1] 0.14
##
## item statistics
## key 0 1 2 3 4 miss r n mean sd
## Q1 3 0.00 0.08 0.23 0.45 0.24 0 0.42 184 0.45 0.50
## Q2 2 0.00 0.31 0.23 0.14 0.32 0 0.40 184 0.23 0.42
## Q3 1 0.00 0.49 0.04 0.38 0.09 0 0.52 184 0.49 0.50
## Q4 1 0.00 0.33 0.09 0.08 0.51 0 -0.12 184 0.33 0.47
## Q5 4 0.01 0.21 0.46 0.13 0.20 0 0.28 184 0.20 0.40
## Q6 1 0.01 0.41 0.16 0.17 0.24 0 0.45 184 0.41 0.49
## Q7 4 0.00 0.10 0.28 0.24 0.38 0 0.47 184 0.38 0.49
## Q8 1 0.00 0.64 0.11 0.13 0.11 0 0.54 184 0.64 0.48
## Q9 4 0.00 0.11 0.19 0.04 0.65 0 0.45 184 0.65 0.48
## Q10 4 0.00 0.08 0.29 0.10 0.53 0 0.40 184 0.53 0.50
## Q11 1 0.00 0.62 0.06 0.14 0.17 0 0.38 184 0.62 0.49
## Q12 2 0.00 0.07 0.66 0.19 0.08 0 0.33 184 0.66 0.47
## Q13 3 0.00 0.34 0.18 0.40 0.08 0 0.55 184 0.40 0.49
## Q14 1 0.00 0.51 0.15 0.17 0.17 0 0.45 184 0.51 0.50
## Q15 3 0.01 0.11 0.18 0.44 0.26 0 0.50 184 0.44 0.50
## Q16 2 0.00 0.15 0.61 0.06 0.18 0 0.58 184 0.61 0.49
## Q17 3 0.01 0.05 0.05 0.64 0.25 0 0.39 184 0.64 0.48
## Q18 1 0.00 0.61 0.05 0.33 0.01 0 0.49 184 0.61 0.49
## Q19 2 0.01 0.09 0.59 0.05 0.27 0 0.49 184 0.59 0.49
## Q20 2 0.00 0.04 0.65 0.07 0.24 0 0.56 184 0.65 0.48
## Q21 2 0.01 0.11 0.61 0.10 0.17 0 0.32 184 0.61 0.49
## Q22 3 0.01 0.24 0.10 0.47 0.17 0 0.40 184 0.47 0.50
## Q23 2 0.00 0.05 0.77 0.06 0.12 0 0.29 184 0.77 0.42This 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.jml estimates a Rasch model using Joint Maximum Likelihood estimation, which is the same technique used in Winsteps.
rasch <- tam.jml(resp = mcrs)## ....................................................
## Iteration 1 2023-02-16 14:42:06
## MLE/WLE estimation |---
## Item parameter estimation |---
## Deviance= 4707.4852
## Maximum MLE/WLE change: 0.219437
## Maximum item parameter change: 0.243583
## ....................................................
## Iteration 2 2023-02-16 14:42:06
## MLE/WLE estimation |---
## Item parameter estimation |---
## Deviance= 4706.4386 | Deviance change: 1.0466
## Maximum MLE/WLE change: 0.08789
## Maximum item parameter change: 0.038265
## ....................................................
## Iteration 3 2023-02-16 14:42:06
## MLE/WLE estimation |---
## Item parameter estimation |--
## Deviance= 4706.3873 | Deviance change: 0.0513
## Maximum MLE/WLE change: 0.021142
## Maximum item parameter change: 0.007122
## ....................................................
## Iteration 4 2023-02-16 14:42:06
## MLE/WLE estimation |--
## Item parameter estimation |--
## Deviance= 4706.3851 | Deviance change: 0.0022
## Maximum MLE/WLE change: 0.004468
## Maximum item parameter change: 0.001401
## ....................................................
## Iteration 5 2023-02-16 14:42:06
## MLE/WLE estimation |--
## Item parameter estimation |--
## Deviance= 4706.385 | Deviance change: 1e-04
## Maximum MLE/WLE change: 0.000913
## Maximum item parameter change: 0.00028
## ....................................................
## Iteration 6 2023-02-16 14:42:06
## MLE/WLE estimation |--
## Item parameter estimation |-
## Deviance= 4706.385 | Deviance change: 0
## Maximum MLE/WLE change: 0.000186
## Maximum item parameter change: 5.7e-05
## ....................................................
## Iteration 7 2023-02-16 14:42:06
## MLE/WLE estimation |-
## Item parameter estimation |-
## Deviance= 4706.385 | Deviance change: 0
## Maximum MLE/WLE change: 3.7e-05
## Maximum item parameter change: 1.1e-05
##
## MLE/WLE estimation |-----
## ....................................................
##
## Start: 2023-02-16 14:42:06
## End: 2023-02-16 14:42:06
## Time difference of 0.03007102 secssummary(rasch)## ------------------------------------------------------------
## TAM 4.1-4 (2022-08-28 16:03:54)
## R version 4.2.2 (2022-10-31 ucrt) x86_64, mingw32 | nodename=LAPTOP-C3SSQ1UB | login=actio
##
## Start of Analysis: 2023-02-16 14:42:06
## End of Analysis: 2023-02-16 14:42:06
## Time difference of 0.03007102 secs
## Computation time: 0.03007102
##
## Joint Maximum Likelihood Estimation in TAM
##
## IRT Model
## Call:
## tam.jml(resp = mcrs)
##
## ------------------------------------------------------------
## Number of iterations = 7
##
## Deviance = 4706.39 | Log Likelihood = -2353.19
## Number of persons = 184
## Number of items = 23
## constraint = cases
## bias = TRUE
## ------------------------------------------------------------
## Person Parameters xsi
## M = 0
## SD = 1.05
## ------------------------------------------------------------
## Item Parameters xsi
## item N M xsi.item AXsi_.Cat1 B.Cat1.Dim1
## Q1 Q1 184 0.446 0.246 0.246 1
## Q2 Q2 184 0.234 1.374 1.374 1
## Q3 Q3 184 0.495 0.015 0.015 1
## Q4 Q4 184 0.326 0.842 0.842 1
## Q5 Q5 184 0.196 1.630 1.630 1
## Q6 Q6 184 0.413 0.403 0.403 1
## Q7 Q7 184 0.380 0.563 0.563 1
## Q8 Q8 184 0.641 -0.691 -0.691 1
## Q9 Q9 184 0.652 -0.747 -0.747 1
## Q10 Q10 184 0.533 -0.165 -0.165 1
## Q11 Q11 184 0.625 -0.610 -0.610 1
## Q12 Q12 184 0.663 -0.802 -0.802 1
## Q13 Q13 184 0.402 0.456 0.456 1
## Q14 Q14 184 0.505 -0.037 -0.037 1
## Q15 Q15 184 0.440 0.272 0.272 1
## Q16 Q16 184 0.614 -0.556 -0.556 1
## Q17 Q17 184 0.641 -0.691 -0.691 1
## Q18 Q18 184 0.609 -0.530 -0.530 1
## Q19 Q19 184 0.587 -0.424 -0.424 1
## Q20 Q20 184 0.647 -0.719 -0.719 1
## Q21 Q21 184 0.614 -0.556 -0.556 1
## Q22 Q22 184 0.473 0.117 0.117 1
## Q23 Q23 184 0.766 -1.383 -1.383 1
## ------------------------------------------------------------
## Item Parameters -A*Xsi
## xsi.label xsi.index xsi se.xsi
## 1 Q1 1 0.246 0.165
## 2 Q2 2 1.374 0.190
## 3 Q3 3 0.015 0.164
## 4 Q4 4 0.842 0.174
## 5 Q5 5 1.630 0.202
## 6 Q6 6 0.403 0.166
## 7 Q7 7 0.563 0.168
## 8 Q8 8 -0.691 0.169
## 9 Q9 9 -0.747 0.170
## 10 Q10 10 -0.165 0.164
## 11 Q11 11 -0.610 0.168
## 12 Q12 12 -0.802 0.171
## 13 Q13 13 0.456 0.167
## 14 Q14 14 -0.037 0.164
## 15 Q15 15 0.272 0.165
## 16 Q16 16 -0.556 0.167
## 17 Q17 17 -0.691 0.169
## 18 Q18 18 -0.530 0.167
## 19 Q19 19 -0.424 0.166
## 20 Q20 20 -0.719 0.170
## 21 Q21 21 -0.556 0.167
## 22 Q22 22 0.117 0.164
## 23 Q23 23 -1.383 0.188You might notice that the item difficulty estimates (xsi.item) are a bit different from Winsteps. This is due to which parameter gets constrained to 0. TAM by default sets the average person ability to 0 instead of average item difficulty. This accounts for the minor discrepancies.
TAM makes it pretty easy to get person reliability:
rasch$WLEreliability # this is person reliability## [1] 0.7583397Gathering Item Stats
Now we’re going to start pulling together item stats from various parts of the TAM model object.
rasch_items <- tam.jml.fit(rasch)$fit.item
rasch_items$beta <- rasch$xsi
rasch_items$se <- rasch$errorPLooking at Item Characteristic Curves
We can also look at ICCs for each item; these show theoretical and empirical curves of successful response probability.
plot(rasch)
## ....................................................
## Plots exported in png format into folder:
## C:/Users/actio/Google Drive/SLS 671 Research in Second Language Testing/Module 5 - Dichotomous Rasch Model/PlotsGathering Person Stats
Now we’ll gather person stats…
rasch_persons <- tam.jml.fit(rasch)$fit.person
rasch_persons$theta <- rasch$theta
rasch_persons$se <- rasch$errorWLEMaking a Wright Map
Wright Maps can also be created in R. We’re trying out a package called WrightMap that makes it pretty easy to get decent-looking maps. By setting item.side = “itemClassic”, we get that kind of old-school text-based histogram that Winsteps does on the item side.
wrightMap(rasch_persons$theta, rasch_items$beta, item.side = "itemClassic")
## [,1]
## [1,] 0.24620769
## [2,] 1.37389024
## [3,] 0.01453566
## [4,] 0.84164877
## [5,] 1.63030766
## [6,] 0.40300551
## [7,] 0.56315646
## [8,] -0.69149535
## [9,] -0.74658129
## [10,] -0.16491828
## [11,] -0.60999061
## [12,] -0.80234450
## [13,] 0.45594884
## [14,] -0.03670778
## [15,] 0.27216379
## [16,] -0.55631358
## [17,] -0.69149535
## [18,] -0.52965119
## [19,] -0.42401861
## [20,] -0.71895797
## [21,] -0.55631358
## [22,] 0.11717860
## [23,] -1.38282877