#### Stochastik IV, Übungsblatt 6

### Aufgabe 3

library(GPArotation)

## Teilaufgabe (a) 

# Test "general" (general intelligence) testet anscheinend
# auf allgemeine Intelligenz, was wohl als eine Art Kombination mehrerer 
# verschiedener Fähigkeits- und Intelligenzdimensionen gedeutet werden 
# könnte. Anscheinend werden hier mitunter auch die beiden nächsten 
# Fähigkeitsdimensionen miteingeschlossen.
  
# Die Tests "picture" ("picture completion"; Bildvervollständigung), 
# "blocks" ("block design"; Blockmotive) und "maze" ("mazes"; Labyrinthe) 
# testen anscheinend auf VISUELLE FÄHIGKEIT.

# Die Tests "reading" ("reading comprehension"; Leseverständnis) und "vocab"
# ("vocabulary"; Sprach-/Wortschatz) testen anscheinend auf VERBALE FÄHIGKEIT.

# Somit sollten, anhand der vorliegenden Indikatoren, zwei 'atomare' 
# Fähigkeitsdimensionen, d.h. zwei (latente) Faktoren, 'gemessen' werden; 
# nämlich: VISUELLE FÄHIGKEIT und VERBALE FÄHIGKEIT.

## Teilaufgabe (b)

# Beachte: (1) In function factanal() wird immer die Maximum-Likelihood 
#              Methode verwendet (Rückgabeobjekt-Komponente 'method = "mle"');
#          (2) Die Varimax-Rotation ist Standardeinstellung in factanal() 

# ML FA mit 1 Faktor

(ability.fa.1f <- factanal(factors = 1, covmat = ability.cov))

# Korrelationsmatrix aus den Daten

(daten.R <- round(ability.fa.1f$correlation, 3))  
# oder: round(cov2cor(ability.cov$cov), 3)

# 'Gefittete' Korrelationsmatrix unter dem 1-Faktormodell 

( fa.1f.R <- 
round( (loadings(ability.fa.1f) %*% t(loadings(ability.fa.1f))) + 
diag(ability.fa.1f$uniquenesses), 3) )

# Residualmatrix (beachte hierbei, dass auch für unzentrierte und unskalierte 
# Ausgangsdaten x, die Korrelationsmatrix R_x, unter dem FA-Modell und den 
# FA-Annahmen, die Bedingung R_x = LL' + Psi erfüllt! Die Residualmatrix soll
# hier definiert sein als: 
#           R_x(Daten) - L(MLE-geschätzt)L(MLE-geschätzt)' + Psi(MLE-geschätzt)

(res.R.1f <- daten.R - fa.1f.R)

# ML FA mit 2 Faktor

(ability.fa.2f <- factanal(factors = 2, covmat = ability.cov))

# Korrelationsmatrix aus den Daten (siehe oben)

daten.R

# 'Gefittete' Korrelationsmatrix unter dem 2-Faktormodell 

( fa.2f.R <- 
round( (loadings(ability.fa.2f) %*% t(loadings(ability.fa.2f))) + 
diag(ability.fa.2f$uniquenesses), 3) )

# Residualmatrix

(res.R.2f <- daten.R - fa.2f.R)

# Was auffallend nicht gut fittet (d.h. auffallend grosses Residuum):
# die Korrelation zwischen picture und maze in den Daten ist 0.193, 
# unter dem 2-Faktormodell gefittet zu 0.308 (Residuum von, im Betrag, 0.115). 

# 3 und 4 Faktoren

(ability.fa.3f <- factanal(factors = 3, covmat = ability.cov))  

factanal(factors = 4, covmat = ability.cov)

# Für 3 Faktoren gilt für die Freiheitsgrade, df_3, der Chi-Quadrat Verteilung 
# im Bartlett-Test (Folie S. 155), df_3 = 0. Für 4 Faktoren gilt df_4 < 0. 
# In der Implementation in R wird im Fall df = 0 MLE-geschätzt, aber ohne einen 
# möglichen Bartlett-Test in diesem Fall; der Fall df < 0 führt zu einer
# Fehlermeldung. 

## Teilaufgabe (c)

# Faktorladungen nach Varimax-Rotation (Achtung: Im Plot Punkte per Maus anklicken, 
# um sie zu benennen!)

par(pty = "s")
plot(loadings(ability.fa.2f), xlim = c(0,1), ylim = c(0,1), xlab = "", ylab = "", pch = 19)
identify(loadings(ability.fa.2f)[, 1], loadings(ability.fa.2f)[, 2], 
         dimnames(loadings(ability.fa.2f))[[1]])  

# Oblimin-Rotation

ability.fa.2f.obl <- oblimin(loadings(ability.fa.2f))

# 'Rotierte' neue Faktormatrix \Lambda_rot = \Lambda_orig * T, wobei T 
# die Transformationsmatrix (siehe unten) ist.

loadings(ability.fa.2f.obl)  # oder: ability.fa.2f.obl$loadings

# Achsen der Oblimin-Rotation, eingezeichnet in der 2D-Ebene der Varimax-Achsen 
# (Achtung: Die Transformationsmatrix T ist hier gegeben durch die Inverse 
# der Transponierten von ability.fa.2f.obl$Th !) 

arrows(rep(0, 2), rep(0, 2), t(ability.fa.2f.obl$Th)[, 1], 
       t(ability.fa.2f.obl$Th)[, 2], col = "lightblue", length = 0.1, lwd = 2)

text(t(ability.fa.2f.obl$Th)[, 1] - 0.085, t(ability.fa.2f.obl$Th)[, 2] + 0.04, 
     c("Faktor1 (Oblimin)", "Faktor2 (Oblimin)"), col = "lightblue")
title(main = "Faktorladungen nach Varimax-Rotation, mit Achsen der Oblimin-Rotation als Pfeile",
      xlab = "Faktor1 (Varimax)", ylab = "Faktor2 (Varimax)")

## Teilaufgabe (d)

# Nein, man kann hier die Faktorwerte nicht berechnen, mittels keiner der beiden Methoden, 
# da nur die Kovarianzmatrix und nicht die Datenmatrix vorliegt, und letztere durch die 
# erstere nicht eindeutig bestimmt ist. In die jeweiligen Schätzgleichungen (Folie Seite 159) 
# gehen benötigte Grössen der Datenmatrix ein (x bzw. x_j), die nicht aus der Kovarianzmatrix 
# eindeutig bestimmt werden können. (Entsprechende Versuche in R mittels der Funktion 
# factanal() führen zu Fehlermeldungen.)