## ----include = FALSE---------------------------------------------------------- knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 3.6, dpi = 96 ) ## ----setup-------------------------------------------------------------------- library(freqTLS) ## ----identifiability---------------------------------------------------------- # Sparse design: few temperatures, so the CTmax/z slope is weakly identified. sparse <- simulate_tls( temps = c(35, 36), times = c(1, 4), reps = 2, n = 20, CTmax = 36, z = 4, family = "binomial", seed = 11 ) fit_sparse <- suppressWarnings(fit_tls( sparse, y = survived, n = total, time = duration, temp = temp, family = "binomial", tref = 1 )) # freqTLS says so: a profile that does not close returns NA on the open side # (here with the prior-free bootstrap fallback turned off to show the raw signal). suppressWarnings(confint(fit_sparse, "z", method = "profile", fallback = FALSE)) ## ----calibration, echo = FALSE------------------------------------------------ cal <- readRDS(system.file("extdata", "calibration_results.rds", package = "freqTLS")) knitr::kable( data.frame( "Sample size" = sub("mid", "medium", cal$cell), "Residual df" = cal$median_df, "Coverage (Wald, z-ref)" = sprintf("%.3f", cal$cov_CTmax_z), "Coverage (Wald, t-ref)" = sprintf("%.3f", cal$cov_CTmax_t), "Extra width" = sprintf("%+.0f%%", 100 * (cal$width_CTmax_t / cal$width_CTmax_z - 1)), check.names = FALSE ), caption = paste( "95% confidence-interval coverage for CTmax (nominal 0.95),", "over ~500 simulated datasets per cell. Monte-Carlo SE ≈ 0.01." ) ) ## ----shrinkage---------------------------------------------------------------- # Data-poor colonies (few assays each) so partial pooling has work to do. d_re <- simulate_tls( family = "binomial", temps = c(34, 36, 38), times = c(1, 4), reps = 1, n = 8, CTmax = 36, z = 4, re_sd = 1.5, n_re_groups = 12, seed = 42 ) # Partial pooling: a random intercept shrinks the colony effects toward 0. fit_re <- suppressWarnings(fit_tls( tls_bf(survived | trials(total) ~ time(duration) + temp(temp), CTmax ~ 1 + (1 | colony)), data = d_re, family = "binomial", tref = 1 )) blup <- ranef(fit_re)$estimate # No pooling: a fixed CTmax per colony (the unshrunk per-group estimates). fit_fix <- suppressWarnings(fit_tls( d_re, y = survived, n = total, time = duration, temp = temp, group = colony, family = "binomial", tref = 1 )) ct_fix <- fit_fix$estimates$estimate[startsWith(fit_fix$estimates$parameter, "CTmax:")] # The Gaussian prior on b_g (scale sigma_CTmax) pulls the spread in. c(no_pooling_sd = round(sd(ct_fix - mean(ct_fix)), 2), partial_pooling_sd = round(sd(blup), 2), sigma_CTmax = round(fit_re$estimates$estimate[fit_re$estimates$parameter == "sigma_CTmax"], 2))