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
b <- create.bspline.basis(c(0, Tmax), nbasis = m, norder = 4)
The optimal encoding is computed using:
fmca <- compute_optimal_encoding(d_JKT, b, nCores = 2, 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.94391 1.826387 1.580183 1.444652 1.343374 1.250424
##
## Explained variance:
## 0.093 0.181 0.256 0.325 0.39 0.45
##
## Optimal encoding:
## A C G T
## [1,] 0.002527816 0.6581468 -0.3285742 -0.002813819
## [2,] -0.017726679 -0.3705776 0.3315846 -0.025042589
## [3,] -0.225499860 -0.6291487 0.6640907 0.944240906
## [4,] -0.399692155 0.5459657 -0.3712249 -0.182539384
## [5,] 0.351928506 0.1955320 -0.7963197 0.291925346
## [6,] 0.226710987 0.5287919 -0.8474339 -0.085294795
##
## Principal components:
## [,1] [,2] [,3] [,4] [,5] [,6]
## 1 0.49121450 -0.7761565 2.2237978 0.8687437 0.20384095 0.4283276
## 2 0.80112812 -0.9389699 1.9379884 0.1778935 0.01221841 -0.7867856
## 3 -4.05767245 0.3735209 -0.6378840 0.4504148 0.45332552 -0.7901623
## 4 0.88436329 0.2856844 0.7993695 -1.6531371 1.23759838 0.9458713
## 5 0.76235790 -2.3636970 1.6430662 0.4241371 -0.56541361 -0.1747624
## 6 -0.06611053 1.1208777 1.5759303 -0.5073735 -0.73253829 -1.7198606
##
## Total elapsed time: 19.112 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.002527816 0.6581468 -0.32857421 -0.002813819
## [2,] -0.017726679 -0.3705776 0.33158464 -0.025042589
## [3,] -0.225499860 -0.6291487 0.66409074 0.944240906
## [4,] -0.399692155 0.5459657 -0.37122487 -0.182539384
## [5,] 0.351928506 0.1955320 -0.79631965 0.291925346
## [6,] 0.226710987 0.5287919 -0.84743389 -0.085294795
## [7,] -0.152761503 0.3105041 -0.31618292 0.277709736
## [8,] 0.336894982 -0.3244686 -0.48491253 0.661191654
## [9,] -0.473977007 0.1768115 -0.02103771 0.305263785
## [10,] -0.128365002 -0.1472297 -0.19896052 0.515596657
The resulting encoding can be plotted:
## Warning: Removed 3 row(s) containing missing values (geom_path).
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)
