### R FAQRandom Coefficient Poisson Models

Version info: Code for this page was tested in R Under development (unstable) (2012-07-05 r59734) On: 2012-07-08 With: knitr 0.6.3 

Please Note: The purpose of this page is to show how to use various data analysis commands. It does not cover all aspects of the research process which researchers are expected to do. In particular, it does not cover data cleaning and checking, verification of assumptions, model diagnostics and potential follow-up analyses.

#### Examples of random coefficient poisson models

Poisson models are useful for count data. Examples could include number of traffic tickets an individual receives in a year, number of tumor sites in cancer patients, or number of awards received by students. In each of these cases, we might expect most people to have very few, with a relatively small number of individuals having higher numbers.

#### Background on the Poisson distribution

Unlike the familiar Gaussian distribution which has two parameters $$\mathcal{N}(\mu, \sigma^{2})$$, the Poisson distribution is described by a single parameter, $$\lambda$$ that is both the mean and variance. That is, $$\lambda = E(x)$$ and $$\lambda = Var(x) = E(x^{2}) - E(x)^{2}$$. In many ways, this can be a strong assumption. In some cases, use of the negative binomial distribution which estimates a second scale parameter to allow for overdispersion. Here is what the Poisson distribution looks like for different values of $$\lambda$$.

par(mfrow = c(3, 2), mar = c(1, 1, 2, 1), oma = c(4, 4, 2, 0))
for (i in c(1, 2, 3, 5, 7, 9)) {
curve(dpois(x, lambda = i), from = 0, to = 18, n = 19, type = "p", pch = 15,
cex = 1.5, xlim = c(0, 19), ylim = c(0, 0.4), xlab = "", ylab = "")
text(x = 15, y = 0.3, substitute(lambda == x, list(x = i)), cex = 2)
}
mtext("Probability Mass Function for Poisson Distribution", line = 0.5,
outer = TRUE, cex = 1.2)
mtext(expression(P(X == x)), side = 2, line = 1.5, outer = TRUE)
mtext("X", side = 1, line = 1.5, outer = TRUE)


