Some useful manuals and tutorials for R (in Czech or English):
- Bína, V., Komárek, A. a Komárková, L.: Jak na jazyk R. (in Czech) | (PDF manual)
- Komárek, A.: Základy práce s R. (in Czech) | (PDF manual)
- Kulich, M.: Velmi stručný úvod do R. (in Czech) | (PDF manual)
- Venables, W.N., Smith, D.M., and the R Core Team: An Introduction to R. | (PDF manual)
- De Vries, A. a Meys, J.: R for Dummies. (ISBN-13: 978-1119055808) | (PDF book)
The covariance/correlation structure, which is of the main interest here (and will be later assessed using the factor analysis method), can be either estimated using a standard variance covariance matrix (command `var(bioData[,2:18])`) or the correlations can be considered instead -- they can be visualized using the `corrplot()` function from the 'corrplot' package instead. ```{r, message=F, fig.width=10, fig.height=8} library(corrplot) corrplot(cor(bioData[,2:18]), method="ellipse") ``` Recall, that unlike the covariance which is not scale invariant the already correlation is. Compare the outcomes for the following commands: ```{r, eval = F} cov(bioData[,2:18]) cov(scale(bioData[,2:18])) cor(bioData[,2:18]) cor(scale(bioData[,2:18])) ``` Can you justfify why the last three commands give the same results?
The idea behind the `factanal()` function is to consider a $p$-dimensional (random) vector $\boldsymbol{X}$ and to express this vector using a set of $k$ factors where we require want that $k \ll p$ (dimensionaly reduction). The model which is fitted by the factor analysis approach when appling the `factanal()' function in R is the following one: $$ \boldsymbol{X} = \Lambda \boldsymbol{F} + \boldsymbol{e}, $$ where $\Lambda$ is a $p\times k$ dimensional (nonrandom) matrix of so called factor loadings, $\boldsymbol{F}$ is a $k$-dimensional (random) vector represending $k$ factors (common factors or scores respectively) and $\boldsymbol{e}$ is an approximation error (or specific factors). The form of the expression above closely reminds the model of the linear regression, where $$ \boldsymbol{Y} = \mathbb{X} \boldsymbol{\beta} + \boldsymbol{\varepsilon}, $$ where $\boldsymbol{Y}$ is the vector of dependent covariates, $\mathbb{X}$ is the given model matrix of the type $n \times p$ and $\boldsymbol{\varepsilon}$ is a random vector -- the errror terms. However, the tricky part in the factor decomposition $\boldsymbol{X} = \Lambda \boldsymbol{F} + \boldsymbol{e}$ is that beside the random vector $\boldsymbol{X}$ no other quantity is directly observed. The problem, as such, is unsolvable in a direct and unique way unless we impose some additional restrictions. This is done by specifying a varinace covariance structure among all quantities which appear in the expression above:
- The (random) common factors $\boldsymbol{F}$ are required to be uncorrelated, with a unit variance, such that $Var (\boldsymbol{F}) = \mathbb{I}_{k};$
- The specific factors (error terms) $\boldsymbol{e}$ are independent with some variance, thus $Var (\boldsymbol{e}) = \boldsymbol{\Psi}$, where $Var (e_{i}) = \psi_{ii}$;
- The correlation matrix of $\boldsymbol{X}$ is decomposed as $\Sigma = \Lambda \Lambda^{T} + \boldsymbol{\Psi}$.
- The mutual covariance between the common factors $\boldsymbol{F}$ and the specific factors $\boldsymbol{e}$ is zero, thus $Cov(\boldsymbol{F}, \boldsymbol{e}) = \boldsymbol{0}$;
A common technique used in the factor analysis approach is to assume that the matrix $\Lambda^\top \boldsymbol{\Psi}^{-1} \Lambda$ is dianogal. This introduces additional equations needed to solve the problem.
- from the factor decomposition itself, thus the equation $\Sigma = \Lambda \Lambda^\top + \boldsymbol{\Psi}$, we have $p (p + 1)/2$ equations in total;
- the total number of unknown parameters to find is $p \times k$ parameters from the matrix $\Lambda$ and, in addition $p$ parameters from the diagonal matrix $\boldsymbol{\Psi}$;
- the additional requirement on the diagonality of $\Lambda^\top \boldsymbol{\Psi}^{-1} \Lambda$ introduces additional $k(k + 1)/2$ equations;
- therefore, the overall degrees of freedom of the factor model can be obtained as $$ \frac{(p - k)^2}{2} - \frac{k + p}{2}; $$
#### Note
- Factor analysis is not unique -- any rotation of $X$ gives the same results;
- Factor analysis is invariant with respect to the scale;
Since the factors are not unique with respect to the rotation, it is usefull to provide some rotation which makes sense. The automatic proceduce which tries to find factors in a way that original covariates can be splitted into disjoint sets is called the varimax procedure.
To use the varimax procedure for detecting the right rotation we can use the additional parameter 'rotation="varimax"' when calling the function `factanal()` in R. ```{r} fa1 <- factanal(bioData[,2:18], factors = 3, rotation="varimax") print(fa1, digits=2, cutoff=.6, sort=TRUE) ``` The blank spaces for some covariates and some factors are not trully missing values but they rather represent some correlation which is below a pre-defined threshold limit (too small correlation between the given covariate and the given latent factor). If there is no reported correlation for some specific covariate, maybe the number of latent factors is not enough (no factor is able to sufficiently correlate with the given covariate). Alternatively, one can also try to lower the threshold value instead.
Is such factor representation enough? Compare the model with some others for the increasing number of factors used: ```{r} fa2 <- factanal(bioData[,2:18], factors = 4, rotation="varimax") fa3 <- factanal(bioData[,2:18], factors = 5, rotation="varimax") fa4 <- factanal(bioData[,2:18], factors = 9, rotation="varimax") ``` There are mainly two different approaches on how to define the number of factors which should be used. For instance, it is easy to see that for each element $X_{j}$ from the random vector $\boldsymbol{X} = (X_{1}, \dots, X_{p})^\top$ it holds that $$ var(X_{j}) =h_{j}^2 + \psi_{j}, \qquad \textrm{where} \qquad h_{j}^2 = \sum_{\ell = 1}^{k} q_{j \ell}^2, $$ and $\boldsymbol{\psi} = Diag\{\psi_{1}, \dots, \psi_{p}\}$, and $\Lambda = \{q_{i j}\}_{i, j = 1}^{p,k}$. The quantity $h_{j}^2$ is the overall variability of $X_{j}$ which is explained by the common factors in $F$ (also called the communality) while the second quantity, $\psi_{j}$, is the variability which is left (also called the specific variability).
It is obvious, that for $p$ factors we obtain that $var(X_{j}) = \sum_{\ell = 1}^{p} q_{j \ell}^2$ and the specific fariability is zero, thus $\psi_j = 0$, for any $j = 1, \dots, p$.
Therefore, to decide on how many factors are needed, one can base the decision on a comparison of the communality and the specific variability.
Alternative approaches to judge the right number of factors:
- Statistical pre-analysis
One usually runs, for instance, a principal components analysis to determine how many factors should be enough (in order to have some reasonable proportion of the overal variability explained); - Expert judgement
Sometimes (under some optimal interpretation options) the estimated factors nicely correspond with some latent variables which are not observed, however, can be indentified by some expert judgement. - Maximum Likelihood Test
Using the likelihood approach we can also test the null hypothesis $H_{0}$ that $\Sigma = \Lambda\Lambda^\top + \boldsymbol{\psi}$ against a general altternative $H_{1}$ specifing no restrictions on the variance covariance matrix $\Sigma$. The likelihood ratio test is given by the test statistic $$ T = - 2 log \Bigg[\frac{\textrm{maximum likelihood under $H_{0}$}}{\textrm{maximum likelihood under $H_{1}$}}\Bigg]. $$ However, this approach assumes that the data follow the normal distribution with some specific parameters.
#### Note
- Another option on how to determine the optimal number of factor to be extracted in the analysis is to use the tools from the library 'nFactors' (the library can be installed by `install.packages("nFactors")`);
```{r, warning = F, message = F, fig.width=10, fig.height=6} library(nFactors) ev <- eigen(cor(bioData[,2:18])) ap <- parallel(subject=nrow(bioData[,2:18]),var=ncol(bioData[,2:18]), rep=100, cent=.05) nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) plotnScree(nS) ``` For different qualitative criteria are used to propoes the right number of factors to be used. ```{r, fig.width=10, fig.height=6} ### colors COLS <- c("red", "green", "blue") M <- apply(abs(fa1$loadings[,1:3]), 1, max) cols <- COLS[matrix(rep(c(1,2,3), 17), nrow = 17, byrow = T)[match(M, abs(fa1$loadings[,1:3]))]] load <- fa1$loadings[,1:2] plot(load,type="n") text(load,labels=names(bioData[,2:18]),cex=.7, col = cols) ``` Can you explain why there is one covariate (Saprind) clearly out of the ren buch of other related covariates? Remember, that there are three factors used in total and the plot obove only uses two of these factors. How would you interpret the fact, that the covariate MDS_2 is in between the red ones?
```{r, fig.width=10, fig.height=6} ### colors par(mfrow = c(1,2)) load1 <- fa1$loadings[,c(1,3)] plot(load1,type="n") text(load1,labels=names(bioData[,2:18]),cex=.7, col = cols) load2 <- fa1$loadings[,2:3] plot(load2,type="n") text(load2,labels=names(bioData[,2:18]),cex=.7, col = cols) ``` Does it make more sense now? Remember, that the original data are with the dimension $p = 17$!
Alternatively, we can also try the library `FactoMineR` and the corresponding function `PCA()` from this library to obtain the complete factor map (covariate map with the mutual correlation). ```{r, fig.width = 10} library(FactoMineR) par(mfrow = c(1,2)) result <- PCA(bioData[,2:18]) ```
#### Note
- The `factanal()` function in R uses scaled variables as the starting point for the factor analysis -- the variance-covariance matrix of $\boldsymbol{X}$ is replaced by the corresponding correlation matrix with ones on its diagonal. Remember, that the factor analysis is, unlike the principal component analysis, scale invariant as we are focussing on the covariance structure (the off diagonal elements of the variance-covariance matrix) rather than the variance structure (diagonal elements of the variance-covariance matrix).
Now, we have an idea that maybe three factors should be ok. In orther to get the corresponding score values (those are the values in $\boldsymbol{F}$ in the expression above given for each observations - in our case $n = 65$) from the `factanal()` function we need an additional parameter to be specified: ```{r} fa1 <- factanal(bioData[,2:18], factors = 3, scores = "regression") fa1$scores ``` and the corresponding loadings (matrix $\Lambda$ in the expression above) is obtained as ```{r} fa1$loadings ```
#### Try by yourself
- Use the factor analysis approach and consider various number of factors used to explain the covariance structure.
- Obtain the corresponding regression scores and reconstruct the original covariates using the factor values and the given factor loadings. Compare the quality of this approximation.
- Use some reasonable error measure and compare the approximation error with respect to the various number of factors.
### 2. Application in Regression or SEM Structural Equation Models (SEM) are statistical modelling techniqes showing great potential especially in causal dependencies between endogenous and exogenous variables, and the measurement model showing the relations between latent variables and their indicators (which is also the case of the latent factors $F$ in the factor analysis).
It this situations is common to firstly draw a schematic diagram with the assumed underlying structure. Later, the diagram is used to specify the final model.
An example of such diagram is below.
In the R software there is package 'sem' available for the structural equation modelling. ```{r, warning=F,message=F} library("sem") ``` and a simple example using the dataset 'PoliticalDemocracy' from the package `library("lavaan")`: ```{r, warning=F,message=F} library("lavaan") data(PoliticalDemocracy) PoliticalDemocracy ``` with the corresponding model structure, which we want to fit:
#### Try by yourself
- Consider the dataset `HolzingerSwineford1939` from the R package `library("lavaan")` which can be loadid in by the command ```{r} data(HolzingerSwineford1939) ```
- Apply the factor analysis approach for the covariates $X_{1}, \dots, X_9$ and consider three different factors when explaining the underlying covariance structure in amoung these covariats: ```{r} fa_Hol <- factanal(HolzingerSwineford1939[,7:15], factors = 3, rotation="varimax") print(fa_Hol, digits=2, cutoff=.45, sort=TRUE) ```
- What is the interpretation of the results from above? Could you propose a linear regresion model (or SEM alternatively) to get some quantitative results?
- The confirmatory factor analysis can do the work for you... The R function `cfa()` from the R packages "lavaan" is very similar to the `sem()` function used before. In fact, these two functions are almost identical. ```{r} HS.model <- ' visual =~ x1 + x2 + x3 textual =~ x4 + x5 + x6 speed =~ x7 + x8 + x9 ' fit <- cfa(HS.model, data=HolzingerSwineford1939) summary(fit, fit.measures=TRUE) ```
- What would be the corresponding diagram capturing the dependence structure from the model above?
## Homework Assignment ### (Deadline: Lab Session 8 / 17.04.2020) Use the command below to instal the `SMSdata` packages with various multivariate datasets: ```{r, eval = F} install.packages(pkgs="http://www.karlin.mff.cuni.cz/~hlavka/sms2/SMSdata_1.0.tar.gz", repos=NULL, type="source") ``` Chose one dataset of your choise from the list of all available datasets in the package: ```{r, eval = F} data(package = "SMSdata") ``` There are 21 different datasets and you can load each of them by typing and running `data(dataset_name)`, for instance: ```{r, eval = F} data(plasma) ```
- Consider the same dataset you applied the principal component analysis to and apply the factor analysis now;
- Compare the outcomes of both, PCA and FA;
- Chose one covariate as the dependent covariate and try to use the regression modelling approach to model the expectation of this covariate using the factors obtained from the factor analysis. Try to find some reasonable interpretation of the final model.