### Install and load the packages:
# eRm:
citation("eRm")
## The original JSS article: Article about CML Estimation in eRm: Article
## about LLRAs in eRm: Book Chapter about LLRAs: Article about the
## performance of Quasi-Exact Tests in eRm: Article about the Performance
## of the nonparametric Q3 Tests in eRm:
##
## Mair P, Hatzinger R (2007). "Extended Rasch modeling: The eRm package
## for the application of IRT models in R." _Journal of Statistical
## Software_, *20*. doi:10.18637/jss.v020.i09
## <https://doi.org/10.18637/jss.v020.i09>.
##
## Mair P, Hatzinger R (2007). "CML based estimation of extended Rasch
## models with the eRm package in R." _Psychology Science_, *49*.
## doi:10.18637/jss.v020.i09 <https://doi.org/10.18637/jss.v020.i09>.
##
## Hatzinger R, Rusch T (2009). "IRT models with relaxed assumptions in
## eRm: A manual-like instruction." _Psychology Science Quarterly_,
## *51*.
##
## Rusch T, Maier M, Hatzinger R (2013). "Linear logistic models with
## relaxed assumptions in R." In Lausen B, van den Poel D, Ultsch A
## (eds.), _Algorithms from and for Nature and Life_, series Studies in
## Classification, Data Analysis, and Knowledge Organization.
## doi:10.1007/978-3-319-00035-0_34
## <https://doi.org/10.1007/978-3-319-00035-0_34>.
##
## Koller I, Maier M, Hatzinger R (2015). "An empirical power analysis
## of quasi-exact tests for the Rasch model: Measurement invariance in
## small samples." _Methodology_, *11*. doi:10.1027/1614-2241/a000090
## <https://doi.org/10.1027/1614-2241/a000090>.
##
## Debelak R, Koller I (2019). "Testing the local independence
## assumption of the Rasch model with Q3-based nonparametric model
## tests." _Applied Psychological Measurement_, *44*.
## doi:10.1177/0146621619835501
## <https://doi.org/10.1177/0146621619835501>.
##
## To see these entries in BibTeX format, use 'print(<citation>,
## bibtex=TRUE)', 'toBibtex(.)', or set
## 'options(citation.bibtex.max=999)'.
#install.packages("eRm")
library("eRm")
## Warning: package 'eRm' was built under R version 4.4.1
transreas <- read.csv("transreas.csv")
summary(transreas)
## Student Grade task_01 task_02
## Min. : 1 Min. :2.000 Min. :0.0000 Min. :0.0000
## 1st Qu.:107 1st Qu.:3.000 1st Qu.:1.0000 1st Qu.:1.0000
## Median :213 Median :4.000 Median :1.0000 Median :1.0000
## Mean :213 Mean :4.005 Mean :0.9412 Mean :0.8094
## 3rd Qu.:319 3rd Qu.:5.000 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :425 Max. :6.000 Max. :1.0000 Max. :1.0000
## task_03 task_04 task_05 task_06
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:1.0000 1st Qu.:1.0000 1st Qu.:1.0000 1st Qu.:1.0000
## Median :1.0000 Median :1.0000 Median :1.0000 Median :1.0000
## Mean :0.8847 Mean :0.7835 Mean :0.8024 Mean :0.9741
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000
## task_07 task_08 task_09 task_10
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.00
## 1st Qu.:1.0000 1st Qu.:1.0000 1st Qu.:0.0000 1st Qu.:0.00
## Median :1.0000 Median :1.0000 Median :0.0000 Median :1.00
## Mean :0.8447 Mean :0.9671 Mean :0.3012 Mean :0.52
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.00
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.00
# Remove the Student ID and Grade level variables:
transreas.responses <- subset(transreas, select = -c(Student, Grade))
# Check descriptive statistics:
summary(transreas.responses)
## task_01 task_02 task_03 task_04
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:1.0000 1st Qu.:1.0000 1st Qu.:1.0000 1st Qu.:1.0000
## Median :1.0000 Median :1.0000 Median :1.0000 Median :1.0000
## Mean :0.9412 Mean :0.8094 Mean :0.8847 Mean :0.7835
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000
## task_05 task_06 task_07 task_08
## Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :0.0000
## 1st Qu.:1.0000 1st Qu.:1.0000 1st Qu.:1.0000 1st Qu.:1.0000
## Median :1.0000 Median :1.0000 Median :1.0000 Median :1.0000
## Mean :0.8024 Mean :0.9741 Mean :0.8447 Mean :0.9671
## 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000
## Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :1.0000
## task_09 task_10
## Min. :0.0000 Min. :0.00
## 1st Qu.:0.0000 1st Qu.:0.00
## Median :0.0000 Median :1.00
## Mean :0.3012 Mean :0.52
## 3rd Qu.:1.0000 3rd Qu.:1.00
## Max. :1.0000 Max. :1.00
# Run the Dichotomous Rasch Model
dichot.transreas <- RM(transreas.responses)
# Request a summary of the model results
summary(dichot.transreas)
##
## Results of RM estimation:
##
## Call: RM(X = transreas.responses)
##
## Conditional log-likelihood: -921.3465
## Number of iterations: 18
## Number of parameters: 9
##
## Item (Category) Difficulty Parameters (eta): with 0.95 CI:
## Estimate Std. Error lower CI upper CI
## task_02 0.258 0.133 -0.003 0.518
## task_03 -0.416 0.157 -0.723 -0.109
## task_04 0.441 0.128 0.190 0.692
## task_05 0.309 0.131 0.052 0.567
## task_06 -2.175 0.292 -2.747 -1.604
## task_07 -0.025 0.141 -0.302 0.252
## task_08 -1.909 0.262 -2.423 -1.395
## task_09 2.923 0.130 2.668 3.179
## task_10 1.836 0.115 1.610 2.062
##
## Item Easiness Parameters (beta) with 0.95 CI:
## Estimate Std. Error lower CI upper CI
## beta task_01 1.243 0.204 0.842 1.643
## beta task_02 -0.258 0.133 -0.518 0.003
## beta task_03 0.416 0.157 0.109 0.723
## beta task_04 -0.441 0.128 -0.692 -0.190
## beta task_05 -0.309 0.131 -0.567 -0.052
## beta task_06 2.175 0.292 1.604 2.747
## beta task_07 0.025 0.141 -0.252 0.302
## beta task_08 1.909 0.262 1.395 2.423
## beta task_09 -2.923 0.130 -3.179 -2.668
## beta task_10 -1.836 0.115 -2.062 -1.610
# Plot the Wright Map
plotPImap(dichot.transreas, main = "Transitive Reasoning Assessment Wright Map", irug = FALSE)
### Examine Item Parameters:
# View individual item parameters:
difficulty <- dichot.transreas$etapar
difficulty
## task_02 task_03 task_04 task_05 task_06 task_07
## 0.25775001 -0.41581384 0.44093780 0.30941151 -2.17531698 -0.02481129
## task_08 task_09 task_10
## -1.90884851 2.92343872 1.83581566
# Calculate the location for item 1:
n.items <- ncol(transreas.responses)
i1 <- 0 - sum(difficulty[1:(n.items - 1)])
difficulty.all <- c(i1, difficulty[c(1:(n.items - 1))])
difficulty.all
## task_02 task_03 task_04 task_05 task_06
## -1.24256308 0.25775001 -0.41581384 0.44093780 0.30941151 -2.17531698
## task_07 task_08 task_09 task_10
## -0.02481129 -1.90884851 2.92343872 1.83581566
# Alternative method to calculate item difficulty (easiness * -1):
difficulty2 <- dichot.transreas$betapar * -1
difficulty2
## beta task_01 beta task_02 beta task_03 beta task_04 beta task_05 beta task_06
## -1.24256308 0.25775001 -0.41581384 0.44093780 0.30941151 -2.17531698
## beta task_07 beta task_08 beta task_09 beta task_10
## -0.02481129 -1.90884851 2.92343872 1.83581566
# View standard errors for item parameters:
dichot.transreas$se.beta
## [1] 0.2043488 0.1327967 0.1566695 0.1282289 0.1314316 0.2916385 0.1413876
## [8] 0.2622114 0.1302767 0.1152173
# Examine descriptive statistics for item parameters:
summary(difficulty.all)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.1753 -1.0359 0.1165 0.0000 0.4081 2.9234
sd(difficulty.all)
## [1] 1.576441
summary(dichot.transreas$se.beta)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1152 0.1306 0.1371 0.1694 0.1924 0.2916
sd(dichot.transreas$se.beta)
## [1] 0.06206381
### This is what the book shows - 0.06206381
hist(difficulty.all, main = "Histogram of Item Difficulty Estimates for the Transitive Reasoning Data",
xlab = "Item Difficulty Estimates in Logits")
### Inserted this from the book page 17
###items.to.plot <- c(1:3)
plotICC(dichot.transreas, ask = FALSE)
### She missed the first line of code in this chunk.
student.locations <- person.parameter(dichot.transreas)
achievement <- student.locations$theta.table
achievement$id <-rownames(achievement)
# Add standard errors to person achievement object:
se <- as.data.frame(student.locations$se.theta$NAgroup1)
se$id <- rownames(se)
names(se) <- c("person_se", "id")
achievement.with.se <- merge(achievement, se, by = "id" )
# Examine descriptive statistics for person parameters:
summary(achievement)
## Person Parameter NAgroup Interpolated id
## Min. :-2.916 Min. :1 Mode :logical Length:425
## 1st Qu.: 1.160 1st Qu.:1 FALSE:374 Class :character
## Median : 1.932 Median :1 TRUE :51 Mode :character
## Mean : 2.010 Mean :1
## 3rd Qu.: 3.028 3rd Qu.:1
## Max. : 4.201 Max. :1
hist(achievement$`Person Parameter`, main = "Histogram of Person Achievement Estimates \nfor the Transitive Reasoning Data",
xlab = "Person Achievement Estimates in Logits")
# Examine item fit statistics:
student.locations <- person.parameter(dichot.transreas)
item.fit <- itemfit(student.locations)
item.fit
##
## Itemfit Statistics:
## Chisq df p-value Outfit MSQ Infit MSQ Outfit t Infit t Discrim
## task_01 245.550 373 1.000 0.657 0.761 -1.327 -1.622 0.455
## task_02 421.904 373 0.041 1.128 1.087 1.071 1.187 0.054
## task_03 218.155 373 1.000 0.583 0.745 -2.742 -2.688 0.583
## task_04 478.894 373 0.000 1.280 1.158 2.442 2.268 -0.014
## task_05 346.015 373 0.839 0.925 0.991 -0.625 -0.109 0.184
## task_06 65.613 373 1.000 0.175 0.611 -2.789 -1.781 0.620
## task_07 244.494 373 1.000 0.654 0.804 -2.768 -2.475 0.493
## task_08 125.478 373 1.000 0.336 0.687 -2.205 -1.569 0.563
## task_09 394.207 373 0.216 1.054 0.904 0.462 -1.727 0.022
## task_10 326.534 373 0.960 0.873 0.898 -1.822 -2.408 0.226
## Calculate reliability of item separation:
# Get Item scores
ItemScores <- colSums(transreas.responses)
ItemScores
## task_01 task_02 task_03 task_04 task_05 task_06 task_07 task_08 task_09 task_10
## 400 344 376 333 341 414 359 411 128 221
# Get Item SD
ItemSD <- apply(transreas.responses,2,sd)
ItemSD
## task_01 task_02 task_03 task_04 task_05 task_06 task_07 task_08
## 0.2355714 0.3932279 0.3197530 0.4123240 0.3986938 0.1589714 0.3626117 0.1786930
## task_09 task_10
## 0.4593099 0.5001886
# Calculate the se of the Item
ItemSE <- ItemSD/sqrt(length(ItemSD))
ItemSE
## task_01 task_02 task_03 task_04 task_05 task_06 task_07
## 0.07449423 0.12434958 0.10111476 0.13038830 0.12607804 0.05027118 0.11466788
## task_08 task_09 task_10
## 0.05650769 0.14524655 0.15817354
# compute the Observed Variance (also known as Total Person Variability or Squared Standard Deviation)
SSD.ItemScores <- var(ItemScores)
SSD.ItemScores
## [1] 8268.011
# compute the Mean Square Measurement error (also known as Model Error variance)
Item.MSE <- sum((ItemSE)^2) / length(ItemSE)
Item.MSE
## [1] 0.01291176
# compute the Item Separation Reliability
item.separation.reliability <- (SSD.ItemScores-Item.MSE) / SSD.ItemScores
item.separation.reliability
## [1] 0.9999984
###Item Separability matches the book.
## Examine person fit statistics:
person.fit <- personfit(student.locations)
summary(person.fit$p.infitMSQ)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.4060 0.4910 0.7487 0.9102 1.1740 2.2483
summary(person.fit$p.outfitMSQ)
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1723 0.2399 0.5468 0.7665 0.9592 7.3060
# Examine person reliability:
person_rel <- SepRel(student.locations)
person_rel$sep.rel
## [1] 0.2061478
### person separation reliability matches the book.
### Summarize the results in tables
summary.table.statistics <- c("Logit Scale Location Mean",
"Logit Scale Location SD",
"Standard Error Mean",
"Standard Error SD",
"Outfit MSE Mean",
"Outfit MSE SD",
"Infit MSE Mean",
"Infit MSE SD",
"Std. Outfit Mean",
"Std. Outfit SD",
"Std. Infit Mean",
"Std. Infit SD",
"Reliability of Separation")
item.summary.results <- rbind(mean(difficulty.all),
sd(difficulty.all),
mean(dichot.transreas$se.beta),
sd(dichot.transreas$se.beta),
mean(item.fit$i.outfitMSQ),
sd(item.fit$i.outfitMSQ),
mean(item.fit$i.infitMSQ),
sd(item.fit$i.infitMSQ),
mean(item.fit$i.outfitZ),
sd(item.fit$i.outfitZ),
mean(item.fit$i.infitMSQ),
sd(item.fit$i.infitZ),
item.separation.reliability)
person.summary.results <- rbind(mean(achievement.with.se$`Person Parameter`),
sd(achievement.with.se$`Person Parameter`),
mean(achievement.with.se$person_se, na.rm = TRUE),
sd(achievement.with.se$person_se, na.rm = TRUE),
mean(person.fit$p.outfitMSQ),
sd(person.fit$p.outfitMSQ),
mean(person.fit$p.infitMSQ),
sd(person.fit$p.infitMSQ),
mean(person.fit$p.outfitZ),
sd(person.fit$p.outfitZ),
mean(person.fit$p.infitZ),
sd(person.fit$p.infitZ),
person_rel$sep.rel)
# Round the values for presentation in a table:
item.summary.results_rounded <- round(item.summary.results, digits = 2)
person.summary.results_rounded <- round(person.summary.results, digits = 2)
Table1 <- cbind.data.frame(summary.table.statistics,
item.summary.results_rounded,
person.summary.results_rounded)
### Tables 1 & 2
# add descriptive column labels:
names(Table1) <- c("Statistic", "Items", "Persons")
print.data.frame(Table1, row.names = FALSE)
## Statistic Items Persons
## Logit Scale Location Mean 0.00 1.71
## Logit Scale Location SD 1.58 1.09
## Standard Error Mean 0.17 0.95
## Standard Error SD 0.06 0.16
## Outfit MSE Mean 0.77 0.77
## Outfit MSE SD 0.35 0.80
## Infit MSE Mean 0.86 0.91
## Infit MSE SD 0.18 0.47
## Std. Outfit Mean -1.03 0.08
## Std. Outfit SD 1.83 0.64
## Std. Infit Mean 0.86 -0.12
## Std. Infit SD 1.67 0.91
## Reliability of Separation 1.00 0.21
## Item calibration table:
# Calculate the proportion correct for each task:
PropCorrect <- apply(transreas.responses, 2, mean)
# Calculate the proportion correct for each task:
PropCorrect <- apply(transreas.responses, 2, mean)
# Combine item calibration results in a table:
Table2 <- cbind.data.frame(colnames(transreas.responses),
PropCorrect,
difficulty.all,
dichot.transreas$se.beta,
item.fit$i.outfitMSQ,
item.fit$i.outfitZ,
item.fit$i.infitMSQ,
item.fit$i.infitZ)
names(Table2) <- c("Task ID", "Proportion Correct", "Item Location","Item SE","Outfit MSE","Std. Outfit", "Infit MSE","Std. Infit")
# Sort Table 2 by Item difficulty:
Table2 <- Table2[order(-Table2$`Item Location`),]
# Round the numeric values (all columns except the first one) to 2 digits:
Table2[, -1] <- round(Table2[,-1], digits = 2)
print.data.frame(Table2, row.names = FALSE)
## Task ID Proportion Correct Item Location Item SE Outfit MSE Std. Outfit
## task_09 0.30 2.92 0.13 1.05 0.46
## task_10 0.52 1.84 0.12 0.87 -1.82
## task_04 0.78 0.44 0.13 1.28 2.44
## task_05 0.80 0.31 0.13 0.93 -0.62
## task_02 0.81 0.26 0.13 1.13 1.07
## task_07 0.84 -0.02 0.14 0.65 -2.77
## task_03 0.88 -0.42 0.16 0.58 -2.74
## task_01 0.94 -1.24 0.20 0.66 -1.33
## task_08 0.97 -1.91 0.26 0.34 -2.21
## task_06 0.97 -2.18 0.29 0.18 -2.79
## Infit MSE Std. Infit
## 0.90 -1.73
## 0.90 -2.41
## 1.16 2.27
## 0.99 -0.11
## 1.09 1.19
## 0.80 -2.47
## 0.75 -2.69
## 0.76 -1.62
## 0.69 -1.57
## 0.61 -1.78
## Person calibration table:
# Calculate proportion correct for persons:
PersonPropCorrect <- apply(student.locations$X.ex, 1, mean)
# Combine person calibration results in a table:
PersonPropCorrect <- apply(student.locations$X.ex, 1, mean)
# Combine person calibration results in a table:
Table3 <- cbind.data.frame(achievement.with.se$id,
PersonPropCorrect,
achievement.with.se$`Person Parameter`,
achievement.with.se$person_se,
person.fit$p.outfitMSQ,
person.fit$p.outfitZ,
person.fit$p.infitMSQ,
person.fit$p.infitZ)
names(Table3) <- c("Person ID", "Proportion Correct", "Person Location","Person SE","Outfit MSE","Std. Outfit", "Infit MSE","Std. Infit")
# Round the numeric values (all columns except the first one) to 2 digits:
Table3[, -1] <- round(Table3[,-1], digits = 2)
print.data.frame(Table3 [1:6, ], row.names = FALSE)
## Person ID Proportion Correct Person Location Person SE Outfit MSE Std. Outfit
## 1 0.8 1.93 0.94 0.24 -0.52
## 10 0.6 1.16 0.83 1.10 0.36
## 100 0.4 1.16 0.83 5.62 3.50
## 101 0.7 3.03 1.19 0.56 -0.37
## 102 0.9 0.52 0.77 0.17 0.12
## 103 0.7 3.03 1.19 0.73 -0.09
## Infit MSE Std. Infit
## 0.41 -1.19
## 1.13 0.47
## 2.10 2.55
## 0.71 -0.57
## 0.49 -0.61
## 0.80 -0.34
write.csv(Table3, file = "Table3.CSV", row.names = FALSE)
write.csv(Table2, file = "Table2.CSV", row.names = FALSE)
write.csv(Table1, file = "Table1.CSV", row.names = FALSE)