Преобразование многоуровневых моделей затяжек от широкого до длинного формата

Question

У меня есть многоуровневая модель. Я пытаюсь преобразовать его из широкого формата в длинный, как описано здесь: 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. 

Ok, так что следующий раздел является «широкий формат» зазубрины модель Thats обеспечивает правильные оценки в объекте 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, но я не могу видеть, куда я иду неправильно. Может ли кто-нибудь помочь? Благодарю.

Answer 1

В вашей матрице df1 есть элементы nData * nVars, но ваша длинная форматная модель использует только первые элементы nData (т. Е. Вы просто используете первый столбец данных). Максимум для основного цикла данных должен быть скорректирован так, чтобы он был равен nData * nVars, а не только nData.

Также вам нужен вектор, представляющий номер строки исходного df1, чтобы вы могли правильно индексировать ваш случайный эффект b0s, например. b0s [dfrow [I]]. Кроме того, трудно выполнить спецификацию данных (например, то, что является длиной (out.names)), поэтому я не уверен, что вы уже это сделали, но либо x1 нужно повторять nVars раз, либо вы должны использовать тот же x1 [dfrow [i]] индексирование как для случайного эффекта (желательно последнее для удобочитаемости вашего кода модели).

Мэтт