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

## ----data---------------------------------------------------------------------
library(PMLE4SCR)
data(BMT, package = "SemiCompRisks")
BMT$g <- factor(BMT$g, levels = c(2, 3, 1),
                labels = c("AML-low", "AML-high", "ALL"))

## ----data-table, echo = FALSE-------------------------------------------------
var_desc <- data.frame(
  Variable = c("T1", "T2", "delta1", "delta2", "g"),
  Description = c(
    "Time to death (terminal event)",
    "Time to leukemia relapse (non-terminal event)",
    "Death indicator (1 = observed, 0 = censored)",
    "Relapse indicator (1 = observed, 0 = censored)",
    "Disease group (AML-low, AML-high, ALL)"
  )
)
knitr::kable(var_desc, caption = "Table 1: Key variables in the BMT dataset")

## ----head---------------------------------------------------------------------
head(BMT[, c("T1", "T2", "delta1", "delta2", "g")])

## ----pmle---------------------------------------------------------------------
myfit <- PMLE4SCR(BMT, time = "T2", death = "T1",
                  status_time = "delta2", status_death = "delta1",
                  T.fmla = ~ g, D.fmla = ~ g,
                  copula.family = "Clayton",
                  copula.control = list(formula = ~ g))

## ----gamma--------------------------------------------------------------------
knitr::kable(myfit$gamma,
             digits = 3,
             caption = "Table 2: Estimated copula parameters (PMLE)")

## ----betaT, echo = FALSE------------------------------------------------------
tab <- myfit$betaT
rownames(tab) <- NULL
tab$para <- c("AML-high", "ALL")
knitr::kable(tab,
             col.names = c("Group", "Estimate", "SE"),
             digits = 3,
             caption = "Table 3: Regression coefficients for relapse time (PMLE)")

## ----betaD, echo = FALSE------------------------------------------------------
tab <- myfit$betaD
rownames(tab) <- NULL
tab$para <- c("AML-high", "ALL")
knitr::kable(tab,
             col.names = c("Group", "Estimate", "SE"),
             digits = 3,
             caption = "Table 4: Regression coefficients for death time (PMLE)")

## ----mle----------------------------------------------------------------------
myfit_mle <- MLE4SCR(BMT, time = "T2", death = "T1",
                     status_time = "delta2", status_death = "delta1",
                     T.fmla = ~ g, D.fmla = ~ g,
                     copula.family = "Clayton",
                     copula.control = list(formula = ~ g))

## ----comparison, echo = FALSE-------------------------------------------------
comp <- data.frame(
  Parameter = c("Copula (Intercept)", "Copula (AML-high)", "Copula (ALL)",
                "Relapse: AML-high", "Relapse: ALL",
                "Death: AML-high", "Death: ALL"),
  PMLE_Est = c(myfit$gamma$est,
               myfit$betaT$est,
               myfit$betaD$est),
  PMLE_SE  = c(myfit$gamma$se,
               myfit$betaT$se,
               myfit$betaD$se),
  MLE_Est  = c(myfit_mle$gamma$est,
               myfit_mle$betaT$est,
               myfit_mle$betaD$est),
  MLE_SE   = c(myfit_mle$gamma$se,
               myfit_mle$betaT$se,
               myfit_mle$betaD$se)
)
knitr::kable(comp,
             digits = 3,
             col.names = c("Parameter",
                           "Estimate", "SE",
                           "Estimate", "SE"),
             caption = "Table 5: Comparison of PMLE and simultaneous MLE estimates") |>
  kableExtra::add_header_above(c(" " = 1, "PMLE" = 2, "Simultaneous MLE" = 2))