使用R2OpenBUGS拟合WinBUGS模型的语法错误排查求助

问题描述 #

使用R2OpenBUGS拟合WinBUGS模型时，首先遇到如下解析错误：

Error in parse(text = x, srcfile = src): <text>:5:44: unexpected symbol
4:         for(j in 1:P){
5:             y[i,j] ~ dnorm(mu[i,j], psi[j])I

显然问题出在第5行的语法，原代码行是：

y[i, j] ~ dnorm(mu[i, j], psi[j])I(thd[j, z[i, j]], thd[j, z[i, j] + 1])

尝试修改为R风格的写法：

y[i, j] ~ dnorm(mu[i, j], psi[j]) * I(thd[j, z[i, j]], thd[j, z[i, j] + 1])

解析错误消失，但模型运行崩溃，日志提示：

"invalid or unexpected token scanned error pos 100 (error on line 5)"

完整模型代码 #

model <- function() {
    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] <- xi[i,1]
        mu[i,4] <- lam[2] * xi[i,1]
        mu[i,5] <- lam[3] * 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] <- lam[7] * xi[i,1]
        mu[i,10] <- xi[i,2]
        mu[i,11] <- lam[8] * xi[i,2]
        mu[i,12] <- lam[9] * xi[i,2]
        mu[i,13] <- lam[10] * xi[i,2]
        mu[i,14] <- lam[11] * xi[i,2]
        mu[i,15] <- lam[12] * xi[i,2]
        mu[i,16] <- xi[i,3]
        mu[i,17] <- lam[13] * xi[i,3]
        mu[i,18] <- lam[14] * xi[i,3]
        mu[i,19] <- xi[i,4]
        mu[i,20] <- lam[15] * xi[i,4]
        mu[i,21] <- lam[16] * xi[i,4]
        mu[i,22] <- lam[17] * xi[i,4]
        mu[i,23] <- lam[18] * xi[i,4]
        mu[i,24] <- lam[19] * xi[i,4]
        mu[i,25] <- lam[20] * xi[i,4]
        mu[i,26] <- lam[21] * 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.0
    }

    #priors on loadings and coefficients
    var.lam[1] <- 4.0 * psi[2]
    var.lam[2] <- 4.0 * psi[4]
    var.lam[3] <- 4.0 * psi[5]
    var.lam[4] <- 4.0 * psi[6]
    var.lam[5] <- 4.0 * psi[7]
    var.lam[6] <- 4.0 * psi[8]
    var.lam[7] <- 4.0 * psi[9]
    var.lam[8] <- 4.0 * psi[11]
    var.lam[9] <- 4.0 * psi[12]
    var.lam[10] <- 4.0 * psi[13]
    var.lam[11] <- 4.0 * psi[14]
    var.lam[12] <- 4.0 * psi[15]
    var.lam[13] <- 4.0 * psi[17]
    var.lam[14] <- 4.0 * psi[18]
    var.lam[15] <- 4.0 * psi[20]
    var.lam[16] <- 4.0 * psi[21]
    var.lam[17] <- 4.0 * psi[22]
    var.lam[18] <- 4.0 * psi[23]
    var.lam[19] <- 4.0 * psi[24]
    var.lam[20] <- 4.0 * psi[25]
    var.lam[21] <- 4.0 * psi[26]
    
    for(i in 1:21){
        lam[i] ~ dnorm(0.8,var.lam[i])
    }

    var.gam <- 4.0 * psd
    gam[1] ~ dnorm(0.6,var.gam)
    gam[2] ~ dnorm(0.6,var.gam)
    gam[3] ~ dnorm(0.4,var.gam)
    gam[4] ~ dnorm(0.4,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

问题根源与解决方法 #

核心问题是WinBUGS/OpenBUGS的截断分布语法与R完全不同：

1. 正确的截断分布写法是在分布后直接加I(下限;上限)，用分号分隔上下限，而非逗号；且分布与I()之间需要加空格。
2. 不能用R的*来连接分布和截断条件，这在BUGS语言中是非法语法。

修正后的第5行代码应为：

y[i,j] ~ dnorm(mu[i,j], psi[j]) I(thd[j, z[i,j]]; thd[j, z[i,j]+1])

额外检查点 #

• 确认传入模型的thd、z等变量维度匹配，z[i,j]的取值需保证z[i,j]+1不超出thd的列数，避免索引越界。
• 检查phi先验的dwish分布中，R矩阵必须是正定矩阵，否则会导致模型运行失败。
• 核对var.lam与psi的索引对应关系，确保符合模型的测量结构设计。

内容的提问来源于stack exchange，提问作者Pål Bjartan