Visualizar resultados de modelos mixtos

15

Uno de los problemas que siempre he tenido con los modelos mixtos es descubrir visualizaciones de datos, del tipo que podría terminar en un papel o póster, una vez que uno tiene los resultados.

En este momento, estoy trabajando en un modelo de efectos mixtos de Poisson con una fórmula similar a la siguiente:

a <- glmer(counts ~ X + Y + Time + (Y + Time | Site) + offset(log(people))

Con algo ajustado en glm (), uno podría usar fácilmente el predic () para obtener predicciones para un nuevo conjunto de datos, y construir algo a partir de eso. Pero con un resultado como este: ¿cómo construiría algo así como una gráfica de la tasa a lo largo del tiempo con los cambios desde X (y probablemente con un valor establecido de Y)? Creo que uno podría predecir el ajuste lo suficientemente bien solo a partir de las estimaciones de efectos fijos, pero ¿qué pasa con el IC del 95%?

¿Hay algo más que alguien pueda pensar que ayude a visualizar los resultados? Los resultados del modelo son los siguientes:

Random effects:
 Groups     Name        Variance   Std.Dev.  Corr          
 Site       (Intercept) 5.3678e-01 0.7326513               
            time        2.4173e-05 0.0049167  0.250        
            Y           4.9378e-05 0.0070270 -0.911  0.172 

Fixed effects:
              Estimate Std. Error z value Pr(>|z|)    
(Intercept) -8.1679391  0.1479849  -55.19  < 2e-16
X            0.4130639  0.1013899    4.07 4.62e-05
time         0.0009053  0.0012980    0.70    0.486    
Y            0.0187977  0.0023531    7.99 1.37e-15

Correlation of Fixed Effects:
         (Intr) Y    time
X      -0.178              
time    0.387 -0.305       
Y      -0.589  0.009  0.085
r data-visualization mixed-model Fomite
fuente
1
(+1) @EpiGrad: ¿Por qué le preocupa el CI (es decir, el error estándar) de las predicciones de la parte de efectos fijos de su modelo?
boscovich
1
@andrea Una respuesta intelectual y práctica: intelectualmente, en general, prefiero cuantificar y visualizar la incertidumbre cuando puedo. Prácticamente, porque estoy bastante seguro de que un revisor lo pedirá.
Fomite
Sí, sí, claro, pero quise decir algo diferente. Mi comentario no fue lo suficientemente claro, lo siento. Escribe en tu pregunta "pero ¿qué pasa con el IC del 95%?". Mi comentario es: ¿por qué no calculan el error estándar de la predicción a partir de la parte de efectos fijos del modelo? Si puede calcular los valores pronosticados de la parte de efecto fijo, también puede calcular el SE y, por lo tanto, el CI. @EpiGrad
boscovich
@andrea Ah. La preocupación es que una de las cosas que me gustaría predecir, el tiempo, también tiene un efecto aleatorio, que no tengo idea de qué hacer.
Fomite
Bueno, quieres predecir counts, no time. Usted fija los valores de X, Yy , timeutilizando la parte de efectos fijos de su modelo que predice counts. Es cierto que timeestá incluido en su modelo también como un efecto aleatorio (al igual que la intercepción y Y), pero no importa aquí porque usar solo la parte de efecto fijo de su modelo para la predicción es como establecer los efectos aleatorios en 0 @EpiGrad
boscovich

Respuestas:

4

Predecir countsusando la parte de efectos fijos de su modelo significa que establece en cero (es decir, su media) los efectos aleatorios. Esto significa que puede "olvidarse" de ellos y utilizar maquinaria estándar para calcular las predicciones y los errores estándar de las predicciones (con los que puede calcular los intervalos de confianza).

Este es un ejemplo usando Stata, pero supongo que se puede "traducir" fácilmente al lenguaje R:

webuse epilepsy, clear
xtmepoisson seizures treat visit || subject: visit
predict log_seiz, xb
gen pred_seiz = exp(log_seiz)
predict std_log_seiz, stdp
gen ub = exp(log_seiz+invnorm(.975)*std_log_seiz)
gen lb = exp(log_seiz-invnorm(.975)*std_log_seiz)

tw (line pred_seiz ub lb visit if treat == 0, sort lc(black black black) ///
 lp(l - -)), scheme(s1mono) legend(off) ytitle("Predicted Seizures") ///
 xtitle("Visit")

El gráfico se refiere treat == 0y pretende ser un ejemplo ( visitno es una variable realmente continua, pero es solo para tener una idea). Las líneas discontinuas son intervalos de confianza del 95%.

ingrese la descripción de la imagen aquí

boscovich
fuente