Poisson rate outcomes

metaGLMM authors

Rates from aggregate counts

For a Poisson analysis, use a rate as the response and the corresponding exposure as ni. A small positive variance approximation is useful for the aggregate-data representation, especially when an observed count is zero.

poisson_dat <- data.frame(
  study = paste0("Study ", 1:12),
  sex = factor(rep(c("Female", "Male"), 6)),
  year = c(0, 0, 2, 2, 6, 6, 6, 6, 7, 7, 8, 8),
  events = c(13, 1, 16, 3, 38, 3, 13, 4, 5, 10, 11, 23),
  exposure = c(9.04, 8.46, 23.75, 23.75, 18.47, 18.47,
               9.38, 9.19, 2.27, 3.34, 2.72, 26.00)
)
poisson_dat$rate <- poisson_dat$events / poisson_dat$exposure
poisson_dat$vi <- 1 / pmax(0.5, poisson_dat$events)
poisson_dat
#>       study    sex year events exposure      rate         vi
#> 1   Study 1 Female    0     13     9.04 1.4380531 0.07692308
#> 2   Study 2   Male    0      1     8.46 0.1182033 1.00000000
#> 3   Study 3 Female    2     16    23.75 0.6736842 0.06250000
#> 4   Study 4   Male    2      3    23.75 0.1263158 0.33333333
#> 5   Study 5 Female    6     38    18.47 2.0573904 0.02631579
#> 6   Study 6   Male    6      3    18.47 0.1624256 0.33333333
#> 7   Study 7 Female    6     13     9.38 1.3859275 0.07692308
#> 8   Study 8   Male    6      4     9.19 0.4352557 0.25000000
#> 9   Study 9 Female    7      5     2.27 2.2026432 0.20000000
#> 10 Study 10   Male    7     10     3.34 2.9940120 0.10000000
#> 11 Study 11 Female    8     11     2.72 4.0441176 0.09090909
#> 12 Study 12   Male    8     23    26.00 0.8846154 0.04347826

Basic single-group meta-analysis and forest plot

An intercept-only model estimates the overall log rate while allowing each study to have its own random effect. This is the direct single-group analogue of a conventional random-effects meta-analysis.

poisson_mean_fit <- metaGLMM(
  rate ~ 1,
  data = poisson_dat,
  vi = poisson_dat$vi,
  ni = poisson_dat$exposure,
  tau2 = NA,
  family = poisson(link = "log"),
  tau2_var = TRUE,
  fast = TRUE,
  ghq_Q = 40L
)

summary(poisson_mean_fit)
#> Aggregate-data generalized linear mixed-effects meta-analysis
#> 
#> Call:
#> metaGLMM(formula = rate ~ 1, data = poisson_dat, vi = poisson_dat$vi, 
#>     ni = poisson_dat$exposure, tau2 = NA, family = poisson(link = "log"), 
#>     tau2_var = TRUE, fast = TRUE, ghq_Q = 40L)
#> 
#> Family:poisson(log)
#> Random-effect structure: row_intercept 
#> Integration: ghq 
#> 
#> Fixed effects:
#>             Estimate Std. Error z value Pr(>|z|)
#> (Intercept)  -0.3773     0.2723  -1.386    0.166
#> 
#> Heterogeneity:
#>   tau^2: 0.9192 
#>   tau:   0.9588 
#> 
#> Plug-in prediction interval for a future underlying true effect:
#>      fit  lower upper
#>  -0.3773 -2.331 1.576
#> 
#> Fit diagnostics:
#>   Optimizer: L-BFGS-B 
#>   Convergence code: 0 
#>   Boundary tau^2: no 
#>   Positive-definite Hessian: yes
coef(poisson_mean_fit)
#> (Intercept) 
#>   -0.377292
confint(poisson_mean_fit, method = "wald")
#>                 lower     upper
#> (Intercept) -0.910935 0.1563511
stopifnot(is.finite(poisson_mean_fit$tau), poisson_mean_fit$tau > 0)

Each forest row is one observed study rate. type = "exp" converts the study estimates, pooled intervals, and prediction interval from the log-rate scale to the rate scale.

poisson_studies <- as_metafor_data(
  poisson_mean_fit, labels = poisson_dat$study
)
head(poisson_studies)
#>           yi         vi slab
#> 1  0.3632902 0.07692308    1
#> 2 -2.1353492 1.00000000    2
#> 3 -0.3949938 0.06250000    3
#> 4 -2.0689702 0.33333333    4
#> 5  0.7214384 0.02631579    5
#> 6 -1.8175355 0.33333333    6

forest(
  poisson_mean_fit,
  labels = poisson_dat$study,
  type = "exp",
  xlab = "Rate",
  ci_methods = c("Wald", "profile", "SBC")
)

Meta-regression with moderators

The formula can include factors, continuous moderators, and interactions. The coefficient names are the names generated by model.matrix().

poisson_fit <- metaGLMM(
  rate ~ sex + year + sex:year,
  data = poisson_dat,
  vi = poisson_dat$vi,
  ni = poisson_dat$exposure,
  tau2 = NA,
  family = poisson(link = "log"),
  tau2_var = TRUE,
  fast = TRUE,
  ghq_Q = 40L
)

summary(poisson_fit)
#> Aggregate-data generalized linear mixed-effects meta-analysis
#> 
#> Call:
#> metaGLMM(formula = rate ~ sex + year + sex:year, data = poisson_dat, 
#>     vi = poisson_dat$vi, ni = poisson_dat$exposure, tau2 = NA, 
#>     family = poisson(link = "log"), tau2_var = TRUE, fast = TRUE, 
#>     ghq_Q = 40L)
#> 
#> Family:poisson(log)
#> Random-effect structure: row_intercept 
#> Integration: ghq 
#> 
#> Fixed effects:
#>              Estimate Std. Error z value Pr(>|z|)   
#> (Intercept)   -0.1737     0.4355  -0.399  0.69000   
#> sexMale       -2.6446     0.8959  -2.952  0.00316 **
#> year           0.1386     0.0792   1.750  0.08012 . 
#> sexMale:year   0.2253     0.1459   1.544  0.12259   
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> Heterogeneity:
#>   tau^2: 0.214 
#>   tau:   0.4626 
#> 
#> Fit diagnostics:
#>   Optimizer: L-BFGS-B 
#>   Convergence code: 0 
#>   Boundary tau^2: no 
#>   Positive-definite Hessian: yes
coef(poisson_fit)
#>  (Intercept)      sexMale         year sexMale:year 
#>   -0.1736831   -2.6445864    0.1385934    0.2253256
exp(coef(poisson_fit))
#>  (Intercept)      sexMale         year sexMale:year 
#>   0.84056321   0.07103473   1.14865692   1.25273057

The exponentiated coefficients are rate ratios relative to the factor reference level and at year = 0. Center year before fitting when another reference year is more meaningful. Because this model has several fixed effects, a pooled forest summary would require an explicitly chosen coefficient or contrast; the basic intercept-only fit above is the clearer study-rate display.

Offsets

An offset is part of the ordinary formula interface. For example, when the response is defined on a count scale, the model matrix contains the exposure offset as follows:

count_formula <- events ~ sex + year + offset(log(exposure))
model.matrix(count_formula, data = poisson_dat)

Choose either a rate response with exposure supplied through ni or a count response with a formula offset according to the sampling convention of the analysis; keep the response and the exposure definition consistent.

mirror server hosted at Truenetwork, Russian Federation.