## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment  = "#>",
  message  = FALSE,
  warning  = FALSE
)
has_survival <- requireNamespace("survival", quietly = TRUE)

## -----------------------------------------------------------------------------
library(mmbcv)

data("msdat3")
dim(msdat3)
length(unique(msdat3$id))
length(unique(msdat3$clus_id))
table(msdat3$event)

## -----------------------------------------------------------------------------
with(msdat3, table(from, to))

## -----------------------------------------------------------------------------
library(survival)

fit <- coxph(
  list(
    Surv(Tstart, Tstop, event) ~ 1,
    state("(s0)"):state("S1") + state("S1"):state("S2") + state("S2"):state("S3") ~
      Z + X,
    state("(s0)"):state("D") + state("S1"):state("D") + state("S2"):state("D") +
      state("S3"):state("D") ~ (Z + X)/common
  ),
  data = msdat3,
  id = id,
  ties = "breslow",
  timefix = FALSE
)

fit

## -----------------------------------------------------------------------------
out <- MMBCV(
  fit, msdat3,
  StartTime = Tstart,
  StopTime  = Tstop,
  ClusterID = clus_id,
  SubjectID = id,
  Event     = event,
  tie       = "breslow",
  details   = FALSE
)

names(out)

## -----------------------------------------------------------------------------
Vlist <- out[c("robust","varMR","varMD","varMDMR","varFG","varFGMR","varKC","varKCMR","varMBN","varMBNMR")]
SE <- sapply(Vlist, function(V) sqrt(diag(V)))
SE

## -----------------------------------------------------------------------------
out_efron <- MMBCV(
  fit, msdat3,
  StartTime = Tstart,
  StopTime  = Tstop,
  ClusterID = clus_id,
  SubjectID = id,
  Event     = event,
  tie       = "efron",
  details   = FALSE
)

sqrt(diag(out_efron$robust))

## ----test-model-inputs--------------------------------------------------------
recurrent_transitions <- c("1:2", "2:3", "3:4")
z_index <- unname(fit$cmap["Z", recurrent_transitions])

beta_Z <- fit$coefficients[z_index]
V_Z_MD <- out$varMD[z_index, z_index, drop = FALSE]

## ----test-model-examples------------------------------------------------------
H_Z <- heterogeneity_test(beta_Z, V_Z_MD)

L_Z <- linear_trend_test(
  beta_Z,
  V_Z_MD,
  scores = 0:2,
  alternative = "greater"
)

O_Z <- order_restricted_test(
  beta_Z,
  V_Z_MD,
  alternative = "increasing"
)

H_Z
L_Z
O_Z

## ----test-model-decisions-----------------------------------------------------
test_p_values <- c(
  heterogeneity    = H_Z$p.value,
  linear_trend     = L_Z$p.value,
  order_restricted = O_Z$p.value
)

test_p_values
test_p_values < 0.05

## ----test-seven-coefficients--------------------------------------------------
beta7 <- c(
  beta1 = -0.31,
  beta2 = -0.21,
  beta3 = -0.12,
  beta4 = -0.02,
  beta5 =  0.09,
  beta6 =  0.21,
  beta7 =  0.34
)

se7 <- c(0.12, 0.11, 0.10, 0.10, 0.11, 0.12, 0.13)
position <- seq_along(beta7)
cor7 <- 0.35 ^ abs(outer(position, position, "-"))
V7 <- diag(se7) %*% cor7 %*% diag(se7)
dimnames(V7) <- list(names(beta7), names(beta7))

H7 <- heterogeneity_test(beta7, V7)

L7 <- linear_trend_test(
  beta7,
  V7,
  scores = 0:6,
  alternative = "greater"
)

O7 <- order_restricted_test(
  beta7,
  V7,
  alternative = "increasing",
  nsim = 5000,
  seed = 20260713
)

H7
L7
O7

## ----test-seven-components----------------------------------------------------
H7$contrast_matrix
H7$differences

L7$slope
L7$stderr
L7$scores

O7$weights
O7$weight_mcse
O7$p.value.mcse

## ----test-subset--------------------------------------------------------------
keep <- c(1, 3, 4, 5)
selected_scores <- c(0, 2, 3, 4)

H_subset <- heterogeneity_test(
  beta7,
  V7,
  index = keep
)

L_subset <- linear_trend_test(
  beta7,
  V7,
  index = keep,
  scores = selected_scores,
  alternative = "greater"
)

O_subset <- order_restricted_test(
  beta7,
  V7,
  index = keep,
  alternative = "increasing",
  nsim = 2000,
  seed = 20260714
)

c(
  heterogeneity    = H_subset$p.value,
  linear_trend     = L_subset$p.value,
  order_restricted = O_subset$p.value
)