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

(長(zhǎng)春工業(yè)大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院， 吉林 長(zhǎng)春 130012)

0 引 言

對(duì)于項(xiàng)目反應(yīng)模型，從連接函數(shù)角度來看，常用的連接函數(shù)有l(wèi)ogistic和probit；從參數(shù)角度來看，有單參數(shù)、雙參數(shù)和三參數(shù)的項(xiàng)目反應(yīng)模型，并且它們均可以用logistic和probit連接函數(shù)。Leventhal Brian C等[1]、張雪[2]、張國(guó)紅[3]、Bolsinova M等[4]項(xiàng)目反應(yīng)模型的連接函數(shù)為logistic連接函數(shù)。Lu J等[5]、Bolsinova M等[6]項(xiàng)目反應(yīng)模型是probit連接函數(shù)。

文中首先介紹了單參數(shù)、雙參數(shù)和三參數(shù)的項(xiàng)目反應(yīng)模型。其次，分別采用連接函數(shù)logistic和probit對(duì)單參數(shù)、雙參數(shù)和三參數(shù)的項(xiàng)目反應(yīng)模型進(jìn)行計(jì)算并比較。應(yīng)用了DIC和LPML兩個(gè)模型方法，DIC小說明模型好，LPML大說明模型好。最后，在實(shí)證研究部分，通過計(jì)算每個(gè)模型的DIC值和LPML值得出結(jié)論，當(dāng)連接函數(shù)為logistic的雙參數(shù)模型最好，因?yàn)镈IC值最小，同時(shí)LPML值最大。由于貝葉斯計(jì)算項(xiàng)目反應(yīng)模型存在一定困難，鏈特別長(zhǎng)，所以針對(duì)實(shí)證研究的程序是由nimble包編寫的，nimble包的語法結(jié)構(gòu)類似于WinBUGS和JAGS，但是比WinBUGS和JAGS更靈活，因?yàn)榫幊虝r(shí)用C++代碼，所以運(yùn)算速度快[7]。

1 模型分析

1.1 模型的建立

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

用雙參數(shù)的logistic連接函數(shù)表示第i個(gè)人答對(duì)第j個(gè)題目的概率為

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

用三參數(shù)的logistic連接函數(shù)表示第i個(gè)人答對(duì)第j個(gè)題目的概率為

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

式中：aj、bj、cj----分別表示項(xiàng)目j的區(qū)分度參數(shù)、難度參數(shù)和偽猜測(cè)參數(shù)。

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

此外，還可用單參數(shù)的probit連接函數(shù)表示第i個(gè)人答對(duì)第j個(gè)題目的概率為[9]

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

用雙參數(shù)的probit連接函數(shù)表示第i個(gè)人答對(duì)第j個(gè)題目的概率為

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

Φ(ajθi-bj)，

用三參數(shù)的probit連接函數(shù)來表示第i個(gè)人答對(duì)第j個(gè)題目的概率為

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

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

式中：Φ----正態(tài)累積分布函數(shù)。

1.2 模型參數(shù)的貝葉斯估計(jì)

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

雙參數(shù)的連接函數(shù)logistic模型的似然函數(shù)為

(1)

三參數(shù)的連接函數(shù)logistic模型的似然函數(shù)為

f(yi|θi)=cj+

當(dāng)單參數(shù)的連接函數(shù)為probit模型時(shí)，第i個(gè)人反應(yīng)部分的似然函數(shù)為

雙參數(shù)的連接函數(shù)為probit模型的似然函數(shù)

三參數(shù)的連接函數(shù)為probit模型的似然函數(shù)

f(yi|θi)=cj+

根據(jù)Linden W[10]的思想進(jìn)行局部獨(dú)立性假設(shè),令

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

我們以連接函數(shù)為logistic的雙參數(shù)項(xiàng)目反應(yīng)模型為例進(jìn)行計(jì)算，利用式(1)得到N個(gè)人的聯(lián)合似然函數(shù)為

γ的后驗(yàn)分布為

(2)

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

π(aj)----aj的先驗(yàn);

π(bj)----bj的先驗(yàn);

π(θi)----θi的先驗(yàn)。

它們具體的先驗(yàn)在前面部分已給出。

1.3 模型選擇方法

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

PD----有效參數(shù)的個(gè)數(shù);

另一個(gè)評(píng)估模型相對(duì)擬合效果的方法是基于Chen M H等[13]、Gelfand A E等[14]的思想，通過采用條件預(yù)測(cè)縱坐標(biāo)(Conditional Predictive Ordinate, CPO)指標(biāo)計(jì)算LPML，設(shè)

則CPO的蒙特卡羅估計(jì)為如下形式

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

R----總迭代次數(shù)。

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

2 實(shí)證研究

2.1 數(shù)據(jù)結(jié)構(gòu)

采用國(guó)際學(xué)生評(píng)估項(xiàng)目2015年數(shù)據(jù)(PISA)。數(shù)據(jù)中有548(N=548)個(gè)學(xué)生，每個(gè)人回答16(J=16)道題的考試，計(jì)算機(jī)會(huì)自動(dòng)記錄并存儲(chǔ)被試者答題正確與否。基于單參數(shù)、雙參數(shù)和三參數(shù)的項(xiàng)目反應(yīng)模型，進(jìn)行30 000次迭代，燒掉前25 000次。被試者答對(duì)題目正確與否的部分?jǐn)?shù)據(jù)見表1。

表1 被試者答對(duì)題目正確與否的數(shù)據(jù)

2.2 項(xiàng)目反應(yīng)模型的nimble實(shí)現(xiàn)(以三參數(shù)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表示考生人數(shù)，J表示題目數(shù)量，為了讓模型的維度一致，方便編程，所以引入了velocity[id[i]+j]，應(yīng)用nimble包，先寫出三參數(shù)logistic模型，然后給出參數(shù)服從的分布以及初值，其中，項(xiàng)目區(qū)分度參數(shù)logaj～N(0,1)，猜測(cè)參數(shù)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)))

第二部分程序是為了計(jì)算三參數(shù)logistic模型的DIC、PD和LPML值，linsam是讀取第一部分產(chǎn)生的項(xiàng)目反應(yīng)模型的MCMC，loglikelir是計(jì)算這個(gè)模型的似然函數(shù)，最后通過for循環(huán)計(jì)算DIC、PD和LPML值。

2.3 結(jié)果分析

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

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

表中通過比較這幾種情況的DIC值和LPML值發(fā)現(xiàn)，連接函數(shù)為logistic的雙參數(shù)項(xiàng)目反應(yīng)模型的DIC值最小，LPML值最大，模型表現(xiàn)最好。單參數(shù)的項(xiàng)目反應(yīng)模型的DIC值都較大，模型表現(xiàn)不好。三參數(shù)的項(xiàng)目反應(yīng)模型比雙參數(shù)項(xiàng)目反應(yīng)模型增加了偽猜測(cè)參數(shù)，正常情況下，DIC值應(yīng)該增加32左右，但結(jié)果增加了100多，所以增加偽猜測(cè)參數(shù)的項(xiàng)目反應(yīng)模型表現(xiàn)也不好。

3 結(jié) 語

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