Messaging conflicts in `estimate_means()` given offset and `estimate = 'population'`

[ejvalencia]

Summary

For models with an offset term, modelbased::estimate_means() will automatically detect the presence of an offset term and provide an informative message that the user can fix the offset term to a single value via the function's offset argument.

• For modelbased::estimate_means(..., estimate = 'population', offset = NULL) the message is printed (as expected).
• For modelbased::estimate_means(..., estimate = 'population', offset = 1) a new message is printed indicating that the offset term is always ignored when estimate = 'population'.

If the offset argument is always ignored under estimate = 'population' then it may be useful to suppress the "hint" that users can specify the offset argument when modelbased::estimate_means(..., estimate = 'population', offset = NULL).

Worked Example

# Installing required packages
if (!require("modelbased")) install.packages("modelbased", quiet = TRUE)

# Generating data
df <-
  data.frame(
    agecat = factor(c(1, 2, 3, 4, 5, 1, 2, 3, 4, 5)),
    smokes = factor(c(1, 1, 1, 1, 1, 0, 0, 0, 0, 0)),
    deaths = c(32, 104, 206, 186, 102, 2, 12, 28, 28, 31),
    pyears = c(52407, 43248, 28612, 12663, 5317, 18790, 10673, 5710, 2585, 1462)
  )

# Estimating age-sex adjusted expected deaths via regression w/ g-computation
# NOTE: The function prints a message suggesting we specify the `offset` arg.
glm(deaths ~ smokes * agecat + offset(log(pyears)), family = poisson, data = df) |>
  modelbased::estimate_means(by = 'smokes', 
                             predict = 'response', 
                             estimate = 'population',
                             offset = NULL)

#> Model contains an offset-term and you average predictions over the
#>   distribution of that offset. If you want to fix the offset to a specific
#>   value, for instance `1`, use `offset = 1` and set `estimate =
#>   "typical"`.
#>   We also found that the model has a log-transformed offset term. If you
#>   use the `offset` argument, the log-transformation will automatically be
#>   applied to the provided offset-value. I.e., consider using, for
#>   instance, `offset = 10` and not `offset = log(10)`.
#> Average Counterfactual Predictions
#> 
#> smokes |  Mean |   SE |         95% CI |     z
#> ----------------------------------------------
#> 0      | 54.54 | 5.47 | [43.83, 65.25] |  9.98
#> 1      | 77.43 | 3.09 | [71.38, 83.48] | 25.09
#> 
#> Variable predicted: deaths
#> Predictors modulated: smokes
#> Predictors averaged: agecat, pyears (1.8e+04)
#> Predictions are on the response-scale.

# The same, but specifying the `offset` arg.
# NOTE: The function prints a message saying the `offset` arg is ignored.
glm(deaths ~ smokes * agecat + offset(log(pyears)), family = poisson, data = df) |>
  modelbased::estimate_means(by = 'smokes', 
                             predict = 'response', 
                             estimate = 'population',
                             offset = 1000)

#> For `estimate = "population"`, predictions are averaged over the
#>   distribution of the offset and the `offset` argument is ignored. If you
#>   want to fix the offset to a specific value, for instance `1`, use
#>   `offset = 1` and set `estimate = "typical"`.
#> Average Counterfactual Predictions
#> 
#> smokes |  Mean |   SE |         95% CI |     z
#> ----------------------------------------------
#> 0      | 54.54 | 5.47 | [43.83, 65.25] |  9.98
#> 1      | 77.43 | 3.09 | [71.38, 83.48] | 25.09
#> 
#> Variable predicted: deaths
#> Predictors modulated: smokes
#> Predictors averaged: agecat, pyears (1e+03)
#> Predictions are on the response-scale.

[strengejacke]

I don't think this is contradicting. The first message (when offset is present , but argument offset is NULL) just informs you about the offset, and suggests that if you want to use it, you have to define a different estimate value.

In the second example, it tells that offset is specified, but ignored for estimate = "population".

