#################################################
#豊田秀樹(2002) 項目反応理論【入門編】朝倉書店 
#第3-6章 尺度値の推定，項目母数の推定，テスト情報関数，テストの構成 p.30-87
#################################################
library(irtoys)
source("http://blue.zero.jp/yokumura/R/testtheory/ctt.txt")


##(a) データの読み込み 
dat1 <- read.csv("http://blue.zero.jp/yokumura/R/testtheory/chap01_01.csv", header=F)
CTT(dat1)


##(b) 項目の削除 p.52
del1 <- mean(dat1)
del2 <- item.total(dat1)$it.cor
del.item <- which(del1 >.9 | del1 < .1 | del2 <.3)
dat2 <- dat1[,-del.item]
dat3 <- dat2[, setdiff(names(dat2), c("V24", "V25", "V33", "V34"))]


##(c) 項目母数の推定 (2母数モデルの周辺最尤推定法) p.52
pl2gaple <- est.gaple(data = dat2 + 1, method="ML", conv=.001, iter=100, ng=1, 
                      path="Z:/data/irt_begginner/gaple")
gaple.parm <- data.frame(a = pl2gaple[[2]][,"a"], b = pl2gaple[[2]][,"b1"])
rownames(gaple.parm) <- names(dat2)
gaple.parm


##(d) 表4.2の結果と比較 p.52
toyoda.parm <- data.frame(a = c(.45,.319,.627,.615,.366,.615,.811,.553,.809,.486,.451,1.225,.922,.607,.666,
                 .710,.650,.569,.560,.388,.712,.703,.858,.511,.612,.407,.646,.332,.630,.595),
                          b = c(-1.468,-1.702,-2.027,-1.735,-.698,-1.461,-1.535,.602,-1.595,-1.507,-1.665,-1.557,-1.751,-.907,-.898,
                                -.228,.367,.851,-1.191,.197,-.535,-.825,-.091,-.156,-2.023,1.878,.494,.535,.768,-.919))
rownames(toyoda.parm) <- names(dat2)
mean(abs(toyoda.parm$a - gaple.parm$a))   #実質的な違いがない 
mean(abs(toyoda.parm$b - gaple.parm$b))   #実質的な違いがない



##(e) 表4.5の結果と比較 p.56
toyoda.parm2 <- data.frame(a = c(.420, .326, .547, .561, .418, 
                                 .580, .719, .573, .731, .505, 
                                 .443, 1.044, .780, .598, .630,
                                 .710, .803, .729, .558, .422, 
                                 .719, .696, .983, .602, .538, 
                                 .397, .635, .353, .623, .580),                                 
                          b = c(-1.584, -1.726, -2.235, -1.861, -.637, 
                                 -1.530, -1.647, .612, -1.690, -1.481, 
                                 -1.719, -1.662, -1.915, -.920, -.932, 
                                 -.211, .339, .739, -1.205, .201, 
                                 -.521, -.827, -.063, -.130, -2.224,
                                 1.962, .524, .535, .800, -.940))
rownames(toyoda.parm2) <- names(dat2)
mean(abs(toyoda.parm2$a - gaple.parm$a))   #実質的な違いがない 
mean(abs(toyoda.parm2$b - gaple.parm$b))   #実質的な違いがない




##(f) 項目特性曲線と項目情報量 (図4.2, 図5.3) p.53 p.70
windows(width=14)
par(mfcol=c(2, 5), ask=T)
for(i in 1:length(gaple.parm$a)){
  logi2.icc(a=gaple.parm$a[i], b = gaple.parm$b[i])
  mtext(names(dat2)[i], 3, font=2, line=1)
  logi2.iteminfo(a=gaple.parm$a[i], b = gaple.parm$b[i], ylim=c(0, .8))
  mtext(names(dat2)[i], 3, font=2, line=1)
}
graphics.off() 


##(g) テスト情報曲線 (図5.2, 表5.1) p.67-68
logi2.testinfo(a = gaple.parm$a, b = gaple.parm$b)
tab5.1 <- logi2.testinfo(a = toyoda.parm2$a, b = toyoda.parm2$b, theta=seq(-3, 3, by=.25))
round(tab5.1, 3)      #表5.1 


##(h) 信頼性係数 (図5.5，表5.2) p.73-74
tab5.2 <- logi2.reliability(a=toyoda.parm2$a, b = toyoda.parm2$b, 
          mu.theta=seq(-4,  2, by=.5), sd.theta=seq(.5,1.5, length=13))
round(tab5.2$rho[,c("0.5","1", "1.5")], 3)
round(tab5.2$residual[,c("0.5","1", "1.5")], 3)
graphics.off()

tab5.2 <- logi2.reliability(a=gaple.parm$a, b = gaple.parm$b, 
          mu.theta=seq(-4,  2, by=.5), sd.theta=seq(.5,1.5, length=13))
round(tab5.2$rho[,c("0.5","1", "1.5")], 3)
round(tab5.2$residual[,c("0.5","1", "1.5")], 3)
graphics.off()


##(i) テスト特性曲線 (図6.1) p.78
a <- c(0.420,0.326,0.547,0.561)
b <- c(-1.584,-1.726,-2.235,-1.861)
weight <- rep(25, 4)

logi2.tcc.mean(a = a, b= b, weight = weight)
graphics.off()


##(j) 尺度のレベルごとの得点の分散 (図6.3) p.79
logi2.tcc.sd(a = a, b= b, weight = weight) 
graphics.off()


##(k) テスト得点の予測 (表6.1，図6.3) p.82
a <- c(0.420,0.326,0.547,0.561)
b <- c(-1.584,-1.726,-2.235,-1.861)
weight <- rep(1, 4)
mu <- -2
sigma.sq <- 1
logi2.predict.mean(a, b, weight, mu=mu, sigma.sq=sigma.sq)

mu <- seq(-2, 0, by=.5)
uy <- numeric(length(mu))
for(i in 1:length(mu)){
  uy[i] <- logi2.predict.mean(a, b, weight, mu=mu[i], sigma.sq=sigma.sq)
}
uy  #表6.1 

a <- toyoda.parm2$a
b <- toyoda.parm2$b
weight <- rep(1, length(a))
mu       <- seq(-4, 2, by=.2)
sigma.sq <- sqrt(seq(.5, 1.5, length=31))
uy <- matrix(nrow=length(mu), ncol=length(sigma.sq))
for(i in 1:length(mu)){
  for(j in 1:length(sigma.sq)){
      uy[i,j] <- logi2.predict.mean(a, b, weight, mu=mu[i], sigma.sq=sigma.sq[j])
  }
}
persp(mu, sigma.sq, uy,
      theta=200, phi=15,expand=.5,shade=1,col="lightblue",
      ticktype="detailed", ylab="sigma", xlab="mu", zlab="predicted mean")
graphics.off()  #図6.3


##(l) 尺度値の最尤推定/ベイズ推定 
est.ability <- logi2.ability(data=dat2, a=gaple.parm$a, b=gaple.parm$b)
mean(abs(pl2gaple[[3]]$MLE - est.ability$MLE), na.rm=T) #GAPLEの尺度値 (最尤推定) とオリジナル関数の相違は小さい 
mean(abs(pl2gaple[[3]]$EAP - est.ability$EAP), na.rm=T) #GAPLEの尺度値 (ベイズ推定) とオリジナル関数の相違は小さい