Mathematical Background
For a fixed \(t\in [0,T_{\max}]\), let \(E^{t}\) the projection operator, \[E^{t}(z) = \mathbb{E}(z|X_{t}) = \sum_{x\in S}{\mathbb{E}(z|X_t=x) \mathbf{1}_{X_t=x}}\] Then, the correlation coefficient between \(z\) and \(X_{t}\) is: \[\eta^{2}(z, X_{t}) = \frac{var(\mathbb{E}(z|X_{t}))}{var(z)}\]
For \(t, s \in [0,T_{\max}]\), the (simple) correspondence analysis looks for \(z\) such that it maximizes \[\displaystyle\frac{1}{2}\left(\eta^{2}(z, X_{t})+\eta^{2}(z, X_{s})\right)\] and the solution is:
\[(E^{t}+E^{s})z = \lambda z\] \(z\) is called a principal component.
As an extension, the functional correspondence analysis is defined as the solution of the optimization problem:
\[ \arg\max_{z}\displaystyle\frac{1}{T_{max}}\int_{0}^{T_{max}}\eta^{2}(z, X_{t})dt\] Solution: \[\int_{0}^{T_{max}}E^{s}z\ ds = \lambda z\]
Denote \(\xi_t = E^{t}z\),
\[\int_0^{T_{\max}}K(t,s)\xi_s\ ds = \lambda \xi_t, \ \ \ \ \mbox{with}\ K(s,t) = E^tE^s. \]
It follows that, \[\xi_{t} = \sum_{x\in S}a_{x}(t)\mathbf{1}_{X_t=x}\] with \(a_x(t) =\mathbb{E}(z|X_t= x), \ \ \forall x \in S\).
For each \(x\in S\), the functions \(a_x: [0,T_{\max}|\rightarrow \mathbb{R}\) are called optimal encoding of the state \(x\).
The eigenvalue problem becomes in terms of encodings:
\[\lambda a_x(t) = \sum_{y\in S}\displaystyle\int_{0}^{T_{\max}}\frac{p_{x,y}(t,s)}{p_x(t)}a_y(s)ds, \ \ \ \ \ \forall x \in S, \forall t\in [0,T_{\max}], \] where\(p_{x,y}(t,s) = \mathbb{P}(X_t=x, X_s=y)\) and \(p_x(t)= \mathbb{P}(X_t=x)\).
Under general conditions (continuity in probability of \(X\)), there exists \(\{\lambda_i\}_{i\geq1}\) positive eigen-values and \(\{a_x^i\}_{i\geq 1}\) eigen-functions.
The following expansion formula holds (Mercer thm.):
\[p_{x,y}(t,s) = p_x(t)p_y(s)\sum_{i\geq 1}\lambda_{i}a^i_x(t)a^i_y(s)\] For \(x=y\), one obtains:
\[p_x(t) =\displaystyle\frac{1}{\displaystyle\sum_{i\geq 1}\lambda_{i}a^i_x(t)a^i_y(s)}\]
We are interested in approximating the encoding functions, \(a_x\), by
\[a_{x} \approx \sum_{i=1}^m \alpha_{x,i}\phi_i, \] where \(\mathbf{\alpha_x} = (\alpha_{x,1}, \ldots, \alpha_{x,m}) \in \mathbb{R}^m\) are the expansion coefficients onto a basis of functions defined on \([0,T]\), \(\{\phi_1, \ldots, \phi_m\}\), \(m \geq 1\).
The main result (Deville, 1982) is that \(\alpha = (\alpha_{1}, \ldots, \alpha_{K}) \in \mathbb{R}^{Km}\) is the solution of the following eigen-value problem:
\[F^{-1}G\alpha = \lambda\alpha, \] where \(G\) and \(F\) the matrix defined by: \[G = cov(\{V_{(x,i)}\}_{x\in S; i=1:m }),\] \[F = \mathbb{E}\left(\{F_{(x,i),(y,j)}\}_{x,y\in S ; i,j=1:m}\right)\] with the random variables \(V_{(x,i)}\) and \(F_{(x,i),(x,j)}\) defined by: \[V_{(x,i)} = \int_{0}^{T_{\max}}\phi_i(t)\mathbf{1}_xdt \ \ \ \ \mbox{and}\ \ \ \
F_{(x,i), (y,j)} = \int_{0}^{T_{\max}}\phi_i(t)\phi_{j}(t)\mathbf{1}_x\mathbf{1}_y dt\]
Application
Firstly, a basis of functions must be defined. We choose a B-splines basis of order 4.
m <- 10
basis <- create.bspline.basis(c(0, Tmax), nbasis = m, norder = 4)
The optimal encoding is computed using:
fmca <- compute_optimal_encoding(d_JKT, basis, verbose = FALSE)
summary(fmca)
## #### FMCA
##
## ## Data
## Number of individuals: 100
## Number of states: 4
## Time Range: 0 to 10
## States: A C G T
##
## ## Basis
## Type: bspline
## Number of basis functions: 10
##
## ## Outputs
## Eigenvalues:
## 1.943909 1.826387 1.580183 1.444652 1.343375 1.250423
##
## Explained variance:
## 0.093 0.181 0.256 0.325 0.39 0.45
##
## Optimal encoding:
## A C G T
## [1,] 0.002526966 0.6581434 -0.3285766 -0.002817961
## [2,] -0.017722972 -0.3705861 0.3315780 -0.025040522
## [3,] -0.225499785 -0.6291435 0.6640877 0.944239116
## [4,] -0.399691921 0.5459625 -0.3712233 -0.182537725
## [5,] 0.351931808 0.1955323 -0.7963273 0.291927372
## [6,] 0.226710207 0.5288021 -0.8474282 -0.085304733
##
## Principal components:
## [,1] [,2] [,3] [,4] [,5] [,6]
## 1 0.49120988 -0.7761616 2.2237977 -0.8687401 0.20384805 0.4283304
## 2 0.80112172 -0.9389770 1.9379875 -0.1778906 0.01221174 -0.7867779
## 3 -4.05766620 0.3735541 -0.6378811 -0.4504154 0.45331863 -0.7901629
## 4 0.88436661 0.2856756 0.7993696 1.6531387 1.23758544 0.9458807
## 5 0.76233663 -2.3637042 1.6430661 -0.4241343 -0.56541813 -0.1747497
## 6 -0.06609743 1.1208784 1.5759327 0.5073759 -0.73255501 -1.7198499
##
## Total elapsed time: 33.121 s
By default this function computes bootstrap estimates of the encoding functions in order to have a confidence interval. This is controlled by the computeCI arguments. The output is a list containing the different elements computed during the process: eigenvalues, pc, alpha, F, G, V and basisobj. Note that alphacontains the coefficients of the different encoding functions and pcthe principal components. These components can be used with classical statistic methods (k-means, regression…).
The eigenvalues can be computed using:
plotEigenvalues(fmca, cumulative = TRUE)
The first encoding function coefficients \(\alpha\)’s:
## A C G T
## [1,] 0.002526966 0.6581434 -0.32857660 -0.002817961
## [2,] -0.017722972 -0.3705861 0.33157800 -0.025040522
## [3,] -0.225499785 -0.6291435 0.66408772 0.944239116
## [4,] -0.399691921 0.5459625 -0.37122331 -0.182537725
## [5,] 0.351931808 0.1955323 -0.79632725 0.291927372
## [6,] 0.226710207 0.5288021 -0.84742825 -0.085304733
## [7,] -0.152766307 0.3105074 -0.31618803 0.277708291
## [8,] 0.336894106 -0.3244741 -0.48489192 0.661187803
## [9,] -0.473978751 0.1768121 -0.02104794 0.305267426
## [10,] -0.128364681 -0.1472333 -0.19895539 0.515597429
The resulting encoding can be plotted:
## Warning: Removed 3 rows containing missing values (`geom_line()`).
or extracted as a fd object (or matrix):
encoding <- get_encoding(fmca, fdObject = TRUE)
By default, it plots (and returns) the first encoding functions. Other encoding can be accessed with the harm argument.
A confidence interval can be plotted using the addCI argument:
plot(fmca, addCI = TRUE, coeff = 2, states = "A")
Plot the two first components:
plotComponent(fmca, comp = c(1, 2), addNames = FALSE)
