1

我有一个多级锯齿模型。我正在尝试将其从宽格式转换为长格式，如下所述：http: //jeromyanglim.tumblr.com/post/37361593128/jags-converting-multilevel-model-from-wide-to但是我的模型比这个例子，所以我在完成这项工作时遇到了一些麻烦。为了说明困难，我做了一个可重复的例子。第一个块创建数据并设置 jags 参数：

library(ecodist)
library(runjags)
set.seed(10)

##### population n
n <- 250
# num outputs
num.ys <- 10

# Vector binary to indicate which domains have correlation with independent variables
corr.vec <- c(0, 0, 0, 1, 1, 0, 0, 1, 1, 1)
correlation = 0.99

# Function to simulate correlated outcome
sim.fn <- function(i, var1, sw1) {
    if(sw1 ==1){
        temp <- corgen(n , var1, correlation )
        temp <- as.numeric(temp$y * attr(temp$y,'scaled:scale') + attr(temp$y,'scaled:center'))
    } else {
        temp <- rnorm(n, 0, 5)
    }
    return(temp)
}

##### Generate data
df0 <- data.frame(var1=rnorm(n, 15, 2))
df1 <- data.frame(df0, sapply(1:num.ys, function(i) sim.fn(i, df0$var1, corr.vec[i])))

out.names <- paste0("y_", 1:num.ys)
names(df1) <- c("var1", out.names)

### Jags parameters
parameters = c("B1O", "b1", "b1o", "nu", "sd")
adaptSteps = 1000             # Number of steps to "tune" the samplers.
burnInSteps = 10000            # Number of steps to "burn-in" the samplers.
nChains = 2                   # Number of chains to run.
numSavedSteps=1000          # Total number of steps in chains to save.
thinSteps=2                   # Number of steps to "thin" (1=keep every step).
nPerChain = ceiling( ( numSavedSteps * thinSteps ) / nChains ) # Steps per chain.

好的，所以下一部分是“宽格式”锯齿模型，它在对象中提供了正确的估计mcmcChain：

modelstring = "
model {
for( i in 1 : nData ) {
    for(np in 1:nVars){
        y[i, np]  ~ dt( mu[i,np], tau, nu)
        mu[i, np] <- b0s[i] + (b1 + b1o[np]) * x1[i]
    }
}

#Random effects
for(i in 1:nData){
    b0s[i] ~ dnorm(0, b0stau)
}

#Outcome level
for (np in 1:nVars){
    b1o[np] ~ dnorm(0, b1otau)
}

##### Priors
#Overarching Level
b1 ~ dnorm(0, 0.0001)

#
b0stau <- pow(b0ssd, -2)
b0ssd  ~ dt(0, 1/625, 1)T(0,)

# tau & nu priors
nuI ~ dunif(0.001,0.5)
nu <- 1/nuI
tau <- pow(sd, -2)
sd ~ dunif(0, 10)

b1otau <- pow(b1osd, -2)
b1osd  ~ dt(0, 1/625, 1)T(0,)
b1dtau <- pow(b1dsd, -2)
b1dsd  ~ dt(0, 1/625, 1)T(0,)

#Transformations
for(np in 1:nVars){
    B1O[np] <- b1 + b1o[np]
}
}
" # close quote for modelstring
writeLines(modelstring,con="model.jags.no_dom.test.txt")

zy <- (df1[, out.names])
sc_ys <- data.frame(lapply(zy, function(x) scale(x)) )

dataList = list( y = as.matrix(sc_ys), x1 = as.numeric(scale(df1$var1,)), 
             nVars = num.ys, nData = nrow(df1))

# Run this model via run.jags
codaSamples <- run.jags(model="model.jags.no_dom.test.txt" , data=dataList , method ="parallel", n.chains=nChains, monitor=parameters,
                    adapt = adaptSteps, burnin = burnInSteps, sample=nPerChain, thin=thinSteps)

mcmcChain <- data.frame(summary( codaSamples ))
mcmcChain

因此，BO 输出接近生成数据的相关性。接下来是我对类似于上面链接中的解释的“长格式”模型的尝试。

modelstring = "
model {
for( i in 1 : nData ) {
    y[i]  ~ dt( mu[i] , tau, nu )
    mu[i] <- b0s[i] + (b1 + b1o[idx[i]]) * x1[i]
}

#Random effects
for(i in 1:nData){
    b0s[i] ~ dnorm(0, b0stau)
}

#Outcome level
for (y in 1:nVars){
    b1o[y] ~ dnorm(0, b1otau[y])
}

##### Priors
#Overarching Level
b1 ~ dnorm(0, 0.0001)

b0stau <- pow(b0ssd, -2)
b0ssd  ~ dt(0, 1/625, 1)T(0,)

for (y in 1:nVars){
    b1otau[y] <- pow(b1osd[y], -2)
    b1osd[y]  ~ dt(0, 1/625, 1)T(0,)
}

tau <- pow(sd, -2)
sd ~ dunif(0, 10)
nuI ~ dunif(0.001,0.5)
nu <- 1/nuI

#Transformations
for(j in 1:nVars){
    B1O[j]  <- b1 + b1o[j]
}
}
" # close quote for modelstring
writeLines(modelstring,con="model.jags.no_dom.long.test.txt")

# Restructure data into long format
dataList2 = list( y = unlist(sc_ys), x1 = rep (as.numeric(scale(df1$var1,)),  length(out.names)),
             idx = rep(1:length(out.names), each=nrow(df1)),
             nVars = length(out.names), nData = nrow(df1))

codaSamples2 <- run.jags(model="model.jags.no_dom.long.test.txt" , data=dataList2 , method ="parallel", n.chains=nChains, monitor=parameters,
                     adapt = adaptSteps, burnin = burnInSteps, sample=nPerChain, thin=thinSteps)

mcmcChain2 <- data.frame(summary( codaSamples2 ))
mcmcChain2

所以结果与 的结果mcmcChain2不匹配mcmcChain，但我看不出哪里出错了。任何人都可以帮忙吗？谢谢。

4

1 回答 1

1

您的矩阵 df1 具有 nData * nVars 元素，但您的长格式模型仅使用前 nData 元素（即实际上您只是使用数据的第一列）。主数据循环的最大值需要调整为等于 nData*nVars 而不仅仅是 nData。

您还需要一个表示原始 df1 行号的向量，以便您可以正确索引随机效应 b0s，例如 b0s[dfrow[i]]。此外，很难遵循数据规范（例如什么是长度（out.names））所以我不确定你是否已经这样做了，但要么 x1 需要重复 nVars 次，要么你应该使用相同的 x1 [dfrow[i]] 索引作为随机效应（为了模型代码的可读性，最好是后者）。

马特

于 2017-02-08T10:02:39.067 回答