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

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

## ----profile-object-----------------------------------------------------------
set.seed(1)
dat <- simulate_tls(family = "beta_binomial", CTmax = 36, z = 4, phi = 50, seed = 1)
fit <- fit_tls(dat, y = survived, n = total, time = duration, temp = temp,
               family = "beta_binomial", tref = 1)

pc <- profile(fit, "CTmax")
c(estimate = pc$estimate, conf.low = pc$conf.low, conf.high = pc$conf.high,
  cutoff = pc$cutoff, min_deviance = min(pc$deviance))

## ----profile-plot, fig.alt = "Profile-likelihood deviance curve for CTmax: a U-shaped curve with a horizontal profile-t cutoff line; the confidence interval is where the curve lies below the cutoff."----
plot(pc)

## ----equivariance-------------------------------------------------------------
ci_z     <- confint(fit, "z",     method = "profile")
ci_log_z <- confint(fit, "log_z", method = "profile")

# z endpoints equal exp() of the log_z endpoints
rbind(
  z          = c(ci_z$conf.low, ci_z$conf.high),
  exp_log_z  = exp(c(ci_log_z$conf.low, ci_log_z$conf.high))
)

# the interval is asymmetric about the estimate (in general)
with(ci_z, c(lower_gap = estimate - conf.low, upper_gap = conf.high - estimate))

## ----profile-vs-wald----------------------------------------------------------
rbind(
  profile = unlist(confint(fit, "CTmax", method = "profile")[c("conf.low", "conf.high")]),
  wald    = unlist(confint(fit, "CTmax", method = "wald")[c("conf.low", "conf.high")])
)

## ----tidy-methods-------------------------------------------------------------
tidy_parameters(fit, method = "profile")[, c("parameter", "estimate", "conf.low", "conf.high", "interval_type")]

## ----non-closing--------------------------------------------------------------
set.seed(7)
sparse <- simulate_tls(
  temps  = c(35, 36),     # only two temperatures
  times  = c(1, 2),       # only two durations
  reps   = 2, n = 8,
  CTmax  = 36, z = 4,
  family = "binomial", seed = 7
)
sparse_fit <- suppressWarnings(
  fit_tls(sparse, y = survived, n = total, time = duration, temp = temp,
          family = "binomial", tref = 1)
)

## ----bootstrap-fallback-recipe, eval = FALSE----------------------------------
# ci_default <- tryCatch(
#   withCallingHandlers(
#     confint(sparse_fit, "CTmax", method = "profile",
#             fallback = TRUE, nboot = 1000, boot_seed = 7),
#     warning = function(w) {
#       message("caught warning: ", conditionMessage(w))
#       invokeRestart("muffleWarning")
#     }
#   ),
#   error = function(e) e
# )
# ci_default[, c("parameter", "conf.low", "conf.high", "estimate", "method",
#                "conf.status")]

## ----strict-non-closing-------------------------------------------------------
ci_strict <- tryCatch(
  withCallingHandlers(
    confint(sparse_fit, "CTmax", method = "profile", fallback = FALSE),
    warning = function(w) {
      message("caught warning: ", conditionMessage(w))
      invokeRestart("muffleWarning")
    }
  ),
  error = function(e) e
)
ci_strict[, c("parameter", "conf.low", "conf.high", "estimate", "method",
              "conf.status")]

## ----non-closing-eye, fig.alt = "Confidence Eye for a weakly identified fit: a hollow point estimate with no confidence lens, signalling that the profile did not close."----
# fallback = FALSE so the eye refuses to draw a lens when the profile does not
# close (a hollow point only), matching the honest non-closing contract above.
suppressWarnings(plot_confidence_eye(sparse_fit, parm = "CTmax", method = "profile",
                                     fallback = FALSE))

## ----coverage-----------------------------------------------------------------
cov_path <- system.file("extdata", "coverage_results.rds", package = "freqTLS")
if (nzchar(cov_path)) {
  cov <- readRDS(cov_path)
  knitr::kable(
    cov$coverage, digits = 3,
    caption = sprintf("Empirical coverage of 95%% profile CIs (nsim = %d; nominal 0.95).",
                      cov$meta$nsim)
  )
}