Poisson regression on aggregated data with an interaction term

Question

For me, Poisson regression has been a nice tool to estimate risk ratios (setting offset to log-number of group size) and rate ratios (setting offset to log-risktime). Recently I came across a situation where the estimates on individual-level data did not agree with the estimates based on aggregated data, when an interaction term was present.

Here is an example in R:

# Full data:
a <- read.table(textConnection("Z X Y t n
m e 0 0.5 1
m e 1 1.5 1
m j 0 1 1
m j 1 0.5 1
n e 0 1 1
n e 0 1 1
n e 0 0.5 1
n j 1 0.5 1"), header = T, sep = " ")
# Aggregated data, grouped by X, Z, and t:
b <- read.table(textConnection("Z X Y t n
m e 0 0.5 1
m e 1 1.5 1
m j 0 1 1
m j 1 0.5 1
n e 0 2 2
n e 0 0.5 1
n j 1 0.5 1"), header = T, sep = " ")
# Aggregated data, grouped by X and Z:
d <- read.table(textConnection("Z X Y t n
m e 1 2 2
m j 1 1.5 2
n e 0 2.5 3
n j 1 0.5 1"), header = T, sep = " ")

Log risk ratios without interaction:

formu <- as.formula("Y ~ X + Z + offset(log(n))")
ma <- glm(formu, poisson, a); mb <- glm(formu, poisson, b); md <- glm(formu, poisson, d)
cbind(ma = coef(ma), mb = coef(mb), md = coef(md))

                    ma         mb         md
(Intercept) -1.3862944 -1.3862944 -1.3862944
Xj           1.0986123  1.0986123  1.0986123
Zn          -0.4054651 -0.4054651 -0.4054651

df <- data.frame(Z = "n", X = "e", n = 1)
c(ma = predict(ma, df), mb = predict(mb, df), md = predict(md, df))

     ma.1      mb.1      md.1 
-1.791759 -1.791759 -1.791759 

Log risk ratios with an interaction:

formu <- as.formula("Y ~ X*Z + offset(log(n))")
ma <- glm(formu, poisson, a); mb <- glm(formu, poisson, b); md <- glm(formu, poisson, d)
cbind(ma = coef(ma), mb = coef(mb), md = coef(md))

                     ma            mb            md
(Intercept)  -0.6931472 -6.931472e-01 -6.931472e-01
Xj            0.0000000  4.351168e-15  4.351168e-15
Zn          -19.6094379 -1.995601e+01 -2.270805e+01
Xj:Zn        20.3025851  2.064916e+01  2.340120e+01

c(ma = predict(ma, df), mb = predict(mb, df), md = predict(md, df))

     ma.1      mb.1      md.1 
-20.30259 -20.64916 -23.40120 

Log rate ratios without interaction:

formu <- as.formula("Y ~ X + Z + offset(log(t))")
ma <- glm(formu, poisson, a); mb <- glm(formu, poisson, b); md <- glm(formu, poisson, d)
cbind(ma = coef(ma), mb = coef(mb), md = coef(md))

                     ma          mb          md
(Intercept) -1.45485965 -1.45485965 -1.45485965
Xj           1.47670839  1.47670839  1.47670839
Zn          -0.09041403 -0.09041403 -0.09041403

df <- data.frame(Z = "n", X = "e", t = 1)
c(ma = predict(ma, df), mb = predict(mb, df), md = predict(md, df))

     ma.1      mb.1      md.1 
-1.545274 -1.545274 -1.545274 

Log rate ratios with interaction:

formu <- as.formula("Y ~ X*Z + offset(log(t))")
ma <- glm(formu, poisson, a); mb <- glm(formu, poisson, b); md <- glm(formu, poisson, d)
cbind(ma = coef(ma), mb = coef(mb), md = coef(md))

                     ma          mb          md
(Intercept)  -0.6931472  -0.6931472  -0.6931472
Xj            0.2876821   0.2876821   0.2876821
Zn          -19.3783888 -19.6094379 -22.5257286
Xj:Zn        20.4770011  20.7080502  23.6243409

c(ma = predict(ma, df), mb = predict(mb, df), md = predict(md, df))

     ma.1      mb.1      md.1 
-20.07154 -20.30259 -23.21888

How to explain the difference of estimates when the regression includes an interaction and how do the interpretations differ?

Answer 1

Consider what $\exp(20)$ is -- about half a billion!

In this example your problem is that the true maximum likelihood estimate for the coefficient of Zn is $-\infty$ and that of Xj:Zn is $+\infty$. You aren't going to be able to get actually infinite estimates out of glm, so the printed results of about $\pm 20$ are just the values where the algorithm failed. There's no particular reason why it should fail at exactly the same point with aggregated and unaggregated data; in the fourth example, glm does 18 iterations for a and b and 20 iterations for d

The reason you get infinite estimates is partly the $Y=0$ entry: observing 0 in a Poisson variable gives you some sort of upper bound on the mean, but no lower bound at all. Since there's no replication, this 0 is all the at that covariate combination and the fitted model has $E[Y]=0$ for that cell, which is a linear predictor of $-\infty$.

So: these results are as equal as they can reasonably be, given the limits of floating point computation. They aren't a counterexample to aggregation.

Answer 2

This is my first response on StackExchange, so my apologies if it's wonky.

It seems as though your problem with the interaction models stems from the fact that it's estimating 3 parameters and the intercept from a data frame with 4 rows; there are 0 residual degrees of freedom.

Using an example where there is another variable, Q, that isn't aggregated (and therefore gives you more than 4 rows of data to estimate from), it appears that the Poisson models using aggregated data are the same as un-aggregated data, regardless of the inclusion of the interaction term.

set.seed(221)
# Semi lazy way of getting data in the same format as yours, with additional variable Q
a <- expand.grid(Z=c("m","n"),
             X=c("e","j"),
             Q=c("l","y","f"))
a <- rbind(a,a,a)
a$Y <- rbinom(n=36,size=1,prob=0.5)
a$t <- rbinom(n=36,size=3,prob=0.5)*0.5 +0.5
a$n <- 1

# Aggregate by Z, X, t and Q
b <- aggregate(a[,c("Y","n")],
           by=list(Z=a$Z,
                   X=a$X,
                   t=a$t,
                   Q=a$Q),
           sum)

# Aggregate by Z, X and Q
d <- aggregate(a[,c("Y","t","n")],
           by=list(Z=a$Z,
                   X=a$X,
                   Q=a$Q),
           sum)

### Log risk ratios no interaction
formu <- as.formula("Y ~ X + Z + offset(log(n))")
ma <- glm(formu, poisson, a); mb <- glm(formu, poisson, b); md <- glm(formu, poisson, d)
cbind(ma = coef(ma), mb = coef(mb), md = coef(md))

             ma         mb         md
(Intercept) -0.9650809 -0.9650809 -0.9650809
Xj           0.4855078  0.4855078  0.4855078
Zn           0.2876821  0.2876821  0.2876821

df <- data.frame(Z = "n", X = "e", n = 1)
c(ma = predict(ma, df), mb = predict(mb, df), md = predict(md, df))

      ma.1       mb.1       md.1 
-0.6773988 -0.6773988 -0.6773988

### Log risk ratios interaction
formu <- as.formula("Y ~ X*Z + offset(log(n))")
ma <- glm(formu, poisson, a); mb <- glm(formu, poisson, b); md <- glm(formu, poisson, d)
cbind(ma = coef(ma), mb = coef(mb), md = coef(md))

                    ma            mb            md
(Intercept) -0.8109302 -8.109302e-01 -8.109302e-01
Xj           0.2231436  2.231436e-01  2.231436e-01
Zn           0.0000000 -7.083498e-13  1.478183e-16
Xj:Zn        0.4700036  4.700036e-01  4.700036e-01
c(ma = predict(ma, df), mb = predict(mb, df), md = predict(md, df))
      ma.1       mb.1       md.1 
-0.8109302 -0.8109302 -0.8109302

### Log rate ratios without interaction
formu <- as.formula("Y ~ X + Z + offset(log(t))")
ma <- glm(formu, poisson, a); mb <- glm(formu, poisson, b); md <- glm(formu, poisson, d)
cbind(ma = coef(ma), mb = coef(mb), md = coef(md))

                    ma         mb         md
(Intercept) -1.1992298 -0.9732687 -1.1992298
Xj           0.4945926  0.6111582  0.4945926
Zn           0.3920733  0.2734027  0.3920733

df <- data.frame(Z = "n", X = "e", t = 1)
c(ma = predict(ma, df), mb = predict(mb, df), md = predict(md, df))

      ma.1       mb.1       md.1 
-0.8071565 -0.6998660 -0.8071565 

### Log rate ratios with interaction
formu <- as.formula("Y ~ X*Z + offset(log(t))")
ma <- glm(formu, poisson, a); mb <- glm(formu, poisson, b); md <- glm(formu, poisson, d)
cbind(ma = coef(ma), mb = coef(mb), md = coef(md))

                     ma         mb          md
(Intercept) -1.01160091 -0.9162907 -1.01160091
Xj           0.17869179  0.5108256  0.17869179
Zn           0.04652002  0.1625189  0.04652001
Xj:Zn        0.56324556  0.1823216  0.56324556

c(ma = predict(ma, df), mb = predict(mb, df), md = predict(md, df))

      ma.1       mb.1       md.1 
-0.9650809 -0.7537718 -0.9650809

I guess mb is different in the rate ratio models because b is aggregated by t, i.e., aggregated by the offset. Since this is unusual practice, perhaps it doesn't warrant further thought here.