Skip to contents

Model selection with AIC and BIC

Source: vignettes/model-selection.Rmd
model-selection.Rmd

Model selection comes after candidate models have been chosen, fitted, and checked with the post-fit workflow in Checking and using fitted models. In drmTMB, AIC and BIC can compare fitted distributional-regression models, but they do not decide whether a candidate set is biologically sensible or whether a weak fit should be interpreted.

AIC and BIC both use the fitted log likelihood and the number of estimated parameters:

AIC = -2 * logLik + 2 * k
BIC = -2 * logLik + log(n) * k

k is the model degrees of freedom and n is the number of observations used by the fitted model. BIC penalizes extra parameters more strongly than AIC once n > exp(2). Compare models only when they were fitted to the same response scale and the same analysis rows. Do not compare a model fit to raw counts with one fit to transformed counts, or two fits that silently dropped different missing rows.

For Gaussian mixed models, use ML (REML = FALSE) when AIC or BIC compares different fixed-effect formulas. REML is useful for Gaussian variance-component estimation inside a fixed mean structure, but the restricted likelihood depends on the fixed-effect design that was integrated out.

library(drmTMB)
#> 
#> Attaching package: 'drmTMB'
#> The following object is masked from 'package:base':
#> 
#>     beta

criterion_converged <- function(fit) {
  isTRUE(fit$opt$convergence == 0)
}

criterion_table <- function(...) {
  models <- list(...)
  out <- data.frame(
    model = names(models),
    AIC = vapply(models, stats::AIC, numeric(1)),
    BIC = vapply(models, stats::BIC, numeric(1)),
    converged = vapply(models, criterion_converged, logical(1)),
    stringsAsFactors = FALSE
  )
  out$delta_AIC <- out$AIC - min(out$AIC)
  out$delta_BIC <- out$BIC - min(out$BIC)
  out
}

Gaussian REML After Mean-Model Selection

For Gaussian mixed models, first compare different fixed-effect formulas with ML. This example includes a random intercept for repeated measurements and a noise predictor z that was not used to generate the response.

set.seed(2404)
n_id <- 12L
n_each <- 6L
id <- factor(rep(seq_len(n_id), each = n_each))
x <- rnorm(n_id * n_each)
z <- rnorm(n_id * n_each)
u <- rnorm(n_id, sd = 0.65)
mixed_dat <- data.frame(
  id = id,
  x = x,
  z = z,
  y = 0.2 + 0.7 * x + u[id] + rnorm(n_id * n_each, sd = 0.45)
)

fit_mixed_x_ml <- drmTMB(
  bf(y ~ x + (1 | id), sigma ~ 1),
  family = gaussian(),
  data = mixed_dat
)

fit_mixed_x_z_ml <- drmTMB(
  bf(y ~ x + z + (1 | id), sigma ~ 1),
  family = gaussian(),
  data = mixed_dat
)

criterion_table(`y ~ x` = fit_mixed_x_ml, `y ~ x + z` = fit_mixed_x_z_ml)
#>               model      AIC      BIC converged delta_AIC delta_BIC
#> y ~ x         y ~ x 111.6812 120.7878      TRUE  0.000000  0.000000
#> y ~ x + z y ~ x + z 113.6786 125.0620      TRUE  1.997482  4.274148

Once the fixed-effect mean structure is chosen, refit that structure with REML = TRUE when the target is Gaussian variance-component estimation. The first REML slice covers ordinary univariate Gaussian random intercepts and slopes in mu with intercept-only sigma.

fit_mixed_x_reml <- drmTMB(
  bf(y ~ x + (1 | id), sigma ~ 1),
  family = gaussian(),
  data = mixed_dat,
  REML = TRUE
)

ll_reml <- logLik(fit_mixed_x_reml)
data.frame(
  estimator = fit_mixed_x_reml$estimator,
  restricted_logLik = as.numeric(ll_reml),
  df = attr(ll_reml, "df"),
  residual_sigma = unname(sigma(fit_mixed_x_reml)[1]),
  id_sd = unname(fit_mixed_x_reml$sdpars$mu["(1 | id)"])
)
#>   estimator restricted_logLik df residual_sigma     id_sd
#> 1      REML         -54.54347  4      0.3840255 0.7881384

When lme4 is installed, the restricted log likelihood for this overlap model matches lme4::lmer(..., REML = TRUE).

if (requireNamespace("lme4", quietly = TRUE)) {
  fit_lme4_reml <- lme4::lmer(
    y ~ x + (1 | id),
    data = mixed_dat,
    REML = TRUE
  )
  data.frame(
    engine = c("drmTMB", "lme4"),
    restricted_logLik = c(
      as.numeric(logLik(fit_mixed_x_reml)),
      as.numeric(logLik(fit_lme4_reml))
    )
  )
}
#>   engine restricted_logLik
#> 1 drmTMB         -54.54347
#> 2   lme4         -54.54347

Tail Assumption

A robust Student-t model is useful when the biological process or measurement process can produce unusually large residuals. Here the data are generated with heavy-tailed residuals, then fitted with a Gaussian model and a Student-t model that uses the same mu and sigma formulas.

