Como especificar o modelo de efeitos mistos bayesiano no BUGS

Postei isso no início da semana e retirei a pergunta quando encontrei uma boa fonte, não querendo perder tempo das pessoas. Não tenho feito muito progresso, temo. Ao tentar ser um bom cidadão aqui, tornarei o problema o mais claro possível. Eu suspeito que haverá poucos compradores.

Eu tenho um dataframe no RI que desejo analisar em BUGS ou R. É em formato longo. Consiste em várias observações em 120 indivíduos, com um total de 885 linhas. Estou examinando a ocorrência de um resultado categórico - mas isso não é realmente relevante aqui. A questão é sobre algo mais profundo.

O modelo que tenho usado até aqui é

  mymodel<-gee(Category ~ Predictor 1 + Predictor 2..family=binomial(link="logit"),
  data=mydata, 
   id=Person)

com um modelo marginal que explica essencialmente o agrupamento de pacientes. Eu então examinei

 mymodel<-gee(Category ~ Predictor 1 + Predictor 2.. , family=binomial(link="logit"), 
  corstr = "AR-M", 
  data=mydata, id=Person)

a fim de levar em consideração o tempo de ordenação das observações sobre as pessoas individualmente.

Isso não mudou muito.

Depois, tentei modelá-los usando o seguinte conjunto de comandos do MCMCPack:

 mymodel<-MCMCglmm(category~  Predictor1 + Predictor2..,  
 data=mydata, family=binomial(link="logit"))

Um exame do resultado foi emocionante, mostrando significância estatística para muitos preditores. Eu me elogiei como um bayesiano recém-convertido, até que percebi que não havia explicado medidas repetidas em pacientes.

Eu entendo que tenho que explicar isso. Entendo que isso pode significar ajustar uma hiperprioridade para cada indivíduo - está certo? Que forma isso terá nos BUGS?

Aqui está um modelo básico de registro de log: (parabéns a Kruschke, J., Indiana)

    model {
  for( i in 1 : nData ) {
    y[i] ~ dbern( mu[i] )
    mu[i] <- 1/(1+exp(-( b0 + inprod( b[] , x[i,] ))))
  }
  b0 ~ dnorm( 0 , 1.0E-12 )
   for ( j in 1 : nPredictors ) {
    b[j] ~ dnorm( 0 , 1.0E-12 )
  }
}

No entanto, não há hiperprior aqui para o indivíduo. Aqui está minha melhor tentativa até agora de um design dentro do indivíduo, respondendo a medidas repetidas dentro das pessoas:

Aqui está o modelo de Jackman para JAGS

1 model{
2 ## loop over data for likelihood
3 for(i in 1:n){
4  y[i] ~ dbern( mu[i] )
    mu[i] <- 1/(1+exp(-( b0 + inprod( b[] , x[i,] ))))
6 }
7 sigma ˜ dunif(0,20) ## prior on standard deviation
8 tau <- pow(sigma,-2) ## convert to precision
9
10 ## hierarchical model for each state’s intercept & slope
11 for(p in 1:50){
12 beta[p,1:2] ˜ dmnorm(mu[1:2],Tau[,]) ## bivariate normal
13 }
14
15 ## means, hyper-parameters
16 for(q in 1:2){
17 mu[q] ˜ dnorm(0,.0016)

}

Aqui está o meu modelo filho bastardo para BUGS

1 model{
2 ## loop over data for likelihood
3 for(i in 1:n){
4 mu.y[i] <- alpha + beta[s[i],1] + beta[s[i],2]*(j[i]-jbar)
5 demVote[i] ˜ dnorm(mu.y[i],tau)
6 }
7 sigma ˜ dunif(0,20) ## prior on standard deviation
8 tau <- pow(sigma,-2) ## convert to precision
9
10 ## hierarchical model for each state’s intercept & slope
11 for(p in 1:120){
12 beta[p,1:2] ˜ dmnorm(mu[1:2],Tau[,]) ## bivariate normal
13 }
14
15 ## means, hyper-parameters
16 for(q in 1:2){
17 mu[q] ˜ dnorm(0,.0016)
  }

Alguém pode me informar se estou indo na direção certa. Minha compreensão disso está crescendo, mas lentamente. Por favor, seja gentil. Eu sou um médico, não uma estatística! Eu usei o R um pouco, mas sou novo no BUGS e novo no Bayes.

Obrigado,

R

Você está (estava) quase lá. Apenas alguns comentários - você não precisa tornar o anterior para os beta[,1:2]parâmetros um MV comum normal; você pode fazer a prévia de tal forma que beta[i,1]e beta[i,2]são independentes, o que simplifica as coisas (por exemplo, não covariância antes precisa ser especificado.) Note que isso não significa que será independente no posterior.

Outros comentários: Como você tem um termo constante - alpha- na regressão, os componentes beta[,1]devem ter média zero no anterior. Além disso, você não tem um prior alphano código.

Aqui está um modelo com termos de interceptação hierárquica e inclinação; Tentei manter seus antecedentes e notação sempre que possível, dadas as alterações:

model {
  for(i in 1:n){
    mu.y[i] <- alpha + beta0[s[i]] + beta1[s[i]]*(j[i]-jbar)
    demVote[i] ~ dnorm(mu.y[i],tau)
  }

  alpha ~ dnorm(0, 0.001) ## prior on alpha; parameters just made up for illustration
  sigma ~ dunif(0,20) ## prior on standard deviation
  tau <- pow(sigma,-2) ## convert to precision

  ## hierarchical model for each state’s intercept & slope
  for (p in 1:120) {
     beta0[p] ~ dnorm(0, tau0)
     beta1[p] ~ dnorm(mu1, tau1)
  }

  ## Priors on hierarchical components; parameters just made up for illustration
  mu1 ~ dnorm(0, 0.001) 
  sigma0 ~ dunif(0,20)
  sigma1 ~ dunif(0,20)
  tau0 <- pow(sigma0,-2)
  tau1 <- pow(sigma1,-2)
}

Um recurso muito útil para modelos hierárquicos, incluindo alguns "truques" para acelerar a convergência, é Gelman e Hill .

(Um pouco tarde com a resposta, mas pode ser útil para alguns futuros questionadores.)