Functional Regression in R: Functional Response Regressed on Scalar Covariate

Question

I have a functional response $Y(t)$ (i.e., a stochastic process) which I regress on a set of scalar explanatory variables $X_1$ ($=1$, i.e., $X_1$ is a constant term), $X_2$, and $X_3$. The equation I am considering is thus $$Y(t) = \alpha(t) + \beta(t)X_2 + \gamma(t)X_3 + \varepsilon(t).$$

I am using the R package fda. I know how to run the regression:

set.seed(313) # seed

t = 100 # time points
n = 50 # sample size

# generate some mock data
Y = sweep(matrix(rnorm(n * t, 0, 1), ncol = 50), 1, sin(1:100)) 
X2 = sample(c(0, 1), size = n, replace = TRUE)
X3 = rnorm(n)

# generate a basis for Y
Yfd = fda::smooth.basis(1:t, Y, fda::create.bspline.basis(c(1, 100), 2, 1))$fd

# run the regression
ffit = fda::fRegress(Yfd ~ X2 + X3)

That is, $Y(t)$ is generated by $\sin(t) + \nu(t)$, where $\nu(t)$ is a Gaussian White Noise process. The process is observed at grid points $\{1,2,\dots, 100\}$. $X_2$ is sampled from a Laplace distribution over $\{0,1\}$, and $X_3$ is sampled from a standard Normal distribution. The interface of fda::fRegress automatically implements a constant, so there is no need to define $X_1$.

Using plot(fda::predict.fRegress(ffit)) gives me plots of the fitted values. However, I want to to compute the predictions by hand. How to do that? The documentation is, unfortunately, not very helpful.

My idea was to compute $$\boldsymbol\theta(t) = (\boldsymbol X'\boldsymbol X)^{-1}\boldsymbol X'\boldsymbol y(t)$$ for each $t$, where $$\boldsymbol\theta(t) = \begin{pmatrix}\alpha(t)\\ \beta(t) \\ \gamma(t) \end{pmatrix}, \qquad\boldsymbol X = \begin{pmatrix} 1 & X_{21} & X_{31} \\ 1 & X_{22} & X_{32} \\ \vdots & \vdots & \vdots \\ 1 & X_{2n} & X_{3n} \end{pmatrix}\qquad\text{and}\qquad\boldsymbol y(t) = \begin{pmatrix} Y_1(t) \\ Y_2(t) \\ \vdots \\ Y_n(t) \end{pmatrix}.$$ Is this the correct way to do it? I don't think so, because I have never used the basis specified for the functional response...

Answer 1

I went through the code of fda::predict.fRegress and found a solution to my problem. There is a multiplication operation fda::times.fd and addition operation fda::plus.fd which allows multiplication and addition of fda::fda objects, respectively. The code below computes predictions manually:

# create data
xfdlist = list(
  X1 = fda::fd(t(1), fda::create.const.basis(c(1, 100))),
  X2 = fda::fd(t(1), fda::create.constant.basis(c(1,100))),
  X3 = fda::fd(t(0.314), fda::create.constant.basis(c(0, 11)))
)

# extract slope functions
bfdlist = ffit$betaestlist

# compute predictions
yhat = purrr::reduce(purrr::map2(xfdlist, bfdlist, \(x, b) fda::times.fd(x, b$fd)), fda::plus.fd)

in Symbols: $$\hat y(t) = \left(\hat\alpha(t)\otimes 1\right) \oplus \left(\hat\beta(t) \otimes 1\right) \oplus \left(\hat\gamma(t)\otimes 0.314\right),$$

where $\oplus$ = fda::plus.fd and $\otimes$ = fda::times.fd. Evaluating this equation at particular values for $t$ can be done by using fda::eval.fd.

It would still be interesting to learn more about what the operations $\oplus$ and $\otimes$ and how they can be spelled out mathematically.