Intro2FDA.Rmd

---
title: "A Gentle Introduction to Functional Data Analysis"
author: "Dominik Liebl"
output:
  ioslides_presentation: default
  slidy_presentation: default
  widescreen: true
  smaller: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
# install.packages(c("fda","png"))
## Load packages
library("fda")
data(growth)
```


## How to Use these Slides

Accompanying Textbook: 
<br><br>
```{r, fig.height=4, fig.align='center'}
library("png")
library("grid")
FDA_Book_Pic <- png::readPNG("./images/IntroFDA_Book.png")
grid::grid.raster(FDA_Book_Pic)
```


## How to Use these Slides

The **HTML-slides** of this short lecture and the **R-codes** for reproducing the examples can be found at my GitHub repository **Intro2FDA**:
<br><br>
<div class="centered">
<bdi style="font-size:150%;"><b>https://github.com/lidom/Intro2FDA</b></bdi>
</div>


##
<br><br><br><br><br>
<div class="red2">
<div class="centered">
<bdi style="font-size:350%;"><b>Functional Data</b></bdi>
</div>
</div>


## Functional Data
```{r FD1, out.width = "100%"}
knitr::include_graphics("./images/boygrowth.png")
```

## Functional Data
```{r FD2, out.width = "100%"}
knitr::include_graphics("./images/Fig_POI.png")
```

## Functional Data
```{r FD4, out.width = "100%"}
knitr::include_graphics("./images/Fig_Electr.png")
```

## Functional Data
```{r FD5, out.width = "100%"}
knitr::include_graphics("./images/Fig_PM10_1.png")
```

## Functional Data
```{r FD6, out.width = "100%"}
knitr::include_graphics("./images/Fig_PM10_2.png")
```

##
<br><br><br><br>
<div class="red2">
<div class="centered">
<bdi style="font-size:350%;"><b>From Raw to Functional Data</b></bdi>
</div>
</div>

## From Raw to Functional Data

The simplest <span style="color:#FF4136">**data**</span> set encountered in functional data analysis is a sample of the form:
$$
\begin{align}
x_i(t_j),\quad t_j\in[a,b],\quad i=1,\dots,n \quad j=1\dots,J.
\end{align}
$$
The <span style="color:#FF4136">**theoretical objects**</span> we wish to study are smooth curves:
$$
\begin{align}
X_i(t),\quad t\in[a,b],\quad i=1,\dots,n \quad j=1\dots,J.
\end{align}
$$
```{r, fig.width=8, fig.height=3, fig.align='center'}
par(mfrow=c(1,2), mar=c(5.1, 4.1, 2.1, 2.1))
matplot(x=growth$age, y=growth$hgtf[,1:5], type="p", lty=1, pch=21, 
        xlab="age", ylab="height (cm)", cex.lab=1.2,col=1:5,bg=1:5,
        main=expression(x[i](t[j])))
matplot(x=growth$age, y=growth$hgtf[,1:5], type="l", lty=1, pch=1, 
        xlab="age", ylab="height (cm)", cex.lab=1.2,
        main=expression(X[i](t)))
par(mfrow=c(1,1), mar=c(5.1, 4.1, 4.1, 2.1))
```


## From Raw to Functional Data

**Basis expansions:**

Typically, the first step in working with functional data is
to express them by means of a <span style="color:#FF4136">**basis expansion**</span>
$$
X_i(t)\approx\sum_{m=1}^M c_{im} B_m(t),\quad t\in[a,b],
$$
where $B_1(t),\dots,B_M(t)$ are some standard collection of basis functions like: 

- splines
- wavelets
- sine and cosine functions
- etc.


## From Raw to Functional Data

<span style="color:#FF4136">**B-spline**</span> basis functions $B_1(t),\dots,B_M(t)$:

```{r, echo=TRUE, out.width = "80%", fig.align='center'}
bspl_bf <- create.bspline.basis(rangeval=c(0,1), norder=3, 
                                     breaks=seq(0,1,len=7))
plot(bspl_bf, main="B-spline Basis Functions", xlab="[a,b]=[0,1]")
```


## From Raw to Functional Data

<span style="color:#FF4136">**Fourier basis**</span> functions $B_1(t),\dots,B_M(t)$:

```{r, echo=TRUE, out.width = "80%", fig.align='center'}
fourier_bf <- create.fourier.basis(rangeval=c(0,1), nbasis=5)
plot(fourier_bf, main="Fourier Basis Functions", xlab="[a,b]=[0,1]")
```


## From Raw to Functional Data

<span style="color:#FF4136">**Example:**</span>  Berkeley growth study **raw data**
```{r, fig.align='center', out.width = "80%", echo=TRUE}
matplot(x=growth$age, y=growth$hgtf[,1:5], type="p", lty=1, pch=21, 
        xlab="age", ylab="height (cm)", cex.lab=1.2,col=1:5,bg=1:5,
        main="5 Girls in Berkeley Growth Study")
```


## From Raw to Functional Data

<span style="color:#FF4136">**Example:**</span>  Berkeley growth study **pre-processed data**
```{r, fig.align='center', out.width = "80%", echo=TRUE}
SmObj <- smooth.basisPar(growth$age,y=growth$hgtf[,1:5],lam=0.1)
plot(SmObj$fd, xlab="age", ylab="height (cm)", cex.lab=1.2,
     main="5 Girls in Berkeley Growth Study", lty=1)
