Конвертировать многоуровневые модели в широкий формат

Question

У меня многоуровневая модель зазубрин. Я пытаюсь преобразовать его из широкого в длинный формат, как описано здесь: http://jeromyanglim.tumblr.com/post/37361593128/jags-converting-multilevel-model-from-wide-to Однако моя модель более сложная, чем пример, поэтому у меня возникли проблемы с выполнением этой работы. Чтобы проиллюстрировать трудности, я сделал повторяемый пример. Этот первый блок создает данные и устанавливает параметры зазубрин:

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, но я не вижу, где я иду не так. Может кто-нибудь помочь, пожалуйста? Благодарю.

Answer 1

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

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

Matt