Hauptkomponentenanalyse (Principal Component Analysis)
Eine vergleichende Abgrenzung der Arten von Faktoranalysen findet sich unter fa
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
Parameter der Vpn. Werte der Vpn in dem jeweiligen Faktor. In R ergebnisobjekt$scores. In Statistica Faktorwerte genannt. Berechnung: 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
Statistica-Ergebnisdialog: Variablen | Eigenvektor 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. 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^2j einer Variablen j bezeichnet.
Mat
Core R includes a maximum likelihood factor analysis function (factanal) and the psych package includes five alternative factor extraction options within one function, fa.
Ask Bartlett test whether correlation matrix is idendity matrix (all elements are 0), should be significant
run Bartlett test on dataframe
cortest.bartlett(Data-Frame) # alternativly on correlation matrix (add n) cortest.bartlett(Korrelationsmatrix, n=...)
Faktoranalysen reproduzieren Korrelationsmatrizen. Ausgangspunkt müssen also nicht unbeding Rohdaten sein, Korrelationsmatrizen reichen.
Principle Component Analysis (PCA) - Hauptkomponentenanalyse
Befehle
princomp()im base packageprcomp()im base packageprincipal()auslibrary(psych)
Bei PCA zur Variablenreduktion ist eine Rotation nicht sinnvoll. Dies kann durch ein Flag gesetzt werden principal(..., rotate="none")
Beispiel: Head Data (Everitt, 2010)
Minimal-Versteh-Beispiel head data Everitt (2010, p. 196)
Kopfgröße (head length in mm) der ersten beiden Söhne von 25 Familien. Zwei Variablen. Wie sehen die Komponenten aus? Was kann aus ihnen gelesen werden?
Bei zwei Variablen erklären zwei Komponenten die gesamte, erklärbare Varianz (Kovarianz).
Hauptkomponentenanalyse kann auf Basis der Rohwertmatrix oder auf Basis der Varianz-Kovarianz-Matrix berechnet werden.
leichte Abweichungen bei der Standardabweichung im Vergleich zur Berechnung mit Var-Covar-Matr
individuelle Vorhersage Familie 3
ddf <- read.delim("http://r.psych.bio.uni-goettingen.de/mv/data/be/head.txt")
head(ddf)## fam hlfs hlss
## 1 1 191 179
## 2 2 195 201
## 3 3 181 185
## 4 4 183 188
## 5 5 176 171
## 6 6 208 192# we are only interested in first and second son's head size
sizes <- c('hlfs', 'hlss')
# get means of size of head of first and second son
apply(ddf[sizes], 2, mean)## hlfs hlss
## 185.7 183.8# covariance matrix
var(ddf[sizes])## hlfs hlss
## hlfs 95.29 69.66
## hlss 69.66 100.81# pca based on raw data
# you may notice small differences in sd compared to pca based on covariance matrix
pca.h <- princomp(ddf[sizes])
print(summary(pca.h),digits=3,loadings=TRUE)## Importance of components:
## Comp.1 Comp.2
## Standard deviation 12.6908 5.2154
## Proportion of Variance 0.8555 0.1445
## Cumulative Proportion 0.8555 1.0000
##
## Loadings:
## Comp.1 Comp.2
## hlfs 0.693 -0.721
## hlss 0.721 0.693# pca based on covariance matrix
pca.cov <- princomp(covmat=var(ddf[sizes]))
print(summary(pca.cov),digits=3,loadings=TRUE)## Importance of components:
## Comp.1 Comp.2
## Standard deviation 12.9525 5.3230
## Proportion of Variance 0.8555 0.1445
## Cumulative Proportion 0.8555 1.0000
##
## Loadings:
## Comp.1 Comp.2
## hlfs 0.693 -0.721
## hlss 0.721 0.693# values of family 3
ddf[3,]## fam hlfs hlss
## 3 3 181 185# factor loadings for each component can be found in resultobject$loadings
g1 <- pca.h$loadings[,1]
g2 <- pca.h$loadings[,2]
# get factor scores (individual value of the two components) for family 3
sum(ddf[3,2:3] * g1)## [1] 258.8sum(ddf[3,2:3] * g2)## [1] -2.29# predict family 3 component 1 manually
# multiply difference from mean by factor loading
sum((ddf[3,sizes] - apply(ddf[sizes], 2, mean)) * g1)## [1] -2.435# check against scores provided by princomp
princomp(ddf[sizes], scores=T)$scores[3,]## Comp.1 Comp.2
## -2.435 4.207# add factor scores to data.frame
ddf <- cbind(ddf, princomp(ddf[sizes], scores=T)$scores)
# a visualization via included plot()
plot(pca.h)
# you may apply what you have found in pca to new data
ddf.new <- data.frame(
hlfs = c(166, 168, 185),
hlss = c(157, 170, 190)
)
predict(pca.h, ddf.new)## Comp.1 Comp.2
## [1,] -33.016 -4.383
## [2,] -22.258 3.184
## [3,] 3.942 4.788rm(ddf.new)
# what is different using prcomp()
# using principal() of library(psych)
require(psych)## Loading required package: psychpca.h2 <- principal(ddf[sizes], nfactors=2, rotate="none")
pca.h2## Principal Components Analysis
## Call: principal(r = ddf[sizes], nfactors = 2, rotate = "none")
## Standardized loadings (pattern matrix) based upon correlation matrix
## PC1 PC2 h2 u2
## hlfs 0.92 -0.38 1 0
## hlss 0.92 0.38 1 0
##
## PC1 PC2
## SS loadings 1.71 0.29
## Proportion Var 0.86 0.14
## Cumulative Var 0.86 1.00
## Proportion Explained 0.86 0.14
## Cumulative Proportion 0.86 1.00
##
## Test of the hypothesis that 2 components are sufficient.
##
## The degrees of freedom for the null model are 1 and the objective function was 0.7
## The degrees of freedom for the model are -2 and the objective function was 0
## The total number of observations was 25 with MLE Chi Square = 0 with prob < NA
##
## Fit based upon off diagonal values = 1# component scores are available
head(pca.h2$scores)## PC1 PC2
## [1,] 0.03180 -1.3449
## [2,] 1.43792 0.9972
## [3,] -0.19894 0.7876
## [4,] 0.07336 0.9111
## [5,] -1.22967 -0.3723
## [6,] 1.67326 -1.9322‘Loadings’ zeigt die Gewichte, mit denen die Komponenten-Werte gebildet werden können. ‘Proportion Var’ ist der Anteil an erklärter Varianz, den die jeweilige Komponente erklärt. ‘Cumulative Var’ ist der kumulierte Anteil an der Varianz, der durch alle höheren Komponenten inkl. der aktuellen erkärt wird.
Ein Komponenten-Wert ist eine Linearkombination aus den Rohwerten der Variablen und den Gewichten (loadings) für eine Komponente. Konkret wird oben die Abweichungen einer Beobachtung (Vp) vom Gruppenmittelwert multipliziert mit dem zugehörigen loading. Dies passiert für alle Variablen in der Komponente und wird aufaddiert.
Die Komponenten-Werte können im Quell-Data-Frame mit abgespeichert werden. Dies ist u. U. sinnvoll, wenn man die Anzahl der Variablen verringern möchte und mit den Komponenten weiterrechnen will.
Komponenten-Gewichte (Koeffizienten) können auch auf neue Daten angewendet werden. Dazu müssen die Gewichte von einem Data-Frame generiert worden sein und die Variablen (Columns) müssen gleich benannt sein.
Beipiel EI: 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 erfundener Datensatz findet sich im üblichen tab-delimited Textformat unter: fa-ei.txt
Der Datensatz enthält Antworten von 100 Beobachtungen auf 12 Fragen. Die Antworten sind skaliert von 1 bis 4 (fast nie/manchmal/oft/fast immer)
Die Fragestellung hier sei, ob sich die Antworten auf wenige Komponenten reduzieren lassen, damit Berechnungen anzustellen.
Das Beispiel EI wird auch bei der explorativen FA (EFA) wieder aufgegriffen und aus der Perspektive der EFA dargestellt Rmd html.
Ansatz mit princom()
# get data
ddf <- read.delim("http://r.psych.bio.uni-goettingen.de/mv/data/virt/v_ei.txt")
# library(psych) is needed
# fit model, run PCA
# use base package princomp()
fit <- princomp(ddf, cor=FALSE)
# simple scree plot
plot(fit,type="lines")
# print variance accounted for
summary(fit)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
## Standard deviation 1.9386 1.2467 0.92256 0.87556 0.7659 0.75406
## Proportion of Variance 0.3645 0.1508 0.08255 0.07435 0.0569 0.05515
## Cumulative Proportion 0.3645 0.5152 0.59778 0.67213 0.7290 0.78417
## Comp.7 Comp.8 Comp.9 Comp.10 Comp.11 Comp.12
## Standard deviation 0.71095 0.68051 0.6382 0.59496 0.50532 0.4902
## Proportion of Variance 0.04902 0.04491 0.0395 0.03433 0.02477 0.0233
## Cumulative Proportion 0.83319 0.87810 0.9176 0.95193 0.97670 1.0000# component scores (principal components) are available and may be used for further computations
head(round(fit$scores, 2))## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7 Comp.8 Comp.9
## [1,] 0.67 1.07 -0.45 -1.03 1.36 -0.75 0.15 0.77 0.24
## [2,] 1.97 1.03 0.05 -0.63 -0.15 -0.83 -0.84 0.36 -0.26
## [3,] 0.22 2.09 -0.16 -0.76 0.81 -0.93 -0.27 -0.05 1.15
## [4,] -1.34 -1.75 0.00 -2.14 0.78 0.14 0.65 1.04 -0.50
## [5,] 1.38 1.96 0.25 -0.37 0.09 1.66 -0.37 0.79 0.33
## [6,] 0.08 -0.32 0.84 1.15 -0.09 0.06 -0.35 -0.04 1.02
## Comp.10 Comp.11 Comp.12
## [1,] -0.35 0.11 -0.14
## [2,] -0.53 0.46 0.20
## [3,] 0.06 0.55 0.30
## [4,] -0.52 -0.21 0.54
## [5,] 0.18 0.09 -1.80
## [6,] 0.26 -0.69 -0.05# 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 fa() aus der library(psych)
Hier mit einigen Erklärungen und Erläuterungen.
ddf <- read.delim("http://r.psych.bio.uni-goettingen.de/mv/data/virt/v_ei.txt")
# library(psych) is needed
library(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
principal(ddf, rotate='none', nfactors=ncol(ddf))## Principal Components Analysis
## Call: principal(r = ddf, nfactors = ncol(ddf), rotate = "none")
## Standardized loadings (pattern matrix) based upon correlation matrix
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11
## i1 0.37 0.56 -0.51 0.18 -0.20 0.06 0.10 0.41 0.05 0.06 -0.14
## i2 0.78 -0.22 0.06 0.12 0.01 0.07 0.31 0.07 0.16 -0.34 0.27
## i3 0.38 0.58 -0.08 -0.18 0.45 -0.42 0.07 0.05 -0.29 -0.09 0.06
## i4 -0.71 0.34 0.01 0.04 0.18 0.08 0.43 0.02 0.28 0.03 -0.02
## i5 0.37 0.48 0.51 -0.24 0.11 0.39 -0.25 0.24 0.05 -0.14 -0.08
## i6 -0.72 0.31 0.32 -0.05 -0.04 0.05 0.33 -0.15 -0.13 -0.05 -0.13
## i7 0.28 0.28 0.40 0.79 -0.12 -0.15 -0.04 -0.04 -0.09 0.02 -0.02
## i8 0.77 -0.09 0.21 -0.05 0.21 0.12 0.18 0.03 0.03 0.50 0.10
## i9 0.42 0.53 -0.34 0.11 0.09 0.38 -0.08 -0.50 -0.02 -0.03 0.04
## i10 -0.78 0.22 -0.03 0.01 -0.21 0.24 0.00 0.17 -0.31 0.07 0.34
## i11 0.36 0.53 0.16 -0.35 -0.53 -0.30 0.01 -0.16 0.17 0.06 0.10
## i12 -0.70 0.32 -0.01 0.16 0.26 -0.16 -0.33 0.03 0.33 0.05 0.17
## PC12 h2 u2
## i1 0.13 1 -2.2e-16
## i2 0.11 1 3.3e-16
## i3 -0.04 1 1.0e-15
## i4 -0.27 1 1.6e-15
## i5 -0.05 1 4.4e-16
## i6 0.32 1 -1.6e-15
## i7 -0.09 1 2.6e-15
## i8 0.09 1 1.4e-15
## i9 0.00 1 8.9e-16
## i10 -0.05 1 1.3e-15
## i11 -0.04 1 1.1e-16
## i12 0.19 1 4.4e-16
##
## PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10
## SS loadings 4.10 1.95 0.98 0.94 0.74 0.70 0.62 0.57 0.45 0.40
## Proportion Var 0.34 0.16 0.08 0.08 0.06 0.06 0.05 0.05 0.04 0.03
## Cumulative Var 0.34 0.50 0.59 0.66 0.73 0.78 0.84 0.88 0.92 0.95
## Proportion Explained 0.34 0.16 0.08 0.08 0.06 0.06 0.05 0.05 0.04 0.03
## Cumulative Proportion 0.34 0.50 0.59 0.66 0.73 0.78 0.84 0.88 0.92 0.95
## PC11 PC12
## SS loadings 0.29 0.27
## Proportion Var 0.02 0.02
## Cumulative Var 0.98 1.00
## Proportion Explained 0.02 0.02
## Cumulative Proportion 0.98 1.00
##
## Test of the hypothesis that 12 components are sufficient.
##
## The degrees of freedom for the null model are 66 and the objective function was 3.98
## The degrees of freedom for the model are -12 and the objective function was 0
## The total number of observations was 100 with MLE Chi Square = 0 with prob < NA
##
## Fit based upon off diagonal values = 1# what suggests fa.parallel() of library(psych)
items.parallel <- fa.parallel(ddf, fa="pc")## Loading required package: parallel
## Loading required package: MASS
## Parallel analysis suggests that the number of factors = NA and the number of components = 2# we go on with 2 components
pc.2 <- principal(ddf, rotate='none', nfactors=2)
# default output
pc.2## Principal Components Analysis
## Call: principal(r = ddf, nfactors = 2, rotate = "none")
## Standardized loadings (pattern matrix) based upon correlation matrix
## PC1 PC2 h2 u2
## i1 0.37 0.56 0.45 0.55
## i2 0.78 -0.22 0.65 0.35
## i3 0.38 0.58 0.49 0.51
## i4 -0.71 0.34 0.61 0.39
## i5 0.37 0.48 0.36 0.64
## i6 -0.72 0.31 0.61 0.39
## i7 0.28 0.28 0.16 0.84
## i8 0.77 -0.09 0.60 0.40
## i9 0.42 0.53 0.46 0.54
## i10 -0.78 0.22 0.65 0.35
## i11 0.36 0.53 0.41 0.59
## i12 -0.70 0.32 0.59 0.41
##
## PC1 PC2
## SS loadings 4.10 1.95
## Proportion Var 0.34 0.16
## Cumulative Var 0.34 0.50
## Proportion Explained 0.68 0.32
## Cumulative Proportion 0.68 1.00
##
## Test of the hypothesis that 2 components are sufficient.
##
## The degrees of freedom for the null model are 66 and the objective function was 3.98
## The degrees of freedom for the model are 43 and the objective function was 0.77
## The total number of observations was 100 with MLE Chi Square = 71.03 with prob < 0.0046
##
## Fit based upon off diagonal values = 0.93# graph it
plot(pc.2)
# get loadings
# alternatively available via fit$loadings
loadings(pc.2)##
## Loadings:
## PC1 PC2
## i1 0.368 0.560
## i2 0.777 -0.217
## i3 0.383 0.583
## i4 -0.705 0.342
## i5 0.366 0.480
## i6 -0.718 0.314
## i7 0.276 0.282
## i8 0.769
## i9 0.417 0.535
## i10 -0.776 0.223
## i11 0.361 0.528
## i12 -0.699 0.318
##
## PC1 PC2
## SS loadings 4.096 1.949
## Proportion Var 0.341 0.162
## Cumulative Var 0.341 0.504# 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= .4)##
## Loadings:
## PC1 PC2
## i1 0.560
## i2 0.777
## i3 0.583
## i4 -0.705
## i5 0.480
## i6 -0.718
## i7
## i8 0.769
## i9 0.417 0.535
## i10 -0.776
## i11 0.528
## i12 -0.699
##
## PC1 PC2
## SS loadings 4.096 1.949
## Proportion Var 0.341 0.162
## Cumulative Var 0.341 0.504# does varimax rotation help?
pc.2.r <- principal(ddf, nfactors=2, rotate='varimax')## Loading required package: GPArotationprint(loadings(pc.2.r), cutoff= .4)##
## Loadings:
## RC1 RC2
## i1 0.667
## i2 -0.789
## i3 0.694
## i4 0.784
## i5 0.596
## i6 0.782
## i7
## i8 -0.721
## i9 0.667
## i10 0.790
## i11 0.636
## i12 0.767
##
## RC1 RC2
## SS loadings 3.631 2.414
## Proportion Var 0.303 0.201
## Cumulative Var 0.303 0.504plot(pc.2.r)
# a quick look at some factor scores
# unrotated
head(round(pc.2$scores, 2))## PC1 PC2
## [1,] 0.37 0.81
## [2,] 1.02 0.63
## [3,] 0.22 1.59
## [4,] -0.85 -1.41
## [5,] 0.78 1.35
## [6,] 0.03 -0.21# component scores are different for rotated components
head(round(pc.2.r$scores, 2))## RC1 RC2
## [1,] 0.05 0.89
## [2,] -0.61 1.03
## [3,] 0.55 1.51
## [4,] 0.10 -1.64
## [5,] -0.06 1.56
## [6,] -0.12 -0.17# component scores may be saved to dataframe for further computations
ddf <- cbind(ddf, pc.2.r$scores)
head(ddf)## i1 i2 i3 i4 i5 i6 i7 i8 i9 i10 i11 i12 RC1 RC2
## 1 4 3 4 2 3 1 4 3 3 1 3 3 0.04795 0.8910
## 2 4 3 3 1 4 1 4 4 4 1 3 2 -0.60774 1.0320
## 3 4 3 4 2 4 2 4 2 4 1 3 3 0.54579 1.5066
## 4 3 3 2 2 1 2 4 3 1 2 2 3 0.09548 -1.6404
## 5 4 4 4 4 4 1 4 4 4 2 3 1 -0.06131 1.5556
## 6 2 4 3 2 4 2 2 3 3 2 3 2 -0.12486 -0.1744in Output unter Standardized loadings based ...
h2 sind die Kommunalitäten (Anteil an Varianz aller Variablen, der durch die beiden Komponenten erklärt/repliziert wird). u2 sind die uniqueness-Werte (siehe EFA).
Eine Erhöhung des cutoffs macht die Tabelle der Factor-Loadings klarer:
Item i7 ist schlecht einer der beiden Komponenten zuordenbar. Item i9 lädt hoch auf beiden Komponenten. Der Rest ist gut bereits unrotiert einer der beiden Komponenten zuordenbar. Zwei Komponenten erklären 50% der Varianz. Andersherum gesehen gehen 50% der Information verloren, wenn man mit 2 Komponenten anstelle von 12 Items in einem multivariaten Verfahren weiterrechnen würde.
Nach Varimax-Rotation findet auch Item i9 ‘seine’ Komponente
Mit den gespeicherten Komponenten-Werten kann man dann weiterrechnen.
Wichtig: Alles hier passiert ohne inhaltlichen Bezug. Die inhaltliche Interpretation der Komponenten ist beim Ziel der Datenreduktion nicht nötig. Hat man aber mit Komponenten gerechnete Ergebnisse, ist schon fraglich, was sie eigentlich aussagen.
Ausflug Reproduzierte Korrelationsmatrix und Residuen, Chi^2 Modelltest
ddf <- read.delim("http://r.psych.bio.uni-goettingen.de/mv/data/virt/v_ei.txt")
library(psych)
pc.2 <- principal(ddf, rotate='none', nfactors=2)
# original correlation table, that will be reproduced more or less by PCA
head(round(cor(ddf), 2))## i1 i2 i3 i4 i5 i6 i7 i8 i9 i10 i11 i12
## i1 1.00 0.18 0.35 -0.06 0.17 -0.23 0.18 0.14 0.43 -0.08 0.30 -0.11
## i2 0.18 1.00 0.14 -0.48 0.18 -0.52 0.21 0.57 0.18 -0.61 0.12 -0.60
## i3 0.35 0.14 1.00 -0.07 0.30 -0.11 0.13 0.25 0.34 -0.25 0.33 -0.05
## i4 -0.06 -0.48 -0.07 1.00 -0.12 0.63 -0.12 -0.45 -0.11 0.53 -0.15 0.54
## i5 0.17 0.18 0.30 -0.12 1.00 -0.04 0.18 0.31 0.26 -0.13 0.33 -0.11
## i6 -0.23 -0.52 -0.11 0.63 -0.04 1.00 -0.04 -0.47 -0.18 0.58 -0.04 0.45# get the reproduced correlation matrix, which is calculated on base of the loadings
# diagonal are the communalities
head(round(factor.model(pc.2$loadings), 2))## i1 i2 i3 i4 i5 i6 i7 i8 i9 i10 i11 i12
## i1 0.45 0.16 0.47 -0.07 0.40 -0.09 0.26 0.23 0.45 -0.16 0.43 -0.08
## i2 0.16 0.65 0.17 -0.62 0.18 -0.63 0.15 0.62 0.21 -0.65 0.17 -0.61
## i3 0.47 0.17 0.49 -0.07 0.42 -0.09 0.27 0.24 0.47 -0.17 0.45 -0.08
## i4 -0.07 -0.62 -0.07 0.61 -0.09 0.61 -0.10 -0.57 -0.11 0.62 -0.07 0.60
## i5 0.40 0.18 0.42 -0.09 0.36 -0.11 0.24 0.24 0.41 -0.18 0.39 -0.10
## i6 -0.09 -0.63 -0.09 0.61 -0.11 0.61 -0.11 -0.58 -0.13 0.63 -0.09 0.60# get the residuals. These are the differences between original and reproduced correlation matrix
pc.diff <- cor(ddf) - factor.model(pc.2$loadings)
head(round(pc.diff, 2))## i1 i2 i3 i4 i5 i6 i7 i8 i9 i10 i11 i12
## i1 0.55 0.02 -0.12 0.01 -0.23 -0.14 -0.08 -0.09 -0.02 0.08 -0.13 -0.03
## i2 0.02 0.35 -0.03 0.15 0.00 0.11 0.06 -0.05 -0.02 0.04 -0.04 0.02
## i3 -0.12 -0.03 0.51 0.00 -0.12 -0.01 -0.14 0.00 -0.14 -0.08 -0.11 0.04
## i4 0.01 0.15 0.00 0.39 -0.03 0.02 -0.02 0.12 0.00 -0.09 -0.07 -0.06
## i5 -0.23 0.00 -0.12 -0.03 0.64 0.07 -0.06 0.07 -0.15 0.05 -0.05 0.00
## i6 -0.14 0.11 -0.01 0.02 0.07 0.39 0.07 0.11 -0.05 -0.04 0.05 -0.15rm(pc.diff)
# get the same via call to factor.residuals()
# diagonal are the uniquenesses.
pc.diff.2 <- factor.residuals(cor(ddf), pc.2$loadings)
head(round(pc.diff.2, 2))## i1 i2 i3 i4 i5 i6 i7 i8 i9 i10 i11 i12
## i1 0.55 0.02 -0.12 0.01 -0.23 -0.14 -0.08 -0.09 -0.02 0.08 -0.13 -0.03
## i2 0.02 0.35 -0.03 0.15 0.00 0.11 0.06 -0.05 -0.02 0.04 -0.04 0.02
## i3 -0.12 -0.03 0.51 0.00 -0.12 -0.01 -0.14 0.00 -0.14 -0.08 -0.11 0.04
## i4 0.01 0.15 0.00 0.39 -0.03 0.02 -0.02 0.12 0.00 -0.09 -0.07 -0.06
## i5 -0.23 0.00 -0.12 -0.03 0.64 0.07 -0.06 0.07 -0.15 0.05 -0.05 0.00
## i6 -0.14 0.11 -0.01 0.02 0.07 0.39 0.07 0.11 -0.05 -0.04 0.05 -0.15rm(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.
Beispiel MSQ
Quelle Inspiration für Beispiel Quelle
A premature consensus: Are happiness and sadness truly opposite affects? Motivation and Emotion, 30, 1, 1-12.). The data set includes EPI and BFI scale scores (see previous examples).
#specify the name and address of the remote file
datafilename="http://personality-project.org/r/datasets/Maps.mixx.msq1.epi.bf.txt"
#note that it is also available as built in example in the psych package named msq
rawmsq =read.table(datafilename,header=TRUE) #read the data file
head(rawmsq)## id delighted sociable jittery hostile sluggish depressed satisfied
## 1 1 3 3 2 0 0 0 3
## 2 2 1 2 0 0 0 0 2
## 3 3 1 0 0 0 2 0 1
## 4 4 1 1 0 0 0 1 1
## 5 5 0 0 1 2 0 2 0
## 6 6 0 1 1 0 0 0 1
## relaxed warmhearted blue intense strong scared enthusiastic proud sad
## 1 3 3 0 1 3 0 3 3 0
## 2 2 2 0 0 2 1 1 1 0
## 3 1 1 0 0 2 0 0 0 0
## 4 2 2 1 1 0 0 0 0 0
## 5 0 0 2 2 2 0 0 2 2
## 6 1 0 0 1 1 0 0 0 0
## active full.of.pep unhappy lively aroused elated wakeful placid
## 1 3 3 0 3 2 2 3 0
## 2 1 0 0 1 2 0 1 0
## 3 0 0 1 0 0 0 0 0
## 4 1 0 0 0 0 0 0 0
## 5 0 1 3 0 2 0 0 0
## 6 1 0 0 0 0 0 0 2
## irritable at.rest clutched.up tired vigorous astonished dull fearful
## 1 0 2 0 1 1 0 0 0
## 2 0 1 0 1 1 0 1 0
## 3 2 1 0 3 0 0 1 0
## 4 0 1 0 2 0 0 1 0
## 5 2 0 2 0 0 9 1 0
## 6 0 0 0 1 0 0 1 0
## at.ease interested afraid bored guilty quiescent scornful determined
## 1 3 2 0 1 0 1 3 3
## 2 2 1 1 1 1 0 2 2
## 3 1 1 0 1 0 0 0 0
## 4 1 0 0 1 0 0 2 0
## 5 0 0 0 1 2 2 0 3
## 6 2 0 0 1 0 1 1 0
## excited gloomy still angry nervous kindly grouchy upset alone calm
## 1 2 0 2 0 0 3 0 0 0 2
## 2 1 0 1 0 1 1 0 0 0 2
## 3 0 0 2 0 0 0 1 0 0 1
## 4 0 0 0 0 0 0 0 0 2 2
## 5 0 0 3 3 0 0 2 2 0 0
## 6 0 0 1 0 0 0 0 0 0 2
## drowsy alert inspired surprised energetic lonely quiet sorry sleepy
## 1 0 3 2 0 2 0 0 0 1
## 2 1 0 0 0 0 0 1 0 2
## 3 2 1 0 0 0 1 1 0 3
## 4 0 0 0 0 0 0 1 0 3
## 5 0 0 0 0 0 1 0 0 0
## 6 0 2 0 0 0 0 1 0 0
## pleased happy distressed attentive wide.awake frustrated serene
## 1 3 3 0 2 2 0 3
## 2 1 1 0 1 0 0 1
## 3 0 0 0 0 0 0 1
## 4 1 1 0 0 0 0 0
## 5 0 0 2 0 2 2 0
## 6 0 0 0 1 0 0 1
## confident content tense ashamed anxious idle epiE epiS epiImp epilie
## 1 3 3 0 0 0 0 18 10 7 3
## 2 2 1 0 0 0 1 16 8 5 1
## 3 1 1 1 0 0 1 6 1 3 2
## 4 0 0 0 0 1 0 12 6 4 3
## 5 2 0 2 2 2 0 14 6 5 3
## 6 2 1 1 0 0 1 6 4 2 5
## epiNeur bfagree bfcon bfext bfneur bfopen bdi traitanx stateanx
## 1 9 138 96 141 51 138 1 24 22
## 2 12 101 99 107 116 132 7 41 40
## 3 5 143 118 38 68 90 4 37 44
## 4 15 104 106 64 114 101 8 54 40
## 5 2 115 102 103 86 118 8 39 67
## 6 15 110 113 61 54 149 5 51 38msq=rawmsq[,2:72] #select the subset of items in the MSQ
msq[msq=="9"] = NA # change all occurences of 9 to be missing values
msq <- data.frame(msq) #convert the input matrix into a data frame for easier manipulation
names(msq) #what are the variables?## [1] "delighted" "sociable" "jittery" "hostile"
## [5] "sluggish" "depressed" "satisfied" "relaxed"
## [9] "warmhearted" "blue" "intense" "strong"
## [13] "scared" "enthusiastic" "proud" "sad"
## [17] "active" "full.of.pep" "unhappy" "lively"
## [21] "aroused" "elated" "wakeful" "placid"
## [25] "irritable" "at.rest" "clutched.up" "tired"
## [29] "vigorous" "astonished" "dull" "fearful"
## [33] "at.ease" "interested" "afraid" "bored"
## [37] "guilty" "quiescent" "scornful" "determined"
## [41] "excited" "gloomy" "still" "angry"
## [45] "nervous" "kindly" "grouchy" "upset"
## [49] "alone" "calm" "drowsy" "alert"
## [53] "inspired" "surprised" "energetic" "lonely"
## [57] "quiet" "sorry" "sleepy" "pleased"
## [61] "happy" "distressed" "attentive" "wide.awake"
## [65] "frustrated" "serene" "confident" "content"
## [69] "tense" "ashamed" "anxious"# how about missings
cat(length(which(is.na(msq)==T)), ' of ', nrow(msq), ' have missings')## 51 of 231 have missingsrequire(psych)
cleaned <- na.omit(msq) #remove the cases with missing values
head(describe(cleaned)) #basic summary statistics -- check for miscodings## vars n mean sd median trimmed mad min max range skew
## delighted 1 197 0.78 0.90 1 0.67 1.48 0 3 3 0.79
## sociable 2 197 1.32 0.96 1 1.28 1.48 0 3 3 0.06
## jittery 3 197 0.53 0.73 0 0.40 0.00 0 3 3 1.21
## hostile 4 197 0.29 0.59 0 0.16 0.00 0 3 3 2.18
## sluggish 5 197 1.26 0.97 1 1.20 1.48 0 3 3 0.41
## depressed 6 197 0.55 0.81 0 0.39 0.00 0 3 3 1.45
## kurtosis se
## delighted -0.53 0.06
## sociable -1.02 0.07
## jittery 0.74 0.05
## hostile 4.75 0.04
## sluggish -0.79 0.07
## depressed 1.39 0.06# show correlation matrix
cor.cleaned <- cor(cleaned)
# take a look at it
head(round(cor.cleaned, 2))## delighted sociable jittery hostile sluggish depressed satisfied
## delighted 1.00 0.63 0.13 -0.10 -0.27 -0.22 0.61
## sociable 0.63 1.00 0.15 -0.18 -0.35 -0.30 0.52
## jittery 0.13 0.15 1.00 0.20 -0.03 0.18 0.00
## hostile -0.10 -0.18 0.20 1.00 0.21 0.41 -0.12
## sluggish -0.27 -0.35 -0.03 0.21 1.00 0.27 -0.24
## depressed -0.22 -0.30 0.18 0.41 0.27 1.00 -0.37
## relaxed warmhearted blue intense strong scared enthusiastic
## delighted 0.41 0.64 -0.20 0.23 0.47 -0.07 0.69
## sociable 0.42 0.65 -0.22 0.14 0.50 -0.01 0.66
## jittery -0.22 0.06 0.21 0.37 0.08 0.41 0.20
## hostile -0.21 -0.21 0.37 0.23 -0.11 0.26 -0.11
## sluggish -0.15 -0.25 0.21 -0.05 -0.24 0.10 -0.37
## depressed -0.39 -0.23 0.76 0.30 -0.16 0.42 -0.24
## proud sad active full.of.pep unhappy lively aroused elated
## delighted 0.49 -0.21 0.58 0.61 -0.26 0.60 0.50 0.66
## sociable 0.43 -0.25 0.61 0.58 -0.31 0.62 0.52 0.56
## jittery 0.10 0.25 0.23 0.15 0.15 0.16 0.22 0.11
## hostile -0.13 0.42 -0.10 -0.13 0.38 -0.12 -0.04 -0.08
## sluggish -0.21 0.30 -0.33 -0.40 0.33 -0.39 -0.33 -0.31
## depressed -0.22 0.81 -0.16 -0.29 0.73 -0.31 -0.14 -0.23
## wakeful placid irritable at.rest clutched.up tired vigorous
## delighted 0.38 -0.19 -0.24 0.20 -0.08 -0.21 0.52
## sociable 0.48 -0.12 -0.30 0.22 -0.20 -0.27 0.49
## jittery 0.12 -0.06 0.11 -0.15 0.37 0.07 0.21
## hostile -0.13 -0.05 0.56 -0.22 0.33 0.26 -0.06
## sluggish -0.46 0.17 0.37 -0.12 0.33 0.66 -0.30
## depressed -0.23 -0.03 0.40 -0.28 0.41 0.24 -0.17
## astonished dull fearful at.ease interested afraid bored guilty
## delighted 0.25 -0.35 -0.08 0.32 0.45 -0.14 -0.20 -0.08
## sociable 0.22 -0.44 -0.01 0.37 0.46 -0.07 -0.20 -0.02
## jittery 0.26 0.02 0.34 -0.25 0.14 0.36 -0.11 0.22
## hostile 0.16 0.31 0.23 -0.32 -0.08 0.27 0.20 0.27
## sluggish -0.12 0.46 0.10 -0.16 -0.28 0.12 0.28 0.09
## depressed 0.09 0.28 0.40 -0.38 -0.21 0.38 -0.01 0.36
## quiescent scornful determined excited gloomy still angry nervous
## delighted 0.07 0.12 0.43 0.60 -0.19 -0.07 -0.11 0.08
## sociable -0.03 0.17 0.42 0.54 -0.27 -0.13 -0.18 0.07
## jittery -0.10 -0.15 0.24 0.19 0.26 -0.05 0.25 0.56
## hostile -0.01 -0.11 0.14 -0.07 0.36 0.03 0.49 0.24
## sluggish 0.08 -0.09 -0.18 -0.30 0.32 0.17 0.11 0.02
## depressed -0.01 -0.24 0.01 -0.19 0.65 -0.02 0.60 0.34
## kindly grouchy upset alone calm drowsy alert inspired surprised
## delighted 0.64 -0.25 -0.24 -0.33 0.15 -0.19 0.46 0.40 0.28
## sociable 0.59 -0.33 -0.24 -0.33 0.21 -0.23 0.44 0.37 0.23
## jittery 0.14 0.10 0.23 -0.11 -0.15 -0.04 0.06 0.30 0.18
## hostile -0.08 0.45 0.33 0.07 -0.29 0.21 -0.11 0.05 0.06
## sluggish -0.36 0.40 0.23 0.48 -0.10 0.62 -0.44 -0.17 -0.12
## depressed -0.30 0.37 0.62 0.12 -0.35 0.12 -0.21 -0.05 -0.02
## energetic lonely quiet sorry sleepy pleased happy distressed
## delighted 0.52 -0.09 -0.26 -0.10 -0.23 0.54 0.64 -0.17
## sociable 0.50 -0.16 -0.34 -0.21 -0.29 0.53 0.60 -0.21
## jittery 0.16 0.20 -0.11 0.07 -0.03 0.06 0.07 0.27
## hostile -0.09 0.42 0.08 0.33 0.18 -0.13 -0.19 0.33
## sluggish -0.40 0.15 0.30 0.20 0.62 -0.29 -0.38 0.20
## depressed -0.18 0.59 0.19 0.44 0.14 -0.28 -0.37 0.56
## attentive wide.awake frustrated serene confident content tense
## delighted 0.37 0.42 -0.17 0.11 0.41 0.50 -0.07
## sociable 0.35 0.43 -0.18 0.13 0.53 0.51 -0.07
## jittery 0.04 0.08 0.25 -0.06 0.03 -0.04 0.48
## hostile -0.11 -0.13 0.44 -0.08 -0.11 -0.19 0.27
## sluggish -0.32 -0.48 0.24 -0.05 -0.29 -0.28 0.13
## depressed -0.22 -0.18 0.59 -0.19 -0.42 -0.42 0.40
## ashamed anxious
## delighted -0.06 0.04
## sociable -0.11 0.02
## jittery 0.27 0.41
## hostile 0.36 0.22
## sluggish 0.12 0.07
## depressed 0.44 0.35# Ask Bartlett test whether correlation matrix is idendity matrix (all elements are 0), should be significant
# run Bartlett test on dataframe
cortest.bartlett(cleaned)## R was not square, finding R from data## $chisq
## [1] 11556
##
## $p.value
## [1] 0
##
## $df
## [1] 2485# alternativly on correlation matrix (add n)
cortest.bartlett(cor.cleaned, n=nrow(cleaned))## $chisq
## [1] 11556
##
## $p.value
## [1] 0
##
## $df
## [1] 2485# KMO
KMO(cleaned)## Kaiser-Meyer-Olkin factor adequacy
## Call: KMO(r = cleaned)
## Overall MSA = 0.92
## MSA for each item =
## delighted sociable jittery hostile sluggish
## 0.93 0.95 0.87 0.83 0.91
## depressed satisfied relaxed warmhearted blue
## 0.91 0.94 0.91 0.94 0.89
## intense strong scared enthusiastic proud
## 0.87 0.93 0.86 0.95 0.93
## sad active full.of.pep unhappy lively
## 0.92 0.95 0.95 0.92 0.95
## aroused elated wakeful placid irritable
## 0.95 0.95 0.95 0.71 0.92
## at.rest clutched.up tired vigorous astonished
## 0.90 0.92 0.87 0.95 0.89
## dull fearful at.ease interested afraid
## 0.94 0.83 0.95 0.94 0.86
## bored guilty quiescent scornful determined
## 0.81 0.87 0.65 0.86 0.94
## excited gloomy still angry nervous
## 0.95 0.93 0.64 0.88 0.86
## kindly grouchy upset alone calm
## 0.96 0.92 0.89 0.87 0.90
## drowsy alert inspired surprised energetic
## 0.88 0.93 0.92 0.84 0.95
## lonely quiet sorry sleepy pleased
## 0.90 0.82 0.88 0.85 0.96
## happy distressed attentive wide.awake frustrated
## 0.96 0.90 0.95 0.94 0.89
## serene confident content tense ashamed
## 0.88 0.95 0.95 0.91 0.88
## anxious
## 0.87# determinant should be higher than 0.00001
det(cor.cleaned)## [1] 5.463e-30pca.cleaned <- principal(cleaned, nfactors=2, rotate = 'varimax', scores=TRUE)
f2 <- fa(cleaned,2,rotate="varimax") #factor analyze the resulting item
#(f2) #show the result
load=loadings(f2)
print(load,sort=TRUE,digits=2,cutoff=0.01) #show the loadings ##
## Loadings:
## MR1 MR2
## delighted 0.70 -0.10
## sociable 0.69 -0.15
## satisfied 0.68 -0.32
## warmhearted 0.70 -0.15
## strong 0.68
## enthusiastic 0.83
## proud 0.70 -0.05
## active 0.79 0.05
## full.of.pep 0.84 -0.09
## lively 0.86 -0.08
## aroused 0.74 0.06
## elated 0.80 -0.05
## wakeful 0.67 -0.05
## vigorous 0.77 0.05
## interested 0.68 0.03
## determined 0.69 0.30
## excited 0.80 0.04
## kindly 0.82 -0.09
## alert 0.71 -0.07
## inspired 0.63 0.24
## energetic 0.84 0.02
## pleased 0.74 -0.17
## happy 0.83 -0.19
## attentive 0.67 -0.06
## wide.awake 0.66 -0.10
## confident 0.66 -0.30
## content 0.72 -0.33
## depressed -0.27 0.71
## blue -0.20 0.67
## intense 0.36 0.51
## scared 0.10 0.70
## sad -0.20 0.76
## unhappy -0.29 0.71
## irritable -0.30 0.57
## clutched.up -0.07 0.64
## fearful 0.10 0.67
## at.ease 0.40 -0.54
## afraid 0.05 0.70
## gloomy -0.25 0.68
## angry -0.11 0.74
## nervous 0.19 0.67
## grouchy -0.25 0.52
## upset -0.16 0.75
## lonely -0.06 0.61
## sorry -0.08 0.60
## distressed -0.12 0.76
## frustrated -0.15 0.77
## tense 0.08 0.75
## ashamed 0.68
## anxious 0.18 0.63
## jittery 0.22 0.44
## hostile -0.11 0.49
## sluggish -0.44 0.22
## relaxed 0.44 -0.48
## placid -0.12 -0.14
## at.rest 0.38 -0.35
## tired -0.39 0.24
## astonished 0.43 0.40
## dull -0.47 0.25
## bored -0.32
## guilty 0.01 0.50
## quiescent 0.06 -0.09
## scornful 0.22 -0.35
## still -0.09 -0.09
## alone -0.45 0.10
## calm 0.26 -0.43
## drowsy -0.35 0.16
## surprised 0.42 0.24
## quiet -0.29 0.10
## sleepy -0.42 0.16
## serene 0.27 -0.27
##
## MR1 MR2
## SS loadings 17.98 12.93
## Proportion Var 0.25 0.18
## Cumulative Var 0.25 0.44plot(load) #plot factor 1 by 2
identify(load,labels=names(msq)) 
## integer(0) #put names of selected points onto the figure -- to stop, click with command key
plot(f2,labels=names(msq))
Beispiel: Crime Data (Everitt, 2010)
Verbrechensraten in den US-Staaten nach Verbrechenstyp getrennt.
Berechnung auf Basis der Varianz-Covarianzmatrix
d.crime <- read.delim(file="http://r.psych.bio.uni-goettingen.de/mv/data/be/pca_crime.txt")
head(d.crime)## state murder rape robbery assault burglary theft vehicules
## 1 ME 2.0 14.8 28 102 803 2347 164
## 2 NH 2.2 21.5 24 92 755 2208 228
## 3 VT 2.0 21.8 22 103 949 2697 181
## 4 MA 3.6 29.7 193 331 1071 2189 906
## 5 RI 3.5 21.4 119 192 1294 2568 705
## 6 CT 4.6 23.8 192 205 1198 2758 447#figure 10.4
pairs(d.crime[2:8])
#PCA results Table 10.3
crime_pc <- princomp(d.crime[,2:8],cor=FALSE)
summary(crime_pc,loadings=TRUE)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5
## Standard deviation 841.1699 290.1020 174.16302 96.85347 61.10167
## Proportion of Variance 0.8471 0.1008 0.03632 0.01123 0.00447
## Cumulative Proportion 0.8471 0.9479 0.98419 0.99542 0.99989
## Comp.6 Comp.7
## Standard deviation 9.1116055 2.294e+00
## Proportion of Variance 0.0000994 6.299e-06
## Cumulative Proportion 0.9999937 1.000e+00
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## murder 0.998
## rape -0.998
## robbery -0.276 0.361 -0.396 -0.794
## assault -0.106 -0.278 0.106 -0.777 0.544
## burglary -0.427 -0.630 -0.629 0.135
## theft -0.885 0.435 0.163
## vehicules -0.128 -0.510 0.661 0.469 0.256Berechnung auf Basis der Korrelationsmatrix
d.crime <- read.delim(file="http://r.psych.bio.uni-goettingen.de/mv/data/be/pca_crime.txt")
#figure 10.4
pairs(d.crime[2:8])
#PCA results Table 10.3
crime_pc <- princomp(d.crime[,2:8],cor=TRUE)
summary(crime_pc,loadings=TRUE)## Importance of components:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6
## Standard deviation 2.1666 0.9928 0.68138 0.58666 0.48886 0.42270
## Proportion of Variance 0.6706 0.1408 0.06633 0.04917 0.03414 0.02552
## Cumulative Proportion 0.6706 0.8114 0.87771 0.92688 0.96102 0.98655
## Comp.7
## Standard deviation 0.30688
## Proportion of Variance 0.01345
## Cumulative Proportion 1.00000
##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## murder -0.381 -0.350 -0.538 -0.274 0.370 0.480
## rape -0.377 0.279 -0.830 -0.250 -0.151
## robbery -0.391 -0.420 0.131 0.275 -0.387 -0.651
## assault -0.410 -0.124 -0.335 0.564 -0.620
## burglary -0.394 0.367 0.162 0.466 0.622 -0.283
## theft -0.321 0.628 0.449 -0.388 -0.282 0.256
## vehicules -0.366 -0.282 0.758 0.163 0.422Everitt (2010) schlägt als Benennung für die Hauptkomponente 1 ‘dangerousness’ und für die zweite vorsichtig ‘property- vs. person’ vor (+/- scheint diese Verbrechen zu kontrastieren).
property-crimes (theft, burglary) - crimes against the person (robbery, murder)
#Figure 10.5
plot(1:7,crime_pc$sdev^2,type="l",xlab="Component number",ylab="Variance")
#fig 10.6
xlim<-range(crime_pc$scores[,1])
plot(crime_pc$scores[,1],crime_pc$scores[,2],xlab="First PC score",ylab="Second PC score",xlim=xlim,ylim=xlim,type="n")
text(crime_pc$scores[,1],crime_pc$scores[,2],d.crime$state,cex=0.5)
# Eigenwerte sind die quadrierten Standard-Deviations der Hauptkomponenten:
print("Eigenwerte:")## [1] "Eigenwerte:"crime_pc$sdev^2## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## 4.69410 0.98562 0.46428 0.34417 0.23898 0.17867 0.09418# Berechnung der Korrelationen zwischen den Variablen und den Hauptkomponenten
# wird in Statistica 'Faktoren/Variablen-Korr. (Faktorladungen)' bzw 'Faktorkoordinaten der Variablen' genannt.
cor(d.crime[,2:8], crime_pc$scores)## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## murder -0.8249 -0.3473 -0.366562 0.01887 -0.13405 0.15639 0.147209
## rape -0.8173 0.2768 -0.012258 -0.48715 -0.12240 -0.02910 -0.046216
## robbery -0.8461 -0.4166 0.089339 0.16107 -0.18915 -0.03169 -0.199710
## assault -0.8874 -0.1231 -0.228432 0.01853 0.27589 -0.26225 0.002292
## burglary -0.8543 0.3643 -0.008004 0.09494 0.22772 0.26284 -0.086755
## theft -0.6951 0.6232 0.052633 0.26330 -0.18963 -0.11935 0.078489
## vehicules -0.7934 -0.2802 0.516489 -0.04319 0.07955 0.01603 0.129595# diese Faktorladungen bekommt man auch, indem man die loadings mit der Standardabweichung der jeweiligen Komponente multipliziert
# für die erste Komponente z. B.
crime_pc$loadings[,1] * crime_pc$sdev[1]## murder rape robbery assault burglary theft vehicules
## -0.8249 -0.8173 -0.8461 -0.8874 -0.8543 -0.6951 -0.7934# und die wiederum quadriert und über die Variablen aufsummiert ergibt den Eigenwert der jeweiligen Komponente
sum((crime_pc$loadings[,1] * crime_pc$sdev[1])^2)## [1] 4.694# quadrierte Korrelationen zwischen den Variablen und den Hauptkomponenten
# Anteil an der Varianz der jeweiligen Variable, der durch die jeweilige Hauptkomponente erklärt wird.
# ist die Basis dessen, was in Statistica Kommunalität genannt wird.
m.cor <- cor(d.crime[,2:8], crime_pc$scores)
m.cor^2## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## murder 0.6805 0.12064 1.344e-01 0.0003561 0.017970 0.0244574 2.167e-02
## rape 0.6679 0.07662 1.503e-04 0.2373195 0.014981 0.0008467 2.136e-03
## robbery 0.7158 0.17358 7.981e-03 0.0259429 0.035777 0.0010041 3.988e-02
## assault 0.7874 0.01516 5.218e-02 0.0003435 0.076114 0.0687773 5.255e-06
## burglary 0.7297 0.13271 6.406e-05 0.0090145 0.051856 0.0690874 7.526e-03
## theft 0.4831 0.38840 2.770e-03 0.0693286 0.035958 0.0142449 6.160e-03
## vehicules 0.6295 0.07851 2.668e-01 0.0018656 0.006328 0.0002569 1.679e-02# Kommunalitäten
# akkumulierte erklärte Varianz
# wird in Statistica 'Kommunalitäten (Kosinus ^ 2)' genannt.
m.cor <- cor(d.crime[,2:8], crime_pc$scores)
q.m.cor <- m.cor^2
# für die erste Komponente stimmen die akkumulierten Werte
# für die drei wichtigsten Komponenten wäre das:
apply(q.m.cor[,1:3], 1, sum)## murder rape robbery assault burglary theft vehicules
## 0.9355 0.7447 0.8974 0.8548 0.8625 0.8743 0.9748#nvar <- length(m.cor[,1])
#for (i in 1:nvar) {
# print(sum(crime_pc$loadings[1,1:i]^2))
# quadrierte loadings
# in Statistica genannt: 'Beiträge der Variablen'
# Anteil der Varianz der Variablen, der durch die jeweilige Componente erklärt wird
crime_pc$loadings^2##
## Loadings:
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## murder 0.145 0.122 0.289 0.137 0.230
## rape 0.142 0.690
## robbery 0.152 0.176 0.150 0.423
## assault 0.168 0.112 0.318 0.385
## burglary 0.155 0.135 0.217 0.387
## theft 0.103 0.394 0.201 0.150
## vehicules 0.134 0.575 0.178
##
## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## SS loadings 0.145 0.232 0.427 0.522 0.204 0.323 0.275
## Proportion Var 0.021 0.033 0.061 0.075 0.029 0.046 0.039
## Cumulative Var 0.021 0.054 0.115 0.190 0.219 0.265 0.304Berechnung der Scores für die Vpn (hier Bundesstaaten) in den Hauptkomponenten auf Basis der Rohwerte und der loadings.
# Wie kann man die Scores aus den loadings errechnen?
# Multiplikation der Abweichung vom Mittelwert mit dem jeweiligen Ladungsgewicht
# und aufaddieren über alle Variablen
#
# Eine Funktion definieren, die die Standardabweichung n-gewichtet
sd0 <- function(dd) {
return(sqrt(sum((dd - mean(dd)) ^ 2) / length(dd)))
}
# die n-gewichteten Standardabweichungen für alle Variablen berechnen
d.sd0s <- apply(d.crime[2:8], 2, sd0)
# die Mittelwerte für alle Variablen
d.means <- apply(d.crime[2:8], 2, mean)
# die Loadings der ersten Hauptkomponente
c1l <- crime_pc$loadings[,1]
# die erste Vp auslesen
vp <- d.crime[1,2:8]
# und den Score auf Hauptkomponente ermitteln
vp1.score <- sum(((vp - d.means) / d.sd0s) * c1l)
vp1.score## [1] 2.808# vp2
vp <- d.crime[2,2:8]
vp2.score <- sum(((vp - d.means) / d.sd0s) * c1l)
vp2.score## [1] 2.654# entspricht den Werten in crime_pc$scores
crime_pc$scores[1:2,]## Comp.1 Comp.2 Comp.3 Comp.4 Comp.5 Comp.6 Comp.7
## [1,] 2.808 -0.00527 0.07301 0.3561 -0.0147 0.1019 -0.09945
## [2,] 2.654 -0.10890 0.26732 -0.1608 -0.1047 0.1213 -0.02318Beispiel UCLA-Video
Daten nicht erhältlich. Daher nur ein paar Bemerkungen zu CFA und FA
Basis: Quelle
pca <- princomp(factor.data)
summary(pca)
loadings(pca)
plot(pca, type="lines")
biplot(pca)
# exploratory fa
# eigenvalues and scree plot
eigen.values <- eigen(cor(factor.data))
plot(eigen.values$values, type="both")
# exploratory fa
efa <- fatanal(factor.data, 3, rotation="promax")
loadings(efa)Übungen / Exercises
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) Kapitel 3: R Environment