```


## From Raw to Functional Data

<span style="color:#FF4136">**Example:**</span>  Berkeley growth study **1st derivative**
```{r, fig.align='center', out.width = "80%", echo=TRUE}
plot(deriv(SmObj$fd, 1), xlab="age", 
     ylab="growth rate (cm / year)", cex.lab=1.2,
     main="5 Girls in Berkeley Growth Study (1st derivative)", lty=1)
```


## From Raw to Functional Data

<span style="color:#FF4136">**Example:**</span>   Berkeley growth study **2nd derivative**
```{r, fig.align='center', out.width = "80%", echo=TRUE}
plot(deriv(SmObj$fd, 2), xlab="age", cex.lab=1.2,
     ylab="growth acceleration (cm / year^2)",
     main="5 Girls in Berkeley Growth Study (2nd derivative)", lty=1)
```

##
<br><br><br><br>
<div class="red2">
<div class="centered">
<bdi style="font-size:350%;"><b>Sample Mean and Covariance</b></bdi>
</div>
</div>


## Sample Mean and Covariance

<!-- <span style="color:#FF4136">**Pointwise mean:**</span>  -->
**Pointwise mean:**
$$
\bar{X}_n(t)=\frac{1}{n}\sum_{i=1}^n X_i(t)
$$
**Pointwise standard deviation:**
$$
S_n(t)=\sqrt{\frac{1}{n-1}\sum_{i=1}^n\Big(X_i(t)-\bar{X}_n(t)\Big)^2}
$$

## Sample Mean and Covariance

**Pointwise mean and standard deviation:**

```{r, eval=FALSE, echo=TRUE}
BM.Mean <- rowMeans(BM.mat)
BM.SD   <- apply(BM.mat, 1, sd)
```

```{r, echo=FALSE, fig.align='center', fig.width=8, fig.height=4}
set.seed(109)
n       <- 50
J       <- 500
BM.mat  <- matrix(0, ncol=n, nrow=J)
for(i in 1:n){BM.mat[,i] <- cumsum(rnorm(J, sd = 1/100))}
BM.Mean <- rowMeans(BM.mat)
BM.SD   <- apply(BM.mat,1,sd)
xx          <- seq(0,1,len=J) 
par(mfrow=c(1,2))
matplot(x=xx, y=BM.mat, xlab="", ylab="", type="l", col=gray(.7), lty=1, main="Brownian Motions")
lines(x=xx, y=BM.Mean)
lines(x=xx, y=BM.SD, lty=2)
legend("topleft", lty = c(1,2), legend = c("Sample Mean","Sample SD"))
matplot(x=xx, y=BM.mat, xlab="", ylab="", type="l", col=gray(.7), lty=1, main="Brownian Motion")
lines(x=c(0,1), y=c(0,0), lty=1)
lines(x=c(0,1), y=c(0,sqrt(J*(1/100)^2)), lty=2)
legend("topleft", lty = c(1,2), legend = c("True Mean","True SD"))
par(mfrow=c(1,1))
```

## Sample Mean and Covariance

**Pointwise covariance function:** 
$$
\hat{c}_n(t,s)=\frac{1}{n-1}\sum_{i=1}^n\Big(X_i(t)-\bar{X}(t)\Big)\Big(X_i(s)-\bar{X}(s)\Big)
$$


- The sample covariance function is extensively used in Functional Data Analysis (FDA). 
- The interpretation of the values of $\hat{c}_n(t,s)$ is the same as for the usual variance-covariance matrix. 


## Sample Mean and Covariance

**Pointwise covariance function:** 
```{r, eval=FALSE, echo=TRUE}
BM.Cov <- var(BM.mat)
```


```{r, echo=FALSE, fig.align='center', out.width="100%"}
BM.cov <- var(t(BM.mat))
slct   <- c(seq.int(1,500,by=20),500)
par(mfrow=c(1,2), mar=c(1, 1.1, 1.2, 1.1))
persp(xx[slct], xx[slct], BM.cov[slct,slct], xlab="s", ylab="t", zlab="",main="Sample Covariance Function",
      theta = -40, phi = 20, expand = .75, col = "blue", shade = 1.05)
x <- seq(0, 1, length= 30); y <- x
f <- function(x, y){min(x,y)}
f <- Vectorize(f)
z <- outer(x, y, f)
persp(x, y, z, xlab="s", ylab="t", zlab="",
      main="True Covariance Function",
      theta = -40, phi = 20, expand = .75, col = "blue", shade = 1.05)
