Handling missing data • drmTMB

Missing data in drmTMB has two different meanings. A response can be missing, so the observation has no measured outcome for one response cell. A predictor can be missing, so the model needs a value of a covariate that was not observed. The current implementation treats those cases differently. Decide this while writing the formula for the family chosen in Choosing response families, before fitting.

For missing responses, drmTMB can keep the original rows and evaluate the observed-response likelihood. Missing response cells contribute no direct response likelihood, but the fitted object still knows where those rows were. For missing predictors, drmTMB requires an explicit predictor model through mi() and impute. The predictor model is part of the same likelihood as the response model; it is not a preprocessing fill-in step.

The broadest missing-predictor module covers one missing predictor at a time in a univariate Gaussian location model. The supported predictor families are family-aware: Gaussian for continuous predictors, Bernoulli/logit for binary predictors, cumulative logit for ordered categories, baseline softmax for unordered categories, beta and zero-one beta for proportions, beta-binomial for successes out of known trials, Poisson/NB2/zero-truncated NB2 for counts, lognormal/Gamma for positive continuous predictors, and Tweedie for semi-continuous predictors with exact zeros. The non-Gaussian response routes support a poisson(), binomial(), nbinom2(), or beta() response with one binary missing predictor. Complete covariates used to explain the missing predictor can still be ordinary model-matrix inputs. Transformed mi() terms and multiple missing predictors are later work.

On the response-masking side, missing = miss_control(response = "include") retains and marginalises missing responses for gaussian(), bivariate Gaussian, student(), skew_normal(), lognormal(), Gamma(link = "log"), tweedie(), zero_one_beta(), binomial(), poisson(), nbinom2(), and beta() fits, plus whole-row beta_binomial() and ordered-factor cumulative_logit() fits, the non-hurdle truncated_nbinom2() route, and fixed-effect zero-inflated Poisson, zero-inflated NB2, and hurdle NB2 routes. Student-t, lognormal, Gamma, beta-binomial, and truncated NB2 include ordinary random-intercept recovery evidence; skew-normal, Tweedie, and zero-one beta are fixed-effect routes. These ticks do not automatically extend to every structured modifier supported by a family.

The response model and the imputation model use two different formula containers. The main drmTMB() model uses drm_formula() or its short alias bf(). Entries inside impute = list(...) are ordinary two-sided predictor formulas, such as body_mass ~ moisture + canopy, or impute_model(...) objects when a predictor-model family must be stated explicitly. Do not wrap the impute entry itself in drm_formula().

Capability matrix

The two missing-data modes are validated per response family. A checkmark means the mode is implemented and recovery-tested; a dash means it is not implemented and the fit rejects loudly (a family-specific error), so an unsupported request never silently degrades to a wrong likelihood.

Response family	`response = "include"` (mask/marginalise)	`predictor = "model"` (`mi()`)
`gaussian()`	✓	✓ (broad predictor-family catalogue)
bivariate Gaussian	✓ (partial `y1`/`y2` pairs)	—
`binomial()`	✓	✓ (one binary predictor)
`poisson()`	✓	✓ (one binary predictor)
`nbinom2()`	✓	✓ (one binary predictor)
`beta()`	✓	✓ (one binary predictor)
`student()`	✓ (ordinary random-intercept evidence)	—
`lognormal()`	✓ (ordinary random-intercept evidence)	—
`Gamma(link = "log")`	✓ (ordinary random-intercept evidence)	—
`skew_normal()`	✓ (fixed effects only)	—
`tweedie()`	✓ (fixed effects; zero atom plus continuous component)	—
`zero_one_beta()`	✓ (fixed effects; zero/one atoms plus interior)	—
`beta_binomial()`	✓ (whole success/failure row; ordinary random-intercept evidence)	—
`cumulative_logit()`	✓ (ordered factors; all declared categories observed)	—
`truncated_nbinom2()`	✓ (non-hurdle; ordinary `mu` random-intercept evidence)	—
zero-inflated Poisson	✓ (fixed effects; zeros and positives tested separately)	—
zero-inflated NB2	✓ (fixed effects; zeros and positives tested separately)	—
hurdle NB2	✓ (fixed effects; hurdle mass and truncated positive component)	—

The package keeps separate single-source-of-truth inventories for response and predictor missingness. A positive runtime test reconciles all 18 fitted response routes with the generated response ledger, while the predictor anti-drift test still requires every family outside drm_missing_predictor_families() to reject. Route-specific response tests retain the unsupported-neighbour checks. Only the Gaussian response carries the broad predictor-family catalogue tabulated below; the non-Gaussian responses currently model one binary (Bernoulli/logit) missing predictor each.

What is fitted now

The default is still conservative complete-case fitting:

miss_control()
#> $response
#> [1] "drop"
#> 
#> $predictor
#> [1] "fail"
#> 
#> $engine
#> [1] "laplace"
#> 
#> attr(,"class")
#> [1] "drm_missing_control"

Use missing = miss_control(response = "include") when you want supported missing responses retained in the fitted object. Use missing = miss_control(predictor = "model") only when the response formula contains one supported mi() term and impute supplies the corresponding predictor model.

