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

## ----load-data, eval=F--------------------------------------------------------
# ## Alternative sample data
# 
# # data.SoFR <- readRDS("data/example_data_sofr.rds")
# 
# ## Load the example data
# set.seed(123)
# n <- 100
# M <- 50
# tgrid <- seq(0, 1, length.out = M)
# dt    <- tgrid[2] - tgrid[1]
# tmat  <- matrix(rep(tgrid, each = n), nrow = n)
# lmat  <- matrix(dt, nrow = n, ncol = M)
# wmat  <- t(apply(matrix(rnorm(n * M), n, M), 1, cumsum)) / sqrt(M)
# beta_true <- sin(2 * pi * tgrid)
# X1 <- rnorm(n)
# eta <- 0.5 * X1 + wmat %*% (beta_true * dt)
# prob <- plogis(eta)
# y <- rbinom(n, 1, prob)
# data.SoFR <- data.frame(y = y, X1 = X1)
# data.SoFR$tmat <- tmat
# data.SoFR$lmat <- lmat
# data.SoFR$wmat <- wmat

## ----fit-model, eval=F--------------------------------------------------------
# library(refundBayes)
# 
# refundBayes_SoFR <- refundBayes::sofr_bayes(
#   y ~ X1 + s(tmat, by = lmat * wmat, bs = "cc", k = 10),
#   data = data.SoFR,
#   family = binomial(),
#   runStan = TRUE,
#   niter = 1500,
#   nwarmup = 500,
#   nchain = 3,
#   ncores = 3
# )

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

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

## ----freq-comparison, eval=F--------------------------------------------------
# library(mgcv)
# 
# ## Fit frequentist SoFR using mgcv
# fit_freq <- gam(
#   y ~ s(tmat, by = lmat * wmat, bs = "cc", k = 10) + X1,
#   data = data.SoFR,
#   family = "binomial"
# )
# 
# ## Extract frequentist estimates
# freq_result <- plot(fit_freq)

## ----inspect-code, eval=F-----------------------------------------------------
# ## Generate Stan code without running the sampler
# sofr_code <- refundBayes::sofr_bayes(
#   y ~ X1 + s(tmat, by = lmat * wmat, bs = "cc", k = 10),
#   data = data.SoFR,
#   family = binomial(),
#   runStan = FALSE
# )
# 
# ## Print the generated Stan code
# cat(sofr_code$stancode)