---
title: "Bayesian Function-on-Scalar Regression with `refundBayes::fosr_bayes`"
output: rmarkdown::html_vignette
date: "2026-03-03"
vignette: >
  %\VignetteIndexEntry{Bayesian Function-on-Scalar Regression}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# Introduction

This vignette provides a detailed guide to the `fosr_bayes()` function in the `refundBayes` package, which fits Bayesian Function-on-Scalar Regression (FoSR) models using Stan. In contrast to scalar-on-function regression (SoFR), where the outcome is scalar and the predictors include functional variables, FoSR reverses this relationship: the outcome is a functional variable (a curve observed over a continuum) and the predictors are scalar.

The methodology follows the tutorial by Jiang, Crainiceanu, and Cui (2025), *Tutorial on Bayesian Functional Regression Using Stan*, published in *Statistics in Medicine*. The procedure for fitting a generalized multilevel Bayesian FoSR was introduced by Goldsmith, Zipunnikov, and Schrack (2015).

# Install the `refundBayes` Package

The `refundBayes` package can be installed from GitHub:

```r
library(remotes)
remotes::install_github("https://github.com/ZirenJiang/refundBayes")
```

# Statistical Model

## The FoSR Model

Function-on-Scalar Regression (FoSR) models the relationship between a functional outcome and one or more scalar predictors. Two key differences from traditional regression are: (1) the outcome $Y_i(t)$ is multivariate and often high-dimensional, and (2) the residuals $e_i(t)$ are correlated across $t$ (points in the domain).

For subject $i = 1, \ldots, n$, let $Y_i(t)$ be the functional response observed at time points $t = t_1, \ldots, t_M$ over a domain $\mathcal{T}$. The domain $\mathcal{T}$ is **not** restricted to $[0,1]$; it is determined by the actual time points in the data (e.g., $\mathcal{T} = [1, 1440]$ for minute-level 24-hour data). The FoSR model assumes:

$$Y_i(t) = \sum_{p=1}^P X_{ip} \beta_p(t) + e_i(t)$$

where:

- $X_{ip}$ with $p = 1, \ldots, P$ are the scalar predictors (the first covariate is often the intercept, $X_{i1} = 1$),
- $\beta_p(t)$ are the corresponding domain-varying (or functional) coefficients,
- $e_i(t)$ is the residual process, which is correlated across $t$.

Each functional coefficient $\beta_p(t)$ describes how the $p$-th scalar predictor affects the outcome at each point $t$ in the domain. This allows the effect of a scalar predictor to vary smoothly over the domain.

## Modeling the Residual Structure

To account for the within-subject correlation in the residuals, the model decomposes $e_i(t)$ using functional principal components:

$$e_i(t) = \sum_{j=1}^J \xi_{ij} \phi_j(t) + \epsilon_i(t)$$

where $\phi_1(t), \ldots, \phi_J(t)$ are the eigenfunctions estimated via FPCA (using `refund::fpca.face`), $\xi_{ij}$ are the subject-specific scores with $\xi_{ij} \sim N(0, \lambda_j)$, and $\epsilon_i(t) \sim N(0, \sigma_\epsilon^2)$ is independent measurement error. The eigenvalues $\lambda_j$ and the error variance $\sigma_\epsilon^2$ are estimated from the data.

## Functional Coefficients via Penalized Splines

Each functional coefficient $\beta_p(t)$ is represented using $K$ spline basis functions $\psi_1(t), \ldots, \psi_K(t)$:

$$\beta_p(t) = \sum_{k=1}^K b_{pk} \psi_k(t)$$

Substituting into the model:

$$Y_i(t) = \sum_{k=1}^K \left(\sum_{p=1}^P X_{ip} b_{pk}\right) \psi_k(t) + \sum_{j=1}^J \xi_{ij} \phi_j(t) + \epsilon_i(t)$$

Let $B_{ik} = \sum_{p=1}^P X_{ip} b_{pk}$ and denote by $\boldsymbol{\Psi}$ the $M \times K$ matrix of spline basis evaluations, by $\boldsymbol{\Phi}$ the $M \times J$ matrix of eigenfunctions, and by $\mathbf{B}_i = (B_{i1}, \ldots, B_{iK})^t$. The model in matrix form becomes:

$$\mathbf{Y}_i = \boldsymbol{\Psi} \mathbf{B}_i + \boldsymbol{\Phi} \boldsymbol{\xi}_i + \boldsymbol{\epsilon}_i$$

where $\mathbf{Y}_i = \{Y_i(t_1), \ldots, Y_i(t_M)\}^t$ is the $M \times 1$ vector of observed functional data for subject $i$.

## Smoothness Penalty

Smoothness of each $\beta_p(t)$ is induced through a quadratic penalty on the spline coefficients $\mathbf{b}_p = (b_{p1}, \ldots, b_{pK})^t$:

$$p(\mathbf{b}_p) \propto \exp\left(-\frac{\mathbf{b}_p^t \mathbf{S} \mathbf{b}_p}{2\sigma_p^2}\right)$$

