貝葉斯項目反應模型及nimble實現

于鈴玉， 曹 蕾， 張婉婷， 李佳偉

(長春工業大學 數學與統計學院， 吉林 長春 130012)

0 引 言

對于項目反應模型，從連接函數角度來看，常用的連接函數有logistic和probit；從參數角度來看，有單參數、雙參數和三參數的項目反應模型，并且它們均可以用logistic和probit連接函數。Leventhal Brian C等[1]、張雪[2]、張國紅[3]、Bolsinova M等[4]項目反應模型的連接函數為logistic連接函數。Lu J等[5]、Bolsinova M等[6]項目反應模型是probit連接函數。

文中首先介紹了單參數、雙參數和三參數的項目反應模型。其次，分別采用連接函數logistic和probit對單參數、雙參數和三參數的項目反應模型進行計算并比較。應用了DIC和LPML兩個模型方法，DIC小說明模型好，LPML大說明模型好。最后，在實證研究部分，通過計算每個模型的DIC值和LPML值得出結論，當連接函數為logistic的雙參數模型最好，因為DIC值最小，同時LPML值最大。由于貝葉斯計算項目反應模型存在一定困難，鏈特別長，所以針對實證研究的程序是由nimble包編寫的，nimble包的語法結構類似于WinBUGS和JAGS，但是比WinBUGS和JAGS更靈活，因為編程時用C++代碼，所以運算速度快[7]。

1 模型分析

1.1 模型的建立

被試者用i(i=1,2,…,N)表示，題目用j(j=1,2,…,J)表示。Yij表示第i個人對第j道題目的回答，Yij=1表示回答正確;Yij=0表示回答錯誤。假定θi為第i個人的能力參數，用單參數的logistic連接函數表示第i個人答對第j個題目的概率為[8]

用雙參數的logistic連接函數表示第i個人答對第j個題目的概率為

pij=P(Yij=1|θi,aj,bj)=

用三參數的logistic連接函數表示第i個人答對第j個題目的概率為

pij=P(Yij=1|θi,aj,bj,cj)=

式中：aj、bj、cj----分別表示項目j的區分度參數、難度參數和偽猜測參數。

假定aj(j=1,2,…,J)相互獨立，bj(j=1,2,…,J)，cj(j=1,2,…,J)相互獨立，且假設它們的先驗logaj～N(0,1),bj～N(0,1),cj～Beta(4,12)。

此外，還可用單參數的probit連接函數表示第i個人答對第j個題目的概率為[9]

pij=P(Yij=1|θi,bj)=Φ(θi-bj)，

用雙參數的probit連接函數表示第i個人答對第j個題目的概率為

pij=P(Yij=1|θi,aj,bj)=

Φ(ajθi-bj)，

用三參數的probit連接函數來表示第i個人答對第j個題目的概率為

pij=P(Yij=1|θi,aj,bj,cj)=

cj+(1-cj)Φ(ajθi-bj),

式中：Φ----正態累積分布函數。

1.2 模型參數的貝葉斯估計

對于單參數的連接函數logistic模型，假設yi=(yi,1,yi,2,…,yi,J)′，第i個人反應部分的似然函數為

雙參數的連接函數logistic模型的似然函數為

(1)

三參數的連接函數logistic模型的似然函數為

f(yi|θi)=cj+

當單參數的連接函數為probit模型時，第i個人反應部分的似然函數為

雙參數的連接函數為probit模型的似然函數

三參數的連接函數為probit模型的似然函數

f(yi|θi)=cj+

根據Linden W[10]的思想進行局部獨立性假設,令

γ=(a1,a2,…,aJ;b1,b2,…,bJ;θ1,θ2,…,θN)′。

我們以連接函數為logistic的雙參數項目反應模型為例進行計算，利用式(1)得到N個人的聯合似然函數為

γ的后驗分布為

(2)

式中：m(y)----正則化常數(normalizing constant)，

π(aj)----aj的先驗;

π(bj)----bj的先驗;

π(θi)----θi的先驗。

它們具體的先驗在前面部分已給出。

1.3 模型選擇方法

文中用到的模型選擇方法為DIC準則和LPML。由于DIC可以在多元模型的固定或隨機部分不同的模型之間進行選擇，而不必指定模型參數的數量。DIC是由Spiegelhalter D J等[11]提出的,它是一種可以用來模型擬合以及模型復雜性測量的方法，DIC準則是在偏差的后驗分布基礎上建立的，其函數表達式為：

PD----有效參數的個數;

另一個評估模型相對擬合效果的方法是基于Chen M H等[13]、Gelfand A E等[14]的思想，通過采用條件預測縱坐標(Conditional Predictive Ordinate, CPO)指標計算LPML，設

則CPO的蒙特卡羅估計為如下形式

式中：r----MCMC算法的第r次迭代，r=1,2,…,R;

R----總迭代次數。

LPML值越大，表明所選擇的模型效果越好。

2 實證研究

2.1 數據結構

采用國際學生評估項目2015年數據(PISA)。數據中有548(N=548)個學生，每個人回答16(J=16)道題的考試，計算機會自動記錄并存儲被試者答題正確與否。基于單參數、雙參數和三參數的項目反應模型，進行30 000次迭代，燒掉前25 000次。被試者答對題目正確與否的部分數據見表1。

表1 被試者答對題目正確與否的數據

2.2 項目反應模型的nimble實現(以三參數logistic模型為例)

res6 <- read.csv("res61.csv", header=FALSE)

res6=as.matrix(res6)

N=length(res6[,1])

J=length(res6[1,])