Without random coefficients, the standard Poisson model is: $log[E(\mathbf{y_{i}})] = \boldsymbol{\alpha + X'_{i}\beta}$ The log link is the canonical link function for the Poisson distribution, and the expected value of the response is modeled.

With random coefficients, for example a random intercept, the model becomes: $log[E(\mathbf{y_{ij}|u_{j}})] = \boldsymbol{\alpha + X'_{ij}\beta + u_{j}}$ Where $$y_{ij}$$ is the observation for individual $$i$$ in group $$j$$ and $$u_{j}$$ is the random effect for group $$j$$. Thus the two distributions are: $y \sim Pois(\lambda)$ and $u \sim \mathcal{N}(0, \sigma^{2})$ The random coefficient model is conditional on the random effect. To show what this means, consider a simple model: $log[E(\mathbf{y_{ij}|u_{j}})] = -.5 + .3x_{ij} + u_{j}$ In the original units, this becomes: $E(\mathbf{y_{ij}|u_{j}}) = e^{-.5 + .3x_{ij} + u_{j}}$ Now look what happens when we graph the estimated change for a 1 unit change in x for values of the random variable (u) ranging from 0 to 4 by increments of .5.

f <- function(x_ij, u_j) {
exp(-0.5 + x_ij * 0.3 + u_j)
}
for (i in seq(0, 4, 0.5)) {
curve(f(x, u_j = i), from = 0, to = 1, n = 200, add = i > 0, ylim = c(0,
45), ylab = "")
}
title(ylab = bquote(E(Y ~ l ~ x[ij], u[j])), main = "Effect of x for different random effects")


Clearly, on the scale of the original units, a 1 unit increase in x has different effects depending on the value of u, hence the conditionalness of the model. Population "average" effects can be obtained by integrating out the random effect or by fitting a marginal model such as using GEEs.

Although the outcome is assumed to have a Poisson distribution, the random effect (in the above example, u) is typically assumed to have a Gaussian distribution.

#### Description of the data

The example data we use comes from a sample of the high school and beyond data set, with a made up variable, number of awards a student receives, awards. Our main predictor will be sex, female, and students are clustered (grouped) within schools, cid. First, we will load all the packages required for these examples.

## load foreign package to read Stata data files
require(foreign)

## ggplot2 package for graphs
require(ggplot2)

## if the package is not installed, install it from R-forge

require(lme4)

## read in data
dat$cid <- factor(dat$cid)
## look at the first few rows of the dataset

##    id female    ses schtyp     prog read write math science socst
## 1  45 female    low public vocation   34    35   41      29    26
## 2 108   male middle public  general   34    33   41      36    36
## 3  15   male   high public vocation   39    39   44      26    42
## 4  67   male    low public vocation   37    37   42      33    32
## 5 153   male middle public vocation   39    31   40      39    51
## 6  51 female   high public  general   42    36   42      31    39
##         honors awards cid
## 1 not enrolled      0   1
## 2 not enrolled      0   1
## 3 not enrolled      0   1
## 4 not enrolled      0   1
## 5 not enrolled      0   1
## 6 not enrolled      0   1


We can get a sense of the distributions in the data using the ggplot2 package. The first plot is just histograms of number of awards for every cid. The second is a filled density plot. The density sums to 1 and the fill shows the distribution of female at every level of awards. If the distribution of female is equal across all awards, they would fall on the horizontal line.

## awards by school
ggplot(dat, aes(awards)) + geom_histogram(binwidth = 0.5) + facet_wrap(~cid)

## density of awards by sex, line at .5 is the null of no sex differences
## in number of awards
ggplot(dat, aes(factor(awards))) + geom_bar(aes(fill = female), position = "fill") +
geom_hline(y = 0.5)


#### Analysis methods you might consider

• Poisson regression with robust standard errors

#### Random coefficient poisson model analysis

Because generalized linear mixed models (GLMMs) such as random coefficient poisson models are rather difficult to fit, there tends to be some variability in parameter estimates between different programs. We will demonstrate the use of two packages in R that are able to fit these models, lme4 and glmmADMB.

## fit a random intercept only model using the Laplace approximation
## (equivalent to 1 point evaluated per axis in Gauss-Hermite
## approximation)
m1a <- glmer(awards ~ 1 + (1 | cid), data = dat, family = poisson(link = "log"))
## fit a random intercept only model using 100 points per axis in the
## adaptive Gauss-Hermite approximation of the log likelihood more points
## improves accuracy but will take longer
m1b <- glmer(awards ~ 1 + (1 | cid), data = dat, family = poisson(link = "log"),
nAGQ = 100)

## compare (only slightly different)
rbind(m1a = coef(summary(m1a)), m1b = coef(summary(m1b)))

##              Estimate Std. Error  z value Pr(>|z|)
## (Intercept) -0.008853     0.2826 -0.03133   0.9750
## (Intercept) -0.009280     0.2834 -0.03275   0.9739

## view summary output from the more accurate model
summary(m1b)

## Generalized linear mixed model fit by the adaptive Gaussian Hermite approximation
## Formula: awards ~ 1 + (1 | cid)
##    Data: dat
##  AIC BIC logLik deviance
##  208 214   -102      204
## Random effects:
##  Groups Name        Variance Std.Dev.
##  cid    (Intercept) 1.45     1.2
## Number of obs: 200, groups: cid, 20
##
## Fixed effects:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.00928    0.28340   -0.03     0.97


We can visually inspect the random effects with a Q-Q plot or a caterpillar plot

## QQ plot
plot(ranef(m1b))

## $cid  ## Caterpillar plot dotplot(ranef(m1b, postVar = TRUE))  ##$cid


The estimate for the intercept is essentially 0, although the random effects variance indicates that there is some variability in the intercepts between schools. Now we will add in female as an explanatory variable.

m2 <- glmer(awards ~ 1 + female + (1 | cid), data = dat, family = poisson(link = "log"),
nAGQ = 100)
summary(m2)

## Generalized linear mixed model fit by the adaptive Gaussian Hermite approximation
## Formula: awards ~ 1 + female + (1 | cid)
##    Data: dat
##  AIC BIC logLik deviance
##  200 210  -97.1      194
## Random effects:
##  Groups Name        Variance Std.Dev.
##  cid    (Intercept) 1.43     1.19
## Number of obs: 200, groups: cid, 20
##
## Fixed effects:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -0.223      0.291   -0.76   0.4443
## femalefemale    0.363      0.119    3.04   0.0024 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
##             (Intr)
## femalefemal -0.260


There appears to be a fairly strong effect of females such that females tend to get more awards than males. Now we will fit the same models using the glmmADMB package

## random intercept only model
m.alt1 <- glmmadmb(awards ~ 1 + (1 | cid), data = dat, family = "poisson",
m.alt2 <- glmmadmb(awards ~ 1 + female + (1 | cid), data = dat, family = "poisson",

summary(m.alt1)

##
## Call:
## glmmadmb(formula = awards ~ 1 + (1 | cid), data = dat, family = "poisson",
##
##
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -0.0088     0.2890   -0.03     0.98
##
## Number of observations: total=200, cid=20
## Random effect variance(s):
## Group=cid
##             Variance StdDev
## (Intercept)    1.452  1.205
##
## Log-likelihood: -285.5

summary(m.alt2)

##
## Call:
## glmmadmb(formula = awards ~ 1 + female + (1 | cid), data = dat,
##     family = "poisson", link = "log")
##
##
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept)    -0.222      0.296   -0.75   0.4533
## femalefemale    0.363      0.119    3.04   0.0023 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of observations: total=200, cid=20
## Random effect variance(s):
## Group=cid
##             Variance StdDev
## (Intercept)    1.426  1.194
##
## Log-likelihood: -280.8


The results from glmmadmb match closely with those from glmer.

#### Things to consider

There have been reports (e.g., Zhang et al, 2009) of inconsistency and substantial alpha inflation between different statistical packages and estimation techniques (e.g., using penalized quasi-likelihood, Laplace approximation, Gauss-Hermite quadrature).

Although not the primary focus of this page, because of these known difficulties with GLMMs, we also show the Stata and SAS code to run the same models using a variety of estimation techniques.

#### SAS


data dat;
set "Z:\data\hsbdemo.sas7bdat";
run;

proc sort data=dat;
by cid descending female;
run;

PROC GLIMMIX data = dat method=LAPLACE noclprint;
class cid;
model awards = / solution dist=poisson;
random intercept / subject=cid;
run;

PROC GLIMMIX data = dat method=QUAD(QPOINTS = 100) noclprint;
class cid;
model awards = / solution dist=poisson;
random intercept / subject=cid;
run;

PROC GLIMMIX data = dat method=LAPLACE noclprint order = data;
class female cid;
model awards = female / solution dist=poisson;
random intercept / subject=cid;
run;

PROC GLIMMIX data = dat method=QUAD(QPOINTS = 100) noclprint order = data;
class female cid;
model awards = female / solution dist=poisson;
random intercept / subject=cid;
run;

/*We can also manually parameterize the model in NLMIXED */
PROC NLMIXED data = dat METHOD=GAUSS QPOINTS=100;
eta = alpha + u;
lambda = exp(eta);
model awards ~ poisson(lambda);
random u ~ normal(0, exp(2*log_sigma)) subject=cid;
estimate 'intercept variance' exp(2*log_sigma);
run;

PROC NLMIXED data = dat METHOD=GAUSS QPOINTS=100;
eta = alpha + beta*female + u;
lambda = exp(eta);
model awards ~ poisson(lambda);
random u ~ normal(0, exp(2*log_sigma)) subject=cid;
estimate 'intercept variance' exp(2*log_sigma);
run;


Random intercept only model, Laplace approximation)


Covariance Parameter Estimates
Standard
Cov Parm     Subject    Estimate       Error
Intercept    CID          1.4429      0.5970
Solutions for Fixed Effects
Standard
Effect       Estimate       Error       DF    t Value    Pr > |t|
Intercept    -0.00877      0.2890       19      -0.03      0.9761


Random intercept only model, adaptive Gauss-Hermite approximation with 100 evaluation points)


Covariance Parameter Estimates
Standard
Cov Parm     Subject    Estimate       Error
Intercept    CID          1.4579      0.6057
Solutions for Fixed Effects
Standard
Effect       Estimate       Error       DF    t Value    Pr > |t|
Intercept    -0.00954      0.2903       19      -0.03      0.9741


Random intercept plus female model, adaptive Gauss-Hermite approximation with 100 evaluation points)


