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)
:
El modelo que claramente parece mal especificada a partir de ese terrible distribución de azar intercepta, pero me gustaría saber:
¿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...
¿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.