[ejvalencia]

I agree that "contradicting" is a mischaracterization, and I take your point that the messages are informative under either setup.

My first reaction is that if you want to intentionally set estimate = 'population' (or estimate = 'average') because you want to average over the distribution of the offset, then the message workflow is a (minor) nuisance. But one can resolve this issue by setting verbose = FALSE so perhaps nothing new is needed here.

Separately, though, I think the messaging workflow exposes a compelling limitation.

# This code gets at expected counts for the "typical" observation
glm(deaths ~ smokes * agecat + offset(log(pyears)), family = poisson, data = df) |>
  modelbased::estimate_means(by = 'smokes', predict = 'response',
                             estimate = 'typical', offset = NULL)

# This code gets at expected rates for the "typical" observation
glm(deaths ~ smokes * agecat + offset(log(pyears)), family = poisson, data = df) |>
  modelbased::estimate_means(by = 'smokes', predict = 'response',
                             estimate = 'typical', offset = 1000)

# This code gets at expected counts averaged over the observed sample distribution
glm(deaths ~ smokes * agecat + offset(log(pyears)), family = poisson, data = df) |>
  modelbased::estimate_means(by = 'smokes', predict = 'response',
                             estimate = 'average', offset = NULL)

# This code gets at expected rates averaged over the observed sample distribution.
# However, to do this, you need to manually override the offset value in the data.
glm(deaths ~ smokes * agecat + offset(log(pyears)), family = poisson, data = df) |>
  modelbased::estimate_means(by = 'smokes', predict = 'response',
                             estimate = 'average', newdata = within(df, pyears <- 1000))

In short, there may be settings where users want to average over the sample distribution (e.g., estimate = 'average' or estimate = 'population') but they also want to fix the offset so as to recover rates instead of counts. The present messaging workflow directs users to set estimate = 'typical' if they want to fix the offset... but this is a fundamentally different task. Users can still manually "override" the offset to set it as a fixed value using newdata but I think that permitting the offset argument for use under estimate = 'population' and estimate = 'average' is the better option?

[strengejacke]

I have revised the handling of the offset argument for estimate = "population", however, it is still not supported for estimate = "average". The reason is that for this option, we first average across the actual data (i.e., predictions for each observation, than average and possibly stratify in by) and do not use a data grid here. This means that if you want predictions at specific values, these values actually have to appear in the data:

library(modelbased)

set.seed(1)
newdata <- data.frame(
  y = c(602, 38, 616, 256, 21, 723, 245, 176, 89, 1614, 31, 27, 313, 251, 345),
  x = as.factor(sample(letters[1:3], 15, replace = TRUE)),
  offset_1 = c(72, 50, 31, 30, 16, 25, 75, 16, 78, 40, 68, 25, 71, 52, 17)
)

m <- MASS::glm.nb(y ~ x, data = newdata)
moff <- MASS::glm.nb(y ~ x + offset(log(offset_1)), data = newdata)

# offset 100 is not in the data
estimate_means(moff, c("x", "offset_1 = 100"), estimate = "average")
#> Error:
#> ! None of the values specified for the predictors `offset_1` are available
#>   in the data. This is required for `estimate="average"`. Either use a
#>   different option for the `estimate` argument, or use values that are
#>   present in the data, such as `72`, `40` or `17`.

# offset 71 is not in the data, but only for level "a"
estimate_means(moff, c("x", "offset_1 = 71"), estimate = "average")
#> Average Predictions
#> 
#> x | offset_1 |   Mean |            95% CI
#> -----------------------------------------
#> a |       71 | 471.93 | [195.47, 1139.36]
#> 
#> Variable predicted: y
#> Predictors modulated: x, offset_1 = 71
#> Predictions are on the response-scale.

^{Created on 2026-06-22 with reprex v2.1.1}

This is a limitation for estimate = "average" because no data grid is used for that option, thus filtering only works for "valid" values.

[strengejacke]

I maybe have found an option how to include offsets for the average-option.