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.

Usage

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

Arguments

Return Value

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

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:

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.

Load and Prepare Data

## 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

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.

Visualization

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

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:

## 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:

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:

## 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)