## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 6,
  fig.height = 3.6,
  dpi = 96
)

## ----setup--------------------------------------------------------------------
library(freqTLS)

## ----ctmax-re-----------------------------------------------------------------
set.seed(1)
d <- simulate_tls(family = "binomial", CTmax = 36, z = 4,
                  re_sd = 1.5, n_re_groups = 14, seed = 1)
fit <- fit_tls(
  tls_bf(survived | trials(total) ~ time(duration) + temp(temp),
         CTmax ~ 1 + (1 | colony)),
  data = d, family = "binomial", tref = 1, quiet = TRUE)

# The estimates table carries the population CTmax / z and the between-colony SD.
est <- fit$estimates
est[est$parameter %in% c("CTmax", "z", "sigma_CTmax"),
    c("parameter", "estimate", "std.error")]

## ----re-recovery, echo = FALSE------------------------------------------------
rr_path <- system.file("extdata", "re_recovery_results.rds", package = "freqTLS")
if (nzchar(rr_path) && file.exists(rr_path)) {
  rr <- readRDS(rr_path)
  knitr::kable(
    rr$re_recovery[, c("n_groups", "mean_sigma", "rel_bias", "CTmax_wald_coverage")],
    col.names = c("n groups", "mean sigma_CTmax", "rel. bias", "CTmax 95% coverage"),
    caption = sprintf(
      paste("Random-intercept recovery vs number of groups (true sigma_CTmax = %.1f;",
            "%d sims/cell). The ML SD is biased low and the fixed-effect CTmax",
            "interval under-covers with few groups; both settle by ~14 groups."),
      rr$meta$truth_sigma, rr$meta$nsim),
    digits = 3)
}

## ----ranef--------------------------------------------------------------------
head(ranef(fit), 4)

## ----re-predict---------------------------------------------------------------
colony_1 <- as.character(ranef(fit)$group[1])
new_colony <- data.frame(temp = 36, duration = 2, colony = colony_1)
data.frame(
  target = c("population", "colony"),
  survival = c(
    predict(fit, new_colony, re.form = "population"),
    predict(fit, new_colony, re.form = "conditional")
  )
)

## ----ctmax-ci-----------------------------------------------------------------
suppressMessages(confint(fit, "sigma_CTmax", method = "wald"))[
  , c("parameter", "conf.low", "conf.high")]
suppressMessages(confint(fit, "CTmax", method = "profile", npoints = 8))[
  , c("parameter", "conf.low", "conf.high", "method")]

## ----eye, fig.alt = "Confidence Eye for the population CTmax of a random-intercept fit: a pale confidence lens with a hollow point estimate near 36 degrees Celsius, drawn with a Wald interval; a freqTLS uncertainty display, not a posterior density."----
suppressMessages(plot_confidence_eye(fit, parm = "CTmax"))

## ----logz-re------------------------------------------------------------------
set.seed(2)
dz <- simulate_tls(family = "binomial", CTmax = 36, z = 4,
                   re_sd_z = 0.3, n_re_groups = 14, seed = 2)
fit_z <- fit_tls(
  tls_bf(survived | trials(total) ~ time(duration) + temp(temp),
         log_z ~ 1 + (1 | colony)),
  data = dz, family = "binomial", tref = 1, quiet = TRUE)
sg <- fit_z$estimates$estimate[fit_z$estimates$parameter == "sigma_logz"]
c(sigma_logz = round(sg, 3),
  approx_fold_spread_in_z = round(exp(sg), 3))

## ----combined-----------------------------------------------------------------
set.seed(3)
db <- simulate_tls(family = "binomial", CTmax = 36, z = 4,
                   re_sd = 1.2, re_sd_z = 0.3, n_re_groups = 14, seed = 3)
# Both random intercepts on the same `colony` grouping. This fits, but emits a
# warning that the two variances are independent (no correlation term); we pass
# `quiet = TRUE` to silence it here and discuss it below.
fit_both <- fit_tls(
  tls_bf(survived | trials(total) ~ time(duration) + temp(temp),
         CTmax ~ 1 + (1 | colony), log_z ~ 1 + (1 | colony)),
  data = db, family = "binomial", tref = 1, quiet = TRUE)
fit_both$estimates$parameter[grepl("^sigma", fit_both$estimates$parameter)]