set.seed(2401)
n <- 220
x <- rnorm(n)
tail_dat <- data.frame(x = x)
tail_dat$y <- 0.2 + 0.7 * x + exp(-0.25) * rt(n, df = 4)

fit_tail_gaussian <- drmTMB(
  bf(y ~ x, sigma ~ 1),
  family = gaussian(),
  data = tail_dat
)

fit_tail_student <- drmTMB(
  bf(y ~ x, sigma ~ 1, nu ~ 1),
  family = student(),
  data = tail_dat
)

criterion_table(
  Gaussian = fit_tail_gaussian,
  `Student-t` = fit_tail_student
)
#>               model      AIC      BIC converged delta_AIC delta_BIC
#> Gaussian   Gaussian 671.0859 681.2668      TRUE  27.05801  23.66438
#> Student-t Student-t 644.0279 657.6024      TRUE   0.00000   0.00000

The lower AIC/BIC model is the better criterion fit inside this two-model candidate set. That result is not a licence to ignore diagnostics: run check_drm() and inspect whether the fitted Student-t nu parameter is near the lower bound or so large that the Student-t model is effectively Gaussian.

check_drm(fit_tail_student)
#> <drm_check: 13 checks>
#> ok: 13; notes: 0; warnings: 0; errors: 0
#>                      check status
#>      optimizer_convergence     ok
#>           optimizer_budget     ok
#>           finite_objective     ok
#>      logsigma_clamp_active     ok
#>             fixed_gradient     ok
#>            sdreport_status     ok
#>  hessian_positive_definite     ok
#>     standard_errors_finite     ok
#>   standard_errors_inflated     ok
#>               dropped_rows     ok
#>             positive_scale     ok
#>                 student_nu     ok
#>   fixed_effect_design_size     ok
#>                                                                                  value
#>                                                                                      0
#>                                                iterations=21; function=32; gradient=22
#>                                                                                  318.0
#>                                                                                   <NA>
#>                                                  max=0.000001686; component=beta_mu[2]
#>                                                                                     ok
#>                                                                                   TRUE
#>                                                                 range=[0.06345,0.5697]
#>                                         n_inflated=0; max_se=0.5697; median_se=0.07471
#>                                                                    nobs=220; dropped=0
#>                                                                             min=0.7981
#>                                                                    range=[4.164,4.164]
#>  total_mb=0.04914; max_cols=2; largest=mu; largest_class=matrix; largest_density=1.000
#>                                                                              message
#>                                                        nlminb convergence code is 0.
#>  Optimizer evaluation counts recorded; no eval.max or iter.max control was supplied.
#>                                             Objective and log-likelihood are finite.
#>                                   The log(sigma) clamp is not active at the optimum.
#>        Maximum absolute fixed gradient is <= 0.001; largest component is beta_mu[2].
#>                                              TMB::sdreport() completed successfully.
#>                                        sdreport reports a positive-definite Hessian.
#>                                         All fixed-effect standard errors are finite.
#>                   No fixed-effect standard error is inflated relative to the others.
#>                   No rows were dropped by model-frame or known-covariance filtering.
#>                                     All fitted scale values are finite and positive.
#>               All fitted Student-t nu values are finite and above the boundary at 2.
#>                          Dense fixed-effect design matrices are modest for this fit.
coef(fit_tail_student, "nu")
#> (Intercept) 
#>   0.7721585

Structural Zeros

For count responses, the extra parameter in a zero-inflated model has a specific interpretation: it is the probability of a separate structural-zero process. This example generates NB2 counts with structural zeros and compares ordinary NB2 against ZINB2.

set.seed(2402)
n <- 260
x <- rnorm(n)
mu <- exp(log(2.3) + 0.5 * x)
sigma <- 0.65
zi <- plogis(-0.8)

count <- rnbinom(n, size = 1 / sigma^2, mu = mu)
structural_zero <- runif(n) < zi
count[structural_zero] <- 0L
count_dat <- data.frame(count = count, x = x)

fit_nb2 <- drmTMB(
  bf(count ~ x, sigma ~ 1),
  family = nbinom2(),
  data = count_dat
)

fit_zinb2 <- drmTMB(
  bf(count ~ x, sigma ~ 1, zi ~ 1),
  family = nbinom2(),
  data = count_dat
)

criterion_table(NB2 = fit_nb2, ZINB2 = fit_zinb2)
#>       model      AIC      BIC converged delta_AIC delta_BIC
#> NB2     NB2 850.0986 860.7807      TRUE  6.434382    2.8737
#> ZINB2 ZINB2 843.6643 857.9070      TRUE  0.000000    0.0000

If the criterion difference is small, keep both explanations in view. NB2 can absorb some extra zeros by increasing overdispersion, while ZINB2 separates structural zeros from count variation. Use zi only when a structural-zero process is plausible for the response and sampling design.

Scale Formula

Model selection is not only family selection. It can also ask whether a distributional parameter needs a predictor. In this Gaussian example the mean changes with x, and the residual scale also changes with x.