where $\mathbf{S}$ is the penalty matrix and $\sigma_p^2$ is the smoothing parameter for the $p$-th functional coefficient. Importantly, each functional coefficient $\beta_p(t)$ has its own smoothing parameter $\sigma_p^2$, allowing different levels of smoothness across predictors.

Unlike the SoFR and FCox models, the FoSR implementation uses the original spline basis without the spectral reparametrisation. This choice follows Goldsmith et al. (2015) and simplifies interpretation by avoiding the need to back-transform the estimated coefficients.

## Full Bayesian Model

The complete Bayesian FoSR model is:

$$\begin{cases}
\mathbf{Y}_i = \boldsymbol{\Psi} \mathbf{B}_i + \boldsymbol{\Phi} \boldsymbol{\xi}_i + \boldsymbol{\epsilon}_i \\
B_{ik} = \sum_{p=1}^P X_{ip} b_{pk}, \quad k = 1, \ldots, K \\
p(\mathbf{b}_p) \propto \exp(-\mathbf{b}_p^t \mathbf{S} \mathbf{b}_p / 2\sigma_p^2), \quad p = 1, \ldots, P \\
\xi_{ij} \sim N(0, \lambda_j), \quad \lambda_j \sim p(\lambda_j), \quad j = 1, \ldots, J \\
\epsilon_i(t_m) \sim N(0, \sigma_\epsilon^2), \quad i = 1, \ldots, n, \; t = t_1, \ldots, t_M \\
\sigma_\epsilon^2 \sim p(\sigma_\epsilon^2), \quad \sigma_p^2 \sim p(\sigma_p^2), \quad p = 1, \ldots, P
\end{cases}$$

The priors $p(\lambda_j)$, $p(\sigma_\epsilon^2)$, and $p(\sigma_p^2)$ are non-informative priors on variance components, using inverse-Gamma priors $IG(0.001, 0.001)$.

The model assumes that the functional responses are observed on a common set of grid points across all subjects. This ensures that the $\boldsymbol{\Psi}$ and $\boldsymbol{\Phi}$ matrices are identical across subjects, simplifying computation in Stan.

# The `fosr_bayes()` Function

## Usage

```r
fosr_bayes(
  formula,
  data,
  joint_FPCA = NULL,
  runStan = TRUE,
  niter = 3000,
  nwarmup = 1000,
  nchain = 3,
  ncores = 1,
  spline_type = "bs",
  spline_df = 10
)
```

## Arguments

| Argument | Description |
|:---------|:------------|
| `formula` | Functional regression formula. The left-hand side should be the name of the functional response variable (an $n \times M$ matrix in `data`). The right-hand side specifies scalar predictors using standard formula syntax. |
| `data` | A data frame containing all variables used in the model. The functional response should be stored as an $n \times M$ matrix, where each row is one subject and each column is one time point. |
| `joint_FPCA` | A logical (`TRUE`/`FALSE`) vector of the same length as the number of functional predictors, indicating whether to jointly model FPCA. Default is `NULL`, which sets all entries to `FALSE`. |
| `runStan` | Logical. Whether to run the Stan program. If `FALSE`, the function only generates the Stan code and data without sampling. Default is `TRUE`. |
| `niter` | Total number of Bayesian posterior sampling iterations (including warmup). Default is `3000`. |
| `nwarmup` | Number of warmup (burn-in) iterations. Default is `1000`. |
| `nchain` | Number of Markov chains for posterior sampling. Default is `3`. |
| `ncores` | Number of CPU cores to use when executing the chains in parallel. Default is `1`. |
| `spline_type` | Type of spline basis used for the functional coefficients. Default is `"bs"` (B-splines). |
| `spline_df` | Number of degrees of freedom (basis functions) for the spline basis. Default is `10`. |

## Return Value

The function returns a list of class `"refundBayes"` containing the following elements:

| Element | Description |
|:--------|:------------|
| `stanfit` | The Stan fit object (class `stanfit`). Can be used for convergence diagnostics, traceplots, and additional summaries via the `rstan` package. |
| `spline_basis` | Basis functions used to reconstruct the functional coefficients from the posterior samples. |
| `stancode` | A character string containing the generated Stan model code. |
| `standata` | A list containing the data passed to the Stan model. |
| `func_coef` | An array of posterior samples for functional coefficients, with dimensions $Q \times P \times M$, where $Q$ is the number of posterior samples, $P$ is the number of scalar predictors, and $M$ is the number of time points on the functional domain. |
| `family` | The family type: `"functional"`. |

## Formula Syntax

The formula syntax for FoSR differs from SoFR and FCox because the outcome is functional and the predictors are scalar. The left-hand side specifies the functional response, and the right-hand side lists scalar predictors:

```r
y ~ X
```

where:

- `y`: the name of the functional response variable in `data`. This should be an $n \times M$ matrix, where each row contains the functional observations for one subject across $M$ time points.
- `X`: scalar predictor(s). Multiple scalar predictors can be included using standard formula syntax (e.g., `y ~ X1 + X2 + X3`).

Note that the FoSR formula does **not** use the `s()` term with `tmat`, `lmat`, and `wmat` that appears in SoFR and FCox models. This is because in FoSR, the functional variable is the *outcome* (not a predictor), and the spline basis for the functional coefficients is constructed internally using the `spline_type` and `spline_df` arguments.

# Example: Bayesian FoSR

We demonstrate the `fosr_bayes()` function using a simulated example dataset with a functional response and one scalar predictor.

## Load and Prepare Data

```{r 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
```

The example dataset contains:

- `MIMS`: an $n \times M$ matrix of functional response values (physical activity data), assigned to `y`,
- `X`: a scalar predictor variable.

## Fit the Bayesian FoSR Model

```{r 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
)
```

In this call:

- The formula specifies the functional response `y` (an $n \times M$ matrix) with one scalar predictor `X`.
- The spline basis for the functional coefficients is constructed internally using B-splines (`spline_type = "bs"`) with 10 degrees of freedom (`spline_df = 10`) by default.
- The eigenfunctions for the residual structure are estimated automatically via FPCA using `refund::fpca.face`.
- The sampler runs 1 chain with 1500 total iterations (500 warmup + 1000 posterior samples). For production analyses, using 3 or more chains is recommended for convergence assessment.

### A Note on Computation

FoSR models are computationally more demanding than SoFR models because the likelihood involves all $M$ time points for each of the $n$ subjects. The example above uses a single chain for demonstration purposes. In practice, consider the trade-off between the number of basis functions, the number of FPCA components, and the computational budget.

## Visualization

The `plot()` method for `refundBayes` objects displays the estimated functional coefficients with pointwise 95% credible intervals:

```{r plot-model, eval=F}
library(ggplot2)
plot(refundBayes_FoSR)
```

## Extracting Posterior Summaries

Posterior summaries of the functional coefficients can be computed from the `func_coef` element. The `func_coef` object is an array with dimensions $Q \times P \times M$, where $Q$ is the number of posterior samples, $P$ is the number of scalar predictors (including the intercept), and $M$ is the number of time points:

```{r 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))
```

## Comparison with Frequentist Results

The Bayesian results can be compared with frequentist estimates obtained via `mgcv::gam`. As described in Crainiceanu et al. (2024), one key distinction is that the frequentist approach fits the model at each time point independently or using fast algorithms, while the Bayesian approach jointly models all time points with a shared smoothness prior and FPCA-based residual structure:

```{r 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)
```

Simulation studies in Jiang et al. (2025) show that the Bayesian FoSR achieves similar estimation accuracy (RISE) to the frequentist approach, while providing superior coverage of pointwise credible intervals.

## Inspecting the Generated Stan Code

Setting `runStan = FALSE` allows you to inspect or modify the Stan code before running the model:

```{r 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)
```

# Practical Recommendations

- **Number of basis functions (`spline_df`)**: The default is `spline_df = 10`. In practice, the appropriate number depends on the complexity of the underlying functional coefficients and the resolution of the data. More basis functions allow greater flexibility but increase computational cost.
- **Spline type (`spline_type`)**: The default `"bs"` uses B-splines. Other basis types supported by `mgcv` may also be used depending on the application.
- **Number of FPCA components**: The number of eigenfunctions $J$ used to model the residual correlation is determined automatically by `refund::fpca.face`. If the residual structure is complex, more components may be needed.
- **Number of iterations and chains**: FoSR models are computationally intensive. Start with fewer iterations and a single chain for exploratory analysis, then increase for final inference. A recommended setup for production is `niter = 3000`, `nwarmup = 1000`, `nchain = 3`.
- **Convergence diagnostics**: After fitting, examine traceplots and $\hat{R}$ statistics using the `rstan` package (e.g., `rstan::traceplot(refundBayes_FoSR$stanfit)`). Warnings about bulk ESS, tail ESS, or $\hat{R}$ indicate that more iterations or chains may be needed.
- **Common grid assumption**: The current implementation assumes that the functional responses are observed on a common set of grid points across all subjects. Subject-specific observation grids are not yet supported.

# References

- Jiang, Z., Crainiceanu, C., and Cui, E. (2025). Tutorial on Bayesian Functional Regression Using Stan. *Statistics in Medicine*, 44(20--22), e70265.
- Crainiceanu, C. M., Goldsmith, J., Leroux, A., and Cui, E. (2024). *Functional Data Analysis with R*. CRC Press.
- Goldsmith, J., Zipunnikov, V., and Schrack, J. (2015). Generalized Multilevel Function-on-Scalar Regression and Principal Component Analysis. *Biometrics*, 71(2), 344--353.
- Carpenter, B., Gelman, A., Hoffman, M. D., et al. (2017). Stan: A Probabilistic Programming Language. *Journal of Statistical Software*, 76(1), 1--32.