Covariance Parameter Estimates
Standard
Cov Parm     Subject    Estimate       Error
Intercept    CID          1.4314      0.5952
Solutions for Fixed Effects
Standard
Effect       FEMALE    Estimate       Error       DF    t Value    Pr > |t|
Intercept               -0.2229      0.2975       19      -0.75      0.4629
FEMALE       1           0.3632      0.1193      179       3.04      0.0027
FEMALE       0                0           .        .        .         .


Random intercept only model, adaptive Gauss-Hermite approximation with 100 evaluation points) using PROC NLMIXED


Parameter Estimates
Standard
Parameter  Estimate     Error    DF  t Value  Pr > |t|   Alpha     Lower     Upper  Gradient
alpha      -0.00957    0.2903    19    -0.03    0.9740    0.05   -0.6172    0.5980  2.881E-6
log_sigma    0.1885    0.2077    19     0.91    0.3756    0.05   -0.2463    0.6233  7.128E-8
Standard
Label                Estimate      Error     DF   t Value   Pr > |t|    Alpha      Lower      Upper
intercept variance     1.4579     0.6057     19      2.41     0.0264     0.05     0.1901     2.7257


Random intercept plus female model, adaptive Gauss-Hermite approximation with 100 evaluation points) using PROC NLMIXED