par(mfrow=c(1,1), mar=c(5.1, 4.1, 4.1, 2.1))
```

##
<br><br><br>
<div class="red2">
<div class="centered">
<bdi style="font-size:350%;"><b>Principal Component Functions
</b></bdi>
</div>
</div>


## Principal Component Functions

- **Idea:** Use the Estimated Functional Principal Components (EFPC's) $\hat{v}_j$ as <span style="color:#FF4136">**basis functions**</span> for $X_i$:
$$
X_i(t)\approx\bar{X}_n(t) + \sum_{j=1}^p\hat{\xi}_{ij}\hat{v}_j(t)
$$

- <span style="color:#FF4136">**Estimated scores**</span>: $\hat{\xi}_{ij}=\int_a^b (X_i(t)-\bar{X}_n(t))\hat{v}_j(t)dt$
- EFPC's $\hat{v}_j$ are <span style="color:#FF4136">**orthonormal**</span>, i.e.
$$\textstyle\int_a^b\hat{v}_j(t)\hat{v}_k(t)dt=\begin{cases}1, &j=k\\0,&j\neq k.\end{cases}$$
- <span style="color:#FF4136">**Best basis property**</span>

## Principal Component Functions 

Eigendecomposition of $\hat{c}_n(t,s)$: 
$$
\begin{align}
\hat{c}_n(t,s)&\approx\sum_{j=1}^p\hat{\lambda}_j\hat{v}_j(t)\hat{v}_j(s),
\end{align}
$$
where $\hat{\lambda}_1>\dots >\hat{\lambda}_p$ denote the estimated eigenvalues.
```{r echo=FALSE, fig.align='center', fig.height=3}
BSPL.basis <- create.bspline.basis(rangeval=c(0,1), nbasis=200)
BM.fd      <- smooth.basis(argvals=xx, y=BM.mat, fdParobj=BSPL.basis)
BM.pca     <- pca.fd(BM.fd$fd, nharm=3)
par(mfrow=c(1,3), mar=c(2.1, 1.1, 4.1, 1.1))
persp(xx[slct], xx[slct], BM.cov[slct,slct], xlab="s", ylab="t", zlab="",
      main="Sample Covariance Function", theta = -40, phi = 20, expand = .75, col = "blue", shade = 1.05)
invisible(plot(BM.pca$values[1:3], type="o", ylab="", main="Estimated Eigenvalues (p=3)"))
invisible(plot(BM.pca$harmonics, lwd=3, ylab="", main="Estimated FPC's (p=3)"))
par(mfrow=c(1,1))
```





## Principal Component Functions

$$
\begin{align}
X_i(t)&\approx\bar{X}_n(t) + \sum_{j=1}^p\hat{\xi}_{ij}\hat{v}_j(t)
\end{align}
$$
```{r echo=TRUE, echo=FALSE, fig.align='center', fig.height=4}
BM.pca <- pca.fd(BM.fd$fd, nharm=15)
v_hat_mat  <- eval.fd(xx, BM.pca$harmonics)
xi_hat_mat <- BM.pca$scores

# FPCA-approx
X_fpca_fit <- BM.Mean + (xi_hat_mat %*% t(v_hat_mat)) 

# plot
par(mfrow=c(1,1), mar=c(5.1, 4.1, 2.1, 2.1))
invisible(plot(BM.fd$fd[1], lwd=1.3, main="Approximation of Brownian Motion (p=15)"))
lines(y=X_fpca_fit[1,], x=xx, col="red", lwd=1.3)
legend("topright", lty = c(1,1), col=c("black","red"), legend = c("True","Approx"))
```




##
<br><br><br>
<div class="red2">
<div class="centered">
<bdi style="font-size:350%;"><b>Case Study:<br> BOA Stock Returns
</b></bdi>
</div>
</div>

## Case Study: BOA Stock Returns

- Bank of America (BOA) stock values $P_i(t_j)$
    - $i=1,\dots,n=2511$ denotes the **trading days** <br>
      (from April 9th, 1997 to April 2nd, 2007).
    - $t_j\in[0,6.5]$ denotes **intra-day trading time** <br>
      (6.5 hours / per trading day). 
    - Minute-wise measurements: $j=1,\dots,J=390$.


We study the daily **cumulative log-return** functions:
$$R_i(t):=\log(P_i(t))-\log(P_i(0))\approx\frac{P_i(t)-P_i(0)}{P_i(0)}$$


## Case Study: BOA Stock Returns

```{r, echo=FALSE, eval=TRUE}
# Data
#BOA        <- read.csv(file="https://raw.githubusercontent.com/lidom/Teaching_Repo/master/stock_prices.csv", header = TRUE, sep = " ", dec = ",")
BOA        <- read.csv(file = "stock_prices.csv", header = TRUE, sep = " ", dec = ",")
Dates      <- BOA$date
BOA        <- BOA[,-1]
BOA        <- data.matrix(BOA)
BOA        <- BOA[-which(Dates=="08/26/2004"),]# Outlier
n          <- dim(BOA)[1]
J          <- dim(BOA)[2]
Times      <- seq(0,6.5,length=J)

# Cumulative log-return functions (raw data)
log_BOA    <- log(BOA) - matrix(log(BOA)[,1], nrow=n, ncol=J)

# B-spline basis functions
bspl_basis <- create.bspline.basis(rangeval=c(0,6.5),nord=4,nb=200)

# Cumulative log-return functions (with basis functions)
log_BOA_fd <- Data2fd(Times,t(log_BOA),basisobj = bspl_basis)
```

```{r, eval=FALSE, echo=TRUE}
# Plot functional data
plot(log_BOA_fd, xlab="Trading Hours",ylab="",lwd=1, col=gray(.5), 
     main="Cumulative Log-Return Functions")
