---
title: "Hauptkomponentenanalyse (Principal Component Analysis)"
author: "Peter Zezula (pzezula@uni-goettingen.de)"
subtitle: "M.Psy.205"
output:
  html_document:
    theme: paper
    toc: yes
    toc_depth: 3
    toc_float: yes
  pdf_document: default
    # latex_engine: xelatex
header-includes:
 - \usepackage{comment}
params:
 soln: TRUE
 todo: FALSE
mainfont: Helvetica
linkcolor: blue
---



### PCA [Rmd](pca.Rmd)

Eine vergleichende Abgrenzung der Arten von Faktoranalysen findet sich unter [fa](fa.html)


# Erläuterungen

## PCA Zweck

Beschreiben (reproduzieren) der Kovarianz einer Menge korrelierter Variablen durch wenige unkorrelierte Variablen (Komponenten). 
Komponenten geordnet nach 'Wichtigkeit' (Anteil an erklärter Varianz).

- Reduktion vieler Maße auf wenige (einen) aussagefähige Werte (Indices).
- Bei sehr vielen Variablen in multivariatem Verfahren: Weiterrechnen mit weniger Variablen die einen substanziellem Anteil der Varianz erklären.
- Häufig konzentriert sich die Aufmerksamkeit auf die erste, wichtigste Komponente, die den höchsten Anteil an Varianz erklärt.
- Teilweise sind aber auch gerade die zweiten Komponenten interessant, weil die Hauptkomponente sowieso klar ist (Taxonomiesysteme, klin. Psych mit Schwere als principal component, aber erst die übrigen Restkomponenten sind für die Intervention interessant)
- Maßnahme, wenn Multikollinearität bei Regressionen Probleme macht.
- Inhaltliche Benennung der neu gefundenen Komponenten.

(Forschungs-)Fragen

- Anzahl der Faktoren
- Natur der Faktoren (interpretierbar?)
- Wichtigkeit von Lösungen und Faktoren (u. a. Grad der Aufklärung im Verhältnis zur Anzahl der Faktoren)
- Theorie in Faktorstruktur (lässt sich eine theoretisch begründete Faktorenstruktur finden?)
- Schätzung der Faktorscores (und ihre Aussagekraft)



# Begriffe und Äquivalenzen


## Faktorladungen

Koeffizienten/Parameter für jede Variable, die ausdrücken, welchen Einfluss die Variable auf die Bildung des Faktors hat.
In PCA wenn von der Korrelationsmatrix aus gerechnet wurde: Korrelation der Variablen mit den Faktorwerten (Werte der Personen in diesem Faktor).
Auch berechenbar über die Multiplikation der Loadings (Eigenvektor) mit der Standardabweichung des jeweiligen Faktors.
<!-- (In Statistica zu finden unter Faktorkoordinaten der Variablen.) -->


## Faktorwerte, Komponentenwerte

- Parameter der Vpn.
- Werte der Vpn in dem jeweiligen Faktor/Komponente.
- In R ergebnisobjekt$scores.
- In andern Paketen (Statistica) Faktorwerte genannt.


## Eigenwerte der Hauptkomponenten

- Eigenwert ist Parameter der Faktoren. 
- Wurzel(Eigenwert) ist Standardabweichung des Faktors.
- Eigenwert ist die Summe der quadrierten Ladungen eines Faktors (Faktorladungen) über alle Variablen.
- Eigenwert ist die Varianz dieses Faktors.
- Da die Gesamtvarianz der Variablen auf 1 gesetzt wird ist der Eigenwert zugleich der Anteil, 
  den der Faktor an der Gesamtvarianz der beobachteten Variablen erklärt.
- Die Eigenwerte sind die quadrierten Standardabweichungen der Hauptkomponenten, die R im Summary zur PCA ausgibt (ergebnisobjekt$sdev).
- In der PCA ist die Gesamtsumme der Eigenwerte = der Menge der Hauptkomponenten. Der Mittelwert der Eigenwerte ist also 1.


## Eigenvektor - Loadings

- entspricht R-Loadings ergebnisobjekt$loadings
- Die Summe der quadrierten Loadings über alle Variablen hinweg ergibt 1.
- Mit Hilfe der Eigenvektoren können die vorhergesagten Werte (Faktorwerte) errechnet werden.
- (Statistica-Ergebnisdialog: Variablen | Eigenvektor)


## Kommunalität

