Fitting a Poisson GLM mixed model with a random slope and intercept

Question

I'm currently working on a series of Poisson time series models trying to estimate the effect of a change in how the counts were obtained (switching from one diagnostic test to another) while controlling for other trends over time (say a general increase in the incidence of disease). I've got data for a number of different sites.

While I've been tinkering with GAMs as well, I've fit a series of pretty basic GLMs with time trends in them, then pooling the results. The code for this would look something like this in SAS:

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

or this in R:

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

Then taking those estimates, and pooling them over the various sites. It's also been suggested to be that I try using a Poisson mixed model with a random slope and intercept for each site, rather than pooling. So essentially you'd have the fixed effect of dependent_variable, then a random effect for the intercept and time (or ideally time and time^2 though I understand that gets a bit hairy).

My issue is I have no idea how to fit one of these models, and it seems that mixed models are where everyone's documentation goes suddenly very opaque. Anyone have a simple explanation (or code) on how to fit what I'm looking to fit, and what to look out for?

Answer 1

In R:

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

In this case $Y_i \sim {\rm Poisson}(\lambda_i)$ and this code fits the model

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

where $X_i$ is dependent_variable, $t$ is time and $i$ is ID. $\beta_0, \beta_1$ are fixed effects and $\eta_{i1}, \eta_{i2}$ are random effects whose variances are estimated by the model.

Here's an example with some quickly simulated data where the random effect variances truly are 0, the covariate has no effect, every outcome is ${\rm Poisson}(1)$, and each individual is seen 10 times at times $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

One point of caution - the Std.Dev. column is just the square root of the Variance column, NOT the standard error of the variance estimate!

Answer 2

In SAS:

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

But then of course there are plenty of options, more or less useful, to play with.