## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")

## ----setup--------------------------------------------------------------------
library(drmTMB)

criterion_converged <- function(fit) {
  is_converged(fit)
}

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-example----------------------------------------------------
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)

## ----gaussian-reml-fit--------------------------------------------------------
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)
mixed_parameters <- summary(fit_mixed_x_reml)$parameters
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 = mixed_parameters[
    mixed_parameters$parm == "sd:mu:(1 | id)", "estimate"
  ]
)

## ----gaussian-reml-lme4-------------------------------------------------------
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))
    )
  )
}

## ----tail-example-------------------------------------------------------------
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
)

## ----tail-check---------------------------------------------------------------
check_drm(fit_tail_student)
coef(fit_tail_student, "nu")

## ----count-example------------------------------------------------------------
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)

## ----scale-example------------------------------------------------------------
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)

## ----scale-ratio--------------------------------------------------------------
exp(coef(fit_sigma_x, "sigma")["x"])

## ----article-summary----------------------------------------------------------
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)