- Parameter der Variablen. 
- Die Summe der quadrierten Ladungen einer Variablen auf allen Faktoren ergibt die Varianz dieser Variablen, die durch die Faktoren gemeinsam erklärt wird. 
- Diese Größe wird als Kommunalität $h^2_j$ einer Variablen j bezeichnet.


# Möglichkeiten der Berechnung einer PCA in R

Befehle 

- `princomp()` im base package
- `prcomp()` im base package
- `principal()`  aus `library(psych)`

Bei PCA zur Variablenreduktion ist eine Rotation nicht sinnvoll. Dies kann durch ein Flag gesetzt werden `principal(..., rotate="none")`



# Beispiel: Wahrnehmung eigener und fremder Gefühle (Emotionale Intelligenzforschung)


Das Beispiel lehnt sich an eine Studie von Tanja Lischetzke et al (2001, Diagnostica Vol 47, No. 4, S. 167 - 177):

Aus dem Abstract: "Das Erkennen der eigenen Gefühle und der Gefühle anderer Menschen ist eine wichtige Kompetenz im Umgang mit Emotionen und Stimmungen. 
Es werden die bisher vor allem im englischen Sprachraum untersuchten Konstrukte der emotionalen Selbstaufmerksamkeit 
und der Klarheit über eigene Gefühle vorgestellt 
und die konzeptuelle Trennung der Konstrukte erstmals auf die Wahrnehmung fremder Gefühle übertragen."

Die Ausführungen hier beziehen sich allerdings nur auf den Teil der Wahrnehmung eigener Gefühle.