Missing-data case	Current route	Main boundary
Missing univariate Gaussian response	`missing = miss_control(response = "include")`	Predictors, grouping variables, offsets, weights, and structured-effect inputs must be complete.
Missing non-Gaussian response (`student()`, `skew_normal()`, `lognormal()`, `Gamma()`, `tweedie()`, `zero_one_beta()`, `binomial()`, `poisson()`, `nbinom2()`, `beta()`)	`missing = miss_control(response = "include")`	Missing response cells are marginalised out; predictors and other inputs must be complete. Lognormal/Gamma observed responses must be positive; Tweedie allows zero and positive values; zero-one beta needs observed interior data plus optional 0/1 atoms; strict `beta()` needs interior `(0, 1)` responses.
Missing bivariate Gaussian response cell	`missing = miss_control(response = "include")`	Partial `y1`/`y2` rows are supported without dense known sampling covariance.
Missing beta-binomial response row	`missing = miss_control(response = "include")`	Missing either success or failure masks the whole encoded row; a partial pair never contributes.
Missing cumulative-logit response	`missing = miss_control(response = "include")`	Use an ordered factor with explicit levels; every declared category must remain represented among observed rows.
Missing zero-truncated NB2 response	`missing = miss_control(response = "include")`	Responses must be positive integers; the non-hurdle route is verified through an ordinary `mu` random intercept.
Missing count-mixture response	Add `zi ~ ...` to Poisson/NB2 or `hu ~ ...` to `truncated_nbinom2()` with `missing = miss_control(response = "include")`	The verified masking slice is fixed-effect only; missing zeros and positive counts are tested separately.
One missing numeric predictor	`bf(y ~ z + mi(x), sigma ~ 1)` with `impute = list(x = x ~ z)`	The missing `mi()` variable must be numeric and modelled with one Gaussian predictor model.
One missing binary predictor	`bf(y ~ z + mi(treatment), sigma ~ 1)` with `impute = list(treatment = impute_model(treatment ~ z, family = binomial()))`	The missing `mi()` variable must be binary and modelled with one fixed-effect Bernoulli/logit predictor model.
Non-Gaussian response with one missing binary predictor	`bf(y ~ z + mi(treatment))`, `family = poisson()`/`binomial()`/`nbinom2()`/`beta()`, with `impute = list(treatment = impute_model(treatment ~ z, family = binomial()))`	The response must be complete; the missing predictor must be binary; zero inflation, random effects, and structured response terms are not implemented for these routes.
One missing ordered predictor	`bf(y ~ z + mi(score), sigma ~ 1)` with `impute = list(score = impute_model(score ~ z, family = cumulative_logit()))`	The missing `mi()` variable must be ordered categorical and modelled with one fixed-effect cumulative-logit predictor model.
One missing unordered predictor	`bf(y ~ z + mi(habitat), sigma ~ 1)` with `impute = list(habitat = impute_model(habitat ~ z, family = categorical()))`	The missing `mi()` variable must be unordered categorical with at least three levels and modelled with one fixed-effect baseline-category softmax predictor model.
One missing strict proportion predictor	`bf(y ~ z + mi(cover), sigma ~ 1)` with `impute = list(cover = impute_model(cover ~ z, family = beta()))`	The missing `mi()` variable must be numeric and strictly inside `(0, 1)`, then modelled with one fixed-effect beta predictor model.
One missing boundary proportion predictor	`bf(y ~ z + mi(cover), sigma ~ 1)` with `impute = list(cover = impute_model(cover ~ z, family = zero_one_beta()))`	The missing `mi()` variable must be numeric in `[0, 1]`, may include exact zeros or ones, and is modelled with one fixed-effect zero-one beta predictor model.
One missing denominator-aware proportion predictor	`bf(y ~ z + mi(cover), sigma ~ 1)` with `impute = list(cover = impute_model(success ~ z, family = beta_binomial(), trials = trials))`	The `mi()` variable is the proportion used by the response model; the predictor model uses success counts and complete known trial denominators.
One missing Poisson count predictor	`bf(y ~ z + mi(abundance), sigma ~ 1)` with `impute = list(abundance = impute_model(abundance ~ z, family = poisson()))`	The missing `mi()` variable must be a non-negative integer count and modelled with one fixed-effect Poisson predictor model.
One missing negative-binomial count predictor	`bf(y ~ z + mi(abundance), sigma ~ 1)` with `impute = list(abundance = impute_model(abundance ~ z, family = nbinom2()))`	The missing `mi()` variable must be a non-negative integer count and modelled with one fixed-effect NB2 predictor model.
One missing zero-truncated count predictor	`bf(y ~ z + mi(abundance), sigma ~ 1)` with `impute = list(abundance = impute_model(abundance ~ z, family = truncated_nbinom2()))`	The missing `mi()` variable must be a positive integer count and modelled with one fixed-effect zero-truncated NB2 predictor model.
One missing lognormal positive predictor	`bf(y ~ z + mi(biomass), sigma ~ 1)` with `impute = list(biomass = impute_model(biomass ~ z, family = lognormal()))`	The missing `mi()` variable must be positive and modelled with one fixed-effect lognormal predictor model.
One missing Gamma positive predictor	`bf(y ~ z + mi(biomass), sigma ~ 1)` with `impute = list(biomass = impute_model(biomass ~ z, family = Gamma(link = "log")))`	The missing `mi()` variable must be positive and modelled with one fixed-effect Gamma mean-CV predictor model.
One missing Tweedie semi-continuous predictor	`bf(y ~ z + mi(biomass), sigma ~ 1)` with `impute = list(biomass = impute_model(biomass ~ z, family = tweedie()))`	The missing `mi()` variable must be non-negative, may include exact zeros, and is modelled with one fixed-effect Tweedie predictor model.
Grouped or structured predictor model	`impute = list(x = x ~ z + (1 \| group))` or `impute = list(x = x ~ z + relmat(1 \| line, Q = Q))`	One grouped random intercept or one intercept-only structured covariate field.
Missing-predictor summaries	`imputed(fit)`	Conditional modes for Gaussian predictors, conditional probabilities for binary predictors, conditional expected scores for ordered predictors, conditional modal categories for unordered predictors, conditional means for strict proportion, boundary proportion, denominator-aware beta-binomial, lognormal, Gamma, and Tweedie predictors, and conditional expected counts for count predictors.

This is likelihood-based missing-data handling. It is not Bayesian posterior imputation, not multiple imputation, and not an EM interface.

Choosing the predictor family

The predictor family should describe the missing covariate, not the response. This table is the practical starting point:

Missing predictor type	Use in `impute_model()`	What is integrated or summed
Continuous on the real line	bare `x ~ z` or `impute_model(x ~ z, family = gaussian())`	a Gaussian latent predictor by Laplace approximation
Binary, logical, or two-level factor	`impute_model(x ~ z, family = binomial())`	the two possible states
Ordered categories	`impute_model(x ~ z, family = cumulative_logit())`	all ordered levels
Unordered categories	`impute_model(x ~ z, family = categorical())`	all unordered levels
Continuous proportion in `(0, 1)`	`impute_model(x ~ z, family = beta())`	deterministic beta quadrature nodes
Continuous proportion in `[0, 1]` with real zeros/ones	`impute_model(x ~ z, family = zero_one_beta())`	exact zero/one mass plus beta quadrature
Successes out of known trials	`impute_model(success ~ z, family = beta_binomial(), trials = trials)`	all integer successes from `0` to `trials`
Non-negative counts	`impute_model(x ~ z, family = poisson())` or `impute_model(x ~ z, family = nbinom2())`	deterministic count states
Positive counts	`impute_model(x ~ z, family = truncated_nbinom2())`	deterministic positive count states
Positive continuous values	`impute_model(x ~ z, family = lognormal())` or `impute_model(x ~ z, family = Gamma(link = "log"))`	deterministic positive quadrature nodes
Non-negative semi-continuous values with exact zeros	`impute_model(x ~ z, family = tweedie())`	exact zero mass plus positive quadrature nodes

Rendering this article is also a smoke test: the worked examples below fit one model for each supported predictor-family route.

Relation to gllvmTMB and glmmTMB

The sibling gllvmTMB missing-data layer uses the same vocabulary: miss_control(), mi(), impute, fit$missing_data, and imputed() for predictors. gllvmTMB adds predict_missing() for masked responses; in drmTMB that information comes from fitted(), residuals(), and fit$missing_data. The gllvmTMB multivariate trait engine needs per-unit products across trait rows, so its missing-predictor implementation is not a drop-in copy of the scalar drmTMB route. The shared contract is the important part: missing responses are masked or marginalised, missing predictors are modelled inside the likelihood, and fitted summaries are not posterior draws.

