4
$\begingroup$

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?

$\endgroup$
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.