Visualisation des résultats de modèles mixtes

15

L'un des problèmes que j'ai toujours eu avec les modèles mixtes est de trouver des visualisations de données - du type qui pourraient se retrouver sur un papier ou une affiche - une fois que l'on a les résultats.

En ce moment, je travaille sur un modèle d'effets mixtes de Poisson avec une formule qui ressemble à ceci:

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

Avec quelque chose installé dans glm (), on pourrait facilement utiliser la fonction Predict () pour obtenir des prédictions pour un nouvel ensemble de données, et construire quelque chose à partir de cela. Mais avec une sortie comme celle-ci - comment pourriez-vous construire quelque chose comme un graphique du taux au fil du temps avec les décalages de X (et probablement avec une valeur définie de Y)? Je pense que l'on pourrait prévoir l'ajustement assez bien uniquement à partir des estimations des effets fixes, mais qu'en est-il de l'IC à 95%?

Y a-t-il autre chose que quelqu'un puisse penser qui puisse aider à visualiser les résultats? Les résultats du modèle sont ci-dessous:

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
la source
1
(+1) @EpiGrad: Pourquoi vous inquiétez-vous de l'IC (c'est-à-dire de l'erreur standard) des prédictions de la partie à effet fixe de votre modèle?
Boscovich
1
@andrea Une réponse intellectuelle et une réponse pratique: intellectuellement, je préfère généralement quantifier et visualiser l'incertitude quand je le peux. Pratiquement, parce que je suis presque sûr qu'un critique le demandera.
Fomite
Ouais ouais, bien sûr, mais je voulais dire quelque chose de différent. Mon commentaire n'était pas assez clair, désolé. Vous écrivez dans votre question "mais qu'en est-il de l'IC à 95%?". Mon commentaire est: pourquoi ne calculez-vous pas l'erreur-type de la prédiction à partir de la partie à effet fixe du modèle? Si vous êtes en mesure de calculer les valeurs prédites à partir de la partie à effet fixe, vous pouvez également calculer le SE et, donc, le CI. @EpiGrad
boscovich
@andrea Ah. Le problème est que l'une des choses que j'aimerais prédire, le temps, a également un effet aléatoire, dont je ne sais pas quoi faire.
Fomite
Eh bien, vous voulez prédire counts, nontime . Vous fixer des valeurs de X, Yet , timeet en utilisant la partie à effets fixes de votre modèle , vous prédisez counts. Il est vrai que timevotre modèle est également inclus en tant qu'effet aléatoire (tout comme l'interception et Y), mais cela n'a pas d'importance ici car utiliser uniquement la partie à effet fixe de votre modèle pour la prédiction, c'est comme mettre les effets aléatoires à 0 @EpiGrad
boscovich

Réponses:

4

Prédire countsà l'aide de la partie effets fixes de votre modèle signifie que vous mettez à zéro (c'est-à-dire leur moyenne) les effets aléatoires. Cela signifie que vous pouvez les "oublier" et utiliser des machines standard pour calculer les prédictions et les erreurs standard des prédictions (avec lesquelles vous pouvez calculer les intervalles de confiance).

Ceci est un exemple utilisant Stata, mais je suppose qu'il peut être facilement "traduit" en langage 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")

Le graphique fait référence treat == 0et est destiné à être un exemple ( visitn'est pas une variable vraiment continue, mais c'est juste pour avoir l'idée). Les lignes en pointillés représentent des intervalles de confiance à 95%.

entrez la description de l'image ici

Boscovich
la source