lines(log_BOA_fd[1:10],lwd=1.5, lty=1)
```

```{r, eval=TRUE, echo=FALSE, fig.align='center', fig.height=4}
par(mfrow=c(1,1), mar=c(5.1, 4.1, 2.1, 2.1))
invisible(plot(log_BOA_fd, xlab="Trading Hours",ylab="",lwd=1, col=gray(.5), 
               main="Cumulative Log-Return Functions"))
lines(log_BOA_fd[1:10],lwd=1.5, lty=1)
```


## Case Study: BOA Stock Returns
```{r echo=FALSE, fig.align='center', fig.cap="Plot of mean function for BOA cumulative returns. Point-wise 95% confidence intervals are included in red."}
muhat          <- mean.fd(log_BOA_fd)
sdhat          <- sd.fd(log_BOA_fd)
SE_hat_U       <- fd(basisobj=bspl_basis) # create upper CI bound
SE_hat_L       <- fd(basisobj=bspl_basis) # create lower CI bound
SE_hat_U$coefs <-  2*sdhat$coefs/sqrt(n) + muhat$coefs
SE_hat_L$coefs <- -2*sdhat$coefs/sqrt(n) + muhat$coefs 
# plot
invisible(plot.fd(SE_hat_U,ylim=c(-0.002,0.002),col='red',lty=2,
                  xlab="Trading Hours",ylab="Returns", main="Mean Function with 95% CI"))
invisible(plot.fd(SE_hat_L,add=TRUE,col='red',lty=2))
invisible(plot.fd(muhat,add=TRUE))
```


## Case Study: BOA Stock Returns
```{r echo=FALSE, fig.cap="EFPC's of BOA cumulative returns versus the theoretical eigenfunctions of Brownian Motions.", fig.align='center'}
log_BOA_fpca  <- pca.fd(log_BOA_fd, nharm=4)
# Eigenfunctions of BM process
efwp_fun      <- function(x, k=1){sqrt(2)*sin((k-1/2)*pi*x)}
# integrate(function(x){(efwp_fun(x,k=1))^2}, lower = 0, upper = 1)
EFWP_mat      <- cbind(efwp_fun(x=xx,k=1), efwp_fun(x=xx,k=2), efwp_fun(x=xx,k=3), efwp_fun(x=xx,k=4))
# plot
par(mfrow=c(1,2))
invisible(plot(log_BOA_fpca$harmonics, main="Estimated Eigenfunctions of\nCumulative Return Functions",
               xlab="Trading Hours", ylab=""))
