Ajuste de un modelo mixto Poisson GLM con pendiente aleatoria e intercepción

Actualmente estoy trabajando en una serie de modelos de series de tiempo de Poisson tratando de estimar el efecto de un cambio en la forma en que se obtuvieron los recuentos (cambio de una prueba de diagnóstico a otra) mientras se controlan otras tendencias a lo largo del tiempo (digamos un aumento general en el incidencia de la enfermedad). Tengo datos para varios sitios diferentes.

Si bien también he estado jugando con los GAM, he ajustado una serie de GLM bastante básicos con tendencias temporales en ellos, luego agrupando los resultados. El código para esto se vería así en SAS:

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

o esto en R:

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

Luego, tome esas estimaciones y agrúpelas en varios sitios. También se ha sugerido que intento usar un modelo mixto de Poisson con una pendiente aleatoria e intercepción para cada sitio, en lugar de agruparlos. Entonces, esencialmente tendría el efecto fijo de variable dependiente, luego un efecto aleatorio para la intercepción y el tiempo (o idealmente tiempo y tiempo ^ 2, aunque entiendo que se vuelve un poco peludo).

Mi problema es que no tengo idea de cómo ajustar uno de estos modelos, y parece que los modelos mixtos son donde la documentación de todos de repente se vuelve muy opaca. ¿Alguien tiene una explicación simple (o código) sobre cómo ajustar lo que estoy buscando y qué buscar?

En R:

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

$Y_i \sim {\rm Poisson}(\lambda_i)$

$X_i$dependent_variable$t$time$i$ID$\beta_0, \beta_1$$\eta_{i1}, \eta_{i2}$ son efectos aleatorios cuyas variaciones son estimadas por el modelo.

${\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

Un punto de precaución: la Std.Dev.columna es solo la raíz cuadrada de la Variancecolumna, ¡ NO el error estándar de la estimación de la varianza!

En SAS:

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

Pero, por supuesto, hay muchas opciones, más o menos útiles, para jugar.