Cómo calcular el estimado de las proporciones y sus intervalos de confianza a partir de un modelo mixto?

Question

Tengo un experimento con dos tratamientos. Es una división de la parcela de experimentación, la estructura de Bloque/Treatment1/Treatment2. Cada tratamiento tiene 2 niveles. La variable dependiente es la presencia/ausencia de datos funcionales speceis grupos. Yo soy el análisis de los datos utilizando modelos mixtos a través de lme4 y afex.

La estructura del modelo es la siguiente:

m <- lmer (DV ~ Treatment 1 * Treatment 2 + (1|Block/Treatment1),
           family = binomial) 

Para la gráfica de los efectos sin embargo, quiero trazar el tratamiento significa para cada una de las 4 combinaciones de tratamiento (Tratamiento 1 nivel y Tratamiento 2 nivel a, el Tratamiento 1 Tratamiento B 2 nivel a, el Tratamiento 1 nivel B y el Tratamiento 2 nivel a, el Tratamiento 1 nivel B y el Tratamiento 2 nivel B). No podía encontrar una manera intuitiva para extraer estos media y el error estándar de los valores de R para los datos con una binomial distribución de error, por lo tanto yo soy el cálculo de estos valores con la mano. Aprecio que el estándar de los errores de calculo no tomar en cuenta el efecto aleatorio de la estructura en mis datos, pero no puedo encontrar una manera de calcular de ellos para incluir este efecto aleatorio de la estructura (mediante una simulación mcmc en lenguaje R para crear HPD intervalos de confianza sólo funciona para Gaussiano de datos).

Como tal, quiero calcular tratamiento significa para la proporción de un determinado grupo funcional en una particular combinación de tratamiento y, a continuación, los errores estándar de estos 4 tratamiento de los medios.

Considere la posibilidad de que las proporciones que tengo para un tratamiento en particular son:

Plot 1: 24/65
Plot 2 26/64
Plot 3: 25/65
Plot 4: 22/62
Plot 5: 30/66
Plot 6: 29/65

Entiendo que para calcular el promedio de la proporción necesito suma el número total de aciertos (suma de todos los numeradores de arriba) y se divide por el total de la suma de aciertos y errores (suma de todos los denominadores de arriba). Esto me da 0.40.

Mi pregunta entonces es, ¿cómo hace uno para calcular el error estándar de la proporción resultante?

Llamar a la proporción de éxitos p y no p, tengo la siguiente fórmula para el error estándar de una proporción:

sqrt((p*q)/n)

Es esto correcto? Además, ¿qué es n? Sea n el número de proporciones que componen la media de la proporción (en este caso 6) o el número total de aciertos y errores (es decir, la suma de todos los denominadores de arriba, que aquí es de 387)?

Siento que esta es una simple pregunta, yo simplemente no puede encontrar una respuesta en cualquier lugar! También pido disculpas si este no es el foro correcto para este post...

Muchas gracias.

Answer 1

Primera nota: no se puede calcular un estándar decente de error en las probabilidades, usted tiene que hacerlo en una escala logit y las usan para construir los intervalos de confianza. Los intervalos alrededor de las probabilidades son casi simétricos, y definitivamente no cuando se utiliza un modelo mixto.

Usted puede fácilmente trama de los efectos mediante el paquete effects:

Con la función Effect() puede especificar los efectos que desea parcela y parcela de inmediato, o extracto de la información que desea.

Algunos al azar los datos falsos:

ndata <- data.frame(
  DV = sample(0:1,200,TRUE),
  Treatment1 = rep(rep(c('A1','B1'),25),4),
  Treatment2 = rep(rep(c('A2','B2'),each=25),4),
  Block = rep(c("Block1","Block2"),each=100)
  )

m <- glmer (DV ~ Treatment1 * Treatment2 + (1|Block/Treatment1),
           data=ndata, family = binomial) 

Para obtener una parcela de los efectos, usted puede simplemente hacer:

plot(allEffects(m))

Para obtener:

El mismo que se obtiene con plot(Effect(c("Treatment1","Treatment2"),m)

Si desea obtener los datos reales, puede guardar el resultado de una llamada a Effect() en un objeto, y extraer los datos necesarios:

est <- Effect(c("Treatment1","Treatment2"),m)
cbind(est$x,est$fit,est$se,est$lower,est$upper)

para obtener:

  Treatment1 Treatment2       est$fit    est$se  est$lower est$upper
1         A1         A2 -7.696104e-02 0.2775554 -0.6209597 0.4670376
2         B1         A2 -1.670541e-01 0.2896820 -0.7348204 0.4007122
3         A1         B2 -3.364722e-01 0.2927591 -0.9102694 0.2373250
4         B1         B2  1.110223e-16 0.2773501 -0.5435962 0.5435962

Tenga en cuenta que estas están en el original (logit) de escala. Calcular un intervalo de confianza implicaría la transformación de este a la escala original, utilizando, por ejemplo. plogis() :

> cbind(est$x,plogis(est$fit),plogis(est$lower),plogis(est$upper))
  Treatment1 Treatment2 plogis(est$fit) plogis(est$lower) plogis(est$upper)
1         A1         A2       0.4807692         0.3495632         0.6146824
2         B1         A2       0.4583333         0.3241378         0.5988588
3         A1         B2       0.4166667         0.2869447         0.5590543
4         B1         B2       0.5000000         0.3673514         0.6326486

PS : esto no es la más limpia de código, es sólo para fines ilustrativos.