######################################################################
# (1) Scored numeracy data
######################################################################
data(data.numeracy)
dat <- data.numeracy$scored
#Run IRT analysis: Rasch model
mod1 <- tam(dat)
#Item difficulties
mod1$xsi
ItemDiff <- mod1$xsi$xsi
ItemDiff
#Ability estimate - Weighted Likelihood Estimate
Abil <- tam.wle(mod1)
Abil
PersonAbility <- Abil$theta
PersonAbility
#Descriptive statistics of item and person parameters
hist(ItemDiff)
hist(PersonAbility)
mean(ItemDiff)
mean(PersonAbility)
sd(ItemDiff)
sd(PersonAbility)
#Extension
#plot histograms of ability and item parameters in the same graph
oldpar <- par(no.readonly = TRUE) # save writable default graphic settings
windows (width=4.45, height=4.45, pointsize=12)
layout(matrix(c(1,1,2),3,byrow=TRUE))
layout.show(2)
hist(PersonAbility,xlim=c(-3,3),breaks=20)
hist(ItemDiff,xlim=c(-3,3),breaks=20)
par( oldpar ) # restore default graphic settings
hist(PersonAbility,xlim=c(-3,3),breaks=20)
######################################################################
# (2) Raw numeracy data
######################################################################
raw_resp <- data.numeracy$raw
#score responses
key <- c(1, 1, 4, 1, 1, 1, 1, 1, 1, 1, 3, 1, 1, 1, 1)
scored <- sapply( seq(1,length(key)) ,
FUN = function(ii){ 1*(raw_resp[,ii] == key[ii]) } )
#run IRT analysis
mod1 <- tam(scored)
#Ability estimate - Weighted Likelihood Estimate
Abil <- tam.wle(mod1)
#CTT statistics
ctt1 <- tam.ctt(raw_resp, Abil$theta)
write.csv(ctt1,"D1_ctt1.csv") # write statistics into a file
# use maybe write.csv2 if ';' should be the column separator
#Fit statistics
Fit <- tam.fit(mod1)
Fit
# plot expected response curves
plot( mod1 , ask=TRUE )
Run the code above in your browser using DataLab