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

## ----load-data, eval=T--------------------------------------------------------
## 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=T--------------------------------------------------
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=T-----------------------------------------------------
## 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)