Подгонка смешанной модели Пуассона GLM со случайным наклоном и пересечением

В настоящее время я работаю над серией моделей временных рядов Пуассона, пытаясь оценить влияние изменения в том, как были получены подсчеты (переключение с одного диагностического теста на другой), в то же время контролируя другие тренды во времени (скажем, общее увеличение заболеваемость). У меня есть данные для разных сайтов.

Хотя я также работал с GAM, я установил ряд довольно простых GLM с временными тенденциями в них, а затем объединил результаты. Код для этого будет выглядеть примерно так в SAS:

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

или это в R:

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

Затем беру эти оценки и объединяю их по различным сайтам. Предполагается также, что я пытаюсь использовать смешанную модель Пуассона со случайным наклоном и перехватом для каждого участка, а не объединять. Так что, по сути, у вас будет фиксированный эффект variable_variable, а затем случайный эффект для перехвата и времени (или, в идеале, времени и времени ^ 2, хотя я понимаю, что это немного затруднительно).

Моя проблема в том, что я понятия не имею, как подойти к одной из этих моделей, и кажется, что смешанные модели, где документация для всех, вдруг оказывается очень непрозрачной. У кого-нибудь есть простое объяснение (или код) о том, как соответствовать тому, что я ищу, и на что обратить внимание?

В R:

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

В этом случае и этот код соответствует модели $Y_i \sim {\rm Poisson}(\lambda_i)$

где есть , есть и есть . - фиксированные эффекты, а - случайные эффекты, дисперсии которых оцениваются моделью. t i β 0 , β 1 η i 1 , η i $X_i$dependent_variable$t$time$i$ID$\beta_0, \beta_1$$\eta_{i1}, \eta_{i2}$

Вот пример с некоторыми быстро смоделированными данными, где случайные отклонения эффекта действительно равны 0, ковариата не имеет эффекта, каждый результат равен , и каждый человек виден 10 раз в моменты времени .т = 1 , . , , , ${\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

Одно предостережение: Std.Dev.столбец - это просто квадратный корень Varianceстолбца, а НЕ стандартная ошибка оценки дисперсии!

В САС:

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

Но, конечно, есть множество вариантов, более или менее полезных, с которыми можно поиграть.