# Poisson regression on aggregated data with an interaction term

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?

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.

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.