## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----load-data, eval=F--------------------------------------------------------
# ## Load the example data
# Func_Cox_Data <- readRDS("data/example_data_Cox.rds")
# 
# ## Prepare the functional predictor and censoring indicator
# Func_Cox_Data$wmat <- Func_Cox_Data$MIMS
# Func_Cox_Data$cens <- 1 - Func_Cox_Data$event

## ----fit-model, eval=F--------------------------------------------------------
# library(refundBayes)
# 
# refundBayes_FCox <- refundBayes::fcox_bayes(
#   survtime ~ X1 + s(tmat, by = lmat * wmat, bs = "cc", k = 10),
#   data = Func_Cox_Data,
#   cens = Func_Cox_Data$cens,
#   runStan = TRUE,
#   niter = 1500,
#   nwarmup = 500,
#   nchain = 3,
#   ncores = 3
# )

## ----plot-model, eval=F-------------------------------------------------------
# library(ggplot2)
# plot(refundBayes_FCox)

## ----posterior-summary, eval=F------------------------------------------------
# ## Posterior mean of the functional coefficient
# mean_curve <- apply(refundBayes_FCox$func_coef[[1]], 2, mean)
# 
# ## Pointwise 95% credible interval
# upper_curve <- apply(refundBayes_FCox$func_coef[[1]], 2,
#                      function(x) quantile(x, prob = 0.975))
# lower_curve <- apply(refundBayes_FCox$func_coef[[1]], 2,
#                      function(x) quantile(x, prob = 0.025))

## ----freq-comparison, eval=F--------------------------------------------------
# library(mgcv)
# 
# ## Fit frequentist functional Cox model using mgcv
# fit_freq <- gam(
#   survtime ~ s(tmat, by = lmat * wmat, bs = "cc", k = 10) + X1,
#   data = Func_Cox_Data,
#   family = cox.ph(),
#   weights = Func_Cox_Data$event
# )
# 
# ## Extract frequentist estimates
# freq_result <- plot(fit_freq)

## ----inspect-code, eval=F-----------------------------------------------------
# ## Generate Stan code without running the sampler
# fcox_code <- refundBayes::fcox_bayes(
#   survtime ~ X1 + s(tmat, by = lmat * wmat, bs = "cc", k = 10),
#   data = Func_Cox_Data,
#   cens = Func_Cox_Data$cens,
#   runStan = FALSE
# )
# 
# ## Print the generated Stan code
# cat(fcox_code$stancode)