Parameter Estimates
Standard
Parameter  Estimate     Error    DF  t Value  Pr > |t|   Alpha     Lower     Upper  Gradient
alpha       -0.2229    0.2975    19    -0.75    0.4629    0.05   -0.8456    0.3998   8.17E-8
beta         0.3632    0.1193    19     3.04    0.0067    0.05    0.1134    0.6129  -0.00004
log_sigma    0.1793    0.2079    19     0.86    0.3991    0.05   -0.2558    0.6145  9.151E-6
Standard
Label                Estimate      Error     DF   t Value   Pr > |t|    Alpha      Lower      Upper
intercept variance     1.4314     0.5952     19      2.40     0.0265     0.05     0.1856     2.6773


#### Stata


use http://www.ats.ucla.edu/stat/data/hsbdemo, clear

*Random intercept model, Laplace
xtmepoisson awards || cid:, laplace

*Random intercept model, adaptive Gauss-Hermite, 100 points
xtmepoisson awards || cid:, intpoints(100)

*Random intercept and female
xtmepoisson awards female || cid:, laplace

*Random intercept and female
xtmepoisson awards female || cid:, intpoints(100)


Random intercept only model, Laplace approximation)


------------------------------------------------------------------------------
awards |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons |     -.0088   .2890208    -0.03   0.976    -.5752703    .5576703
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
cid: Identity                |
sd(_cons) |    1.20123   .2485066      .8008156    1.801856
------------------------------------------------------------------------------


Random intercept only model, adaptive Gauss-Hermite approximation with 100 evaluation points)


------------------------------------------------------------------------------
awards |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
_cons |  -.0095717   .2902932    -0.03   0.974    -.5785359    .5593924
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
cid: Identity                |
sd(_cons) |   1.207439   .2508294      .8035994    1.814225
------------------------------------------------------------------------------


Random intercept plus female model, adaptive Gauss-Hermite approximation with 100 evaluation points)


------------------------------------------------------------------------------
awards |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
female |   .3631578   .1193071     3.04   0.002     .1293202    .5969954
_cons |  -.2229269     .29752    -0.75   0.454    -.8060553    .3602015
------------------------------------------------------------------------------
------------------------------------------------------------------------------
Random-effects Parameters  |   Estimate   Std. Err.     [95% Conf. Interval]
-----------------------------+------------------------------------------------
cid: Identity                |
sd(_cons) |    1.19643   .2487479      .7960039    1.798288
------------------------------------------------------------------------------


#### References

• Categorical Data Analysis 2nd Ed. by Alan Agresti
• Regression Analysis of Count Data by Cameron and Trivedi
• Zhang, Lu, Feng, Thurston, Xia, Zhu, & Tu, 2009. On fitting generalized linear mixed-effects models for binary responses using different statistical packages. Statistics in Medicine. DOI: 10.1002/sim.4265.

The content of this web site should not be construed as an endorsement of any particular web site, book, or software product by the University of California.