8 votos

Valores desconcertantes predichos en el modelo multinivel generalizado

Me encontré con un multinivel modelo con una variable dicotómica como el resultado de usar la logística función de enlace. Hice esto con la lme4 paquete en R. Una aproximación de mis datos se pueden encontrar en este GitHub Esencial: https://gist.github.com/markhwhiteii/3a954c4da82efdad5d8abee601b8a6a8

El diseño del proceso de recolección de datos es bastante simple: los Participantes son asignados al azar a un tratamiento (t) o control (c) condición. Luego se hicieron tres preguntas; pueden responder a uno, dos o tres de ellos. Para cada una de las preguntas (q1, q2, q3), pueden responder de manera positiva (1) o negativamente (0). No necesariamente la atención sobre el efecto de la variable-estoy considerando cada pregunta como indicadores de más o menos la misma construcción.

La pregunta es simple: ¿el tratamiento es aumentar la probabilidad de responder positivamente (1)? Es importante destacar que, el tamaño del efecto que quiero usar es un modelo implícito punto porcentual de diferencia entre los dos grupos. Esto requiere de la transformación de los valores predichos de logits a la probabilidad.

El modelo que especifica una intersección-sólo en el modelo, con la respuesta en el Nivel 1 (denotado por $i$) y la persona en el Nivel 2 (donado por $j$):

$\text{logit}(\pi_{ij}) = \beta_{0j}$ $\beta_{0j} = \gamma_{00} + \gamma_{01}X_j + u_{0j}$

En la sintaxis de lme4, esto es:

mod_1 <- glmer(y ~ x + (1 | id), dat, binomial)

(Si read.csv() en el .archivo csv en la Esencia enlace de arriba, este código va a estimar el modelo.)

Aquí es donde la confusión aumenta. Puedo transformar los valores de predicción para las respuestas en la condición de control y aquellos en la condición de tratamiento en la probabilidad de espacio. Que me dan:

inv_logit <- function(x) exp(x) / (1 + exp(x))
(control <- inv_logit(fixef(mod_1)[[1]]))
(treatment <- inv_logit(sum(fixef(mod_1))))
treatment - control

Esta se muestran sobre un 0.0003 probabilidad de responder positivamente en la condición de control, así como la condición de tratamiento, o un .03% de probabilidad de ser positivo. La diferencia entre las condiciones es bastante pequeña.

Estas probabilidades condicionales en cada condición parece demasiado pequeño, sobre todo cuando nos fijamos en el ingenuo porcentajes de respuestas positivas en ambas condiciones:

with(dat, tapply(y, x, mean))

Esta se muestran que alrededor de un 22% de las respuestas en la condición de control son positivos, mientras que el 23% en la condición de tratamiento son positivos. Mi pregunta fundamental es: ¿Cómo es que hay una gran diferencia entre el 22% de los resultados positivos empíricamente y sólo .03% implícitas en el modelo?

Podemos conseguir lo mismo si lo primero que promedio por persona (a no permitir que aquellas que responden más a tener más peso) y, a continuación, por la condición:

library(dplyr)
dat %>% 
  group_by(id) %>% 
  mutate(p = mean(y)) %>% 
  slice(1) %>% 
  group_by(x) %>% 
  summarise(p = mean(p))

Tanto en condiciones de demostrar el 24% de las respuestas positivas. De nuevo, este es drásticamente diferente que el .03% implícitas mod_1.

Alguna otra información que puede ser útil:

  • 604 la gente tenía 1 respuesta, 301 había 2 respuestas, y 95 tenía 3 respuestas. Creo que parte del problema es que la mayoría de la gente sólo había una observación? Como lo que yo puedo decir, sin embargo, no necesitamos estimar una diferente distribución normal de las respuestas para cada grupo (en este caso, una persona), así que no entiendo por qué sólo 1 respuesta iba a ser que me da hang-ups.

  • La distribución de azar intercepta fue decididamente no es normal , como puede verse en hist(ranef(mod_1)$id[[1]], breaks = 50):

enter image description here

El modelo que claramente parece mal especificada a partir de ese terrible distribución de azar intercepta, pero me gustaría saber:

  1. ¿Por qué mi ingenuo mira que cerca de 22 a 24% de las respuestas positivas, difieren tanto de la modelo-implícita .03%? Estoy en busca de una respuesta más allá de "el modelo está mal especificada," porque soy curioso en cuanto a por qué esto está ocurriendo, además de que no me acerque a la especificación de un modelo adecuado. Lo que me lleva a...

  2. ¿Cómo puedo especificar un modelo adecuado que responda a mi pregunta sobre la condición experimental? Supongo que también podrían tratar de utilizar algún tipo de sándwich de covarianza de la estimación que maneja los que dependen de las respuestas, tales como sandwich::vcovCL()? Ese enfoque no cambia el SEs demasiado en la circunstancia actual. Esto también no satisfacer mi curiosidad de por qué el modelo mixto es que me extrañas predicciones.

4voto

Bill Puntos 16

La respuesta podría al menos parcialmente oculto en su respuesta:

La predicción de los efectos aleatorios son todos no negativos. La adición de ellos a la lineal predictor de los efectos fijos a empujar hacia arriba el promedio de la predicción para el nivel deseado.

Así que lo que podría ser la razón por la que no centrada en la distribución de la eBLUPs? Creo que está oculto en el extremadamente desequilibrada datos de situación: 604 de los 1000 identificadores de proporcionar una sola medición, lo que puede llevar fácilmente a numérico conflictos entre aleatorios y de efectos fijos.

En tales situaciones, que a menudo corren como un análisis de sensibilidad normal de modelo mixto. Si hacemos esto con sus datos, obtenemos:

mod_2 <- lmer(y ~ x + (1 | id), dat)
summary(mod_2)

# Output
Random effects:
 Groups   Name        Variance Std.Dev.
 id       (Intercept) 0.06961  0.2638  
 Residual             0.10825  0.3290  
Number of obs: 1491, groups:  id, 1000

Fixed effects:
            Estimate Std. Error t value
(Intercept) 0.228315   0.022716  10.051
xt          0.007224   0.026922   0.268

Correlation of Fixed Effects:
   (Intr)
xt -0.844

Los efectos fijos ahora mira como se esperaba por el análisis descriptivo y muy diferente de lo que el GLMER encontrado. Ahora, el eBLUPs son incluso centrado (aunque, por supuesto, todavía lejos de la normal):

summary(unlist(ranef(mod_2, drop = TRUE)))

# Output
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.15513 -0.12844 -0.09218  0.00000  0.14878  0.50347 

Dependiendo de su exacta de las preguntas de investigación y objetivos del análisis, se podría tratar de generalizar la estimación de ecuaciones de GEE, en R:

library(gee)
fit_gee <- gee(y ~ x, 
               id = id,         
               data = dat %>% arrange(id), 
               family = binomial, 
               corstr = "exchangeable")
summary(fit_gee)

# Output
Coefficients:
               Estimate Naive S.E.    Naive z Robust S.E.   Robust z
(Intercept) -1.22561668  0.1257257 -9.7483411   0.1235148 -9.9228358
xt           0.04325299  0.1484306  0.2914022   0.1479156  0.2924167

Estimated Scale Parameter:  0.9859774
Number of Iterations:  3

Working Correlation
          [,1]      [,2]      [,3]
[1,] 1.0000000 0.3133096 0.3133096
[2,] 0.3133096 1.0000000 0.3133096
[3,] 0.3133096 0.3133096 1.0000000

# Distribution of predictions
summary(fitted(fit_gee))

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
 0.2269  0.2269  0.2346  0.2324  0.2346  0.2346 

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X