¿Qué son la estructura R y la estructura G en un glmm?

16

He estado usando el MCMCglmmpaquete recientemente. Estoy confundido por lo que se refiere en la documentación como estructura R y estructura G. Estos parecen relacionarse con los efectos aleatorios, en particular especificando los parámetros para la distribución previa sobre ellos, pero la discusión en la documentación parece suponer que el lector sabe cuáles son estos términos. Por ejemplo:

lista opcional de especificaciones previas que tienen 3 elementos posibles: R (estructura R) G (estructura G) y B (efectos fijos) ............ Los antecedentes de las estructuras de varianza (R y G ) son listas con las (co) varianzas esperadas (V) y el parámetro de grado de creencia (nu) para el inverso-Wishart

... tomado de aquí .

EDITAR: Tenga en cuenta que he reescrito el resto de la pregunta después de los comentarios de Stephane.

¿Alguien puede arrojar luz sobre lo que son la estructura R y la estructura G, en el contexto de un modelo de componentes de varianza simple donde el predictor lineal es con y

β0+e0ij+u0j
e0ijN(0,σ0e2)u0jN(0,σ0u2)

Hice el siguiente ejemplo con algunos datos que vienen con MCMCglmm

> require(MCMCglmm)
> require(lme4)
> data(PlodiaRB)
> prior1 = list(R = list(V = 1, fix=1), G = list(G1 = list(V = 1, nu = 0.002)))
> m1 <- MCMCglmm(Pupated ~1, random = ~FSfamily, family = "categorical", 
+ data = PlodiaRB, prior = prior1, verbose = FALSE)
> summary(m1)


 G-structure:  ~FSfamily

         post.mean l-95% CI u-95% CI eff.samp
FSfamily    0.8529   0.2951    1.455      160

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units         1        1        1        0

 Location effects: Pupated ~ 1 

            post.mean l-95% CI u-95% CI eff.samp  pMCMC    
(Intercept)   -1.1630  -1.4558  -0.8119    463.1 <0.001 ***
---

> prior2 = list(R = list(V = 1, nu = 0), G = list(G1 = list(V = 1, nu = 0.002)))
> m2 <- MCMCglmm(Pupated ~1, random = ~FSfamily, family = "categorical", 
+ data = PlodiaRB, prior = prior2, verbose = FALSE)
> summary(m2)


 G-structure:  ~FSfamily

         post.mean l-95% CI u-95% CI eff.samp
FSfamily    0.8325   0.3101    1.438    79.25

 R-structure:  ~units

      post.mean l-95% CI u-95% CI eff.samp
units    0.7212  0.04808    2.427    3.125

 Location effects: Pupated ~ 1 

            post.mean l-95% CI u-95% CI eff.samp  pMCMC    
(Intercept)   -1.1042  -1.5191  -0.7078    20.99 <0.001 ***
---

> m2 <- glmer(Pupated ~ 1+ (1|FSfamily), family="binomial",data=PlodiaRB)
> summary(m2)
Generalized linear mixed model fit by the Laplace approximation 
Formula: Pupated ~ 1 + (1 | FSfamily) 
   Data: PlodiaRB 
  AIC  BIC logLik deviance
 1020 1029   -508     1016
Random effects:
 Groups   Name        Variance Std.Dev.
 FSfamily (Intercept) 0.56023  0.74849 
Number of obs: 874, groups: FSfamily, 49

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.9861     0.1344  -7.336  2.2e-13 ***

Entonces, según los comentarios de Stephane, creo que la estructura G es para . Pero los comentarios también dicen que la estructura R es para pero esto no parece aparecer en la salida.σ0u2σ0e2lme4

Tenga en cuenta que los resultados lme4/glmer()son consistentes con ambos ejemplos de MCMC MCMCglmm.

Entonces, ¿es la estructura R para y por qué no aparece esto en la salida para ?σ0e2lme4/glmer()

Joe King
fuente
1
Con la terminología SAS (pero posiblemente es una terminología más común), la matriz G es la matriz de varianza de los efectos aleatorios y la matriz R es la matriz de varianza de los "términos de error" (en su caso, tal vez sea el residual estimado varianza )?σ0e2
Stéphane Laurent
@ StéphaneLaurent gracias. Me preguntaba si podría estimarse pero cuando supe por primera vez sobre el modelo lineal generalizado, recuerdo que no se estima, solo se calcula la "desviación" (como con ). Tal vez me estoy perdiendo algo? σ0e2σ0e2lme4
Joe King
1
quizás el sentido de la varianza residual no está claro cuando la familia de distribución no es la gaussiana
Stéphane Laurent
1
@ Stéphane Laurent ¡Sí! Por favor, vea mi comentario a la respuesta de Michael hace un minuto - para el resultado binario, debe ser reparado (como en mis modelos en mi OP)
Joe King
1
Cuando tiene un modelo ME / Multinivel, hay varias variaciones. Imagine el caso más simple: . Hay variación en las intersecciones , y en el término de error . se usa a menudo para la matriz var-covar de los efectos aleatorios (en este caso, un escalar, ) y es para la matriz var-covar de las varianzas residuales después de tener en cuenta los efectos aleatorios fijos y de ese grupo efectos Generalmente se concibe como una matriz diagonal de 's. Además, ambos discos se consideran multivariados normales con media = 0. b i ε i G σ 2 b R i ε i σ 2Yi=β0+β1X+bi+εibiεiGσb2Riεiσ2
gung - Restablece a Monica