nitem=rep(J,N)

cuni=cumsum(nitem)

id=cuni-nitem

item_index=matrix(,N,J)

for (i in 1:J)

{item_index[,i]=i}

y=rep(0,J*N)

for (i in 1:N){

y[((i-1)*J+1):(i*J)]=res6[i,]

}

velocity=rep(1,N*J)

irt<-nimbleCode({

for (j in 1:J){

for(i in 1:N){

y[id[i]+j]～dbern(prob[id[i]+j])

prob[id[i]+j]<-c[j]*velocity[id[i]+j]+(1-c[j]*velocity[id[i]+j])*(exp(a[j]*velocity[id[i]+j]*theta[i]-b[j]*velocity[id[i]+j])/(1+exp(a[j]*velocity[id[i]+j]*theta[i]-b[j]*velocity[id[i]+j])))

}

a[j]～dlnorm(0, sdlog=1)

b[j]～dnorm(0, sd=1)

c[j]～dbeta(4,shape1 = 12)

}

for (i in 1:N){

theta[i]～dnorm(0,1)

}

})

data <- list( y =y)

constants <- list(J=J,N = N,id=id,nitem=nitem,velocity=velocity)

inits <- list(theta=rep(0,N),a=rep(0.5,J),b=rep(0,J),c=rep(0.1,J))

jointmodel <- nimbleModel(irt,

data = data,

constants = constants,

inits = inits,

check = FALSE)

mcmc.out2<- nimbleMCMC(model = jointmodel,

niter = 30000, nchains = 1, nburnin =25000,thin=1,monitors = c('a','b','c','theta'), summary = TRUE,WAIC = TRUE,setSeed = 10)

第一部分程序中的N表示考生人數，J表示題目數量，為了讓模型的維度一致，方便編程，所以引入了velocity[id[i]+j]，應用nimble包，先寫出三參數logistic模型，然后給出參數服從的分布以及初值，其中，項目區分度參數logaj～N(0,1)，猜測參數bj～N(0,1)，迭代了30 000次，燒掉了前25 000次。

linsam=read.csv(file='res_only_logistic_sam_lin1.csv')[,-1]

loglisampler=rep(0,nusample)

prob<- cij+(1-cij)*(exp(aijtheta-bij)/(1+exp(aijtheta-bij)))

coran=which(y==1)

loglikelir= sum(log(prob[coran]))+sum(log(1-prob)[-coran])

loglikelir=(-2)*loglikelir

bigtheta3= bigzeta3=rep(0,cuni[N])

cpoi=matrix(rep(0,N*5000),5000)

probb=matrix(rep(0,nusample*cuni[N]),nusample)

loglisampler=rep(0,nusample)

for (k in 1:5000){

aj<-as.numeric(linsam[k,1:16])

bj<-as.numeric(linsam[k,17:32])

cj<-as.numeric(linsam[k,33:48])

theta<- as.numeric(linsam[k,49:(49+N-1)])

bij=rep(bj,N)

aijtheta=c(t(outer(theta,aj,FUN="*")))

cij=rep(cj,N)

for (qq in 1:(N-1)){

beg1=cuni[qq]+1

end1=cuni[qq+1]

bigtheta3[beg1:end1]=rep(theta[qq+1],nitem[qq+1])

}

beg1=1

end1=nitem[1]

bigtheta3[beg1:end1]=rep(theta[1],nitem[1])

probb[k,]<- cij+(1-cij)*(exp(aijtheta-bij)/(1+exp(aijtheta-bij)))

for (i in 1:(N-1))

{ beg1=cuni[i]+1

end1=cuni[i+1]

aa=which(y[beg1:end1]==1)

bb=probb[k,beg1:end1]

cpoi[k,i+1]=1/(prod(bb[aa])*prod(1-bb[-aa]))

}

beg1=1

end1=nitem[1]

aa=which(y[beg1:end1]==1)

bb=probb[k,beg1:end1]

cpoi[k,1]=1/(prod(bb[aa])*prod(1-bb[-aa]))

loglisampler[k]= sum(log(probb[k,][coran]))+sum(log(1-probb[k,])[-coran]) }#k

loglisampler=(-2)*loglisampler

EDEr1=mean(loglisampler);

PDICr1=EDEr1-loglikelir

DICr1=loglikelir+2*PDICr1

LPML<- sum(log(1/apply(cpoi,2,mean)))

第二部分程序是為了計算三參數logistic模型的DIC、PD和LPML值，linsam是讀取第一部分產生的項目反應模型的MCMC，loglikelir是計算這個模型的似然函數，最后通過for循環計算DIC、PD和LPML值。

2.3 結果分析

不同情況下的DIC、PD和LPML值見表2。

表2 不同情況下的DIC、PD和LPML值

表中通過比較這幾種情況的DIC值和LPML值發現，連接函數為logistic的雙參數項目反應模型的DIC值最小，LPML值最大，模型表現最好。單參數的項目反應模型的DIC值都較大，模型表現不好。三參數的項目反應模型比雙參數項目反應模型增加了偽猜測參數，正常情況下，DIC值應該增加32左右，但結果增加了100多，所以增加偽猜測參數的項目反應模型表現也不好。

3 結 語

基于國際學生評估項目2015年數據(PISA)，采用單參數、雙參數和三參數的項目反應模型分別進行建模，同時在單參數模型、雙參數模型和三參數模型下，分別考慮logistic連接函數和probit連接函數，通過計算6種模型的DIC值和LPML值得出結論，連接函數為logistic雙參數項目反應模型的DIC值最小，并且LPML值最大，此模型表現最好。