Ошибка в jags.model

Question

Я новичок в R и JAGS, и даже не очень опытный в программировании.

Я пытаюсь настроить иерархическую модель для некоторых данных, но я получаю эту ошибку:

####error in line above:
Error in jags.model(file = "modelControl.txt", data = dataList, inits = 
initsList,  : 
  RUNTIME ERROR:
Cannot insert node into beta1[1:158]. Dimension mismatch

В конце выполнения кода ниже. Что я делаю неправильно? Как я могу избежать ошибки?

#declaration
y=as.numeric(PANSS[treat == 0]) #treat == 0 indicates control group
x=as.numeric(time[treat == 0])
meanYcontrol = mean(PANSS[treat == 0])
sdYcontrol = sd(PANSS[treat == 0])
s = pp[treat==0]

#data list
dataList = list(
  y = y,
  x = x,
  Ntotal = length(y),
  Nsubj = length(y)/6 , #each subject had 6 test moments
  s = s
)
#model
modelString = "
model {
for ( i in 1:Ntotal ) {
y[i] ~ dnorm( beta0[s[i]] + beta1[s[i]] * x[i,1], 1/sigma^2 ) 
}
for ( j in 1:Nsubj ) {
beta0[j] ~ dnorm( beta0mu , 1/(beta0sigma)^2 ) 
beta1[j] ~ dnorm( beta1mu , 1/(beta1sigma)^2 )
}
#vague priors
beta0 ~ dnorm( 0, 1/(10)^5 )  
beta1 ~ dt( 0, 1, 1 ) #Cauchy distribution
beta0sigma ~ dunif( 1.0E-5, 1.0E+5 )
beta1sigma ~ dunif( 1.0E-5, 1.0E+5 )
sigma ~ dunif( 1.0E-5, 1.0E+5 )
nu = nuMinusOne+1
nuMinusOne ~ dexp(1/29)
}
"

 #write model to text file
writeLines(modelString, con="modelControl.txt")

#initialization chains
beta0Init = meanYcontrol
beta1Init = 0
sigmaInit = sdYcontrol
initsList = list(beta0=beta0Init, beta1=beta1Init, sigma=sigmaInit)

#run chains
parameters = c("beta0", "beta1", "sigma") #parameters to be monitored
numSavedSteps = 7500 #number of steps in chain to save
adaptSteps = 1000  #number of steps to tune the samplers
burnInSteps = 500 #number of steps to burn-in the samplers
thinSteps = 1 #number of steps to keep (1=keep every step)
nChains = 3 #number of chains to run
nIter = ceiling(numSavedSteps / nChains) #number of steps per chain

jagsModel = jags.model(file="modelControl.txt", data=dataList, inits = 
initsList, n.chains=nChains, n.adapt=adaptSteps)

Answer 1

Ваши комментарии подразумевают, что вы изменили несколько вещей в вашей модели, но определение модели в вашем вопросе не изменилось (или, возможно, было частично отредактировано). Вот что я думаю, что вы сейчас используете:

model {
for ( i in 1:Ntotal ) {
y[i] ~ dnorm( beta0[s[i]] + beta1[s[i]] * x[i,1], 1/sigma^2 ) 
}
for ( j in 1:Nsubj ) {
beta0[j] ~ dnorm( beta0mu , 1/(beta0sigma)^2 ) 
beta1[j] ~ dnorm( beta1mu , 1/(beta1sigma)^2 )
}
#vague priors
beta0mu ~ dnorm( 0, 1/(10)^5 )  
beta1mu ~ dt( 0, 1, 1 ) #Cauchy distribution
beta0sigma ~ dunif( 1.0E-5, 1.0E+5 )
beta1sigma ~ dunif( 1.0E-5, 1.0E+5 )
sigma ~ dunif( 1.0E-5, 1.0E+5 )
nu = nuMinusOne+1
nuMinusOne ~ dexp(1/29)
}

# In R:
meanYcontrol = mean(PANSS[treat == 0])
sdYcontrol = sd(PANSS[treat == 0])
beta0Init = meanYcontrol
beta1Init = 0
sigmaInit = sdYcontrol
initsList = list(beta0=beta0Init, beta1=beta1Init, sigma=sigmaInit)

Таким образом, вы даете начальные значения для beta0 и beta1, но я думаю, что они предназначены для beta0mu и beta1mu. Впрочем, возможны и другие ошибки - это сложно проверить, потому что в данный момент мы не можем запустить модель, поскольку у нас нет ваших данных. В будущем было бы неплохо представить минимальный воспроизводимый пример, включая все необходимые данные и т. Д., Поскольку это поможет вам быстрее и полнее генерировать ответы.

Один из способов избежать подобных ошибок - использовать теги # inits # и # data # с пакетом runjags, что может сделать данные и начальные значения более понятными, например:

model {
    for ( i in 1:Ntotal ) {   #data# Ntotal
        y[i] ~ dnorm( beta0[s[i]] + beta1[s[i]] * x[i,1], 1/sigma^2 ) 
            #data# y, x, s
    }
    for ( j in 1:Nsubj ) {  #data# Nsubj
        beta0[j] ~ dnorm( beta0mu , 1/(beta0sigma)^2 ) 
        beta1[j] ~ dnorm( beta1mu , 1/(beta1sigma)^2 )
            #inits# beta0mu, beta1mu
    }
    #vague priors
    beta0mu ~ dnorm( 0, 1/(10)^5 )  
    beta1mu ~ dt( 0, 1, 1 ) #Cauchy distribution
    beta0sigma ~ dunif( 1.0E-5, 1.0E+5 )
    beta1sigma ~ dunif( 1.0E-5, 1.0E+5 )
    sigma ~ dunif( 1.0E-5, 1.0E+5 )
        #inits# sigma
    nu = nuMinusOne+1
    nuMinusOne ~ dexp(1/29)
}

[Обратите внимание, что я также использовал отступы, чтобы облегчить чтение модели - это часто помогает выявлять ошибки, даже не публикуя модель]

Затем в R вам просто нужно указать данные и единицы измерения, которые должны быть в вашей рабочей среде, а затем использовать:

results <- runjags::run.jags("modelControl.txt", monitor=...)

Надеюсь, это поможет.