## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")

## ----analytical---------------------------------------------------------------
library(likelihood.contr)
library(likelihood.model)

# Exponential exact: log f(t; lambda) = log(lambda) - lambda * t
exp_exact <- contr_fn(
  loglik = function(df, par, ...) {
    sum(dexp(df$t, rate = par[1], log = TRUE))
  },
  score = function(df, par, ...) {
    c(rate = nrow(df) / par[1] - sum(df$t))
  },
  hess = function(df, par, ...) {
    matrix(-nrow(df) / par[1]^2, 1, 1)
  }
)

## ----mixed-model--------------------------------------------------------------
model <- likelihood_contr(
  obs_type = "status",
  exact = exp_exact,
  right = contr_name("exp", "right", ob_col = "t")
)

set.seed(42)
raw <- rexp(200, rate = 2)
df <- data.frame(
  t      = pmin(raw, 1.0),
  status = ifelse(raw <= 1.0, "exact", "right")
)

result <- fit(model)(df, par = c(rate = 1))
summary(result)

## ----lrt----------------------------------------------------------------------
set.seed(99)
df_test <- data.frame(
  t      = rweibull(200, shape = 2, scale = 3),
  status = "exact"
)

null_model <- likelihood_contr(
  obs_type = "status",
  exact = contr_name("exp", "exact", ob_col = "t")
)
alt_model <- likelihood_contr(
  obs_type = "status",
  exact = contr_name("weibull", "exact", ob_col = "t")
)

null_fit <- suppressWarnings(fit(null_model)(df_test, par = c(rate = 0.5)))
alt_fit  <- suppressWarnings(fit(alt_model)(df_test, par = c(shape = 1.5, scale = 2)))

lrt(
  null     = null_model,
  alt      = alt_model,
  data     = df_test,
  null_par = coef(null_fit),
  alt_par  = coef(alt_fit),
  dof      = 1
)

## ----fim----------------------------------------------------------------------
model_with_rdata <- likelihood_contr(
  obs_type = "status",
  exact = contr_name("exp", "exact", ob_col = "t"),
  rdata_fn = function(theta, n, ...) {
    data.frame(t = rexp(n, rate = theta[1]), status = "exact")
  }
)

set.seed(1)
I <- fim(model_with_rdata)(theta = c(rate = 2), n_obs = 100, n_samples = 500)
I