Scalar-on-Function Linear Regression Models
The general form of the scalar-on-function regression model is given by \[T(F_{Y_i|X_i,Z_i}) = \sum_{l=1}^{L} \int_{\Omega_l} \beta_l(t) X_{li}(t) dt + (1,Z_i^T)\gamma\] where
- \(Y_i\) is the scalar response variable;
- \(X_{li}(t)\) denotes the \(l\)th function-valued covariate defined on the domain \(\Omega_l\);
- \(\beta_l(t)\) is the corresponding function-valued regression coefficient, and both \(\beta_l\) and \(X_{li}\) belong to \(L^2(\Omega_l)\);
- \(Z_i = (Z_{1i}, \dots, Z_{qi})^T\) represents the scalar covariates;
- \(\gamma = (\gamma_0, \gamma_1, \dots, \gamma_q)^T\) is the parameter vector for the scalar covariates;
- \(F_{Y_i\mid X_i,Z_i}\) is the cumulative distribution function of \(Y_i\) conditional on \(X_i\) and \(Z_i\);
- \(T(\cdot)\) is a statistical functional.
A statistical functional is defined as a mapping from a probability distribution to a real number. Statistical functionals represent quantifiable properties of probability distributions, (e.g., an expectation, a variance or a quantile).
Two common model types supported in this package are:
Multi-level Generalized Scalar-on-Function Regression:
This model is formulated as \[ T\bigl(F_{Y_i\mid X_i,Z_i}\bigr) = g\Bigl\{E(Y_i\mid X_i,Z_i)\Bigr\} = g\left\{\int_\mathbb{R}ydF_{Y_i|X_i,Z_i}(y)\right\}, \] where \(g(\cdot)\) is a strictly monotonic link function analogous to that used in generalized linear models.Scalar-on-Function Quantile Regression:
This model is expressed as \[ T\bigl(F_{Y_i\mid X_i,Z_i}\bigr) = Q_{Y_i\mid X_i,Z_i}(\tau) = F_{Y_i\mid X_i,Z_i}^{-1}(\tau), \] where \(\tau \in (0,1)\) denotes the quantile of interest.
Representation of Functional Data
Function-valued variables (functional variables) are typically recorded as measurements at specific time points, resulting in an \(n \times m\) data matrix \((x_{ij})\), where each row corresponds to an observation (subject) and each column to a measurement time.
Functional Data
Function-valued variables (functional variables) are typically recorded as the measurements of the values of the functions at specific (time) points in the domain. The data of a function-valued variable are often presented in the form of a matrix, \((x_{ij})_{n\times m}\), where \(x_{ij} = f_i(t_j)\), represents the value of \(f_i(t_j)\), each row corresponds to an observation (subject) and each column to a measurement time.
S4 Class functional_variable
In the MECfda package, the s4 class
functional_variable represents function-valued data stored
in matrix form. The main slots include:
X: The data matrix \((x_{ij})\);t_0andperiod: The left endpoint and length of the function’s domain, respectively (The class only support the domain in the format of interval);t_points: A vector of the time points at which the measurements are recorded.
A dedicated method, dim, returns the number of subjects
and the number of measurement points for a
functional_variable object.
Basis Expansion
Let \(\{\rho_{k}\}_{k=1}^\infty\) be a complete basis for \(L^2(\Omega)\). For an arbitrary function \(\beta(t) \in L^2(\Omega)\), there exists a sequence of (real) number \(\{c_k\}\) such that \[ \beta(t) = \sum_{k=1}^\infty c_k \, \rho_k(t). \]
For the function-valued variable \(X_i(t)\), we have \[ \int_\Omega \beta(t) X_i(t) \, dt = \int_\Omega X_i(t) \sum_{k=1}^\infty c_{k}\rho_{k}(t) dt = \sum_{k=1}^\infty c_k \int_\Omega \rho_k(t) X_i(t) \, dt. \]
In practice, for \(\int_\Omega\beta(t) X_i(t) dt\) in statistical models, we truncate the infinite series to \(p\) terms (\(p<\infty\)), using \(c_{k} \int_\Omega X_i(t) \rho_{k}(t) dt\) to approximate \(\int_\Omega\beta(t) X_i(t) dt\) so that the problem reduces to estimating a finite number of parameters \(c_1, \dots, c_p\), and treating \(\int_\Omega \rho_k(t) X_i(t) \, dt\) as new covariates. The number of basis, \(p\), can be related to the sample size, \(n\), and the \(p\) is in a form of \(p(n)\).
In practice, we may not necessarily use the truncated complete basis of \(L^2\) function space. Instead, we can just use a finite sequence of linearly independent functions as the basis function. Common choices for the basis include the Fourier basis, b-spline basis, and eigenfunction basis.
In numerical computing level, performing basis expansion methods on function-valued variable data in matrix form as we mentioned is to compute \((b_{ik})_{n\times p}\), where \(b_{ik} = \int_\Omega f(t)\rho_k(t) dt\).
Fourier Basis
On the interval \([t_0, t_0+T]\), the Fourier basis consists of
\[ \frac{1}{2}, \quad \cos\left(\frac{2\pi}{T}k(x-t_0)\right), \quad \sin\left(\frac{2\pi}{T}k(x-t_0)\right), \quad k = 1, 2, \dots \]
S4 Class Fourier_series
In the MECfda package, s4 class
Fourier_series represents a Fourier series of the form
\[ \frac{a_0}{2} + \sum_{k=1}^{p_a} a_k \cos\left(\frac{2\pi k (x-t_0)}{T}\right) + \sum_{k=1}^{p_b} b_k \sin\left(\frac{2\pi k (x-t_0)}{T}\right),\quad x \in [t_0, t_0+T]. \]
Its main slots include:
double_constant: The value of \(a_0\);cosandsin: The coefficients \(a_k\) and \(b_k\) for the cosine and sine terms, respectively;k_cosandk_sin: The frequency values corresponding to the cosine and sine coefficients;t_0andperiod: The left endpoint, \(t_0\), and length, \(T\), of the domain (interval \([t_0, t_0+T]\)), respectively.
Additional methods such as plot,
FourierSeries2fun, and extractCoef are
provided for visualization, function evaluation, and coefficient
extraction.
fsc = Fourier_series(
double_constant = 3,
cos = c(0,2/3),
sin = c(1,7/5),
k_cos = 1:2,
k_sin = 1:2,
t_0 = 0,
period = 1
)
plot(fsc)
FourierSeries2fun(fsc,seq(0,1,0.05))
#> [1] 2.1666667 3.1712610 3.6252757 3.4344848 2.7346112 1.8333333
#> [7] 1.0888125 0.7715265 0.9623175 1.5254623 2.1666667 2.5532270
#> [13] 2.4497052 1.8164508 0.8324982 -0.1666667 -0.8133005 -0.8465074
#> [19] -0.2132530 0.9074283 2.1666667
extractCoef(fsc)
#> $a_0
#> 0
#> 3
#>
#> $a_k
#> 1 2
#> 0.0000000 0.6666667
#>
#> $b_k
#> 1 2
#> 1.0 1.4The object fsc represents the summation \[\frac32 + \frac23 \cos(2\pi\cdot2x) + \sin(2\pi
x) + \frac75\sin(2\pi\cdot2x).\]
B-spline Basis
A b-spline basis \(\{B_{i,p}(x)\}_{i=-p}^{k}\) on the interval \([t_0, t_{k+1}]\) is defined recursively as follows: - For \(p = 0\), \[ B_{i,0}(x) = \begin{cases} I_{(t_i, t_{i+1}]}(x), & i = 0, 1, \dots, k, \\ 0, & \text{otherwise}. \end{cases} \] - For \(p > 0\), the recursion is given by \[ B_{i,p}(x) = \frac{x-t_i}{t_{i+p}-t_i} B_{i,p-1}(x) + \frac{t_{i+p+1}-x}{t_{i+p+1}-t_{i+1}} B_{i+1,p-1}(x). \]
At discontinuities in the interval \((t_0,t_k)\), the basis function is defined as its limit, \(B_{i,p}(x) = \lim_{t\to x} B_{i,p}(t)\), to ensure continuity.
S4 Classes bspline_basis and
bspline_series
The s4 class bspline_basis represents a b-spline basis,
\(\{B_{i,p}(x)\}_{i=-p}^{k}\), on the
interval \([t_0, t_{k+1}]\). Its key
slots include:
Boundary.knots: The boundary \([t_0, t_{k+1}]\), (default is \([0,1]\));knots: The internal spline knots, \((t_1,\dots,t_k)\), (automatically generated as equally spaced, \(t_j = t_0 + j\cdot\frac{t_{k+1}-t_0}{k+1}\), if not assigned);intercept: Indicates whether an intercept is included (default isTRUE);df: Degrees of freedom of the basis, which is the number of the splines, equal to \(p + k + 1\) (by default \(k =0\) and \(\text{df} = p+1\));degree: The degree of the spline, which is the degree of piecewise polynomials, \(p\), (default is 3).
The s4 class bspline_series represents the linear
combination
\[ \sum_{i=0}^{k} b_i B_{i,p}(x), \]
where the slot coef stores the
coefficients \(b_i\) (\(i = 0,\dots,k\)) and
bspline_basis specifies the associated
basis (represented by a bspline_basis object).
Methods such as plot and bsplineSeries2fun
are provided for visualization and function evaluation.
bsb = bspline_basis(
Boundary.knots = c(0,24),
df = 7,
degree = 3
)
bss = bspline_series(
coef = c(2,1,3/4,2/3,7/8,5/2,19/10),
bspline_basis = bsb
)
plot(bss)
bsplineSeries2fun(bss,seq(0,24,0.5))
#> [1] 2.0000000 1.7677509 1.5690908 1.4011502 1.2610597 1.1459499 1.0529514
#> [8] 0.9791948 0.9218107 0.8779297 0.8446824 0.8191993 0.7986111 0.7805266
#> [15] 0.7644676 0.7504340 0.7384259 0.7284433 0.7204861 0.7145544 0.7106481
#> [22] 0.7087674 0.7089120 0.7110822 0.7152778 0.7216857 0.7312404 0.7450629
#> [29] 0.7642747 0.7899969 0.8233507 0.8654574 0.9174383 0.9804145 1.0555073
#> [36] 1.1438380 1.2465278 1.3634786 1.4897151 1.6190430 1.7452675 1.8621942
#> [43] 1.9636285 2.0433759 2.0952418 2.1130317 2.0905511 2.0216053 1.9000000The object bsb represents \(\{B_{i,3}(x)\}_{i=-3}^{0}\), and the object
bss represents
\[2B_{i,-3}(x)+B_{i,-2}(x)+\frac34B_{i,-1}(x)+\frac23B_{i,0}(x)
+ \frac78B_{i,1}(x) +
\frac52B_{i,2}(x) +\frac{19}{10}B_{i,3}(x),\] where \(x\in[t_0,t_4]\) and \(t_0=0\), \(t_k =
t_{k-1}+6\) ,\(k=1,2,3,4\).
Eigenfunction basis
Suppose \(K(s,t)\in L^2(\Omega\times \Omega)\), \(f(t)\in L^2(\Omega)\). Then \(K\) induces an linear operator \(\mathcal{K}\) by \[(\mathcal{K}f)(x) = \int_{\Omega} K(t,x)f(t)dt\] If \(\xi \in L^2(\Omega)\) such that \[\mathcal{K}\xi = \lambda \xi\] where \(\lambda\in \mathbb{C}\), we call \(\xi\) an eigenfunction/eigenvector of \(\mathcal{K}\), and \(\lambda\) an eigenvalue associated with \(\xi\).
For a stochastic process \(\{X(t),t\in\Omega\}\) we call the orthogonal basis, \(\{\xi_k\}_{k=1}^\infty\) corresponding to eigenvalues \(\{\lambda_k\}_{k=1}^\infty\) (\(\lambda_1\geq\lambda_2\geq\dots\)), induced by \[K(s,t)=\operatorname{Cov}(X(t),X(s))\] a functional principal component (FPC) basis for \(L^2(\Omega)\).
The eigenfunction basis relies on the covariance function \(K(s,t)\). And we usually cannot analytically express the equation of the basis funcitons in eigenfunction basis. In practice, the covariance function is estimated from the data and the eigenfunctions are computed numerically.
Numerical Basis
Sometimes, analytical expressions for the basis functions are not available. In these cases, they can be represented numerically. Let \(\{\rho_k\}_{k=1}^\infty\) denote a basis of the function space. We can numerically approximate the basis by storing the values of a finite subset of these functions at a finite set of points in the domain, i.e., \(\rho_k(t_j)\) for \(k=1,\dots,p\) and \(j=1,\dots,m\).
S4 Class numeric_basis and
numericBasis_series
In this package, the s4 class numeric_basis is designed
to represent a finite sequence of functions \(\{\rho_k\}_{k=1}^p\) by their values at a
finite set of points within their domain. In this context, all functions
share the same domain, which is assumed to be an interval. The key slots
include:
basis_function: The matrix \((\zeta_{jk})_{m \times p}\), where each element is defined as \(\zeta_{jk} = \rho_k(t_j)\) for \(j = 1, \dots, m\) and \(k = 1, \dots, p\). In this matrix, each row corresponds to a point in the domain, and each column corresponds to a specific basis function.t_points: A numeric vector that represents the points in the domain at which the basis functions are evaluated. The \(j\)-th element of this vector corresponds to the \(j\)-th row of thebasis_functionmatrix.t_0andperiod: The left endpoint and length of the domain, respectively.
Additionally, the package provides an s4 class
numericBasis_series, which represents a linear combination
of the basis functions represented by a numeric_basis
object. The key slots include:
coef: The linear coefficients for the summation series.numeric_basis: Function basis as represented by anumeric_basisobject.
Function basis2fun
The generic function basis2fun is provided for
Fourier_series, bspline_series, and
numericBasis_series objects. For a
Fourier_series object, it is equivalent to
FourierSeries2fun; for a bspline_series
object, it is equivalent to bsplineSeries2fun; For a
numericBasis_series object, it is equivalent to
numericBasisSeries2fun.
basis2fun(fsc,seq(0,1,0.05))
#> [1] 2.1666667 3.1712610 3.6252757 3.4344848 2.7346112 1.8333333
#> [7] 1.0888125 0.7715265 0.9623175 1.5254623 2.1666667 2.5532270
#> [13] 2.4497052 1.8164508 0.8324982 -0.1666667 -0.8133005 -0.8465074
#> [19] -0.2132530 0.9074283 2.1666667
basis2fun(bss,seq(0,24,0.5))
#> [1] 2.0000000 1.7677509 1.5690908 1.4011502 1.2610597 1.1459499 1.0529514
#> [8] 0.9791948 0.9218107 0.8779297 0.8446824 0.8191993 0.7986111 0.7805266
#> [15] 0.7644676 0.7504340 0.7384259 0.7284433 0.7204861 0.7145544 0.7106481
#> [22] 0.7087674 0.7089120 0.7110822 0.7152778 0.7216857 0.7312404 0.7450629
#> [29] 0.7642747 0.7899969 0.8233507 0.8654574 0.9174383 0.9804145 1.0555073
#> [36] 1.1438380 1.2465278 1.3634786 1.4897151 1.6190430 1.7452675 1.8621942
#> [43] 1.9636285 2.0433759 2.0952418 2.1130317 2.0905511 2.0216053 1.9000000Basis expansion methods for the functional_variable
class
The MECfda pakcage provide the methods
fourier_basis_expansion,
bspline_basis_expansion, FPC_basis_expansion,
and numeric_basis_expansion for the class
functional_variable, which perform basis expansion using
Fourier basis, b-spline basis, functional principal component basis, and
numerical basis respectively.
Numerical Computation of Integrals
The package approximates the integral
\[ \int_{\Omega} \rho_{k}(t) X_{i}(t) \, dt \approx \frac{\mu(\Omega)}{|T|} \sum_{t \in T} \rho_{k}(t) X_{i}(t), \]
where \(\mu(\cdot)\) denote the Lebesgue measure. \(T\) denotes the set of measurement time points for \(X_i(t)\) and \(|T|\) represents the cardinal number of \(T\).