glmmTMB has a different missing-data surface. Its documented na.action argument follows ordinary R model-frame handling: the default drops observations containing NA in model variables, while na.exclude still drops those observations for fitting but returns NA predictions and residuals for excluded cases (glmmTMB reference). The drmTMB route described here is therefore an in-model missing-data layer, not a wrapper around na.action.

Missing responses

For a univariate Gaussian response, the included-row likelihood is:

$\ell_i = I_i \log p(y_i \mid \mu_i, \sigma_i),$

where $I_i = 1$ when the response is observed and $I_i = 0$ when the response is missing. A missing response row is retained for row accounting, but it adds zero direct response likelihood.

This example creates a simple growth data set with three missing response values:

set.seed(20260531)
n <- 40
growth_data <- data.frame(
  temperature = seq(-1.6, 1.6, length.out = n)
)
growth_data$growth <- 0.5 + 0.8 * growth_data$temperature +
  rnorm(n, sd = 0.25)
growth_data$growth[c(6, 17, 32)] <- NA_real_

fit_response <- drmTMB(
  bf(growth ~ temperature, sigma ~ 1),
  family = gaussian(),
  data = growth_data,
  missing = miss_control(response = "include"),
  control = drm_control(se = FALSE)
)

nobs() reports the number of likelihood-contributing responses, not the number of retained rows:

nobs(fit_response)
#> [1] 37
fit_response$missing_data$counts
#> $retained_rows
#> [1] 40
#> 
#> $observed_response
#> [1] 37
#> 
#> $missing_response
#> [1] 3
#> 
#> $likelihood_rows
#> [1] 37

The fitted object preserves original-row accounting:

head(data.frame(
  original_row = fit_response$missing_data$original_row,
  observed_y = fit_response$missing_data$observed_y,
  fitted = fitted(fit_response),
  residual = residuals(fit_response)
), 10)
#>    original_row observed_y     fitted    residual
#> 1             1       TRUE -0.8553559 -0.04142081
#> 2             2       TRUE -0.7879308  0.14755851
#> 3             3       TRUE -0.7205056  0.16248531
#> 4             4       TRUE -0.6530804 -0.19202655
#> 5             5       TRUE -0.5856552 -0.35960220
#> 6             6      FALSE -0.5182301          NA
#> 7             7       TRUE -0.4508049  0.12518721
#> 8             8       TRUE -0.3833797  0.12351978
#> 9             9       TRUE -0.3159546  0.05197083
#> 10           10       TRUE -0.2485294  0.08098015

Rows with missing responses have fitted values because their predictors are known. Their residuals are NA because no response was observed.

Partial bivariate response pairs

For bivariate Gaussian models, one response can be observed while the other is missing. A complete pair contributes the bivariate Gaussian density with rho12. A row with only y1 observed contributes the marginal Gaussian density for y1; a row with only y2 observed contributes the marginal density for y2. A row with both responses missing contributes zero direct response likelihood but keeps its row identity.

Symbolically, the row contribution is:

both observed     -> log p(y1_i, y2_i | mu1_i, mu2_i, sigma1_i, sigma2_i, rho12_i)
only y1 observed  -> log p(y1_i | mu1_i, sigma1_i)
only y2 observed  -> log p(y2_i | mu2_i, sigma2_i)
both missing      -> 0

The residual correlation rho12 is learned directly from complete response pairs. Rows with one observed response still help the marginal location and scale models, but they do not directly identify rho12.

trait_data <- data.frame(
  x = seq(-1.4, 1.4, length.out = 36)
)
trait_data$activity <- 0.3 + 0.7 * trait_data$x + rnorm(36, sd = 0.25)
trait_data$boldness <- -0.1 - 0.4 * trait_data$x + rnorm(36, sd = 0.30)
trait_data$activity[c(5, 16, 29)] <- NA_real_
trait_data$boldness[c(8, 16, 31)] <- NA_real_

fit_pair <- drmTMB(
  bf(
    mu1 = activity ~ x,
    mu2 = boldness ~ x,
    sigma1 = ~1,
    sigma2 = ~1,
    rho12 = ~1
  ),
  family = c(gaussian(), gaussian()),
  data = trait_data,
  missing = miss_control(response = "include"),
  control = drm_control(se = FALSE)
)

fit_pair$missing_data$counts
#> $retained_rows
#> [1] 36
#> 
#> $observed_y1
#> [1] 33
#> 
#> $observed_y2
#> [1] 33
#> 
#> $complete_response
#> [1] 31
#> 
#> $y1_only
#> [1] 2
#> 
#> $y2_only
#> [1] 2
#> 
#> $both_missing
#> [1] 1
#> 
#> $likelihood_rows
#> [1] 35
table(fit_pair$missing_data$response_pattern)
#> 
#>  both_missing both_observed   y1_observed   y2_observed 
#>             1            31             2             2

Use the response-pattern counts before interpreting rho12. If complete pairs are scarce, drmTMB warns that the residual correlation is weakly identified.

Missing predictors

Missing predictors are not handled by dropping a row into the response likelihood with an arbitrary filled-in covariate. The response model needs the predictor value, so drmTMB asks for an explicit model for that predictor.

For the first implemented route, the response model is a univariate Gaussian location model:

$\begin{aligned} y_i \mid x_i, \mu_i, \sigma &\sim \operatorname{Normal}(\mu_i, \sigma^2),\\ \mu_i &= \beta_0 + \beta_1 z_i + \beta_2 x_i. \end{aligned}$

The impute formula defines a Gaussian model for the missing predictor:

$x_i \sim \operatorname{Normal}(W_i \alpha, \sigma_x^2).$

Observed x values contribute their predictor-model density and enter the response model directly. Missing x values are latent TMB random effects. TMB’s Laplace approximation integrates them from the joint likelihood.

set.seed(20260532)
n <- 48
predictor_data <- data.frame(
  moisture = seq(-1.5, 1.5, length.out = n),
  canopy = cos(seq_len(n) / 5)
)
predictor_data$body_mass_full <- 0.2 + 0.7 * predictor_data$moisture -
  0.2 * predictor_data$canopy + rnorm(n, sd = 0.08)
predictor_data$growth <- 0.6 + 1.1 * predictor_data$body_mass_full -
  0.3 * predictor_data$moisture + rnorm(n, sd = 0.20)
predictor_data$body_mass <- predictor_data$body_mass_full
predictor_data$body_mass[c(7, 19, 34, 43)] <- NA_real_

fit_predictor <- drmTMB(
  bf(growth ~ moisture + mi(body_mass), sigma ~ 1),
  family = gaussian(),
  data = predictor_data,
  impute = list(body_mass = body_mass ~ moisture + canopy),
  missing = miss_control(predictor = "model")
)

coef(fit_predictor, "mu")
#>   (Intercept)      moisture mi(body_mass) 
#>     0.6499983    -0.1405535     0.8939445
imputed(fit_predictor)
#>    variable original_row model_row observed   estimate  std_error
#> 1 body_mass            7         7    FALSE -0.5977839 0.06694745
#> 2 body_mass           19        19    FALSE  0.1230898 0.06692573
#> 3 body_mass           34        34    FALSE  0.4920569 0.06762264
#> 4 body_mass           43        43    FALSE  1.1676530 0.06729264
#>             source uncertainty_status
#> 1 conditional_mode                 ok
#> 2 conditional_mode                 ok
#> 3 conditional_mode                 ok
#> 4 conditional_mode                 ok

