Dopasowanie mieszanego modelu Poissona GLM z losowym nachyleniem i punktem przecięcia

Obecnie pracuję nad serią modeli szeregów czasowych Poissona, próbując oszacować efekt zmiany sposobu zliczania (przejście z jednego testu diagnostycznego do drugiego), jednocześnie kontrolując inne trendy w czasie (powiedzmy ogólny wzrost w występowanie choroby). Mam dane dla wielu różnych stron.

Chociaż majstrowałem również przy GAM, dopasowałem serię dość podstawowych GLM z trendami czasowymi, a następnie zsumowałem wyniki. Kod tego mógłby wyglądać mniej więcej tak w SAS:

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

lub to w R:

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

Następnie biorąc te szacunki i łącząc je w różnych witrynach. Sugeruje się również, że próbuję używać mieszanego modelu Poissona z losowym nachyleniem i przechwytywaniem dla każdego miejsca, zamiast łączenia. Zasadniczo miałbyś więc stały efekt zmiennej zależnej, a następnie losowy efekt przechwytywania i czasu (lub idealnie czasu i czasu ^ 2, choć rozumiem, że robi się trochę owłosiony).

Mój problem polega na tym, że nie mam pojęcia, jak dopasować jeden z tych modeli, i wydaje się, że modele mieszane są tam, gdzie dokumentacja każdego staje się nagle bardzo nieprzejrzysta. Czy ktoś ma proste wyjaśnienie (lub kod), jak dopasować to, co chcę dopasować i na co uważać?

— Fomite
źródło

Odpowiedzi:

W R:

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

W tym przypadku i ten kod pasuje do modelu $Y_i \sim {\rm Poisson}(\lambda_i)$

\log (λ_{i}) = β_{0} + β_{1} X_{i} + η_{i 1} + η_{i 2} t

$\log(\lambda_i) = \beta_0 + \beta_1 X_i + \eta_{i1} + \eta_{i2}t$

gdzie jest , jest, a jest . są efektami stałymi, a są efektami losowymi, których wariancje są szacowane przez model. $X_i$ dependent_variable $t$ time $i$ ID $\beta_0, \beta_1$ $\eta_{i1}, \eta_{i2}$

Oto przykład z niektórymi szybko symulowanymi danymi, w których wariancje efektu losowego naprawdę wynoszą 0, zmienna towarzysząca nie ma żadnego efektu, każdy wynik to , a każda osoba jest widziana 10 razy w czasie . ${\rm 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

Jeden punkt ostrożności - Std.Dev.kolumna jest tylko pierwiastkiem kwadratowym Variancekolumny, a NIE standardowym błędem oszacowania wariancji!

— Makro
źródło

A jego ηi1 powoduje losowe przechwytywanie?

— Fomite

η_{i 1}

$\eta_{i1}$ jest losowym przechwytywaniem, tak.

— Makro

Dzięki za odpowiedź - jeszcze jedno pytanie. W niektórych GLM niektóre strony korzystają z terminów ^ 2. Wydaje się, że większość ludzi zaznacza jeden lub dwa losowe efekty - jak paskudny będzie trzeci pod względem obliczeń?

— Fomite

Cóż, w powyższym symulowanym przykładzie, który miał tylko 100 obserwacji i 10 grup, dodanie trzeciego losowego efektu (przez wpisanie g <- lmer(y ~ x + (1+t+I(t^2)|ID), family="poisson")) zwiększyło czas obliczeń z około 0,75 sekundy do około 11 sekund. Wraz ze wzrostem wielkości próbki prawdopodobnie zwiększa się również czas obliczeń.

— Makro

@ andrea, co odzwierciedlałoby przekonanie, że w zbiorze danych istnieje świecki trend - nie zakładałem z góry, że ma to sens. Gdyby jednostkami były osoby wiek, to z pewnością zgadzam się, że ustalony efekt w czasie miał wiele sensu, ale w innych sytuacjach nachylenie w czasie byłoby bardziej losowym kierunkiem, w którym zmierza każda jednostka. Dlatego nie uwzględniłem tego efektu (a epigrad nie zapytał, jak uwzględnić efekt ustalonego czasu)

t

$t$

— Makro

W SAS:

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

Ale oczywiście jest wiele opcji, mniej lub bardziej użytecznych, do zabawy.

— Boskowicz
źródło

Dzięki :) Niestety, wydaje mi się, że wpadłem na ziemię w kwestiach konwergencji, ale majstruję przy nich.

— Fomite