matplot(x = xx, y=EFWP_mat, type="l",  main="True Eigenfunctions of\nBrownian Motion",ylab="",xlab="")
```


##
<br><br><br>
<div class="centered">
<div class="red2">
<bdi style="font-size:290%;">
<b>Mathematical Framework:<br>1. Square Integrable Functions</b>
</bdi>
</div>
</div>

## Square Integrable Functions

<!-- - Let's consider standardized intervals $[a,b]=[0,1]$.<br> -->
<!--   Standardization: $t^{\text{new}}=(t^{\text{old}}-a)/(b-a)$. -->

- A function $f$ is said to be <span style="color:#FF4136">**square integrable**</span>, i.e., $f\in L^2([a,b])$ if
$$
\textstyle\int_a^b \big(f(t)\big)^2 dt<\infty.
$$
- Square integrable functions form a <span style="color:#FF4136">**vector space**</span>, i.e., for $c,d\in\mathbb{R}$: 
$$g,f\in L^2([a,b])\quad\Rightarrow\quad cf+dg\in L^2([a,b]),$$ 
where addition is pointwise, i.e. <br>
$(cf+dg)(t)=cf(t)+dg(t)$ for *almost all* $t\in[0,1]$. 
<br><br>
(We will ignore measure-theoretic considerations.)


##  Square Integrable Functions 

What makes the space $L^2([a,b])$ so convenient is its structure:

- <span style="color:#FF4136">**Inner product**</span> of two functions is given by
$$
\textstyle\langle f,g\rangle=\int_a^bf(t)g(t)dt
$$
-  Two functions $f$ and $g$ are <span style="color:#FF4136">**orthogonal**</span> if 
$$\langle f,g\rangle=0$$

- The <span style="color:#FF4136">**norm**</span> $||f||=\sqrt{\langle f,f\rangle}$, gives us a notion for the distance between two functions: $$d(f,g)=||f-g||$$


##  Square Integrable Functions 

As we have already seen, <span style="color:#FF4136">**basis expansions**</span> play an important role in
methodology for functional data. 

- We say that a set of functions $\{e_1 ,e_2,e_3,\dots\}$
is a <span style="color:#FF4136">**basis**</span> in $L^2([a,b])$ if every function $f\in L^2([a,b])$ admits a *unique* expansion
$$
\textstyle f(t)=\sum_{j=1}^\infty a_j e_j(t)
$$
- We say that $\{e_1 ,e_2,e_3,\dots\}$ is an <span style="color:#FF4136">**orthonormal basis**</span> if, in addition
$$\langle e_j , e_k\rangle=\begin{cases}1, &j=k\\0&j\neq k.\end{cases}
$$ 


##
<br><br><br>
<div class="centered">
<div class="red2">
<bdi style="font-size:290%;">
<b>Mathematical Framework:<br>2. Random Functions</b>
</bdi>
</div>
</div>


## Random Functions

Let $X$ denote a <span style="color:#FF4136">**random function**</span> defined on a probability space, say
$\Omega$.


- We assume that *all* <span style="color:#FF4136">**realizations**</span> $X(\omega)$, $\omega\in\Omega$, are elements of the space $L^2([a,b])$ of square integrable functions, i.e <br><br>
$||X(\omega)||=\sqrt{\int_a^b \big(X(\omega)(t)\big)^2dt} < \infty$ for all $\omega\in\Omega$.

- $||X||\in\mathbb{R}$ is thus a <span style="color:#FF4136">**real random variable**</span>. 

- If $E(||X||^2)<\infty$, we say that the <span style="color:#FF4136">**random function $X$ is square integrable**</span>. 
**Caution:** Integration is here with respect to $\omega$, not $t$. It might be more pedagogical to say that the **random function $X$ has a finite second moment**. 


## Random Functions

- <span style="color:#FF4136">**Mean function:**</span>
$$
\mu(t)=E(X(t))
$$
- <span style="color:#FF4136">**Covariance function:**</span>
$$
c(t,s)=E\Big((X(t)-\mu(t))(X(s)-\mu(s))\Big)
$$
- The sample mean function $\bar{X}_n(t)$ and the sample covariance function $\hat{c}_n(t,s)$ are viewed as <span style="color:#FF4136">**estimators**</span> of the population parameters $\mu(t)$ and $c(t,s)$.


## Random Functions 

The population covariance function leads to <span style="color:#FF4136">**Functional Principal Component Analysis (FPCA)**</span>, which allows us to represent a square integrable random function $X$ as:
$$
\textstyle X(t)=\mu(t)+\sum_{j=1}^\infty\xi_j v_j(t)
$$

- The <span style="color:#FF4136">**eigenfunctions**</span> $v_j$ are the solutions of the eigen-equation
$\int_a^bc(t,s)v_j(s)ds=\lambda_jv_j(t)$. 

- $\lambda_1\geq\lambda_2\geq\dots$ denote the <span style="color:#FF4136">**eigenvalues**</span>. 

- The *random variables* $\xi_j$ are called the <span style="color:#FF4136">**scores**</span>
$$
\textstyle \xi_j=\langle X-\mu,v_j\rangle=\int_a^b(X(t)-\mu(t))v_j(t)dt 
$$
<!-- - Typically, the eigenvalues, and so the corresponding eigenfunctions $v_j$ are arranged in nonincreasing order: $\lambda_1 \geq \lambda_2 \geq\dots$.  -->

## Random Functions 

<span style="color:#FF4136">**FPCA**</span> allows us to represent a square integrable random function $X$ as (<span style="color:#FF4136">**Karhunen-Loéve expansion**</span>):
$$
\textstyle X(t)=\mu(t)+\sum_{j=1}^\infty\xi_j v_j(t)
$$

It can be shown that:
$$
E(\xi_j)=0,\quad V(\xi_j)=\lambda_j,\quad \operatorname{Cov}(\xi_j,\xi_k)=0,\,j\neq k,
$$
and $E\int_a^b(X(t)-\mu(t))^2dt=E||X-\mu||^2=\sum_{j=1}^\infty\lambda_j$.

That is, $\lambda_j$ is the variance of the random function $X$ <span style="color:#FF4136">**in the principal direction**</span> $v_j$ and the sum of these variances is the <span style="color:#FF4136">**total variance**</span> of $X$. 


##
<br><br><br>
<div class="red2">
<div class="centered">
<bdi style="font-size:350%;"><b>Functional Regression Models</b>
</bdi>
</div>
</div>


## Functional Regression Models

**Scalar–on–function regression:**
$$Y_i=\int_a^b\beta(s)X_i(s)ds+\varepsilon_i$$
**Function–on-scalar regression:**
$$Y_i(t)=\sum_{k=1}^p\beta_k(t)x_{ik}+\varepsilon_i(t)$$
**Function–on-function regression:**
$$Y_i(t)=\int_a^b\beta(t,s)X_{i}(s)ds+\varepsilon_i(t)$$


## Functional Regression Models

**Nonparametric scalar–on–function regression:**
$$Y_i=m(X_i)+\varepsilon_i,$$
where $Y_i\in\mathbb{R}$ and $X_i\in L^2([a,b])$.

**Nonparametric function–on–function regression:**
$$Y_i=m(X_i)+\varepsilon_i,$$
where $Y_i,\,X_i\in L^2([a,b])$.

The above models are just prototypes illustrating the general idea. The main point is that the functional regression models are **infinite dimensional objects** which must be estimated from a **finite sample**.



##
<br><br><br>
<div class="red2">
<div class="centered">
<bdi style="font-size:350%;"><b>Difficulties in Functional Regression</b>
</bdi>
</div>
</div>


## Difficulties in Functional Regression

**Remember: The usual regression model**
$$
\textstyle Y_i=X\boldsymbol{\beta}+\varepsilon
$$
**Normal equations and solutions (population & sample):**
$$
\textstyle \mathbf{C}_{XY}=\mathbf{C}_{XX}\boldsymbol{\beta}\quad\Rightarrow\quad  \boldsymbol{\beta}=\mathbf{C}_{XX}^{-1}\mathbf{C}_{XY}
$$
$$
\textstyle \widehat{\mathbf{C}}_{XY}=\widehat{\mathbf{C}}_{XX}\widehat{\boldsymbol{\beta}}_n\quad\Rightarrow\quad
\widehat{\boldsymbol{\beta}}_n=\widehat{\mathbf{C}}_{XX}^{-1}\widehat{\mathbf{C}}_{XY}
$$
where <br>
$[\mathbf{C}_{XX}]_{r,\ell}=E(X_rX_\ell)$, <br>
$[\mathbf{C}_{XY}]_r=E(X_rY)$, <br>
$\big[\widehat{\mathbf{C}}_{XX}\big]_{r,\ell}=\frac{1}{n}\sum_{i=1}^n X_{ir}X_{i\ell}$, and <br> 
$\big[\widehat{\mathbf{C}}_{XY}\big]_r=\frac{1}{n}\sum_{i=1}^n X_{ir}Y_{i}$.


## Difficulties in Functional Regression

**Normal equations and solutions (population & sample):**
$$
\begin{align}
\textstyle \mathbf{C}_{XY}=\mathbf{C}_{XX}\boldsymbol{\beta}\quad&\Rightarrow\quad  \boldsymbol{\beta}=\mathbf{C}_{XX}^{-1}\mathbf{C}_{XY}\\
\textstyle \widehat{\mathbf{C}}_{XY}=\widehat{\mathbf{C}}_{XX}\widehat{\boldsymbol{\beta}}_n\quad&\Rightarrow\quad
\widehat{\boldsymbol{\beta}}_n=\widehat{\mathbf{C}}_{XX}^{-1}\widehat{\mathbf{C}}_{XY}
\end{align}
$$

To obtain $\boldsymbol{\beta}$ and $\widehat{\boldsymbol{\beta}}_n$, we need <span style="color:#FF4136">**invertibility**</span> of the $p\times p$ matrices $\mathbf{C}_{XX}$ and $\widehat{\mathbf{C}}_{XX}$. 


## Difficulties in Functional Regression 

**Functional regression**
$$
\textstyle Y_i=\int_a^b\beta(t)X_i(t)+\varepsilon_i
$$
**Normal equations (population & sample):**
$$
\textstyle c_{XY}(t)=\int_a^b c_{XX}(t,s)\beta(s)ds
$$
$$
\textstyle \hat{c}_{XY}(t)=\int_a^b \hat{c}_{XX}(t,s)\beta(s)ds
$$
where <br>
$c_{XX}(t,s)=E(X(t)X(s))$, 
$c_{XY}(t)=E(X(t)Y)$, <br>
$\hat{c}_{XX}(t,s)=\frac{1}{n}\sum_{i=1}^n(X_i(t)X_i(s))$, and<br>
$\hat{c}_{XY}(t)=\frac{1}{n}\sum_{i=1}^n(X_i(t)Y_i)$.



## Difficulties in Functional Regression 

**Functional regression**
$$
\textstyle Y_i=\int_a^b\beta(t)X_i(t)+\varepsilon_i
$$
**Normal equations (sample):**
<!-- $$ -->
<!-- \textstyle c_{XY}(t)=\int_a^b c_{XX}(t,s)\beta(s)ds -->
<!-- $$ -->
$$
\textstyle \hat{c}_{XY}(t)=\int_a^b \hat{c}_{XX}(t,s)\beta(s)ds
$$

**Problem:** 

- Ill-posed inversion problem 
- Eigenvalues $\approx 0$

**Classical solutions:** 

- Schrinkage methods (ridge regression, principal component regression, etc)


##
<br><br><br>
<div class="red2">
<div class="centered">
<bdi style="font-size:350%;"><b>Estimation Through Basis Expansion</b>
</bdi>
</div>
</div>


## Estimation Through Basis Expansion

Let us consider the scalar-on-function model
$$
\textstyle Y_i=\alpha + \int_a^b \beta(s)X_i(s)ds+\varepsilon_i,\quad i=1,\dots,n
$$

Simplest approach: Expand $\beta$ using deterministic basis functions (splines, fourier, etc) $$
\textstyle \beta(t)=\sum_{k=1}^Kc_kB_k(t)
$$


Rewriting: 
<!-- the functional regression model into an approximative multivariate reggression problem -->
$$
\begin{align}
\textstyle\int_a^b\beta(s)X_i(s)ds
&\textstyle =\sum_{k=1}^Kc_k\underbrace{\int_a^bB_k(s)X_i(s)ds}_{=x_{ik}}
%&\textstyle =\sum_{k=1}^Kc_k
\end{align}
$$



## Estimation Through Basis Expansion

<!-- This reduces the the scalar-on-function model to a classical regression model: -->
$$
\begin{align}
\textstyle Y_i        &\textstyle =       \alpha + \int_a^b \beta(s)X_i(s)ds  +\varepsilon_i\\
\textstyle Y_i        &\textstyle = \alpha + \sum_{k=1}^Kc_k x_{ik}     +\varepsilon_i\\
\textstyle \mathbf{Y} &\textstyle = \mathbf{X}\mathbf{c}+\varepsilon_i,\quad\quad
\end{align}
$$

The estimator of the parameter vector $\mathbf{c}=(\alpha,c_1,\dots,c_{K})'$ 
$$
\textstyle \hat{\mathbf{c}}=(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\mathbf{Y},
$$
leads to an estimator of the parameter function $\beta(t)$
$$
\textstyle \widehat{\beta}(t)=\sum_{k=1}^K\hat{c}_kB_k(t).
$$


## Estimation Through Basis Expansion

**Biased** estimation due to the **truncation error:** $\delta_K(t)$
$$
\begin{align}
\textstyle\beta(t)
&\textstyle =\sum_{k=1}^\infty c_k B_k(t)\\
&\textstyle =\sum_{k=1}^K c_k B_k(t) + \underbrace{\sum_{\ell=K+1}^\infty c_\ell B_\ell(t)}_{=\delta_K(t)}
\end{align}
$$

It can be shown that (see pp. 55-56 of the accompanying textbook):
$$
\textstyle\hat{\mathbf{c}}=\mathbf{c}+(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\boldsymbol{\delta}_K+(\mathbf{X}'\mathbf{X})^{-1}\mathbf{X}'\boldsymbol{\varepsilon},
$$
where $[\boldsymbol{\delta}_K]_i=\int_a^b\delta_K(s)X_i(s)ds$.



##
<br><br><br>
<div class="red2">
<div class="centered">
<bdi style="font-size:350%;"><b>Estimation with a Roughness Penalty</b>
</bdi>
</div>
</div>

## Estimation with a Roughness Penalty

$$
\begin{align}
\textstyle P_{\lambda}(\alpha,\beta)=
&\textstyle \sum_{i=1}^n\Big(Y_i-\alpha-\int_a^b\beta(s)X_i(s)ds\Big)^2\\
&\textstyle +\lambda\int_a^b\big[(L\beta)(s)\big]^2ds
\end{align}
$$

**Idea:** 

- Use (again) the expansion $\beta(s)=\sum_{k=1}^K c_kB_k(t)$, but with a very high $K\gg 0$.
- Regularization is controlled by $\lambda$ and $L$. <br>
  Usually: $(L\beta)(s)=\beta''(s)$.

It can be shown that (see pp. 56-57 of the accompanying textbook):
$$
\textstyle \mathbf{c}_{\lambda}=(\mathbf{X}'\mathbf{X}+\lambda \mathbf{R})^{-1}\mathbf{X}'\mathbf{Y}.
$$


## Estimation with a roughness penalty

Continued:

$$
\textstyle \mathbf{c}_{\lambda,K}=(\mathbf{X}'\mathbf{X}+\lambda \mathbf{R})^{-1}\mathbf{X}'\mathbf{Y}
$$
where
$$
\mathbf{R}=
\left(\begin{matrix}
0      &0     &\dots &0     \\
0      &R_{11}&\dots &R_{1K}\\
0      &R_{21}&\dots &R_{2K}\\
\vdots &\vdots&\vdots&\vdots\\
0      &R_{K1}&\dots &R_{KK}\\
\end{matrix}\right)
$$
with $R_{rk}=\int_a^b(LB_r)(s)(LB_k)(s)ds$.


##
<br><br><br><br>
<div class="red2">
<div class="centered">
<bdi style="font-size:350%;"><b>Regression on FPCs</b>
</bdi>
</div>
</div>


## Regression on FPCs

$$
\begin{align}
\textstyle Y_i        &\textstyle =       \alpha + \int_a^b \beta(s)X_i(s)ds +\varepsilon_i\\
\textstyle Y_i        &\textstyle =       \alpha + \int_a^b \beta(s)\big(\hat{\mu}(s)+\sum_{j=1}^p\hat{\xi}_{ij}\hat{v}_j(s)\big)ds      +\varepsilon_i\\
\textstyle Y_i        &\textstyle =\beta_0+\sum_{j=1}^p\beta_j\hat{\xi}_{ij} +\varepsilon_i,
\end{align}
$$
where 
$$\textstyle\beta_0=\alpha+\int_a^b\beta(s)\hat{\mu}(s)ds$$ 
and 
$$\textstyle\beta_j=\int_a^b\beta(s)\hat{v}_j(s)ds$$
are treated as unknown parameters. 



## Regression on FPCs

The parameter vector $\mathbf{\beta}=(\beta_0,\beta_1,\dots,\beta_p)'$ can be estimated by (see page 59 of the accompanying textbook)
$$
\textstyle\mathbf{\beta}=(\mathbb{X}'\mathbb{X})^{-1}\mathbb{X}'Y
$$

where
$$
\mathbb{X}=
\left(\begin{matrix}
1     &\hat{\xi}_{11}&\dots &\hat{\xi}_{1K}\\
\vdots&\vdots        &\vdots&\vdots        \\
1     &\hat{\xi}_{K1}&\dots &\hat{\xi}_{KK}\\
\end{matrix}\right)
$$



## Regression on FPCs

The parameter vector $\mathbf{\beta}=(\beta_0,\beta_1,\dots,\beta_p)'$ can be estimated by (see page 59 of the accompanying textbook)
$$
\textstyle\mathbf{\beta}=(\mathbb{X}'\mathbb{X})^{-1}\mathbb{X}'Y
$$

The estimators $\hat\beta_0,\hat\beta_1,\dots,\hat\beta_p$ lead to the estimators 
$$\textstyle \hat{\beta}(t)=\sum_{j=1}^p\beta_j\hat{v}_j(t)$$
and
$$\textstyle \hat{\alpha}  =\hat\beta_0-\sum_{j=1}^p\beta_j\int_a^b\hat{v}_j(s)\hat\mu(s)ds$$
of our actual interest.



##
<br><br><br>
<div class="red2">
<div class="centered">
<bdi style="font-size:350%;"><b>Functional Regression:<br>Hands On</b>
</bdi>
</div>
</div>


## Functional Regression: Hands On

```{r, echo=TRUE}
set.seed(9000)
n    <- 1000
grid <- seq(0, 1, length = 101)
##
beta <- sin(grid * 2 * pi)
X    <- matrix(0, nrow=n, ncol=length(grid))
for(i2 in 1:n){
 X[i2,] <- X[i2,]+rnorm(length(grid), 0, 1)
 X[i2,] <- X[i2,]+runif(1, 0, 5)
 X[i2,] <- X[i2,]+rnorm(1, 1, 0.2)*grid
 for(j2 in 1:10){
  e      <- rnorm(2, 0, 1/j2^(2))
  X[i2,] <- X[i2,]+e[1]*sin((2*pi)*grid*j2)
  X[i2,] <- X[i2,]+e[2]*cos((2*pi)*grid*j2)
}}
## Note the Integral-Approximation using '* .01':
Y = X %*% beta * .01 + rnorm(n, 0, .4)
```


## Functional Regression: Hands On
```{r, echo=TRUE}
matplot(x = grid, y=t(X), type = "l", lty=1)
```


## Functional Regression: Hands On

Estimating with the roughness penalty approach and the FPCA-appoach:
```{r, echo=TRUE}
library("refund")

