Cómo especificar el modelo de efectos mixtos bayesianos en BUGS

Publiqué esto a principios de semana y luego retracté la pregunta cuando encontré una buena fuente, no queriendo perder el tiempo de las personas. No he progresado mucho, me temo. Al tratar de ser un buen ciudadano aquí, dejaré el problema lo más claro posible. Sospecho que habrá pocos tomadores.

Tengo un marco de datos en RI que deseo analizar en ERRORES o R. Está en formato largo. Consiste en múltiples observaciones en 120 individuos, con un total de 885 filas. Estoy examinando la aparición de un resultado categórico, pero eso no es realmente relevante aquí. La pregunta es sobre algo más profundo.

El modelo que he estado usando hasta aquí es

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

con un modelo marginal que explica esencialmente la agrupación de pacientes. Entonces examiné

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

para dar cuenta del tiempo ordenado de las observaciones en las personas individuales.

Esto no cambió mucho.

Luego intenté modelarlos usando el siguiente conjunto de comandos MCMCPack:

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

Un examen de la salida fue emocionante, mostrando significancia estadística para muchos predictores. Me alabé como un bayesiano recién convertido, hasta que me di cuenta de que no había tenido en cuenta las medidas repetidas dentro de los pacientes.

Entiendo que tengo que dar cuenta de eso. Entiendo que esto puede significar ajustar un hiperprior para cada individuo, ¿es así? ¿Qué forma tendrá esto en ERRORES?

Aquí hay un modelo básico de registro: (felicitaciones 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 )
  }
}

Sin embargo, no hay hiperprior aquí para el individuo. Aquí está mi mejor intento hasta ahora en un diseño individual, que tenga en cuenta las medidas repetidas dentro de las personas:

Aquí está el 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)

}

Aquí está mi modelo bastardo-niño 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)
  }

¿Alguien puede decirme si me dirijo en la dirección correcta? Mi comprensión de esto está creciendo, pero lentamente. Por favor se gentil. ¡Soy un médico, no una estadística! He usado R bastante, pero soy nuevo en BUGS y nuevo en Bayes.

Gracias,

R

Estás (estabas) casi allí. Solo unos pocos comentarios: no es necesario que el previo para los beta[,1:2]parámetros sea una MV conjunta normal; puede hacer que el anterior, tal que beta[i,1]y beta[i,2]es independiente, lo que simplifica las cosas (por ejemplo, hay covarianza antes es necesario especificar.) Tenga en cuenta que al hacerlo no significa que van a ser independiente en la parte posterior.

Otros comentarios: Dado que tiene un término constante alpha- en la regresión, los componentes beta[,1]deben tener una media cero en el anterior. Además, no tiene un previo alphaen el código.

Aquí hay un modelo con intersección jerárquica y términos de pendiente; He tratado de mantener sus antecedentes y notación cuando sea posible, dados los cambios:

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)
}

Un recurso muy útil para los modelos jerárquicos, incluidos algunos "trucos" para acelerar la convergencia, es Gelman y Hill .

(Un poco tarde con la respuesta, pero puede ser útil para algún futuro interrogador).