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

## ----load-data, eval=F--------------------------------------------------------
# ## Load the example data
# FoSR_exp_data <- readRDS("data/example_data_FoSR.rds")
# 
# ## Set the functional response
# FoSR_exp_data$y <- FoSR_exp_data$MIMS

## ----fit-model, eval=F--------------------------------------------------------
# library(refundBayes)
# 
# refundBayes_FoSR <- refundBayes::fosr_bayes(
#   y ~ X,
#   data = FoSR_exp_data,
#   runStan = TRUE,
#   niter = 1500,
#   nwarmup = 500,
#   nchain = 1,
#   ncores = 1
# )

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

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

## ----freq-comparison, eval=F--------------------------------------------------
# library(refund)
# 
# ## The frequentist approach can be implemented using mgcv or refund::pffr
# ## See Crainiceanu et al. (2024) for details
# 
# fit.freq = pffr(y~-1+X,data=FoSR_exp_data,bs.yindex=list(bs="cc", k=10))
# plotfot = plot(fit.freq)

## ----inspect-code, eval=F-----------------------------------------------------
# ## Generate Stan code without running the sampler
# fosr_code <- refundBayes::fosr_bayes(
#   y ~ X,
#   data = FoSR_exp_data,
#   runStan = FALSE
# )
# 
# ## Print the generated Stan code
# cat(fosr_code$stancode)