data(data.timssAusTwn)
raw_resp <- data.timssAusTwn
#Recode data
resp <- raw_resp[,1:11]
#Column 12 is country code. Column 13 is gender code. Column 14 is Book ID.
all.na <- rowMeans( is.na(resp) )==1
#Find records where all responses are missing.
resp <- resp[!all.na,] #Delete records with all missing responses
resp[resp==20 | resp==21] <- 2 #TIMSS double-digit coding: "20" or "21" is a score of 2
resp[resp==10 | resp==11] <- 1 #TIMSS double-digit coding: "10" or "11" is a score of 1
resp[resp==70 | resp==79] <- 0 #TIMSS double-digit coding: "70" or "79" is a score of 0
resp[resp==99] <- 0 #"99" is omitted responses. Score it as wrong here.
resp[resp==96 | resp==6] <- NA #"96" and "6" are not-reached items. Treat these as missing.
#Score multiple-choice items #"resp" contains raw responses for MC items.
Scored <- resp
Scored[,9] <- (resp[,9]==4)*1 #Key for item 9 is D.
Scored[,c(1,2)] <- (resp[,c(1,2)]==2)*1 #Key for items 1 and 2 is B.
Scored[,c(10,11)] <- (resp[,c(10,11)]==3)*1 #Key for items 10 and 11 is C.
#Run IRT analysis for partial credit model (MML estimation)
mod1 <- TAM::tam.mml(Scored)
#Item parameters
mod1$xsi
#Thurstonian thresholds
tthresh <- TAM::tam.threshold(mod1)
tthresh
if (FALSE) {
#Plot Thurstonian thresholds
windows (width=8, height=7)
par(ps=9)
dotchart(t(tthresh), pch=19)
# plot expected response curves
plot( mod1, ask=TRUE)
#Re-run IRT analysis in JML
mod1.2 <- TAM::tam.jml(Scored)
stats::var(mod1.2$WLE)
#Re-run the model with "not-reached" coded as incorrect.
Scored2 <- Scored
Scored2[is.na(Scored2)] <- 0
#Prepare anchor parameter values
nparam <- length(mod1$xsi$xsi)
xsi <- mod1$xsi$xsi
anchor <- matrix(c(seq(1,nparam),xsi), ncol=2)
#Run IRT with item parameters anchored on mod1 values
mod2 <- TAM::tam.mml(Scored2, xsi.fixed=anchor)
#WLE ability estimates
ability <- TAM::tam.wle(mod2)
ability
#CTT statistics
ctt <- TAM::tam.ctt(resp, ability$theta)
write.csv(ctt,"TIMSS_CTT.csv")
#plot histograms of ability and item parameters in the same graph
windows(width=4.45, height=4.45, pointsize=12)
layout(matrix(c(1,1,2),3,byrow=TRUE))
layout.show(2)
hist(ability$theta,xlim=c(-3,3),breaks=20)
hist(tthresh,xlim=c(-3,3),breaks=20)
#Extension
#Score equivalence table
dummy <- matrix(0,nrow=16,ncol=11)
dummy[lower.tri(dummy)] <- 1
dummy[12:16,c(3,4,7,8)][lower.tri(dummy[12:16,c(3,4,7,8)])]<-2
mod3 <- TAM::tam.mml(dummy, xsi.fixed=anchor)
wle3 <- TAM::tam.wle(mod3)
}
Run the code above in your browser using DataLab