set.seed(2403)
n <- 220
x <- rnorm(n)
sigma <- exp(-0.45 + 0.55 * x)
scale_dat <- data.frame(
  x = x,
  y = 0.3 + 0.55 * x + rnorm(n, sd = sigma)
)

fit_sigma_constant <- drmTMB(
  bf(y ~ x, sigma ~ 1),
  family = gaussian(),
  data = scale_dat
)

fit_sigma_x <- drmTMB(
  bf(y ~ x, sigma ~ x),
  family = gaussian(),
  data = scale_dat
)

criterion_table(`sigma ~ 1` = fit_sigma_constant, `sigma ~ x` = fit_sigma_x)
#>               model      AIC      BIC converged delta_AIC delta_BIC
#> sigma ~ 1 sigma ~ 1 634.2400 644.4209      TRUE  181.3262  177.9326
#> sigma ~ x sigma ~ x 452.9138 466.4883      TRUE    0.0000    0.0000

The sigma ~ x model estimates a log-scale slope. A positive slope means the modelled residual scale increases as x increases. After selecting this candidate, interpret the sigma coefficient on the ratio scale:

exp(coef(fit_sigma_x, "sigma")["x"])
#>        x 
#> 1.859795

A 200-Replicate Article-Support Simulation

The package includes a seeded article-support simulation. It uses 200 replicates per scenario, which is enough to make the table more stable than a plumbing smoke run while still being small enough to ship as documentation evidence. It is not a formal operating-characteristic study over sample sizes, effect sizes, and candidate-set misspecification.

summary_path <- system.file(
  "sim/reports/model-selection-article-summary.csv",
  package = "drmTMB"
)
if (!nzchar(summary_path)) {
  candidates <- c(
    "../inst/sim/reports/model-selection-article-summary.csv",
    "inst/sim/reports/model-selection-article-summary.csv"
  )
  summary_path <- candidates[file.exists(candidates)][1L]
}

model_selection_article <- read.csv(summary_path)
display_article <- model_selection_article[, c(
  "scenario",
  "selection_target",
  "n_replicate",
  "aic_truth_selection_rate",
  "aic_truth_selection_mcse",
  "bic_truth_selection_rate",
  "bic_truth_selection_mcse",
  "candidate_convergence_rate",
  "candidate_pdHess_rate",
  "candidate_warning_rate"
)]
names(display_article) <- c(
  "scenario",
  "target",
  "replicates",
  "AIC selected target",
  "AIC MCSE",
  "BIC selected target",
  "BIC MCSE",
  "candidate convergence",
  "candidate pdHess",
  "candidate warning"
)
knitr::kable(display_article, digits = 3)
scenario target replicates AIC selected target AIC MCSE BIC selected target BIC MCSE candidate convergence candidate pdHess candidate warning
constant_sigma sigma ~ 1 200 0.825 0.027 0.975 0.011 1.000 1.000 0.000
extra_zeros ZINB2 200 0.920 0.019 0.730 0.031 1.000 1.000 0.000
heavy_tail Student-t 200 0.840 0.026 0.665 0.033 0.998 0.995 0.005
nb2_counts NB2 200 0.935 0.017 0.990 0.007 1.000 1.000 0.000
normal_tail Gaussian 200 0.945 0.016 0.980 0.010 0.760 0.818 0.180
sigma_signal sigma ~ x 200 1.000 0.000 1.000 0.000 1.000 1.000 0.000

The normal_tail row is deliberately useful: the Gaussian model is selected, but the unnecessary Student-t candidate often carries warning or weak-Hessian status because the extra tail parameter is unnecessary for normal data. The heavy_tail and extra_zeros rows show the criterion tradeoff: AIC more often keeps the extra Student-t or zero-inflation parameter, while BIC often prefers the simpler candidate under this sample size because its penalty is stronger. That is the workflow lesson. Selection criteria answer “which candidate has the smaller penalized likelihood score?”, while diagnostics answer “is this candidate fit stable enough to use?”.

Practical Checklist

Use this sequence when comparing fitted drmTMB models:

  1. Define the candidate set before looking at AIC/BIC.
  2. Fit all candidates to the same response, row set, offsets, and weights. For Gaussian mixed models with different fixed-effect formulas, keep REML = FALSE.
  3. Run check_drm() and keep convergence, Hessian status, warnings, and boundary diagnostics beside the criterion table.
  4. Exclude errored fits from selection and treat nonconverged or weak-Hessian fits as diagnostic findings, not as ordinary winners.
  5. Report AIC and BIC differences, not only the winning model.
  6. Interpret the selected model’s mu, sigma, nu, zi, or rho12 parameters in scientific units or ratios.

When AIC and BIC disagree, describe the tradeoff. AIC is more willing to keep extra parameters when they improve fit; BIC asks for stronger evidence as sample size grows. In applied drmTMB work, that disagreement is often a sign to inspect predictions, residuals, and the scientific meaning of the extra distributional parameter rather than to declare one criterion universally correct.