0

我目前正在使用 WinBUGS 研究具有有序分类变量的 SEM（结构方程建模）。我对 WinBUGS 比较陌生，所以我按照Lee 和 Song 的书5.2.4 节中的代码进行操作。我的目标是将数据拟合到指定的模型（参见下面的 WinBUGS 代码）并获得参数的估计值。

下面提供了我的 WinBUGS 代码。数据集太大（Sample size 1986 with 19 variables）无法直接包含在下面的代码中，所以我在这里附上了两个数据集z.data和xi.data。

model{
    for(i in 1:N){
    #measurement equation model
        for(j in 1:P){
        y[i,j]~dnorm(mu[i,j],psi[j])I(thd[j,z[i,j]],thd[j,z[i,j]+1])
        ephat[i,j]<-y[i,j]-mu[i,j]
        }
    mu[i,1]<-eta[i]
    mu[i,2]<-lam[1]*eta[i]
    mu[i,3]<-lam[2]*eta[i]
    mu[i,4]<-lam[3]*eta[i]
    mu[i,5]<-xi[i,1]
    mu[i,6]<-lam[4]*xi[i,1]
    mu[i,7]<-lam[5]*xi[i,1]
    mu[i,8]<-lam[6]*xi[i,1]
    mu[i,9]<-xi[i,2]
    mu[i,10]<-lam[7]*xi[i,2]
    mu[i,11]<-xi[i,3]
    mu[i,12]<-lam[8]*xi[i,3]
    mu[i,13]<-lam[9]*xi[i,3]
    mu[i,14]<-xi[i,4]
    mu[i,15]<-lam[10]*xi[i,4]
    mu[i,16]<-lam[11]*xi[i,4]
    mu[i,17]<-lam[12]*xi[i,4]
    mu[i,18]<-lam[13]*xi[i,4]
    mu[i,19]<-lam[14]*xi[i,4]

    #structural equation model
    xi[i,1:4]~dmnorm(u[1:4],phi[1:4,1:4])
    eta[i]~dnorm(nu[i],psd)
    nu[i]<-gam[1]*xi[i,1]+gam[2]*xi[i,2]+gam[3]*xi[i,3]+gam[4]*xi[i,4]
    dthat[i]<-eta[i]-nu[i]
}#End of i

    for(i in 1:4){u[i]<-0}

    #priors on loadings and coefficients
    var.lam[1]<-4*psi[2]    var.lam[2]<-4*psi[3]    var.lam[3]<-4*psi[4]
    var.lam[4]<-4*psi[6]    var.lam[5]<-4*psi[7]    var.lam[6]<-4*psi[8]
    var.lam[7]<-4*psi[10]   var.lam[8]<-4*psi[12]   var.lam[9]<-4*psi[13]
    var.lam[10]<-4*psi[15]  var.lam[11]<-4*psi[16]  var.lam[12]<-4*psi[17]
    var.lam[13]<-4*psi[18]  var.lam[14]<-4*psi[19]

    for(i in 1:14){lam[i]~dnorm(0.8,var.lam[i])}

    var.gam<-4*psd
    gam[1]~dnorm(0.5,var.gam)   gam[2]~dnorm(0.5,var.gam)
    gam[3]~dnorm(0.5,var.gam)   gam[4]~dnorm(0.5,var.gam)

    #priors on precisions
    for(j in 1:P){
        psi[j]~dgamma(10,8)
        sgm[j]<-1/psi[j]
    }
    psd~dgamma(10,8)
    sgd<-1/psd
    phi[1:4,1:4]~dwish(R[1:4,1:4],30)
    phx[1:4,1:4]<-inverse(phi[1:4,1:4])

}#End of model

