Poisson rate outcomes

The hardware and bandwidth for this mirror is donated by METANET, the Webhosting and Full Service-Cloud Provider.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]metanet.ch.

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.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.