---
title: 'Testing for the number of interesting components in ICA'
author: "Klaus Nordhausen, Hannu Oja, David E. Tyler, Joni Virta"
date: "`r Sys.Date()`"
output:
html_document:
toc: yes
toc_depth: 2
rmarkdown::html_vignette:
fig_caption: yes
toc: yes
toc_depth: 2
pdf_document:
fig_caption: yes
toc: yes
toc_depth: 2
vignette: |
%\VignetteIndexEntry{ICA}
%\usepackage[utf8]{inputenc}
%\VignetteEngine{knitr::rmarkdown}
---
## Fourth order blind identification
In independent component analysis (ICA) it is usually assumed that the observable $p$-vector $x$ is a linear mixture of $p$ mutually independent components, where the vector of independent components is $s$. Hence the model can be written as
$$
x = As + \mu,
$$
where the full rank $A$ represents the linear mixing and $\mu$ the location.
It is then usually assumed that $E(s)=0$ and $Cov(s)=I_p$. The goal is then to find
an unmixing matrix $W$ such that $Wx$ has independent components.
In ICA gaussian components are considered the uninteresting ones whereas the non-gaussian ones are considered interesting. Actually, if there is more than one gaussian component only the non-gaussian components can be recovered but not the gaussian ones.
Therefore it is naturally of interest to test how many interesting components there are. This is done by testing if a subvector of $Wx$ is gaussian.
FOBI is one of the first ICA methods - it jointly diagonalizes the regular covariance matrix $Cov$ and the so-called matrix of fourth moments
$$
Cov_4(x) = E(r^2(x-\mu)(x-\mu)^T),
$$
where $\mu = E(x)$ and $r^2=(x-\mu)^T Cov(x)^{-1} (x-\mu)$.
Note that under this definition $Cov_4$ is not consistent under the normal model and the scale of a normal component would be $\sqrt{p+2}$. (Therefore for example the function `cov4` from the `ICS` package divides $Cov_4$ by $p+2$.)
Hence FOBI is the solution to the generalized eigenvector-eigenvalue problem which solves
$$
W Cov(x) W^T = I_p \ and \ W Cov_4(x) W^T = D,
$$
where D is a diagonal matrix.
For the context under consideration it is natural to order the general eigenvectors in $W$ such
that the eigenvalues, which are the diagonal elements of $D$, $d_1,...,d_p$, fulfill
$$
(d_1 - (p+2))^2 \geq ... \geq (d_p - (p+2))^2,
$$
which means gaussian components would be put at the end of $s$ and values of $d_i$'s which correspond to a gaussian component should be $p+2$.
Hence, assuming that there are $k$ interesting components the hypothesis is that
$$
H_0: d_{k+1} = ... = d_p = p+2.
$$
The package `ICtest` provides different asymptotic and bootstrapping based tests to test this hypothesis
### Asymptotic tests for a specific value of $k$
The function for the asymptotic test is `FOBIasymp` and it provides three different tests,
denoted by `S1`, `S2` and `S3`, and specified using the `type` argument:
1. `type="S1"`: The test statistic $T$ is the variance of the last $p-k$ eigenvalues around p+2:
$$
T = n \sum_{i=k+1}^p (d_i-(p+2))^2
$$
the limiting distribution of $T$ under the null is the sum of two weighted chisquare distributions with weights $w_1 = 2 \sigma_1 / (p-k)$ and $w_2 = 2 \sigma_1 / (p-k) + \sigma_2$ and degrees of freedom $df_1 = (p-k-1)(p-k+2)/2$ and $df_2 = 1$.
2. `type="S2"`: Another possible version for the test statistic is a scaled sum of the variance of the eigenvalues around the mean plus the variance around
the expected value under normality ($p+2$). Denote $VAR_{dpk}$ as the variance of the last $p-k$ eigenvalues and $VAR2_{dpk}$ as the variance of these eigenvalues around $p+2$.
Then the test statistic is:
$$
T = (n (p-k) VAR_{dpk}) / (2 \sigma_1) + (n VAR2_{dpk}) / (2 \sigma_1 / (p-k) + \sigma_2)
$$
This test statistic has a limiting chisquare distribution with $(p-k-1)(p-q+2)/2 + 1$ degrees of freedom.
3. `type="S3"`: The third possible test statistic just checks the equality of the last $p-k$ eigenvalues using only the first part of the test statistic of `type="S2"`.
The test statistic is then:
$$
T = (n (p-k) VAR_{dpk}) / (2 \sigma_1)
$$
and has a limiting chisquare distribution with $(p-k-1)(p-q+2)/2$ degrees of freedom.
The constants $\sigma_1$ and $\sigma_2$ depend on the kurtosis values of the independent components and on the dimension $p$ and are estimated from the data.
To demonstrate the function `FOBIasymp` consider the following artificial data
```{r, message = FALSE, warning = FALSE, fig.show = "hold"}
library(ICtest)
set.seed(1)
n <- 1500
S <- cbind(runif(n), rchisq(n, 2), rexp(n), rnorm(n), rnorm(n), rnorm(n))
A <- matrix(rnorm(36), ncol = 6)
X <- S %*% t(A)
```
Therefore there are three interesting components and three noise components.
```{r, message = FALSE, warning = FALSE, fig.show = "hold"}
FOBI1 <- FOBIasymp(X, k = 3, model="ICA")
screeplot(FOBI1)
abline(h = 8) # p=6, i.e. p + 2 = 8
```
The screeplot shows that the components are in the right order and 3 components have an eigenvalue close to the value of interest.
The test decision is
```{r, message = FALSE, warning = FALSE, fig.show = "hold"}
FOBI1
```
Which we can compare to the other test version:
```{r, message = FALSE, warning = FALSE, fig.show = "hold"}
FOBIasymp(X, k = 3, type = "S1", model = "ICA")
FOBIasymp(X, k = 3, type = "S2", model = "ICA")
```
which means all three tests do not reject that the last three components are gaussian.
Using the default, `model = "NGCA"` should be also consistent but does not explicitely assume an ICA model but the more general non-gaussian component analysis (NGCA) model with possible dependent interesting signals.
To see if there are 4 gaussian components using the default `model` argument one could use for example:
```{r, message = FALSE, warning = FALSE, fig.show = "hold"}
FOBIasymp(X, k = 2)
```
which would be rejected at the $0.05$ level.
### Bootstrapping tests for a specific value of $k$
The test statistic used here is again the variance of the last $p-k$ components around the value of interest $p+2$:
$$
T = n \sum_{i=k+1}^p (d_i-(p+2))^2
$$
Two possible bootstrap testing strategies, denoted `B1` and `B2`, are provided for testing the null hypothesis of interest.
Let $X$ be the centered data matrix of the $n$ $p$-variate observations and $S$ the corresponding estimated independent components which can be decomposed into $S_1$, the non-gaussian components and $S_2$ the gaussian components. Let $W$ be the estimated unmixing matrix.
1. `s.boot="B1"`:
The first strategy has the following steps:
a. Take a bootstrap sample $S_1^*$ of size $n$ from $S_1$.
b. Take a bootstrap sample $S_2^*$ consisting of a matrix of standard normally distributed elements.
c. Combine $S^*=(S_1^*, S_2^*)$ and create $X^*= S^* W^{-1}$.
d. Compute the test statistic $T^*$ based on $X^*$.
e. Repeat the previous steps `n.boot` times to get the test statistics $T^*_1,..., T^*_{n.boot}$.
Note that in this bootstrapping test the assumption of ''independent components'' is not used, it is only used that the last $p-k$ components are gaussian and independent from the first $k$ components. Therefore this strategy can be applied in an independent component analysis (ICA) framework and in a non-gaussian components analysis (NGCA) framework.
2. `s.boot="B2"`:
The second strategy has the following steps:
a. Take a bootstrap sample $S_1^*$ of size $n$ from $S_1$ where the subsampling is done separately for each independent component.
b. Take a bootstrap sample $S_2^*$ consisting of a matrix of standard normally distributed elements.
c. Combine $S^*=(S_1^*, S_2^*)$ and create $X^*= S^* W^{-1}$.
d. Compute the test statistic $T^*$ based on $X^*$.
e. Repeat the previous steps `n.boot` times to get the test statistics $T^*_1,..., T^*_{n.boot}$.
The p-value is then in both cases
$$
\frac{\#(T_i^* \geq T)+1}{n.boot+1}.
$$
This bootstrapping strategy assumes a full ICA model and cannot be used in an NGCA framework.
To test for the previous data whether there are three gaussian components, the bootstrap `FOBIboot` can be used as follows:
```{r, message = FALSE, warning = FALSE, fig.show = "hold"}
FOBIboot(X, k = 3, s.boot = "B1")
FOBIboot(X, k = 3, s.boot = "B2")
```
## Non-Gaussian projection pursuit
Non-Gaussian projection pursuit (NGPP) assumes that the observed multivariate data is generated by a linear mixture of $k$ signals and $p - k$ channels of Gaussian noise. The package provides both methods for estimating the signals when $k$ is known *a priori* and a technique for testing and estimating the value of $k$ itself. The estimation is based on the maximization of convex combinations of squared objective functions, of which four popular ones are included in the package:
- Third cumulant, `skew`
- Fourth cumulant, `pow3`
- Logarithmic hyperbolic cosine, `tanh`
- Gaussian function, `gauss`
The names (`skew` etc.) of the functions in the package are based on the corresponding derivatives, or *non-linearities*.
### Estimating the signal components
The simplest way to use NGPP is to just call the function using default settings and indicating the number of signal components `k` you want to extract. The default objective function is a certain linear combination of squared third and fourth cumulants, the *Jarque-Bera* test statistic for normality. To use some other objective function you can specify a vector of objective function names in `nl` and the corresponding weights in `alpha`.
```{r, message = FALSE, warning = FALSE, fig.show = "hold"}
library(ICtest)
X <- as.matrix(iris[, 1:4])
res1 <- NGPP(X, k = 2)
plot(res1$S, col = as.factor(iris[, 5]), xlab = "signal_1", ylab = "signal_2")
res2 <- NGPP(X, k = 2, nl=c("tanh", "pow3"), alpha = c(0.5, 0.5))
plot(res2$S, col = as.factor(iris[, 5]), xlab = "signal_1", ylab = "signal_2")
```
### Testing for a specific value of $k$
The above example used the default `method="symm"` to specify that all `k` signals should be estimated simultaneously and the alternative `method="defl"` could have been used to estimate the signals one-by-one, leading into different solutions. But perhaps more interestingly, the one-by-one estimation can also be used to test the null hypothesis $H_0: true\_k \leq k$ with the function `NGPPsim`. Let's create next some mixed toy data and try to infer its true signal dimension.
```{r, message = FALSE, warning = FALSE}
set.seed(2016)
S <- cbind(rexp(100), rnorm(100), rnorm(100))
A <- matrix(rnorm(9), 3, 3)
X <- S%*%A
res0 <- NGPPsim(X, k = 0, N = 200)
res0$p.value
res1 <- NGPPsim(X, k = 1, N = 200)
res1$p.value
```
The obtained $p$-values provide sufficient evidence to conclude that $k\_true > 0$ and $k\_true \leq 1$, correctly implying that the signal dimension is equal to one. The parameter `N` controls the number of repetitions used to simulate the distribution of the test statistic under the null hypothesis and increasing it gives more accurate results but increases the computation time.
### Estimation of the true dimension
In the above example we had to separately call `NGPPsim` for each value of `k`. This process is automated by the function `NGPPest` which performs the hypothesis test for all $k = 0, \ldots , p-1$. Let's create another toy data and try to estimate its dimension, using this time only the non-linearity `pow3`.
```{r, message = FALSE, warning = FALSE}
set.seed(2016)
S <- cbind(rexp(100), runif(100), rnorm(100))
A <- matrix(rnorm(9), 3, 3)
X <- S%*%A
res <- NGPPest(X, nl = "pow3", N = 200)
res$p.value
```
The resulting vector of $p$-values can now be directly interpreted as testing the null hypotheses whether each particular component is just noise. The correct conclusion of two signal components is again reached.