#Data set
list(N=1986,P=19,
    R=structure(
    .Data=c(8,0,0,0,
        0,8,0,0,
        0,0,8,0,
        0,0,0,8),
    .Dim=c(4,4)),
    thd=structure(
    .Data=c(-200,-1.24951,0.860111,2.07197,200,200,200,200,200,200,200,
    -200,-2.25443,-2.00170,-1.59472,-1.19878,-0.661082,-0.194325,0.423813,1.13399,1.51446,200,
    -200,-2.23869,-2.05085,-1.58136,-1.26340,-0.691222,-0.171223,0.457192,1.21446,1.63661,200,
    -200,-1.28040,-0.849199,-0.381381,-0.606203,0.330280,0.640802,1.05806,1.59472,1.97438,200,
    -200,-2.09407,-1.74742,-1.38135,-1.04489,-0.499650,-0.00504862,0.578479,1.27470,1.68680,200,
    -200,-1.83527,-1.09179,-0.448801,200,200,200,200,200,200,200,
    -200,-1.96559,-1.68680,-1.27470,-0.973485,-0.498221,-0.00883516,0.550357,1.30960,1.62233,200,
    -200,-1.10797,-0.193039,1.10797,200,200,200,200,200,200,200,
    -200,-1.49488,-0.0795994,1.39793,2.34291,200,200,200,200,200,200,
    -200,-1.98332,-1.18085,-0.450198,200,200,200,200,200,200,200,
    -200,-2.30533,-1.49488,-0.872971,200,200,200,200,200,200,200,
    -200,-1.20658,0.370544,1.51050,200,200,200,200,200,200,200,
    -200,-0.993961,0.860111,2.03060,200,200,200,200,200,200,200,
    -200,-1.32462,0.0353475,1.14608,200,200,200,200,200,200,200,
    -200,-1.61299,-1.28905,-0.858285,-0.472661,-0.0467165,0.320964,0.694431,1.35877,1.78962,200,
    -200,-1.74164,-0.743465,0.493939,1.48721,200,200,200,200,200,200,
    -200,-1.78962,-1.48341,-1.04707,-0.705720,-0.230473,0.298453,0.931761,1.57696,1.92379,200,
    -200,-1.46841,-1.12444,-0.662653,-0.349000,0.113839,0.611626,1.18593,1.77724,2.07197,200,
    -200,-1.44650,0.565103,1.91582,200,200,200,200,200,200,200),
    .Dim=c(19,11)),
    z=structure(
    .Data=c(paste z.data here),
    .Dim=c(1986,19)))

#Two different initial values
#Initials 1
list(
    lam=c(0,0,0,0,0,0,0,0,0,0,0,0,0,0),
    gam=c(0,0,0,0),
    psi=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1),
    psd=1,
    phi=structure(
    .Data=c(1,0,0,0,
        0,1,0,0,
        0,0,1,0,
        0,0,0,1),
    .Dim=c(4,4)),
    xi=structure(
    .Data=c(paste xi.data here),
    .Dim=c(1986,4)))

#Initials 2
list(
    lam=c(1,1,1,1,1,1,1,1,1,1,1,1,1,1),
    gam=c(1,1,1,1),
    psi=c(2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2),
    psd=2,
    phi=structure(
    .Data=c(2,0,0,0,
        0,2,0,0,
        0,0,2,0,
        0,0,0,2),
    .Dim=c(4,4)),
    xi=structure(
    .Data=c(paste xi.data here),
    .Dim=c(1986,4)))

所有 19 个变量都是在各种尺度上测量的有序分类变量，范围从四点量表到十点量表。我遵循书中所述的常见做法，使用 R 找到阈值并将它们作为“thd”包含在数据集中。

我希望上述代码应该在 WinBUGS 中成功运行。模型检查、数据加载和模型编译步骤没有问题。但是，当我继续在 WinBUGS 中加载两组初始值时，两条链都显示“此链包含未初始化的变量”消息，但我不知道哪些变量未初始化。当我尝试使用“gen inits”选项生成初始值时，在初始化模型时，会显示一个陷阱日志消息“未定义的实际结果”和一长串错误。我还用Lee 和 Song 的书5.2.4 节中的代码进行了试运行，但也出现了类似的错误（代码和数据集可以在同一个站点下载）。

我想到的问题是：

  1. 模型中的哪些变量未初始化？我以为我已经指定了所有必要的参数。
  2. 当我尝试使用“gen inits”选项生成初始值时，是什么导致了“未定义的实际结果”问题以及如何解决这个问题？

您能否提一些建议？非常感谢！

4

0 回答 0