項目反応理論（段階反応モデル，ベイズ推定）
データの読み込み
library(readr) SNFC_1 <- read_csv("SNFC_6.csv") View(SNFC_1) dat <- data.frame(SNFC_1[,8:22]) #特定の項目だけ取り出し dat=dat-1 dat2=subset(dat,select=c("Q1","Q2","Q3","Q4","Q6","Q7","Q8","Q9","Q10","Q11","Q12","Q13","Q14","Q15"))
IRTによる尺度構成
局所独立の検討
library(mirt) x <- mirt(dat2, 1) residuals(x) Theta <- fscores(x) round(residuals(x, type = 'Q3', Theta=Theta),2)
単調性の確認
library(mokken) mono.list<-check.monotonicity(dat2) summary(mono.list) (plot(mono.list,curves="ISRF", ci=TRUE,alpha = .05)) coefH(dat2)
段階反応モデルの実行
基本情報
library(ltm) descript(dat2,chi.squared = T)
特性反応曲線
par(mfrow=c(1,1)) par(ps=24) par(mar=c(5, 5, 5, 2)) plot(fitirt, legend=T, cx="left",main = "",xlab="特性値",ylab="確率",lwd=3)
テスト情報関数
par(mfrow=c(1,1)) par(ps=24) par(mar=c(5, 5, 5, 2)) plot(fitirt, type = "IIC", items = 0, plot = T,main="",lwd=3,xlab="特性値",ylab="情報量") par(mfrow=c(1,1)) vals <- plot(fitirt, type = "IIC", items = 0, plot = T) plot(vals[, "z"], 1 / sqrt(vals[, "test.info"]), type = "l", lwd = 3,xlab="特性値",ylab="標準誤差", main = "")
項目情報関数
par(mfrow=c(1,1)) #グラフ出力デフォルトに戻す plot(fitirt, type = "IIC", plot = T, legend=T, cx="left") information(fitirt,c(-3.1,3.1))
個人の能力値
fitt=factor.scores(fitirt, resp.pattern=dat2) SNFC_1$theta=fitt$score.dat$z1 View(SNFC_1)
因子得点と能力値θのプロット
plot(SNFC_1$SNFC,SNFC_1$theta) cor(SNFC_1[,23:29],use="complete.obs") #リストワイズ削除
水平テストの構築
#データの修正 typeA dat4=subset(dat,select=c("Q8","Q4","Q2","Q1","Q15","Q6","Q3")) #段階反応モデルの実行 (fitirtA=grm(dat4,Hessian = T)) #テスト情報関数 par(mfrow=c(1,1)) #グラフ出力デフォルトに戻す plot(fitirtA, type = "IIC", items = 0, plot = T,main="テスト情報曲線",xlim=c(-4,4), ylim=c(0,14), ylab="", lty=2) #データの修正 typeB dat5=subset(dat,select=c("Q14","Q13","Q9","Q11","Q7","Q10","Q12")) #段階反応モデルの実行 (fitirtB=grm(dat5,Hessian = T)) #テスト情報関数 par(new=T) plot(fitirtB, type = "IIC", items = 0, plot = T,main="",xlim=c(-4,4), ylim=c(0,14))
段階反応モデル verベイズ
呼び出し
require(rstan) options(max.print=400000000)
stanコード格納
stancode <- "
data{
int ni;
int nj;
int nc;
real D;
int<lower=1,upper=6> y[ni,nj];
}
parameters{
vector<lower=0,upper=5>[nj] a;
ordered[nc-1] ba[nj];
vector<lower=-4,upper=4>[ni] theta;
}
transformed parameters{
real b[nj,nc];
vector<lower=0,upper=1>[nc-1] pa[ni,nj];
simplex[nc] p[ni,nj];
for (j in 1:nj){
for (c in 1:nc){
if (c ==1){
b[j,c]=ba[j,c];
}else if (c ==nc){
b[j,c]=ba[j,c-1];
}else{
b[j,c]=(ba[j,c-1]+ba[j,c])/2;
}
}
}
for (i in 1:ni){
for (j in 1:nj){
for (c in 1:nc-1){
pa[i,j,c]= 1/(1+exp(-D*a[j]*(theta[i]-ba[j,c])));
}
}
}
for (i in 1:ni){
for (j in 1:nj){
for(c in 1:nc){
if (c==1){
p[i,j,c]=1-pa[i,j,c];
}else if(c==nc){
p[i,j,c]= pa[i,j,c-1];
}else{
p[i,j,c]= pa[i,j,c-1]-pa[i,j,c];
}
}
}
}
}
model{
for (i in 1:ni){
theta[i]~normal(0,1);
for (j in 1:nj){
y[i,j]~categorical(p[i,j]);
}
}
for (j in 1:nj){
a[j]~lognormal(0,sqrt(0.5));
for (c in 1:nc-1){
ba[j,c]~normal(0,2);
}
}
}
generated quantities{
real bg[nj,nc];
bg=b;
}このstanコードは豊田(2017)を引用
実行
u_ori<-dat2[1:14]+1 ni<-1508 nj<-14 nc<-5 D<-1 data_ori<-list(y=u_ori,nj=nj,ni=ni,nc=nc,D=D) par<-c("theta","ba","a","b") war<-1000 ite<-21000 see<-1234 dig<-3 cha<-5 fit_ori<-stan(file=stancode,model_name=stancode,data=data_ori,pars=par,verbose=F,seed=see,chains=cha,warmup=war,iter=ite) print(fit_ori,pars=par,digits_summary=dig) sink("./result1.txt") #保存 print(fit_ori,pars=par,digits_summary=dig) sink()