Read imputed(fit_predictor) as a fitted missing-predictor summary. The estimate column is the conditional mode used by the fitted likelihood. The std_error column is a likelihood-based conditional standard error when TMB::sdreport() is available. It is not a posterior interval and it is not a set of multiple-imputation draws.

Binary missing predictors

Binary predictors are a separate likelihood problem. They should not be forced through a Gaussian predictor model just because some values are missing. For a binary missing predictor, put the predictor formula in the impute list and state the predictor-model family with impute_model():

set.seed(20260533)
n <- 50
binary_data <- data.frame(
  z = seq(-1.6, 1.6, length.out = n)
)
binary_data$treatment_full <- as.numeric(
  sin(seq_len(n) * 1.7) + 0.35 * binary_data$z > 0
)
binary_data$growth <- 0.45 + 0.55 * binary_data$z +
  1.25 * binary_data$treatment_full + rnorm(n, sd = 0.08)
binary_data$treatment <- factor(binary_data$treatment_full, levels = c(0, 1))
binary_data$treatment[c(6, 15, 28, 43)] <- NA

fit_binary <- drmTMB(
  bf(growth ~ z + mi(treatment), sigma ~ 1),
  family = gaussian(),
  data = binary_data,
  impute = list(
    treatment = impute_model(treatment ~ z, family = binomial())
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_binary, "mi_treatment")
#> (Intercept)           z 
#>   0.1630819   0.2781205
imputed(fit_binary)
#>    variable original_row model_row observed     estimate std_error
#> 1 treatment            6         6    FALSE 9.476287e-49        NA
#> 2 treatment           15        15    FALSE 1.000000e+00        NA
#> 3 treatment           28        28    FALSE 1.611721e-45        NA
#> 4 treatment           43        43    FALSE 7.082334e-40        NA
#>                    source uncertainty_status
#> 1 conditional_probability   sdreport_skipped
#> 2 conditional_probability   sdreport_skipped
#> 3 conditional_probability   sdreport_skipped
#> 4 conditional_probability   sdreport_skipped

Here treatment ~ z is the formula for the missing predictor itself. It is not wrapped in drm_formula(). Internally, drmTMB sums over the two possible states for each missing treatment value. imputed(fit_binary) reports the fitted conditional probability that the missing predictor is in the second factor level.

Non-Gaussian responses with binary missing predictors

The non-Gaussian response routes use the same binary predictor model, but the response likelihood changes. This is useful when the outcome is a count, proportion, or Bernoulli/binomial outcome and a two-level covariate is partly missing. The first example uses a Poisson count response:

set.seed(20260543)
n <- 64
poisson_binary_data <- data.frame(
  z = seq(-1.5, 1.8, length.out = n)
)
poisson_binary_data$treatment_full <- as.numeric(
  sin(seq_len(n) * 1.3) + 0.3 * poisson_binary_data$z > 0
)
lambda_count <- exp(
  0.25 + 0.45 * poisson_binary_data$z +
    0.75 * poisson_binary_data$treatment_full
)
poisson_binary_data$count <- pmax(
  0,
  round(lambda_count + sqrt(lambda_count) * cos(seq_len(n) / 4))
)
poisson_binary_data$treatment <- factor(
  poisson_binary_data$treatment_full,
  levels = c(0, 1)
)
poisson_binary_data$treatment[c(8, 19, 31, 46, 57)] <- NA

fit_poisson_binary <- drmTMB(
  bf(count ~ z + mi(treatment)),
  family = poisson(),
  data = poisson_binary_data,
  impute = list(
    treatment = impute_model(treatment ~ z, family = binomial())
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_poisson_binary, "mu")
#>   (Intercept)             z mi(treatment) 
#>     0.2347720     0.3805174     0.7857868
coef(fit_poisson_binary, "mi_treatment")
#> (Intercept)           z 
#>   0.0381178   0.4203158
imputed(fit_poisson_binary)
#>    variable original_row model_row observed  estimate std_error
#> 1 treatment            8         8    FALSE 0.1947620        NA
#> 2 treatment           19        19    FALSE 0.3470786        NA
#> 3 treatment           31        31    FALSE 0.7055292        NA
#> 4 treatment           46        46    FALSE 0.9033351        NA
#> 5 treatment           57        57    FALSE 0.4029249        NA
#>                    source uncertainty_status
#> 1 conditional_probability   sdreport_skipped
#> 2 conditional_probability   sdreport_skipped
#> 3 conditional_probability   sdreport_skipped
#> 4 conditional_probability   sdreport_skipped
#> 5 conditional_probability   sdreport_skipped

For a missing treatment value, drmTMB sums over treatment 0 and treatment 1, combining the Bernoulli predictor-model probability with the Poisson response likelihood for the observed count. This route is intentionally narrow: the count response must be observed, and zero-inflated, grouped, or structured Poisson response models with mi() remain later work.

The overdispersed-count (nbinom2()) and proportion (beta()) responses reuse the identical 2-point sum; only the response density and its dispersion change. The negative-binomial response carries size = exp(-2 * log_sigma):

fit_nbinom2_binary <- drmTMB(
  bf(count ~ z + mi(treatment), sigma ~ 1),
  family = nbinom2(),
  data = poisson_binary_data,
  impute = list(
    treatment = impute_model(treatment ~ z, family = binomial())
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_nbinom2_binary, "mu")
#>   (Intercept)             z mi(treatment) 
#>     0.2347714     0.3805169     0.7857859
coef(fit_nbinom2_binary, "mi_treatment")
#> (Intercept)           z 
#>   0.0381136   0.4203145

The beta response models an interior proportion with a logit mean and precision phi = exp(-2 * log_sigma):

beta_binary_data <- poisson_binary_data
beta_binary_data$cover <- plogis(
  -0.2 + 0.5 * beta_binary_data$z +
    0.7 * beta_binary_data$treatment_full + 0.15 * cos(seq_len(nrow(beta_binary_data)))
)

fit_beta_binary <- drmTMB(
  bf(cover ~ z + mi(treatment), sigma ~ 1),
  family = beta(),
  data = beta_binary_data,
  impute = list(
    treatment = impute_model(treatment ~ z, family = binomial())
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_beta_binary, "mu")
#>   (Intercept)             z mi(treatment) 
#>    -0.1915409     0.5060492     0.6835154
coef(fit_beta_binary, "mi_treatment")
#> (Intercept)           z 
#> 0.001890477 0.421393905

A binomial() response works the same way with a cbind(successes, failures) or 0/1 outcome. Each route is intentionally narrow in the same way as the Poisson one: the response must be observed, and zero inflation, random effects, and structured response terms with mi() remain later work.

Ordered missing predictors

Ordered predictors are also finite-state missing predictors. They need their own ordinal predictor model rather than a Gaussian model on arbitrary category codes:

set.seed(20260535)
n <- 60
ordered_data <- data.frame(
  z = seq(-1.5, 1.5, length.out = n)
)
ordered_levels <- c("low", "medium", "high")
ordered_data$score_full <- cut(
  ordered_data$z + sin(seq_len(n) * 1.4),
  breaks = c(-Inf, -0.45, 0.55, Inf),
  labels = ordered_levels,
  ordered_result = TRUE
)
ordered_data$growth <- 0.25 + 0.5 * ordered_data$z +
  0.45 * as.numeric(ordered_data$score_full) + rnorm(n, sd = 0.18)
ordered_data$score <- ordered_data$score_full
ordered_data$score[c(8, 18, 36, 51)] <- NA

fit_ordered <- drmTMB(
  bf(growth ~ z + mi(score), sigma ~ 1),
  family = gaussian(),
  data = ordered_data,
  impute = list(
    score = impute_model(score ~ z, family = cumulative_logit())
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_ordered, "mi_score")
#>        z 
#> 2.086121
imputed(fit_ordered)
#>   variable original_row model_row observed estimate std_error
#> 1    score            8         8    FALSE 1.244165        NA
#> 2    score           18        18    FALSE 1.012190        NA
#> 3    score           36        36    FALSE 2.956817        NA
#> 4    score           51        51    FALSE 2.931307        NA
#>                       source uncertainty_status
#> 1 conditional_expected_score   sdreport_skipped
#> 2 conditional_expected_score   sdreport_skipped
#> 3 conditional_expected_score   sdreport_skipped
#> 4 conditional_expected_score   sdreport_skipped

Here the predictor model is a cumulative-logit model for the ordered score. For a missing score, drmTMB sums over all ordered categories. The reported estimate is the fitted conditional expected category score, and the missing-data metadata stores the conditional probability for each ordered level.

Unordered missing predictors

Unordered categorical predictors use a baseline-category softmax model. This is for nominal categories where the order of the levels is not meaningful:

set.seed(20260536)
n <- 66
categorical_data <- data.frame(
  z = seq(-1.7, 1.7, length.out = n)
)
categorical_score <- sin(seq_len(n) / 4) + 0.45 * categorical_data$z
categorical_data$habitat_full <- factor(
  ifelse(
    categorical_score < -0.35,
    "forest",
    ifelse(categorical_score < 0.55, "grass", "wetland")
  ),
  levels = c("forest", "grass", "wetland")
)
habitat_effect <- c(forest = -0.35, grass = 0.2, wetland = 0.75)
categorical_data$growth <- 0.15 + 0.5 * categorical_data$z +
  unname(habitat_effect[as.character(categorical_data$habitat_full)]) +
  rnorm(n, sd = 0.12)
categorical_data$habitat <- categorical_data$habitat_full
categorical_data$habitat[c(7, 16, 31, 48, 60)] <- NA

fit_categorical <- drmTMB(
  bf(growth ~ z + mi(habitat), sigma ~ 1),
  family = gaussian(),
  data = categorical_data,
  impute = list(
    habitat = impute_model(habitat ~ z, family = categorical())
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_categorical, "mi_habitat")
#>   grass:(Intercept)             grass:z wetland:(Intercept)           wetland:z 
#>          0.11840042          0.07411497         -0.46305607          1.48837772
imputed(fit_categorical)
#>   variable original_row model_row observed estimate std_error
#> 1  habitat            7         7    FALSE        2        NA
#> 2  habitat           16        16    FALSE        1        NA
#> 3  habitat           31        31    FALSE        3        NA
#> 4  habitat           48        48    FALSE        2        NA
#> 5  habitat           60        60    FALSE        3        NA
#>                       source uncertainty_status
#> 1 conditional_modal_category   sdreport_skipped
#> 2 conditional_modal_category   sdreport_skipped
#> 3 conditional_modal_category   sdreport_skipped
#> 4 conditional_modal_category   sdreport_skipped
#> 5 conditional_modal_category   sdreport_skipped

Here the first factor level is the baseline category in the predictor model. For a missing habitat value, drmTMB sums over all habitat levels. imputed(fit_categorical) reports the fitted modal category code, and the missing-data metadata stores the fitted conditional probability for each unordered level.

This is the pattern for broader non-Gaussian missing predictors too: put the predictor model in the impute list, and make the family explicit. The next sections use the same idea for strict proportion, boundary proportion, denominator-aware proportion, count, and positive continuous predictors.

Strict proportion missing predictors

Strict proportion predictors use a beta predictor model. This route is for continuous proportions inside (0, 1), not exact zeros, exact ones, or counted successes out of a known denominator:

set.seed(20260537)
n <- 72
proportion_data <- data.frame(
  z = seq(-1.8, 1.8, length.out = n)
)
proportion_data$cover_full <- plogis(
  -0.25 + 0.9 * proportion_data$z + 0.18 * sin(seq_len(n) / 5)
)
proportion_data$growth <- 0.35 + 1.25 * proportion_data$cover_full -
  0.30 * proportion_data$z + rnorm(n, sd = 0.10)
proportion_data$cover <- proportion_data$cover_full
proportion_data$cover[c(8, 19, 34, 51, 67)] <- NA_real_

fit_proportion <- drmTMB(
  bf(growth ~ z + mi(cover), sigma ~ 1),
  family = gaussian(),
  data = proportion_data,
  impute = list(
    cover = impute_model(cover ~ z, family = beta())
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_proportion, "mi_cover")
#> (Intercept)           z 
#>  -0.2354154   0.8854963
coef(fit_proportion, "sigma_mi_cover")
#>      cover 
#> 0.05899296
imputed(fit_proportion)
#>   variable original_row model_row observed  estimate std_error
#> 1    cover            8         8    FALSE 0.2073250        NA
#> 2    cover           19        19    FALSE 0.2934076        NA
#> 3    cover           34        34    FALSE 0.4005025        NA
#> 4    cover           51        51    FALSE 0.5943262        NA
#> 5    cover           67        67    FALSE 0.7827985        NA
#>                        source uncertainty_status
#> 1 conditional_quadrature_mean   sdreport_skipped
#> 2 conditional_quadrature_mean   sdreport_skipped
#> 3 conditional_quadrature_mean   sdreport_skipped
#> 4 conditional_quadrature_mean   sdreport_skipped
#> 5 conditional_quadrature_mean   sdreport_skipped

For a missing cover value, drmTMB integrates over the possible proportion values using deterministic quadrature under the fitted beta predictor model and the Gaussian response likelihood. imputed(fit_proportion) reports the fitted conditional quadrature mean for each missing cover value.

Boundary proportion missing predictors

Use a zero-one beta predictor model when a proportion covariate is continuous on [0, 1] and exact zeros or ones are real observations, not coding mistakes. The first route uses cover ~ z for the predictor-model mean and estimates constant interior scale, exact-boundary probability zoi, and conditional exact-one probability coi:

set.seed(20260541)
n <- 78
boundary_prop_data <- data.frame(
  z = seq(-1.8, 1.8, length.out = n)
)
boundary_prop_data$cover_full <- plogis(
  -0.15 + 0.8 * boundary_prop_data$z + 0.12 * sin(seq_len(n) / 4)
)
boundary_prop_data$cover_full[seq(6, n, by = 18)] <- 0
boundary_prop_data$cover_full[seq(13, n, by = 19)] <- 1
boundary_prop_data$growth <- 0.30 + 1.10 * boundary_prop_data$cover_full -
  0.25 * boundary_prop_data$z + rnorm(n, sd = 0.10)
boundary_prop_data$cover <- boundary_prop_data$cover_full
boundary_prop_data$cover[c(8, 21, 39, 58, 73)] <- NA_real_

fit_boundary_proportion <- drmTMB(
  bf(growth ~ z + mi(cover), sigma ~ 1),
  family = gaussian(),
  data = boundary_prop_data,
  impute = list(
    cover = impute_model(cover ~ z, family = zero_one_beta())
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_boundary_proportion, "mi_cover")
#> (Intercept)           z 
#>  -0.1497048   0.7815367
coef(fit_boundary_proportion, "sigma_mi_cover")
#>      cover 
#> 0.04055868
coef(fit_boundary_proportion, "zoi_mi_cover")
#>     cover 
#> 0.1154065
coef(fit_boundary_proportion, "coi_mi_cover")
#>     cover 
#> 0.4443852
imputed(fit_boundary_proportion)
#>   variable original_row model_row observed  estimate std_error
#> 1    cover            8         8    FALSE 0.2145062        NA
#> 2    cover           21        21    FALSE 0.3024846        NA
#> 3    cover           39        39    FALSE 0.4935249        NA
#> 4    cover           58        58    FALSE 0.6006695        NA
#> 5    cover           73        73    FALSE 0.7472980        NA
#>                        source uncertainty_status
#> 1 conditional_quadrature_mean   sdreport_skipped
#> 2 conditional_quadrature_mean   sdreport_skipped
#> 3 conditional_quadrature_mean   sdreport_skipped
#> 4 conditional_quadrature_mean   sdreport_skipped
#> 5 conditional_quadrature_mean   sdreport_skipped

For a missing boundary-proportion value, drmTMB integrates over exact zero mass, exact one mass, and interior beta quadrature under the fitted predictor model and the Gaussian response likelihood. imputed(fit_boundary_proportion) reports the fitted conditional quadrature mean for each missing cover value.

Denominator-aware proportion missing predictors

Use a beta-binomial predictor model when the covariate used in the response model is a proportion, but the observed data are success counts out of known trial totals. The main formula still uses the proportion through mi(cover). The impute_model() formula uses the success-count column on the left-hand side and trials for the known denominator:

set.seed(20260542)
n <- 70
denom_prop_data <- data.frame(
  z = seq(-1.7, 1.7, length.out = n),
  trials = rep(8:16, length.out = n)
)
p_cover <- plogis(-0.20 + 0.80 * denom_prop_data$z)
denom_prop_data$success_full <- qbinom(
  ppoints(n),
  size = denom_prop_data$trials,
  prob = p_cover
)
denom_prop_data$cover_full <-
  denom_prop_data$success_full / denom_prop_data$trials
denom_prop_data$growth <- 0.25 + 1.10 * denom_prop_data$cover_full -
  0.22 * denom_prop_data$z + rnorm(n, sd = 0.10)
denom_prop_data$success <- denom_prop_data$success_full
denom_prop_data$cover <- denom_prop_data$cover_full
denom_prop_data$success[c(8, 21, 39, 58)] <- NA_real_
denom_prop_data$cover[c(8, 21, 39, 58)] <- NA_real_

fit_denominator_proportion <- drmTMB(
  bf(growth ~ z + mi(cover), sigma ~ 1),
  family = gaussian(),
  data = denom_prop_data,
  impute = list(
    cover = impute_model(
      success ~ z,
      family = beta_binomial(),
      trials = trials
    )
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)
#> Warning: `drmTMB()`: optimizer reported non-convergence (code 1: false convergence (8)).
#> ℹ Treat the estimates and standard errors with caution. The optimizer escalates
#>   its preset ladder automatically, so inspect `fit$optimizer_attempts` and run
#>   `check_drm()` to diagnose; consider rescaling the data, within-group
#>   replication, or a penalized/MAP fit.

coef(fit_denominator_proportion, "mi_cover")
#> (Intercept)           z 
#>  -0.2295829   1.5665969
coef(fit_denominator_proportion, "sigma_mi_cover")
#>        cover 
#> 0.0001725686
imputed(fit_denominator_proportion)
#>   variable original_row model_row observed   estimate std_error
#> 1    cover            8         8    FALSE 0.09237150        NA
#> 2    cover           21        21    FALSE 0.08036363        NA
#> 3    cover           39        39    FALSE 0.55044334        NA
#> 4    cover           58        58    FALSE 0.77509226        NA
#>                        source uncertainty_status
#> 1 conditional_proportion_mean   sdreport_skipped
#> 2 conditional_proportion_mean   sdreport_skipped
#> 3 conditional_proportion_mean   sdreport_skipped
#> 4 conditional_proportion_mean   sdreport_skipped

For a missing success count, drmTMB sums over all possible successes from zero to the known trial total. imputed(fit_denominator_proportion) reports the fitted conditional proportion mean, not a filled-in integer count.

Count missing predictors

Count predictors use a Poisson predictor model when the mean-variance relationship is close to Poisson. This route is for non-negative integer counts, not continuous biomass:

set.seed(20260538)
n <- 76
count_data <- data.frame(
  z = seq(-1.6, 1.7, length.out = n)
)
lambda_count <- exp(0.35 + 0.55 * count_data$z)
count_data$abundance_full <- pmax(
  0,
  round(lambda_count + sqrt(lambda_count) * sin(seq_len(n) / 3))
)
count_data$growth <- 0.20 + 0.42 * count_data$abundance_full -
  0.25 * count_data$z + rnorm(n, sd = 0.10)
count_data$abundance <- count_data$abundance_full
count_data$abundance[c(7, 18, 33, 49, 68)] <- NA_real_

fit_count <- drmTMB(
  bf(growth ~ z + mi(abundance), sigma ~ 1),
  family = gaussian(),
  data = count_data,
  impute = list(
    abundance = impute_model(abundance ~ z, family = poisson())
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_count, "mi_abundance")
#> (Intercept)           z 
#>   0.3748932   0.4665136
imputed(fit_count)
#>    variable original_row model_row observed    estimate std_error
#> 1 abundance            7         7    FALSE 1.039293126        NA
#> 2 abundance           18        18    FALSE 0.978534027        NA
#> 3 abundance           33        33    FALSE 0.005923382        NA
#> 4 abundance           49        49    FALSE 0.980589533        NA
#> 5 abundance           68        68    FALSE 2.001062614        NA
#>                       source uncertainty_status
#> 1 conditional_expected_count   sdreport_skipped
#> 2 conditional_expected_count   sdreport_skipped
#> 3 conditional_expected_count   sdreport_skipped
#> 4 conditional_expected_count   sdreport_skipped
#> 5 conditional_expected_count   sdreport_skipped

For a missing abundance value, drmTMB sums over a deterministic count support under the fitted Poisson predictor model and the Gaussian response likelihood. imputed(fit_count) reports the fitted conditional expected count for each missing abundance value.

Use nbinom2() instead when the missing count predictor is overdispersed:

fit_count_nb <- drmTMB(
  bf(growth ~ z + mi(abundance), sigma ~ 1),
  family = gaussian(),
  data = count_data,
  impute = list(
    abundance = impute_model(abundance ~ z, family = nbinom2())
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_count_nb, "mi_abundance")
#> (Intercept)           z 
#>   0.3748932   0.4665137
coef(fit_count_nb, "sigma_mi_abundance")
#>    abundance 
#> 1.143018e-05
imputed(fit_count_nb)
#>    variable original_row model_row observed    estimate std_error
#> 1 abundance            7         7    FALSE 1.039293135        NA
#> 2 abundance           18        18    FALSE 0.978533969        NA
#> 3 abundance           33        33    FALSE 0.005923377        NA
#> 4 abundance           49        49    FALSE 0.980589475        NA
#> 5 abundance           68        68    FALSE 2.001062615        NA
#>                       source uncertainty_status
#> 1 conditional_expected_count   sdreport_skipped
#> 2 conditional_expected_count   sdreport_skipped
#> 3 conditional_expected_count   sdreport_skipped
#> 4 conditional_expected_count   sdreport_skipped
#> 5 conditional_expected_count   sdreport_skipped

The negative-binomial route uses the same NB2 parameterization as a count response model, $Var(x_i) = \mu_i + \sigma_x^2 \mu_i^2$ . Missing counts are summed over a deterministic count support under the fitted NB2 predictor model.

Use truncated_nbinom2() when the predictor is a positive count by design, such as days until first event among records where the event occurred:

positive_count_data <- count_data
positive_count_data$abundance_full <- count_data$abundance_full + 1
positive_count_data$abundance <- positive_count_data$abundance_full
positive_count_data$abundance[c(7, 18, 33, 49, 68)] <- NA_real_
positive_count_data$growth <- 0.20 + 0.42 * positive_count_data$abundance_full -
  0.25 * positive_count_data$z + rnorm(nrow(positive_count_data), sd = 0.10)

fit_count_trunc_nb <- drmTMB(
  bf(growth ~ z + mi(abundance), sigma ~ 1),
  family = gaussian(),
  data = positive_count_data,
  impute = list(
    abundance = impute_model(abundance ~ z, family = truncated_nbinom2())
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_count_trunc_nb, "mi_abundance")
#> (Intercept)           z 
#>   0.7737085   0.3838097
coef(fit_count_trunc_nb, "sigma_mi_abundance")
#>    abundance 
#> 1.031833e-05
imputed(fit_count_trunc_nb)
#>    variable original_row model_row observed estimate std_error
#> 1 abundance            7         7    FALSE 2.000065        NA
#> 2 abundance           18        18    FALSE 1.999675        NA
#> 3 abundance           33        33    FALSE 1.000162        NA
#> 4 abundance           49        49    FALSE 1.999248        NA
#> 5 abundance           68        68    FALSE 2.999541        NA
#>                       source uncertainty_status
#> 1 conditional_expected_count   sdreport_skipped
#> 2 conditional_expected_count   sdreport_skipped
#> 3 conditional_expected_count   sdreport_skipped
#> 4 conditional_expected_count   sdreport_skipped
#> 5 conditional_expected_count   sdreport_skipped

The zero-truncated route excludes zero from the missing-count support. Use nbinom2() instead when zero is a possible observed predictor value.

Positive continuous missing predictors

Positive continuous predictors can use a lognormal or Gamma predictor model. These routes are for strictly positive biomass, concentration, duration, or size predictors, not exact zeros or semi-continuous variables:

set.seed(20260539)
n <- 76
positive_data <- data.frame(
  z = seq(-1.6, 1.6, length.out = n)
)
positive_data$biomass_full <- exp(
  0.15 + 0.55 * positive_data$z + 0.10 * sin(seq_len(n) / 4)
)
positive_data$growth <- 0.40 + 0.85 * positive_data$biomass_full -
  0.25 * positive_data$z + rnorm(n, sd = 0.10)
positive_data$biomass <- positive_data$biomass_full
positive_data$biomass[c(7, 18, 33, 49, 65)] <- NA_real_

fit_positive <- drmTMB(
  bf(growth ~ z + mi(biomass), sigma ~ 1),
  family = gaussian(),
  data = positive_data,
  impute = list(
    biomass = impute_model(biomass ~ z, family = lognormal())
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_positive, "mi_biomass")
#> (Intercept)           z 
#>   0.1497577   0.5317314
coef(fit_positive, "sigma_mi_biomass")
#>    biomass 
#> 0.06855065
imputed(fit_positive)
#>   variable original_row model_row observed  estimate std_error
#> 1  biomass            7         7    FALSE 0.5638418        NA
#> 2  biomass           18        18    FALSE 0.6975571        NA
#> 3  biomass           33        33    FALSE 1.1319169        NA
#> 4  biomass           49        49    FALSE 1.4092666        NA
#> 5  biomass           65        65    FALSE 2.0610705        NA
#>                        source uncertainty_status
#> 1 conditional_quadrature_mean   sdreport_skipped
#> 2 conditional_quadrature_mean   sdreport_skipped
#> 3 conditional_quadrature_mean   sdreport_skipped
#> 4 conditional_quadrature_mean   sdreport_skipped
#> 5 conditional_quadrature_mean   sdreport_skipped

For a missing biomass value, drmTMB integrates over log-scale predictor states under the fitted lognormal predictor model and the Gaussian response likelihood. imputed(fit_positive) reports the fitted conditional quadrature mean for each missing biomass value.

Use a Gamma predictor model when the positive covariate is better described by a mean-CV model on the original scale:

fit_positive_gamma <- drmTMB(
  bf(growth ~ z + mi(biomass), sigma ~ 1),
  family = gaussian(),
  data = positive_data,
  impute = list(
    biomass = impute_model(biomass ~ z, family = Gamma(link = "log"))
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_positive_gamma, "mi_biomass")
#> (Intercept)           z 
#>   0.1535083   0.5306886
coef(fit_positive_gamma, "sigma_mi_biomass")
#>   biomass 
#> 0.1206742
imputed(fit_positive_gamma)
#>   variable original_row model_row observed  estimate std_error
#> 1  biomass            7         7    FALSE 0.5358504        NA
#> 2  biomass           18        18    FALSE 0.6547371        NA
#> 3  biomass           33        33    FALSE 0.9969145        NA
#> 4  biomass           49        49    FALSE 1.3790015        NA
#> 5  biomass           65        65    FALSE 2.0517131        NA
#>                        source uncertainty_status
#> 1 conditional_quadrature_mean   sdreport_skipped
#> 2 conditional_quadrature_mean   sdreport_skipped
#> 3 conditional_quadrature_mean   sdreport_skipped
#> 4 conditional_quadrature_mean   sdreport_skipped
#> 5 conditional_quadrature_mean   sdreport_skipped

For a missing biomass value under the Gamma route, drmTMB integrates over positive predictor states under the fitted Gamma mean-CV predictor model and the Gaussian response likelihood.

Semi-continuous missing predictors with exact zeros

Use a Tweedie predictor model when the missing covariate is non-negative and can have exact zeros as part of the same semi-continuous process. The first route fits the predictor mean and scale, fixes the predictor-model Tweedie power at 1.5, and integrates missing predictor values over the exact zero mass plus a deterministic positive quadrature grid:

set.seed(20260540)
n <- 82
semi_data <- data.frame(
  z = seq(-1.7, 1.8, length.out = n)
)
semi_data$biomass_full <- ifelse(
  seq_len(n) %% 6 %in% c(0, 1),
  0,
  exp(0.10 + 0.45 * semi_data$z) *
    (1 + 0.16 * sin(seq_len(n) / 4))
)
semi_data$growth <- 0.25 + 0.62 * semi_data$biomass_full -
  0.18 * semi_data$z + rnorm(n, sd = 0.10)
semi_data$biomass <- semi_data$biomass_full
semi_data$biomass[c(7, 18, 31, 44, 58, 76)] <- NA_real_

fit_tweedie_predictor <- drmTMB(
  bf(growth ~ z + mi(biomass), sigma ~ 1),
  family = gaussian(),
  data = semi_data,
  impute = list(
    biomass = impute_model(biomass ~ z, family = tweedie())
  ),
  missing = miss_control(predictor = "model"),
  control = drm_control(se = FALSE)
)

coef(fit_tweedie_predictor, "mi_biomass")
#> (Intercept)           z 
#>  -0.2861113   0.4578044
coef(fit_tweedie_predictor, "sigma_mi_biomass")
#> biomass 
#> 1.15473
imputed(fit_tweedie_predictor)
#>   variable original_row model_row observed   estimate std_error
#> 1  biomass            7         7    FALSE 0.02112717        NA
#> 2  biomass           18        18    FALSE 0.04817500        NA
#> 3  biomass           31        31    FALSE 0.05043276        NA
#> 4  biomass           44        44    FALSE 1.06593905        NA
#> 5  biomass           58        58    FALSE 1.85991131        NA
#> 6  biomass           76        76    FALSE 2.17920872        NA
#>                        source uncertainty_status
#> 1 conditional_quadrature_mean   sdreport_skipped
#> 2 conditional_quadrature_mean   sdreport_skipped
#> 3 conditional_quadrature_mean   sdreport_skipped
#> 4 conditional_quadrature_mean   sdreport_skipped
#> 5 conditional_quadrature_mean   sdreport_skipped
#> 6 conditional_quadrature_mean   sdreport_skipped

imputed(fit_tweedie_predictor) reports the fitted conditional quadrature mean for each missing biomass value. This is still a likelihood-based summary, not a draw from a posterior distribution or a multiple-imputation data set.

Grouped and structured predictor models

The same one-predictor route can include one grouped random intercept in the predictor model:

drmTMB(
  bf(growth ~ moisture + mi(body_mass), sigma ~ 1),
  family = gaussian(),
  data = predictor_data,
  impute = list(body_mass = body_mass ~ moisture + (1 | site)),
  missing = miss_control(predictor = "model")
)

Use this when missingness in the predictor is tied to a grouping level and a separate covariate model is scientifically defensible. The grouped effect is a covariate-model random intercept, not a response-model random effect.

Structured covariate models are also explicit. For example, a missing species-level trait can use a known relatedness matrix:

drmTMB(
  bf(growth ~ moisture + mi(body_mass), sigma ~ 1),
  family = gaussian(),
  data = predictor_data,
  impute = list(body_mass = body_mass ~ moisture + relmat(1 | line, Q = Q)),
  missing = miss_control(predictor = "model")
)

The current structured route supports one intercept-only phylo(), coordinate spatial(), animal(), or relmat() covariate field. It does not automatically inherit structure from the response model. That explicitness is intentional: a model for a missing trait is a scientific model for the trait, not just a technical patch for a response regression.

Inspect the missing-data metadata

Every fitted route records a fit$missing_data object. The most useful fields are:

Field	Meaning
`version`	The fitted implementation slice, such as `"MD1"`, `"MD2"`, `"MD3a"`, `"MD3b"`, `"MD4"`, `"MD6a"`, `"MD6b"`, `"MD6c"`, `"MD7a"`, `"MD7b"`, `"MD7c"`, `"MD7d"`, `"MD7e"`, `"MD7f"`, `"MD8a"`, `"MD8b"`, `"MD8c"`, `"MD9a"` (Poisson-response `mi()`), `"MD-binomial-mi"`, `"MD-nbinom2-mi"`, or `"MD-beta-mi"`.
`response_policy`	Whether response rows were dropped or included.
`predictor_policy`	Whether missing predictors fail or are modelled.
`original_row`	Row indices in the user’s original data.
`model_row`	Row indices in the retained model frame.
`observed_y`	Logical response mask for univariate fits, or a two-column mask for bivariate fits.
`counts`	Retained rows, observed responses, missing responses, and likelihood-contributing rows.
`predictors`	Metadata and fitted values for explicitly modelled `mi()` predictors.

For ordinary user work, prefer nobs(), fitted(), residuals(), and imputed() over hand-parsing internals. The metadata is there so row accounting stays auditable.

Current boundaries

The current implementation deliberately keeps several features out of scope:

EM, profile, and REML missing-data engines are not implemented.
Multiple missing predictors are not implemented.
Missing predictors in non-Gaussian response models are implemented for the poisson(), binomial(), nbinom2(), and beta() responses, each with one binary (Bernoulli/logit) missing predictor. Other non-Gaussian responses, and non-binary missing predictors on non-Gaussian responses, remain later work.
Hurdle count predictor models, transformed mi() terms, interactions inside mi(), splines inside mi(), and predictor-model random slopes are not implemented.
Grouped and structured binary, ordered, unordered, beta/proportion, zero-one beta, Poisson/NB2, lognormal, Gamma, or Tweedie predictor models are not implemented.
Dense known sampling-covariance matrices with bivariate partial-response rows are not implemented.
imputed() does not summarize missing responses.
Simulated imputations, posterior means, credible intervals, and multiple-imputation pooling are not implemented.
Measurement-error models are a later modelling lane, not part of this missing-data interface.

When you need one of those routes, treat the current error as a design boundary rather than a failed fit. For now, fit the supported likelihood-based route, use complete-case fitting, or handle the unsupported missingness outside drmTMB with a method whose assumptions you can state.