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.mml(Scored)
#Item parameters
mod1$xsi
#Thurstonian thresholds
tthresh <- tam.threshold(mod1)
tthresh
## Not run:
# #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.jml2(Scored)
# 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.mml(Scored2, xsi.fixed=anchor)
#
# #WLE ability estimates
# ability <- tam.wle(mod2)
# ability
#
# #CTT statistics
# ctt <- 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.mml(dummy, xsi.fixed=anchor)
# wle3 <- tam.wle(mod3)
# ## End(Not run)
Run the code above in your browser using DataLab