From a course by Prof. Camelia Goga
April 24, 2024
library(dda)
library(tidyverse)
library(ggplot2)
data(vietnam_temperature)
vntemp <- vietnam_temperature |>
filter(year == 2016) |>
mutate(t_max = lapply(t_max, unclass), id = seq_len(n())) |>
unnest_longer(t_max, indices_to = "day")
vntemp |> ggplot(aes(day, t_max, group = id, color = region)) +
geom_point() +
geom_line() +
xlab("Day number") +
ylab("Temperature (°C)")Figure 1: There are more discretisation points (365 or 366) than observations (63 provinces).
Figure 2: The discretisation points are not equally spaced and different for each curve. Source: CONTENT study, Sixth Framework Programme of the European Union.
Figure 3: The x-axis is not time.
We observe a sample of random functions \(Y_1, \ldots, Y_n:\) \[ Y_i:\Omega\longrightarrow \mathcal F \] and we consider that \(Y_i, i=1, \ldots, n\) are iid.
We always observe the functions \(Y_i, i=1, \ldots, n\) in a finite (often large) number of discretisation points (which may be different from one curve to another): \[t_{i1}, \ldots, t_{iD_i}\subset [0, T]\] so the data is presented in vector form.
To overcome these difficulties,
In functional statistics, the objects of interest are functions: \[ Y=(Y(t))_{t \in [0, T]} \]
Therefore, in theory, \(Y\) could be measured at any time \(t \in [0, T]\).
Nevertheless, in practice:
To overcome these difficulties, the true curve \(Y=(Y(t))_{t\in [0,T]}\) can be “reconstructed” from the discretised data using non-parametric smoothing methods.
library(fda)
y <- filter(vietnam_temperature, year == 2016 & province == "HANOI") |>
pull(t_max) |>
unlist()
breaks <- seq(1, length(y), length.out = 10)
bsp <- create.bspline.basis(breaks = breaks)
yhat <- smooth.basis(seq_len(length(y)), y, bsp)
invisible(
plot(yhat, xlab = "Day number", ylab = "Temperature (°C)")
)
points(seq_len(length(y)), y)Figure 4: Smoothing of daily maximum temperature in 2016 in Hanoi.
We can use non-parametric methods like
to give an estimate \(\widehat f(t_d)\) of \(f(t)\);
The re-constructed curve is \(Y(t)=\widehat f(t)\) for all \(t \in [0, T]\).
Let \(b_1(\cdot), \ldots, b_q(\cdot)\) be a basis of spline functions (\(B\)-splines or truncated polynomials) of dimension \(q=K+m\) with \(K\) equidistant interior knots and piecewise polynomials of degree \(m-1\);
We denote by \(\mathbf{b}^T(t)=(\textit{b}_1(t), \ldots, b_q(t))\) the function vectors of the basis calculated at a point \(t\in [0, T]\) and by \(\mathbf{B}\) the matrix containing the values of \(\mathbf{b}^T(t)\) for the discretisation points \(t_d, d=1, \ldots, D\): \[ \mathbf{B}=\left(\begin{array}{c} \mathbf{b}^T(t_1)\\ \mathbf{b}^T(t_2)\\ \ldots\\ \mathbf{b}^T(t_D) \end{array}\right)=\left(\begin{array}{cccc} b_1(t_1) \quad \ldots \quad b_q(t_1)\\ b_1(t_2) \quad \ldots \quad b_q(t_2)\\ \ldots \quad \ldots \quad \ldots\\ b_1(t_D) \quad \ldots \quad b_q(t_D) \end{array}\right) \]
The estimator \(\widehat f(t_d)\) of \(f(t_d)\) is given by \[\begin{aligned} \widehat f(t_d) &= \mathbf{b}^T(t_d)\boldsymbol{\widehat{\theta}} \quad d=1, \ldots, D\quad \mbox{with}\\ \boldsymbol{\widehat{\theta}} &= \underset{\boldsymbol{\theta}\in \mathbb R^q}{\operatorname{argmin}} \sum_{d=1}^D(y_d-\mathbf{b}^T(t_d)\boldsymbol{\theta})^2=(\mathbf{B}^T\mathbf{B})^{-1}\mathbf{B}^T\mathbf{Y}_d \end{aligned}\] where \(\mathbf{Y}_d=(Y_d)_{d=1}^D\).
In the end, we have the curve reconstructed by non-parametric regression using spline functions: \[ \widehat Y(t)=\widehat f(t)=\mathbf{b}^T(t)\boldsymbol{\widehat\theta}=\mathbf{a}^T(t)\mathbf{Y}_d \quad \text{for all}\quad t\in [0, T]; \]
The number of knots \(K\) can be chosen using the cross-validation or generalised cross-validation criterion and \(m\) is small (3 or 4). The knots can be positioned at the quantiles of the instants \(t_1, \dots, t_D\).
We have a sample of \(n\) curves, \(Y_1, \ldots, Y_n\) which are, for simplicity, measured at the same discretisation times \(t_1, \ldots, t_D;\)
Let \(Y_i(t_d)=Y_{id}\) for \(i=1, \ldots, n\) and \(d=1, \ldots, D\).
Let \(n\) be non-parametric models: \[\begin{aligned} Y_{1d} & = f_1(t_d)+\varepsilon_{1d}, \quad d=1, \ldots, D \\ Y_{2d} & = f_2(t_d)+\varepsilon_{2d}, \quad d=1, \ldots, D \\ &\vdots \\ Y_{nd} & = f_n(t_d)+\varepsilon_{nd}, \quad d=1, \ldots, D \end{aligned}\]
The objective is to estimate \(f_i, i=1, \ldots, n\);
Let \(b_1(\cdot), \ldots, b_q(\cdot)\) be the basis of spline functions of order \(q=K+m\): \(K\) equidistant interior knots of degree \(m-1;\)
This same basis is used to re-construct all the \(n\) curves;
Using the method described for smoothing a single curve, we obtain: \[\begin{aligned} \widehat Y_i(t)= \widehat f_i(t) &= \mathbf{b}^T(t)\widehat{\boldsymbol{\theta}}_i \\ &= \mathbf{b}^T(t)(\mathbf{B}^T\mathbf{B})^{-1}\mathbf{B}^T\mathbf{Y}_i, \quad i=1, \ldots, n \end{aligned}\] where \(\mathbf{Y}_i=(Y_{id})_{d=1}^D;\)
The same method can be used in the case where the discretisation times are not the same for all the curves. In this situation, we consider the union of all the measurement instants and the knots are placed at the quantiles of this new set.
The variables \(Y_1, \ldots, Y_n\) are now functions of \(t\in [0, T]:\) \[ Y_i=(Y_i(t))_{t\in [0, T]}, \] and we consider that they are independent and identically distributed;
We consider that the \(Y_i, i=1, \ldots, n\) belong to the Hilbert space \(L^2[0, T]\), the space of integrable square functions with the usual scalar product and the associated norm: \[ \langle f,g \rangle = \int_0^Tf(t)g(t)dt, \quad \|f \| =\sqrt{\int_0^Tf^2(t) \mathrm dt} \]
The processes (functions) \(Y_i\) are defined on a probability space \((\Omega, \mathcal K, \mathbb P)\) with values in \(L^2[0, T]:\) \[ Y_i(\omega)=(Y_i(t,\omega))_{t\in [0, T]}, \quad i=1, \ldots, n \]
Remark. This is not the Lebesgue integral but rather the Bochner integral (its generalisation to Banach spaces).
Let \(Y_i, i=1, \ldots, n\) be independent and identically distributed with \(Y\).
Figure 5: Empirical mean of the NIR curves, computed pointwise.
Under the preceding conditions: \[ \mathbb E(\langle Y, f \rangle)=\langle \mathbb E(Y), f\rangle \quad \mbox{for all}\quad f\in L^2[0, T]; \] so the integral and the scalar product commute.
“Pythagorean theorem”: \[ \mathbb E\|Y-\mu\|^2=\mathbb E\|Y\|^2-\mu^2 \]
The \(\Gamma\) operator has the following properties (assuming that \(\mathbb E(\|Y\|^2)<\infty\)):
\(\displaystyle \langle \Gamma f, g\rangle=\mathbb E\left(\langle Y-\mu, f\rangle\langle Y-\mu, g\rangle\right)\) for all \(f, g\in L^2[0, T];\)
\(\Gamma=\mathbb E (Y\otimes Y)-\mu\otimes\mu\)
it is a self-adjoint operator (\(\Gamma=\Gamma^*\), where \(\Gamma^*\) is the adjoint operator) because the kernel \(\gamma(\cdot, \cdot)\) is symmetric and non-negative (\(\langle \Gamma u, u \rangle \geq 0\));
it has finite trace: \(\|\Gamma\|_{TR}=\sum_{j\geq 1}\langle (\Gamma^*\Gamma)^{1/2} e_j, e_j\rangle=\mathbb E\|Y\|^2<\infty\);
is a Hilbert-Schmidt operator: \[|\Gamma\|^2_{HS}=\sum_{j\geq 1}\|\Gamma e_j\|^2<\infty,\]
where \(\{e_j\}_{j\geq 1}\) is any orthonormal basis of \(L^2[0, T]\).
Let \((\lambda_j, v_j)_{j\leq 1}\) with \(\lambda_j\in \mathbb R\) and \(v_j=(v_j(t))_{t\in [0, T]}\in L^2[0, T]\) be called eigenvalues and associated eigenvectors (functions) of \(\Gamma\): \[\begin{aligned} \Gamma v_j &= \lambda_j \ v_j\\ \Gamma v_j (t)& = \lambda_j \ v_j(t), \quad t \in [0,T], \quad j\geq 1 \end{aligned}\]
\(\Gamma\) is non-negative so \(\lambda_j\geq 0\) for all \(j\geq 1;\)
\(\Gamma\) is self-adjoint, so the eigenvectors associated with the distinct eigenvalues are linearly independent and orthogonal;
The operator \(\Gamma\) is compact and self-adjoint so the set of non-zero eigenvalues is finite or is a sequence which converges to zero. Let \((\lambda_j)_{j\geq 1}\) be the ordered eigenvalues: \(\lambda_1\geq\lambda_2\geq\ldots>0\). The associated eigenvectors \((v_j)_{j\geq 1}\) form an orthonormal basis in \(L^2[0,T]\) for \(\overline{\operatorname{Im} (\gamma)}\) i.e \(\langle v_j, v_{j'} \rangle = \delta_{jj'}\). The covariance then has the following decomposition: \[ \Gamma=\sum_{j\geq 1}\lambda_jv_j\otimes v_j \quad \mbox{et}\quad \Gamma f=\sum_{j\geq 1}\lambda_j \langle f,v_j\rangle v_j; \]
The covariance function \(\gamma\) associated with \(\Gamma\) has the following decomposition: \[ \gamma(s,t)=\sum_{j\geq 1}\lambda_jv_j(s) v_j(t) \quad \mbox{for any}\quad s,t\in [0, T]; \]
\(Y-\mu=\sum_{j\geq 1}\langle Y-\mu, v_j\rangle v_j\) and \(\mathbb E(\langle Y-\mu, v_j\rangle v_j)=0;\)
The trace norm of \(\Gamma\) is given by \[ \|\Gamma\|_{TR}=\sum_{j\geq 1}\langle \Gamma v_j, v_j\rangle=\sum_{j\geq 1} \lambda_j<\infty \] and the Hilbert-Schmidt norm by \[ \|\Gamma\|_{HS}=\sum_{j\geq 1} \lambda^2_j<\infty. \]
\(\Gamma^{1/2}=\sum_{j\geq 1}\lambda_j^{1/2}v_j\otimes v_j\) and \(\Gamma^p=\sum_{j\geq 1}\lambda_j^pv_j\otimes v_j.\)
Therefore, the optimal subspace of dimension \(q,\) (unique if \(\lambda_j\) is distinct), is the space generated by the \(q\) eigenvectors associated to the \(q\) largest eigenvalues of \(\Gamma;\)%Thus the optimal subspace of dimension \(q,\) which is unique if \(\lambda_q > \lambda_{q+1},\) is the space generated by the \(q\) eigenfunctions of \(\Gamma\) associated to the \(q\) largest eigenvalues. \[ Y_k(t) = \mu(t) + \sum_{j =1}^q \langle Y_k - \mu,v_j \rangle v_j(t) + R_{q,k}(t), \quad t \in [0,T] \]
The space generated by \(v_1, \cdots, v_q\) gives a representation of the variation over time of the data around the mean.
The variance of the projection on \(v_j\) is given by: \[ \lambda_j = \frac{1}{n} \sum_{k=1}^n \langle Y_k - \mu,v_j \rangle^2 \ \] because \(\frac{1}{n} \sum_{k=1}^n \langle Y_k - \mu,v_j \rangle=0.\)
The notion of quantile is not generalisable to \(\mathbb R^p\) because we cannot give a total order relation in \(\mathbb R^p\); the same thing in a more general space like \(L^2[0, T];\)
Several definitions have been proposed of a concept that could be considered a quantile but none has the basic property of a quantile of order \(\alpha\) in \(R,\) having \(\alpha\%\) of the distribution less than this value;
A definition for the median \[ m = \underset{y\in L^2[0,T]}{\operatorname{argmin}} \sum _{i=1}^n \|Y_i - y\|; \]
Let \(X_1, \ldots, X_p\) be real variables and \(\mathbf{x}'_i=(1,X_{i1}, \ldots, X_{ip});\)
The data are \((\mathbf{x}'_i, Y_i(t))\) for \(i=1, \ldots, n;\)
The following linear model (Faraway 1997) \[ M1:\quad Y_i(t)=\mathbf{x}'_i\boldsymbol{\beta}(t)+\varepsilon_{i}(t), \quad t\in [0,T] \] where \(\boldsymbol{\beta}(t)=(\beta_0(t), \beta_1(t), \ldots, \beta_p(t))'\) is the functional regression coefficient and the functional residuals \(\varepsilon_{k}(t)\) are centred : \[ \mathbb E(\varepsilon_{i}(t))=0 \] and decorrelated with each the same covariance function \(\gamma\): \[ \operatorname{cov} (\varepsilon_{i}(t),\varepsilon_{j}(t))=0, \quad i\neq j, t\in [0,T] \] \[ \operatorname{cov} (\varepsilon_{i}(t),\varepsilon_{i}(r))=\gamma(t,r), \quad t,r\in [0,T] \]
It is assumed that the matrix \(\mathbf{X}=(\mathbf{x}'_i)_{i=1}^n\) of size \((n,p+1)\) is of rank \(p+1;\)
The pointwise estimator of \(\boldsymbol{\beta}\) (for \(t\) fixed) is obtained by ordinary least squares: \[\begin{aligned} \widehat{\boldsymbol{\beta}}(t) & = \underset{\beta(t)}{\operatorname{argmin}} \sum_{i=1}^n(Y_i(t)-\mathbf{x}'_i\boldsymbol{\beta}(t))^2=\|\mathbf{Y}(t)-\mathbf{X}\boldsymbol{\beta}(t)\|^2 \\ & = (\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}\mathbf{Y}(t) \end{aligned}\] where \(\mathbf{Y}(t)=(Y_i(t))_{i=1}^n;\)
It can be shown that the estimator \(\widehat{\boldsymbol{\beta}}=(\widehat{\boldsymbol{\beta}}(t))_{t\in [0, T]}\) satisfies the global criterion: \[ \widehat{\boldsymbol{\beta}}=\mbox{argmin}_{\beta}\int_{0}^T\|\mathbf{Y}(t)-\mathbf{X}\boldsymbol{\beta}(t)\|^2dt \]
The predictions are \(\widehat{\mathbf{Y}}(t)=\mathbf{X}\widehat{\boldsymbol{\beta}}(t)=\mathbf{P}_X\mathbf{Y}(t);\)
The estimated residuals are \(\widehat{\boldsymbol{\varepsilon}}(t)=\mathbf{Y}(t)-\widehat{\mathbf{Y}}(t)=\left(\mathbf{I}_n-\mathbf{P}_X\right)\mathbf{Y}(t):\) \[ \mathbb E(\widehat{\boldsymbol{\varepsilon}}(t))=0 \] \[ \operatorname{cov}(\widehat{\boldsymbol{\varepsilon}}(t), \widehat{\boldsymbol{\varepsilon}}(r))=\gamma(t,r)\left(\mathbf{I}_n-\mathbf{P}_X\right) \] \[ \mathbb E(\widehat{\boldsymbol{\varepsilon}}(t)'\widehat{\boldsymbol{\varepsilon}}(r))=(n-p-1)\gamma(t,r) \] so an unbiased estimator of \(\gamma(t,r)\) is given by \[ \widehat{\gamma}(t,r)=\frac{1}{n-p-1}\sum_{i=1}^n(Y_i(t)-\widehat Y_i(t))(Y_i(r)-\widehat Y_i(r)) \]
R
Figure 7: Implementation of scalar-on-function regression with refund::pfr.
Let \(X_1, \ldots, X_p\) be functional variables, \(X_j=(X_j(t))_{t\in [0, T]};\) let \(\mathbf{x}'_i(t)=(1,X_{i1}(t), \ldots, X_{ip}(t));\)
The data are \((\mathbf{x}'_i(t), Y_i(t))\) for \(i=1, \ldots, n;\) \[ M1:\quad Y_i(t)=\mathbf{x}'_i(t)\boldsymbol{\beta}(t)+\varepsilon_{i}(t), \quad t\in [0,T] \] where \(\boldsymbol{\beta}(t)=(\beta_0(t), \beta_1(t), \ldots, \beta_p(t))'\) is the functional regression coefficient and the functional residuals \(\varepsilon_{k}(t)\) are centred: \[ \mathbb E(\varepsilon_{i}(t))=0 \] and decorrelated with each the same covariance function \(\gamma\): \[ \operatorname{cov} (\varepsilon_{i}(t),\varepsilon_{j}(t))=0, \quad i\neq j, t\in [0,T] \] \[ \operatorname{cov} (\varepsilon_{i}(t),\varepsilon_{i}(r))=\gamma(t,r), \quad t,r\in [0,T] \]
It is assumed that the matrix \(\mathbf{X}(t)=(\mathbf{x}'_i(t))_{i=1}^n\) of size \((n,p+1)\) is of rank \(p+1\) for each \(t\in [0, T];\)
The point estimator of \(\boldsymbol{\beta}\) (for \(t\) fixed) is obtained by ordinary least squares: \[\begin{aligned} \widehat{\boldsymbol{\beta}}(t) & = \underset{\beta(t)}{\operatorname{argmin}} \sum_{i=1}^n(Y_i(t)-\mathbf{x}'_i(t)\boldsymbol{\beta}(t))^2=\|\mathbf{Y}(t)-\mathbf{X}(t)\boldsymbol{\beta}(t)\|^2\\ & = (\mathbf{X(t)}'\mathbf{X}(t))^{-1}\mathbf{X}(t)\mathbf{Y}(t) \end{aligned}\] where \(\mathbf{Y}(t)=(Y_i(t))_{i=1}^n;\)
It can be shown that the estimator \(\widehat{\boldsymbol{\beta}}=(\widehat{\boldsymbol{\beta}}(t))_{t\in [0, T]}\) satisfies the global criterion: \[ \widehat{\boldsymbol{\beta}}= \underset{\beta}{\operatorname{argmin}} \int_{0}^T\|\mathbf{Y}(t)-\mathbf{X}(t)\boldsymbol{\beta}(t)\|^2 \mathrm dt \]