Ein (simulierter) Datensatz findet sich im üblichen tab-delimited Textformat unter: 
[fa-ei.txt](http://md.psych.bio.uni-goettingen.de/data/virt/v_ei.txt)

Die Fragen/Items sind:

- i1	Ich denke über meine Gefühle nach.
- i2	Ich kann meine Gefühle benennen.
- i3 	Ich schenke meinen Gefühlen Aufmerksamkeit.
- i4	Ich bin mir im unklaren darüber, was ich fühle.
- i5	Ich beschäftige mich mit meinen Gefühlen.
- i6	Ich habe Schwierigkeiten, meine Gefühle zu beschreiben.
- i7	Ich denke darüber nach, wie ich mich fühle.
- i8	Ich weiß, was ich fühle.
- i9	Ich beobachte meine Gefühle.
- i10	Ich habe Schwierigkeiten, meinen Gefühlen einen Namen zu geben.
- i11	Ich ache darauf, wie ich mich fühle.
- i12	Ich bin mir unsicher, was ich eigentlich fühle.

Die Items sind skaliert von 1 bis 4 (fast nie/manchmal/oft/fast immer)

Hier soll die Frage geklärt werden, ob sich die obigen Items auf ein paar Komponenten reduzieren lassen, wenn ja, wie viele sind sinnvoll.

## Pretests

- Bartlett Test, ob es Korrelationen (!= 0) gibt
- Kaiser-Meyer-Olkin-Kriterium zur Beurteilung der Eignung der Daten zur Durchführung einer FA.

``` {r}
# get data
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_ei.txt")
require(psych)
# Ask Bartlett test whether correlation matrix is idendity matrix (all elements are 0), should be significant
# run Bartlett test on dataframe
cortest.bartlett(dd)

# alternativly on correlation matrix (add n), here calculated on the fly using cor()
cortest.bartlett(cor(dd), n=length(dd[,1]))

# get the KMO coefficients
kmo <- KMO(dd)
# show it
kmo
# get ordered MSAi
kmo$MSAi[order(kmo$MSAi)]

```


- Die Korrelationstabelle ist sig. verschieden von der Identitätsmatrix.
- Nach dem Kaiser-Meyer-Olkin Kriterium .81 eignen sich die Items recht gut zur Faktorisierung.
- Kein MSAi ist kleiner als .64 und sie steigen bis auf .88.
  Es gibt zwei 'mäßige' MSAi-Werte unter .70 aber kein Item unter .50. 
  Damit brauchen wir kein Item auszuschließen, was bei MSAi unter .5 empfohlen wird.
  Allerdings sollten wir die Items i1 und i7 im Auge behalten.



## Anzahl der Komponenten


Via Parallelanalyse `psych::fa.parallel()`

``` {r}
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_ei.txt")
require(psych)
# get parallel analysis for principal component analysis
# plot is generated automagically
psych::fa.parallel(dd, fa="pc")
# parallel analysis suggests a one component solution
```

Die Parallelanalyse empfiehlt zwei Komponenten.



## Ansatz mit `princom()`

`princomp()` ist bereits im base System von R enthalten.

```{r}
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_ei.txt")

# fit model, run PCA
# use base package princomp()
fit <- princomp(dd, cor=FALSE)

# simple scree plot would be
# plot(fit,type="lines")
# but we already did fa.parallel

# print variance accounted for
summary(fit)

# component scores (principal components) are available and may be used for further computations
head(round(fit$scores, 2))
# biplot(fit)
```

Gut zu sehen ist (Cumulative Var), dass dieselbe Anzahl an Komponenten wie Variablen die gesamte Varianz erklärt, 
die Original Korrelationsmatrix also repliziert.

Das Problem der subjektiven Festlegung der Anzahl der Komponenten bleibt, obwohl der Screeplot hilft.

`princomp()` erlaubt keine Anpassung auf eine bestimmte Anzahl von Komponenten und keine Rotation.

`psych::principal()` und noch mehr `psych::fa()` sind flexibler und damit häufig auch anschaulicher.



## Ansatz mit `principal()` bzw. `fa()` aus der `library(psych)`

Hier mit einigen Erklärungen und Erläuterungen.


```{r}
# get data
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_ei.txt")
# take a look at the data
head(dd)
# we need library(psych)
require(psych)

# some descriptives
psych::describe(dd)

# Ask Bartlett test whether correlation matrix is idendity matrix (all elements are 0), should be significant
cortest.bartlett(dd)
# alternativly on correlation matrix (add n) 
cortest.bartlett(dd, n=nrow(dd))

# get kmo: are data good for factorization
psych::KMO(dd)

# determinant ('area' of data) should be higher than 0.00001, singularity might wait ...
det(cor(dd))

# number of factors suggested by fa.parallel() is 2
items.parallel <- fa.parallel(dd, fa="pc")


# original correlation table, that will be reproduced more or less by PCA
round(cor(dd), 2)

require(psych)
# equivalence to base package's princomp()
# full model: as much components as items, all variance explained by components.
# Cumulative Proportion adds up to 1, communalities h2 are set to 1
psych::principal(dd, rotate='none', nfactors=ncol(dd))

# what suggests fa.parallel() of library(psych)
# already done above
# items.parallel <- fa.parallel(dd, fa="pc")

# we go on with 2 components
pc.2 <- psych::principal(dd, rotate='none', nfactors=2)
# default output
# h2 are communalities
# u2 are residuals
# h2 and u2 sum up to 1
pc.2

# graph it
plot(pc.2)

# get loadings
# alternatively available via fit$loadings
loadings(pc.2)
# loadings understands parameter cutoff. only loadings above cutoff will be printed. default is .1
# a higher cutoff might result in a better understanding of the structure
print(loadings(pc.2), cutoff= .45)

# does varimax rotation help?
pc.2.r <- principal(dd, nfactors=2, rotate='varimax')
print(loadings(pc.2.r), cutoff= .3)
# even better
print.psych(pc.2.r, cut=0.3, sort=T)

# a quick look at some factor scores
# unrotated
head(round(pc.2$scores, 2))
# component scores are different for rotated components
head(round(pc.2.r$scores, 2))

# component scores may be saved to dataframe for further computations
dd <- cbind(dd, pc.2.r$scores)
head(dd)




# finally we might want to visualize the models
# we can use the usual plot() method
plot(pc.2)
plot(pc.2.r)
# or the fa.diagram() of library(psych), parameter cut is equivalent to looking at scores, parameter simple suppresses inclusion of cross-loadings
fa.diagram(pc.2, simple=TRUE, cut=.2, digits=2)
fa.diagram(pc.2.r, simple=TRUE, cut=.2, digits=2)

```


in `Output unter Standardized loadings based ...`

`h2` sind die Kommunalitäten. Anteil an Varianz der jeweiligen Variable, der durch die (hier beiden) Komponenten erklärt/repliziert wird.

`u2` sind die uniqueness-Werte (siehe [efa](efa.html)).

`com` Complexity index, wieviele Faktoren brauchen wir, um eine Variable 'hinreichend gut' zu erklären (siehe [efa](efa.html))
      Je kleiner der Wert, desto weniger Faktoren erklären schon einen großen Anteil der Varianz einer Variablen.
      The complexity index is a positive number indicating on the average how many factors are used to explain each variable in a factor solution.

Eine Erhöhung des cutoffs auf `.45` macht die Tabelle der Factor-Loadings klarer:

Die Variablen sind gut bereits unrotiert gut einer der beiden Komponenten zuordenbar.
Wir sehen sehr schön die Struktur, dass die gerade bzw. die ungeraden Items auf je einer Komponente laden.
Das Item i7 kann nicht eindeutig einer der beiden Komponenten zugeordnet werden.
Cross-loadings, also "Falschzuordnungen" sind nicht zu beobachten.
Zwei Komponenten erklären 50% der Varianz. Andersherum gesehen gehen 50% der Information verloren, 
wenn wir mit 2 Komponenten anstelle von 12 Items in einem multivariaten Verfahren weiterrechnen würde.

Nach Varimax-Rotation findet sich dieselbe Struktur.
Hier wird das kritische Item i7 zumindest knapp (mit niedrigem Ladungsgewicht) der zweiten Komponente zugeordnet, wie alle ungeraden Items. 

Mit den gespeicherten Komponenten-Werten könnten wir dann weiterrechnen.

Alles hier passiert zunächst ohne inhaltlichen Bezug. 
Hat man aber mit Komponenten gerechnete Ergebnisse, ist u. U. fraglich, was sie eigentlich aussagen.
In diesem Beispiel lassen sich die beiden inhaltlichen Itemgruppen aber klar in den beiden Komponenten identifizieren.









## Ausflug Reproduzierte Korrelationsmatrix und Residuen, Chi^2 Modelltest

todo: adapt to EI data

```{r}
dd.raw <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_bmi_preprocessed.txt", fileEncoding = "UTF-8")
require(dplyr)
#dd.f <- dplyr::select(dd, c(c_phys_app:c_dress, i_c_bad_way, c_figure))
dd <- dplyr::select(dd.raw, c(height, weight, c_phys_app:c_dress, i_c_bad_way, c_figure))
require(psych)

pc.2 <- principal(dd, rotate='none', nfactors=2)
# model test whether model with 2 components is sig. different from model with all components
pc.2

# original correlation table, that will be reproduced more or less by PCA
round(cor(dd), 2)

# get the reproduced correlation matrix, which is calculated on base of the loadings
# diagonal are the communalities
round(factor.model(pc.2$loadings), 2)
# get the residuals. These are the differences between original and reproduced correlation matrix
pc.diff <- cor(dd) - factor.model(pc.2$loadings)
round(pc.diff, 2)
rm(pc.diff)
# get the same via call to factor.residuals()
# diagonal are the uniquenesses.
pc.diff.2 <- factor.residuals(cor(dd), pc.2$loadings)
round(pc.diff.2, 2)
rm(pc.diff.2)

```

In der Diagonalen der reproduzierten Korrelationstabelle stehen die Kommunalitäten.
In der Diagonalen der Differenzen (Residuen) stehen die spezifischen Anteile (Uniquenesses).

Ein Maß für die Qualität des Fits kann auf diesen Residuen aufbauen.
Die Residuen im Verhältnis zur Höhe der Korrelation sind eine Grundidee für einen Fit-Index.
Die Summe der quadrierten Residuen im Verhältnis zur Summe der quadrierten Korrelationen, von 1 abgezogen, 
ist der Wert hinter
'Fit based upon off diagonal values = ...'
Werte über .95 sprechen für einen guten Fit.

Der Likelihood-Quotienten-Test prüft auf Signifikanz ob die Anzahl der Faktoren ausreicht, um die Daten genüugend zu erklären. 
Er setzt Normalverteilung voraus.
Geprüft wird die 0-Hypothese, dass die Daten und das Modell mit k Faktoren statistisch gleich sind. 
Dabei hat die Aufnahme eines Faktors ihren 'Preis'.
Die Anzahl der Elementgleichungen geht in die Modellprüfung mit ein.




## Ausflug Matrizen

[from](https://www.itl.nist.gov/div898/handbook/pmc/section5/pmc552.htm)

https://aaronschlegel.me/principal-component-analysis-r-example.html

```{r}
# data, 3 vars, 10 subjects
XX <- matrix(c(7,4,6,8,8,7,5,9,7,8,4,1,3,6,5,2,3,5,4,2,3,8,5,1,7,9,3,8,5,2), ncol=3)
colnames(XX) <- c("v1", "v2", "v3")
# center
cXX <- scale(XX, scale=FALSE)
# z transform
sXX <- scale(XX)
# correlation matrix
CC <- cor(XX)
VC <- cov(XX)
#solve(CC)
iCC <- solve(CC)
# Solve for the roots of R
# i.e. get the eigenvalues
eigen(CC)
# Notice that:
# - Each eigenvalue satisfies |R−λI|=0. # ??
# diag(3)  # identity matrix for three
GG <- CC - eigen(CC)$values[1] * diag(3)  # should become 0
rowSums(GG)
colSums(GG)
# - The sum of the eigenvalues =3=p, which is equal to the trace of R (i.e., the sum of the main diagonal elements).
# - The determinant of R is the product of the eigenvalues.
# - The product is λ1×λ2×λ3=0.499.
det(CC)
# the V matrix
VV <- eigen(CC)$vectors
# Notice that if you multiply V by its transpose, the result is an identity matrix, V′V=I.
VV %*% t(VV)


# check with various PCA possibilities
r.psych.nr <- psych::principal(XX, rotate='none', nfactors=ncol(XX))
r.psych    <- psych::principal(XX,                nfactors=ncol(XX))
# prcomp() and princomp()
r.prcomp    <- prcomp(XX, scale. = TRUE)
r.prcomp.nr <- prcomp(XX) 
r.princomp  <- princomp(XX, cor = TRUE)

r.psych <- psych::principal(XX, rotate='none', nfactors=ncol(XX))
head(XX)
head(psych::predict.psych(r.psych, XX))
head(r.psych$scores)
XX[1,1] * eigen(CC)$vectors[1,1] + XX[1,2] * eigen(CC)$vectors[2,1] + XX[1,3] * eigen(CC)$vectors[3,1]
XX[1,1] * eigen(VC)$vectors[1,1] + XX[1,2] * eigen(VC)$vectors[2,1] + XX[1,3] * eigen(VC)$vectors[3,1]

sXX[1,1] * eigen(CC)$vectors[1,1] + sXX[1,2] * eigen(CC)$vectors[2,1] + sXX[1,3] * eigen(CC)$vectors[3,1]

XX[1,1] * r.psych$loadings[1,1] + XX[1,2] * r.psych$loadings[2,1] + XX[1,3] * r.psych$loadings[3,1]

# XX %*% eigen(CC)$values

```


<!--
ignore

```{r eval=FALSE}
r.psych["Structure"]
r.psych$weights  # coefficien
t matrix

r.princomp$loadings
r.prcomp$rotation

# defaults
head(r.psych$scores, n=3) # = scale(XX) %*% r.psych$weights
head(predict(r.psych, data=XX), n=3)
head(predict(r.prcomp), n=3)
head(predict(r.princomp), n=3)

```




```{r}
r.pca <- prcomp(XX, scale.=TRUE) 
(variance <- (r.pca$sdev)^2)  # Eigenvalues
( varPercent <- variance/sum(variance) * 100 ) # explained variance
(loadings <- r.pca$rotation) 
(rownames(loadings) <- colnames(XX) )
(scores <- r.pca$x)

```

```{r}
data(iris)
head(iris)
# do pca
res <- prcomp(x = iris[,1:4], center = FALSE, scale. = FALSE)
res

# factor scores for first six observations
head( res$x )

# as an example: lets take a look at the factor score
# of the first observation at the first principal component
res$x[1,1]

# by hand
# multiply loadings with data and sum together to get that value:
sum( res$rotation[,1] * iris[1,1:4] )

```


-->




# Beispiele, Übungen / Exercises


Beispiel [Body-Data](example_bmi.html)

Beispiel [head data](example_head.html) aus dem Everitt-Buch.

Beispiel [crime data](example_crime.html) aus dem Everitt-Buch.

In den [Beispielen](/mv/example) finden sich zu PCA:

EI PCA

MSQ



# Links und Referenzen


R-Bloggers Tutorial John Quick
[http://www.r-bloggers.com/r-tutorial-series-exploratory-factor-analysis/]

Tutorial Wollschläger:
[http://www.uni-kiel.de/psychologie/rexrepos/posts/multFA.html]


[Field (2012, p. 112)][field] **Kapitel 3: R Environment**

<!-- this does not appear in html-page. It's used only to arm links and references. -->
[field]: http://www.uk.sagepub.com/dsur/main.htm "Field et al (2012). Discovering Statistics using R"

[everitt]: http://www.crcpress.com/product/isbn/9781439807699 "Brian S. Everitt (2010). Multivariable Modeling and Multivariate Analysis for the Behavioral Sciences"



Version: `r format(Sys.time(), '%d %B, %Y %H:%M')`