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

## ----define-states------------------------------------------------------------
library(RFmstate)

# Use the built-in clinical trial structure
ms <- clinical_states()
print(ms)

## ----custom-states, eval=FALSE------------------------------------------------
# # A simple 3-state illness-death model
# ms_simple <- define_multistate(
#   state_names = c("Healthy", "Sick", "Dead"),
#   absorbing = "Dead",
#   transitions = list(
#     Healthy = c("Sick", "Dead"),
#     Sick = c("Dead")
#   )
# )
# 
# # A 4-state model with recovery
# ms_recovery <- define_multistate(
#   state_names = c("Healthy", "Sick", "Recovered", "Dead"),
#   absorbing = "Dead",
#   transitions = list(
#     Healthy = c("Sick", "Dead"),
#     Sick = c("Recovered", "Dead"),
#     Recovered = c("Dead")
#   )
# )

## ----simulate-----------------------------------------------------------------
set.seed(42)
dat <- sim_clinical_data(n = 300, structure = ms)
head(dat)

## ----prepare------------------------------------------------------------------
msdata <- prepare_data(
  data = dat, id = "ID", structure = ms,
  time_map = list(
    Responded = "time_Responded",
    Unresponded = "time_Unresponded",
    Stabilized = "time_Stabilized",
    Progressed = "time_Progressed",
    Death = "time_Death"
  ),
  censor_col = "time_censored",
  covariates = c("age", "sex", "BMI", "treatment")
)
print(msdata)

## ----aj-----------------------------------------------------------------------
aj <- aalen_johansen(msdata)
print(aj)

## ----aj-plot, fig.cap="State occupation probabilities from Aalen-Johansen estimator"----
plot(aj, type = "state_occupation")

## ----aj-hazard, fig.cap="Nelson-Aalen cumulative hazards"---------------------
plot(aj, type = "cumulative_hazard")

## ----fit----------------------------------------------------------------------
fit <- rfmstate(
  msdata,
  covariates = c("age", "sex", "BMI", "treatment"),
  num.trees = 200,
  seed = 42
)
print(fit)

## ----summary------------------------------------------------------------------
s <- summary(fit)

## ----importance, fig.cap="Feature importance per transition"------------------
imp <- importance(fit)
print(imp)
plot(imp, type = "barplot")

## ----importance-heat, fig.cap="Feature importance heatmap"--------------------
plot(imp, type = "heatmap")

## ----predict, fig.cap="Predicted state occupation for two patient profiles"----
newdata <- data.frame(
  age = c(50, 70),
  sex = c(0, 1),
  BMI = c(24, 32),
  treatment = c(1, 0)
)

pred <- predict(fit, newdata = newdata, times = seq(10, 365, by = 10))

# Plot for patient 1 (young, treated)
plot(pred, type = "state_occupation", subject = 1)

# Plot for patient 2 (older, untreated)
plot(pred, type = "state_occupation", subject = 2)

## ----diagnostics--------------------------------------------------------------
diag <- diagnose(fit)
print(diag)

## ----diag-brier, fig.cap="Time-dependent Brier score"-------------------------
plot(diag, type = "brier")

## ----diag-concordance, fig.cap="Concordance index per transition"-------------
plot(diag, type = "concordance")

## ----diag-bv, fig.cap="Bias-variance decomposition"---------------------------
plot(diag, type = "bias_variance")

## ----diagram, fig.cap="Transition diagram with event counts"------------------
plot_transition_diagram(ms, msdata)