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

## ----libs---------------------------------------------------------------------
library(maskedhaz)
library(hypothesize)
library(algebraic.mle)

## ----data---------------------------------------------------------------------
set.seed(42)

true_model <- dfr_series_md(components = list(
  dfr_weibull(shape = 2, scale = 5),
  dfr_weibull(shape = 2, scale = 8)
))
df <- rdata(true_model)(theta = c(2, 5, 2, 8), n = 300, tau = 10, p = 0.3)
table(df$omega)

## ----fit-both-----------------------------------------------------------------
null_model <- dfr_series_md(
  components = list(dfr_exponential(), dfr_exponential()),
  n_par = c(1L, 1L)
)
alt_model <- dfr_series_md(
  components = list(dfr_weibull(), dfr_weibull()),
  n_par = c(2L, 2L)
)

fit_null <- suppressWarnings(fit(null_model)(df, par = c(0.1, 0.1)))
fit_alt  <- suppressWarnings(fit(alt_model)(df, par = c(1.5, 4, 1.5, 6)))

coef(fit_null)
coef(fit_alt)

## ----lrt----------------------------------------------------------------------
# lrt() accepts logLik objects directly and derives the dof difference
# from the `df` attribute attached by likelihood.model::logLik.fisher_mle
lrt_result <- lrt(logLik(fit_null), logLik(fit_alt))
lrt_result

## ----wald-per-param-----------------------------------------------------------
est <- coef(fit_alt)
se_vec <- se(fit_alt)         # from algebraic.mle

# Parameter layout: (shape1, scale1, shape2, scale2) = positions 1, 2, 3, 4
w_shape1 <- wald_test(est[1], se = se_vec[1], null_value = 1)
w_shape2 <- wald_test(est[3], se = se_vec[3], null_value = 1)

c(shape1_pval = pval(w_shape1),
  shape2_pval = pval(w_shape2))

## ----adjust-------------------------------------------------------------------
adjusted <- adjust_pval(list(w_shape1, w_shape2), method = "bonferroni")
c(shape1_adj = pval(adjusted[[1]]),
  shape2_adj = pval(adjusted[[2]]))

## ----intersection-------------------------------------------------------------
comp_intersection <- intersection_test(w_shape1, w_shape2)
pval(comp_intersection)

## ----fisher-combine-----------------------------------------------------------
comp_fisher <- fisher_combine(w_shape1, w_shape2)
comp_fisher

## ----score-test---------------------------------------------------------------
# Map the exp-series MLE into the Weibull parameter space:
# (shape = 1, scale = 1/rate) for each component
null_in_alt <- c(1, 1/coef(fit_null)[1],
                 1, 1/coef(fit_null)[2])

# Score and observed Fisher info of the ALT model at the null point
s_at_null <- score(alt_model)(df, par = null_in_alt)
I_at_null <- -hess_loglik(alt_model)(df, par = null_in_alt)

s_at_null

sc <- score_test(s_at_null, I_at_null, null_value = null_in_alt)
sc

## ----invert-test--------------------------------------------------------------
# Build a Wald test parameterised by shape1's null value
wald_for_shape1 <- function(theta0) {
  wald_test(est[1], se = se_vec[1], null_value = theta0)
}

ci_shape1 <- invert_test(wald_for_shape1, grid = seq(1.0, 3.0, by = 0.01), alpha = 0.05)
c(lower = lower(ci_shape1), upper = upper(ci_shape1))

## ----direct-ci----------------------------------------------------------------
confint(wald_test(est[1], se = se_vec[1], null_value = est[1]))

## ----cross-check, eval = requireNamespace("maskedcauses", quietly = TRUE)-----
library(maskedcauses)

fit_null_mc <- suppressWarnings(
  likelihood.model::fit(exp_series_md_c1_c2_c3())(df, par = c(0.1, 0.1))
)

# LRT using maskedcauses null fit + maskedhaz alt fit, across two packages
lrt_cross <- lrt(logLik(fit_null_mc), logLik(fit_alt))
c(
  lrt_same_package = pval(lrt_result),
  lrt_cross_package = pval(lrt_cross)
)