# Researchers plan to assign diets to 16 pens, each holding six grower pigs.
# The numbers below are assumptions for a new study.
library("dplyr")
library("simglm")
library("lme4")
library("performance")

settings <- list(
  formula = gain ~ diet + (1 | pen),
  fixed = list(
    diet = list(
      var_type = "factor",
      levels = c("Control", "Supplement"),
      var_level = 2,
      force_equal = TRUE
    )
  ),
  sample_size = list(level1 = 6, level2 = 16),
  reg_weights = c(700, 35),
  randomeffect = list(int_pen = list(variance = 40^2, var_level = 2)),
  error = list(variance = 35^2),
  model_fit = list(model_function = lme4::lmer, REML = TRUE),
  power = list(alpha = 0.05, dist = "qnorm", direction = "two-tailed"),
  replications = 100
)
# level1 counts pigs per pen; level2 counts pens receiving diets.
# var_level = 2 gives every animal in a pen the same assigned diet.
# force_equal = TRUE assigns half the pens to each diet.
# 700 g/day is the assumed Control mean; 35 g/day is the worthwhile change.
# 40 and 35 g/day are shared pen and individual SDs; variance is SD squared.
# This simplified future experiment has no room blocks or section treatments.
# REML estimates Gaussian group and residual variation for the planned analysis.

set.seed(821029)
one_trial <- simglm(settings)
print(one_trial |> select(pen, diet, gain) |> head(3))
print(one_trial |> distinct(pen, diet) |> count(diet))

# Preserve the original plan order: 16/40, 32/40, 16/60, 32/60.
plans <- tidyr::expand_grid(pen_sd = c(40, 60), pens = c(16, 32)) |>
  select(pens, pen_sd) |>
  mutate(
    attempted = settings$replications,
    usable = 0L,
    detected = 0L,
    power_among_usable = NA_real_
  )
for (i in seq_len(nrow(plans))) {
  settings$sample_size$level2 <- plans$pens[i]
  settings$randomeffect$int_pen$variance <- plans$pen_sd[i]^2
  # One fresh trial per call counts an error without losing the batch.
  trial_settings <- settings
  trial_settings$replications <- 1
  trials <- vector("list", plans$attempted[i])
  estimates <- list()
  for (j in seq_along(trials)) {
    result <- tryCatch(replicate_simulation(trial_settings), error = identity)
    if (inherits(result, "error")) {
      message("Plan ", i, ", trial ", j, ": ", conditionMessage(result))
      next
    }
    if (length(result) != 1L || !inherits(result[[1]], "lmerMod")) {
      message("Plan ", i, ", trial ", j, ": no lmer fit returned")
      next
    }
    fit <- result[[1]]
    trials[[j]] <- fit
    convergence <- tryCatch(check_convergence(fit), error = identity)
    singular <- tryCatch(check_singularity(fit), error = identity)
    if (!isTRUE(convergence) || !identical(singular, FALSE)) next
    coefficients <- tryCatch(extract_coefficients(fit), error = identity)
    if (inherits(coefficients, "error")) {
      message("Plan ", i, ", trial ", j, ": ", conditionMessage(coefficients))
      next
    }
    # simglm's lmer extractor returns fixed terms, without an effect column.
    required <- c("term", "estimate", "std.error", "statistic")
    if (!all(required %in% names(coefficients))) next
    usable <- sum(coefficients$term == "dietSupplement") == 1L &&
      all(is.finite(coefficients$estimate)) &&
      all(is.finite(coefficients$std.error) & coefficients$std.error > 0) &&
      all(is.finite(coefficients$statistic))
    if (!usable) next
    # simglm applies its two-sided normal critical value to the statistic.
    # Version 1.0.0 halves alpha again if supplied p-values are two-sided.
    coefficients <- coefficients |> select(-any_of("p.value"))
    estimates[[length(estimates) + 1L]] <- coefficients
  }
  plans <- plans |>
    mutate(usable = if_else(row_number() == i, length(estimates), usable))
  if (length(estimates) > 0) {
    power <- compute_statistics(
      estimates,
      settings,
      power = TRUE,
      type_1_error = FALSE,
      precision = FALSE
    )
    diet_power <- power |>
      filter(term == "dietSupplement") |>
      pull(power)
    stopifnot(length(diet_power) == 1L, is.finite(diet_power))
    plans <- plans |>
      mutate(
        power_among_usable = if_else(
          row_number() == i, diet_power, power_among_usable
        ),
        detected = if_else(
          row_number() == i, as.integer(round(diet_power * usable)), detected
        )
      )
  }
}
plans <- plans |>
  mutate(
    failed = attempted - usable,
    simulation_se = sqrt(
      power_among_usable * (1 - power_among_usable) / usable
    ),
    lower_if_failures_miss = detected / attempted,
    upper_if_failures_detect = (detected + failed) / attempted
  )
# All-failed plans retain detected = 0, bounds 0 to 1, and NA power and SE.
print(plans, width = Inf)
# Power counts two-sided approximate normal Wald detections below 0.05 for
# the assumed diet change. It is conditional on this generating process and
# analysis; it is not a small-sample t-test or a test of a fitted lmer p-value.
# Errors, nonconvergence, singularity and unusable coefficients are failures.
# Investigate failures. The bounds count each failure first as a missed
# detection, then as a detection. Do not replace failures with fresh draws.
# 100 trials give a rough demonstration: around power 0.5 the simulation SE is
# about 0.05. Use more trials and justified assumptions for an actual decision.
dir.create("outputs", showWarnings = FALSE)
readr::write_csv(plans, "outputs/power-plans.csv")