## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 4.1,
  dpi = 144
)
library(drmTMB)

## -----------------------------------------------------------------------------
miss_control()

## -----------------------------------------------------------------------------
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(fit_response)
length(fitted(fit_response))
sum(is.na(residuals(fit_response)))

## -----------------------------------------------------------------------------
head(data.frame(
  original_row = seq_along(fitted(fit_response)),
  observed_y = !is.na(residuals(fit_response)),
  fitted = fitted(fit_response),
  residual = residuals(fit_response)
), 10)

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

pair_residuals <- residuals(fit_pair)
table(
  activity_observed = !is.na(pair_residuals[, 1]),
  boldness_observed = !is.na(pair_residuals[, 2])
)

## -----------------------------------------------------------------------------
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")
imputed(fit_predictor)

## -----------------------------------------------------------------------------
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")
imputed(fit_binary)

## -----------------------------------------------------------------------------
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")
coef(fit_poisson_binary, "mi_treatment")
imputed(fit_poisson_binary)

## -----------------------------------------------------------------------------
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")
coef(fit_nbinom2_binary, "mi_treatment")

## -----------------------------------------------------------------------------
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")
coef(fit_beta_binary, "mi_treatment")

## -----------------------------------------------------------------------------
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")
imputed(fit_ordered)

## -----------------------------------------------------------------------------
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")
imputed(fit_categorical)

## -----------------------------------------------------------------------------
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")
coef(fit_proportion, "sigma_mi_cover")
imputed(fit_proportion)

## -----------------------------------------------------------------------------
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")
coef(fit_boundary_proportion, "sigma_mi_cover")
coef(fit_boundary_proportion, "zoi_mi_cover")
coef(fit_boundary_proportion, "coi_mi_cover")
imputed(fit_boundary_proportion)

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

coef(fit_denominator_proportion, "mi_cover")
coef(fit_denominator_proportion, "sigma_mi_cover")
imputed(fit_denominator_proportion)

## -----------------------------------------------------------------------------
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")
imputed(fit_count)

## -----------------------------------------------------------------------------
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")
coef(fit_count_nb, "sigma_mi_abundance")
imputed(fit_count_nb)

## -----------------------------------------------------------------------------
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")
coef(fit_count_trunc_nb, "sigma_mi_abundance")
imputed(fit_count_trunc_nb)

## -----------------------------------------------------------------------------
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")
coef(fit_positive, "sigma_mi_biomass")
imputed(fit_positive)

## -----------------------------------------------------------------------------
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")
coef(fit_positive_gamma, "sigma_mi_biomass")
imputed(fit_positive_gamma)

## -----------------------------------------------------------------------------
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")
coef(fit_tweedie_predictor, "sigma_mi_biomass")
imputed(fit_tweedie_predictor)