Ajustement d'un modèle mixte de Poisson GLM avec une pente et une interception aléatoires

9

Je travaille actuellement sur une série de modèles de séries chronologiques de Poisson essayant d'estimer l'effet d'un changement dans la façon dont les chiffres ont été obtenus (passage d'un test de diagnostic à un autre) tout en contrôlant d'autres tendances au fil du temps (par exemple, une augmentation générale de la incidence de la maladie). J'ai des données pour un certain nombre de sites différents.

Bien que j'aie aussi bricolé avec les GAM, j'ai adapté une série de GLM assez basiques avec les tendances temporelles, puis j'ai regroupé les résultats. Le code pour cela ressemblerait à ceci dans SAS:

PROC GENMOD data=work.data descending;
  model counts = dependent_variable time time*time / link=log dist = poisson;
run;

ou ceci dans R:

glm(counts ~ dependent_variable + time + time*time, family="poisson")

Prendre ensuite ces estimations et les regrouper sur les différents sites. Il a également été suggéré que j'essaie d'utiliser un modèle mixte de Poisson avec une pente aléatoire et une interception pour chaque site, plutôt que de mettre en commun. Donc, essentiellement, vous auriez l'effet fixe de variable_dépendante, puis un effet aléatoire pour l'interception et le temps (ou idéalement le temps et le temps ^ 2, même si je comprends que cela devient un peu velu).

Mon problème est que je n'ai aucune idée de comment adapter l'un de ces modèles, et il semble que les modèles mixtes sont où la documentation de tout le monde devient soudainement très opaque. Quelqu'un a-t-il une explication (ou un code) simple sur la façon de s'adapter à ce que je cherche à adapter, et à quoi faire attention?

Fomite
la source

Réponses:

14

Dans R:

library(lme4)
lmer(counts ~ dependent_variable + (1+time|ID), family="poisson")

YiPoisson(λi)

log(λi)=β0+β1Xi+ηi1+ηi2t

Xidependent_variablettimeiIDβ0,β1ηi1,ηi2

Poisson(1)t=1,...,10

x = rnorm(100)
t = rep(1:10,each=10)
ID = rep(1:10,10)
y = rpois(100,1)
g <- lmer(y ~ x + (1+t|ID), family="poisson")
summary(g)
Generalized linear mixed model fit by the Laplace approximation 
Formula: y ~ x + (1 + t | ID) 
   AIC   BIC logLik deviance
 108.8 121.9 -49.42    98.85
Random effects:
 Groups Name        Variance  Std.Dev. Corr   
 ID     (Intercept) 0.0285038 0.168831        
        t           0.0027741 0.052669 -1.000 
Number of obs: 100, groups: ID, 10

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.09078    0.11808  -0.769    0.442
x            0.13670    0.08845   1.546    0.122

Correlation of Fixed Effects:
  (Intr)
x -0.127

Un point de prudence - la Std.Dev.colonne est juste la racine carrée de la Variancecolonne, PAS l'erreur standard de l'estimation de la variance!

Macro
la source
Et son ηi1 qui se traduit par l'ordonnée à l'origine aléatoire?
Fomite
ηi1
Merci pour la réponse - encore une question. Dans certains GLM, certains sites bénéficient grandement à ^ 2 terme. Il semble que la plupart des gens identifient un ou deux effets aléatoires - à quel point un troisième sera-t-il méchant en termes de calcul?
Fomite
Eh bien, dans l'exemple simulé ci-dessus, qui n'avait que 100 observations et 10 groupes, en ajoutant le troisième effet aléatoire (en tapant g <- lmer(y ~ x + (1+t+I(t^2)|ID), family="poisson")), le temps de calcul est passé d'environ 0,75 seconde à environ 11 secondes. À mesure que la taille de l'échantillon augmente, l'augmentation du temps de calcul augmente probablement également.
Macro
1
t
2

En SAS:

proc glimmix data = yourdata ic = q;
    class id;
    model y = x / dist = poisson solution;
    random intercept t / subject = id;
run;

Mais bien sûr, il existe de nombreuses options, plus ou moins utiles, pour jouer avec.

Boscovich
la source
Merci :) Malheureusement, je semble avoir échoué sur les problèmes de convergence, mais je vais bricoler avec ceux-ci.
Fomite