## Roughness Penalty
fit.RP   <- pfr(Y ~ lf(X, bs = "ps", k = 50))

## FPCA
fit.FPCA <- pfr(Y ~ fpc(X, pve = 0.99))
```



## Functional Regression: Hands On

```{r, echo=FALSE, fig.align="center"}
suppressMessages(library(ggplot2))
suppressMessages(library(dplyr))
suppressMessages(library(reshape2))
##
coefs <-  data.frame(grid      = grid, 
                     FPCA      = coef(fit.FPCA)$value, 
                     Penalized = coef(fit.RP)$value,
                     Truth     = beta)
coefs.m           <- melt(coefs, id = "grid")
colnames(coefs.m) <- c("grid", "Method", "Value")
##
# dev.new(width=6/2,height=3/2)
ggplot(coefs.m, aes(x=grid,y=Value,color=Method,group=Method), 
       width=12,height=6) + geom_path() + theme_bw()
```



## Functional Regression: Hands On

A more complex $\beta(t)$ function:

```{r, echo=TRUE}
##
beta <- -1 * dnorm(grid, mean=.20, sd=.03) + 
         3 * dnorm(grid, mean=.50, sd=.04) + 
         1 * dnorm(grid, mean=.75, sd=.05)

## Note the Integral-Approximation using '* .01':
Y = X %*% beta * .01 + rnorm(n, 0, .4)
```


## Functional Regression: Hands On

```{r, echo=FALSE, fig.align="center"}
## Roughness Penalty
fit.RP   <- pfr(Y ~ lf(X, bs = "ps", k = 50))
## FPCA
fit.FPCA <- pfr(Y ~ fpc(X, pve = 0.99))
##
coefs <-  data.frame(grid      = grid, 
                     FPCA      = coef(fit.FPCA)$value, 
                     Penalized = coef(fit.RP)$value,
                     Truth     = beta)
