段階反応モデルにおける特性値θのベイズ推定

困難度・識別力パラメータの読み込み

#Q7,9,10,11,12,13,14 雲財中村（2018）に基づくEAP推定
Extrmt1=c(-2.09,-1.55,-1.99,-1.86,-2.35,-1.48,-1.38)
Extrmt2=c(-0.79,-0.55,-0.90,-0.78,-1.35,-0.49,-0.42)
Extrmt3=c(0.21,0.43,0.07,0.36,-0.43,0.60,0.64)
Extrmt4=c(1.51,1.37,0.96,1.47,0.74,1.47,1.54)
Dscrmn=c(1.35,2.49,1.95,1.72,1.38,2.46,2.71)
para2=data.frame(Extrmt1=Extrmt1,Extrmt2=Extrmt2,Extrmt3=Extrmt3,Extrmt4=Extrmt4,Dscrmn=Dscrmn)

特性値θの推定

library(irtoys)
library(ltm)
grm.theta(dat2, a=para2[,4],bc=para2[,c(1:4)],  D=1.0, method = "ML")

#EAP推定用の関数を定義
grm.eap <- function(X, a, bc, D=1.7, mu=0, sigma=1, n=50){
  X<-as.matrix(X);     a<-as.vector(a);   bc<-as.matrix(bc)
  ns  <- nrow(X);      ni<- ncol(X);      nc  <- ncol(bc)
  qp = normal.qu(n,mu=mu,sigma=sigma)$quad.points
  qw = normal.qu(n,mu=mu,sigma=sigma)$quad.weights
  o=sapply(1:ns, function(i) eap.inside(x=X[i, ],qp,qw,a,bc,D,nc,ni,n),USE.NAMES=F)
  rownames(o) = c("est", "sem")
  return(t(o))
}
eap.inside <- function (x,qp,qw,a,bc,D,nc,ni,n) {
  grm.LL <- function(theta){  
    LL01<-numeric(ni)
    p01<-1/(1+exp(-D*(a %*% t(rep(1,nc)))*(matrix(theta,ni,nc)-bc)))  #ni*nc
    p02<- cbind( rep(1,ni), p01, rep(0,ni))                       #ni*(nc+2) 
    for (j in 1:ni){ LL01[j]<-p02[j,x[j]+1]-p02[j,x[j]+2]}
    LL02<-sum(log(LL01))
    return(LL02)
  }
  ll<-numeric(n)
  for (i in 1:n) {
    ll[i]<-grm.LL(qp[i])
  }
  wl = exp(ll) * qw
  swl = sum(wl)
  xx = sum(wl * qp)/swl
  dev = qp - xx
  sem = sqrt(sum(wl * dev * dev)/swl)
  return(c(xx, sem))
}

cc3<-grm.eap(dat2,a=para2[,5],bc=para2[,c(1:4)],D=1.0,mu=0,sigma=1,n=50) #ロジスティック計量
dat3=cbind(dat2,cc3)
write.table(dat3, file="eap.txt")  　#書き出し