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

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?

-

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!

-
And its ηi1 that results in the random intercept? –  EpiGrad May 6 '12 at 2:24
$\eta_{i1}$ is the random intercept, yes. –  Macro May 6 '12 at 2:27
Thanks for the answer - one more question. In some of the GLMs, some of the sites benefit greatly from a t^2 term. It seems that most people tag out at one or two random effects - how nasty is a third one going to be in terms of computation? –  EpiGrad May 6 '12 at 2:39
Well, in the simulated example above, which only had 100 observations and 10 groups, adding the third random effect (by typing g <- lmer(y ~ x + (1+t+I(t^2)|ID), family="poisson")), increased computation time from about .75 seconds to about 11 seconds. As the sample size grows the increase in computing time probably also increases. –  Macro May 6 '12 at 2:46
@andrea, that would reflect a belief that there is a secular trend in the data set - I didn't assume a priori that this made sense. If the units were people and $t$ were age, then I certainly agree that a fixed effect in time made a lot of sense, but in other situations, the slope in time would be more of a random direction each individual is headed. That's why I didn't include that effect (and epigrad didn't ask how to include the fixed time effect) –  Macro May 6 '12 at 13:50

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.

-
Thanks :) Sadly, I seem to have run aground on convergence issues, but I'll tinker with those. –  EpiGrad May 6 '12 at 7:49