coefs.m           <- melt(coefs, id = "grid")
colnames(coefs.m) <- c("grid", "Method", "Value")
##
ggplot(coefs.m, aes(x=grid,y=Value,color=Method,group=Method), 
       width=12,height=6) + geom_path() + theme_bw()
```


##
<br><br><br>
<div class="red2">
<div class="centered">
<bdi style="font-size:350%;"><b>Nonparametric Functional Regression</b>
</bdi>
</div>
</div>


## Nonparametric Functional Regression

**Nonparametric scalar–on–function regression:**
$$Y_i=m(X_i)+\varepsilon_i,$$
where $Y_i\in\mathbb{R}$ and $X_i\in L^2([a,b])$.


Functional version of the **Nadaraya-Watson** estimator:

$$
\widehat{m}(X)=\frac{\sum_{i=1}^n K(h^{-1}d(X,X_i))Y_i}{\sum_{j=1}^nK(h^{-1}d(X,X_j))},
$$
where $d$ is a **semi-metric** such as, for instance, the FPCA-based semi-metric:
$$
\textstyle d(X_i,X_j)=\sqrt{\sum_{k=1}^K(\hat{\xi}_{ik}-\hat{\xi}_{jk})^2}
$$



## Nonparametric Functional Regression

```{r, echo=TRUE, fig.align="center"}
suppressMessages(library("fda.usc")); data(tecator)
absorp <- tecator$absorp.fdata
x      <- absorp[1:129,]
y      <- tecator$y$Fat[1:129]
```
```{r, echo=FALSE}
plot(x)
```


## Nonparametric Functional Regression
```{r, echo=TRUE, fig.align="center"}
res.pca1 <- fregre.np(x,y,Ker=AKer.tri,metri=semimetric.pca,q=1)
summary.fregre.fd(res.pca1, draw = FALSE)
```