Respuestas:

8

Preferiría publicar mis comentarios a continuación como un comentario, pero esto no sería suficiente. Estas son preguntas en lugar de una respuesta (de manera similar a @gung, no me siento lo suficientemente fuerte sobre el tema).

Tengo la impresión de que MCMCglmm no implementa un "mmm" bayesiano "verdadero". El verdadero modelo bayesiano se describe en la sección 2 de este documento . De manera similar al modelo frecuentista, uno tiene y hay un requisito previo en el parámetro de dispersión ϕ 1 además de los parámetros fijos β y la varianza "G" de la efecto aleatorio u .g(E(yu))=Xβ+Zuϕ1βu

Pero de acuerdo con esta viñeta MCMCglmm , el modelo implementado en MCMCglmm viene dado por , y no involucra el parámetro de dispersión ϕ 1 . No es similar al modelo frecuentista clásico.g(E(yu,e))=Xβ+Zu+eϕ1

Por lo tanto, no me sorprendería que no haya un análogo de con glmer.σe

Por favor, discúlpate por estos comentarios, solo eché un vistazo rápido al respecto.

Stéphane Laurent
fuente
Gracias. ¿Se supone que este tema es difícil porque lo encuentro bastante difícil? Creo que estoy satisfecho con el significado de la estructura R y G ahora. Todavía estoy confundido por la falta de con y soy muy curioso acerca de su comentario de que no es verdaderamente bayesiano. Honestamente, no puedo decir que entiendo todo el documento al que se vinculó y también estoy luchando con partes de la viñeta, pero solo desde la perspectiva de mi ejemplo, creo que el parámetro de dispersión ϕ 1 debería ser constante (porque el ejemplo es binomial). Qué me estoy perdiendo ? σeglmerMCMCglmmMCMCglmmϕ1
Joe King
Lo siento, mis palabras no fueron totalmente apropiadas. MCMCglmm es verdaderamente bayesiano, pero no implementa exactamente el glmm clásico (creo). Además, debe tener en cuenta que es difícil establecer prioridades que generen una inferencia sobre los componentes de la varianza cercanos a la inferencia frecuentista.
Stéphane Laurent
Gracias de nuevo. En mis estudios, descubrí que puedo usar la distribución inversa inversa de Wishart para los componentes de varianza al MCMCglmmusar una variedad de parámetros, y los intervalos creíbles del 95% siempre contienen el valor de varianza para la estimación de efectos aleatorios por glmerlo que sentí que esto era razonable , pero ¿cómo debo interpretar este caso, que puede no ser típico, donde el resultado es que los MCMCglmmintervalos no son muy sensibles a la elección de antes? ¿Tal vez debería hacer una nueva pregunta sobre esto?
Joe King
¿Quizás tiene un gran tamaño de muestra? Con respecto a su pregunta inicial, tengo la impresión de que, al menos para el caso binomial, el modelo glmer es equivalente al modelo MCMCglmm con . ¿Qué sucede si establece un previo en σ e altamente concentrado en 0 ? σe=0σe0
Stéphane Laurent
σeσeσuglmer
11

Rσe21

GG

Una nota final, debido a que la varianza residual no está fijada en cero, las estimaciones no coincidirán con las de glmer. Necesitas reescalarlos. Aquí hay un pequeño ejemplo (no usa efectos aleatorios, pero generaliza). Observe cómo la varianza de la estructura R se fija en 1.

# example showing how close the match is to ML without separation
m2 <- MCMCglmm(vs ~ mpg, data = mtcars, family = "categorical",
  prior = list(
    B = list(mu = c(0, 0), V = diag(2) * 1e10),
    R = list(V = 1, fix = 1)),
  nitt = 1e6, thin = 500, burnin = 10000)
summary(m2)

Aquí está la constante de reescalado para la familia binomial:

k <- ((16*sqrt(3))/(15*pi))^2

Ahora divide la solución entre ella y obtén los modos posteriores

posterior.mode(m2$Sol/(sqrt(1 + k)))

Que debería estar bastante cerca de lo que obtenemos de glm

summary(glm(vs ~mpg, data = mtcars, family = binomial))
Joshua
fuente
¿Sabrías cómo especificar la heterocedasticidad en el nivel uno en MCMCglmm? ¿Es esa la estructura R? ¿Cuál es la sintaxis entonces?
Maxim.K
@Joshua, ¿puedes explicar la "constante de reescalado para la familia binomial"? PD: Para semilla 123, obtengo (con la corrección) de m2los valores -8.164y 0.421; y de glmlos valores -8.833y 0.430.
Qaswed
La constante de reescalado se puede encontrar en Diggle et. Alabama. ( amazon.de/Analysis-Longitudinal-Oxford-Statistical-Science/dp/… ) - de acuerdo con cran.r-project.org/web/packages/MCMCglmm/vignettes/… eq. 2.14 en la página 47.
Qaswed