# Faktoren- und Reliabilitätsanalyse
#
# Erweitertes Beispiel aus Wirtz & Nachtigall (2006, S. 203 ff.)

# ==============================================================================

BaseURL="http://www2.jura.uni-hamburg.de/instkrim/kriminologie/Mitarbeiter/Enzmann/Lehre/StatIIKrim/Uebungen/"
Pfad = ""

# Gegebenfalls die folgende Anweisung entkommentieren und den Pfadnamen (mit ab-
# schließendem Schrägstrich-Vorwärts) eintragen, z.B.:

# Pfad = 'd:/lehre/statistik2/ss10/Uebungen/'

# ------------------------------------------------------------------------------
# Nur bei Internet-Anschluss verwenden:
download.file(paste(BaseURL,'desc.r',sep=""),paste(Pfad,'desc.r',sep=""),mode='wb')
download.file(paste(BaseURL,'meanscore.r',sep=""),paste(Pfad,'meanscore.r',sep=""),mode='wb')
download.file(paste(BaseURL,'recode.r',sep=""),paste(Pfad,'recode.r',sep=""),mode='wb')
download.file(paste(BaseURL,'pca.r',sep=""),paste(Pfad,'pca.r',sep=""),mode='wb')
download.file(paste(BaseURL,'fa.promax.r',sep=""),paste(Pfad,'fa.promax.r',sep=""),mode='wb')
download.file(paste(BaseURL,'reliability.r',sep=""),paste(Pfad,'reliability.r',sep=""),mode='wb')
download.file(paste(BaseURL,'factor1.dat',sep=""),paste(Pfad,'factor1.dat',sep=""),mode='wb')

# ------------------------------------------------------------------------------

# Die Dateien müssen im aktuellen Arbeitsverzeichnis vorliegen!
source(paste(Pfad,'desc.r',sep=""))
source(paste(Pfad,'meanscore.r',sep=""))
source(paste(Pfad,'recode.r',sep=""))
source(paste(Pfad,'pca.r',sep=""))
source(paste(Pfad,'fa.promax.r',sep=""))
source(paste(Pfad,'reliability.r',sep=""))
ErwStud = read.table(paste(Pfad,'factor1.dat',sep=""),header=T,sep='\t')

# ==============================================================================

# Deskriptive Statistiken:
desc(ErwStud)

# Korrelationsmatrix der Variablen:
round(cor(ErwStud),3)

# ------------------------------------------------------------------------------

# Hauptkomponentenanalyse (Kriterium für die Anzahl der Faktoren: Eigenwert >=
# 1, Output: sortierte Ladungen):
erg1 = pca(ErwStud,sort=T)

# Varimaxrotierte Ladungsmatrix:
erg1$varimax.loadings

# HINWEIS: Der folgende Teil der Syntax dient dazu, eine Grafik der unrotierten
# und Varimax rotierten Lösung zu erzeugen, die Syntax dazu ist komplex (und
# sicher nicht optimal) - sie muss nicht nachvollzogen werden!

windows(record=T)

dcos.inv = function(cosinus) acos(cosinus)*180/pi
degree.rot = function(rotmat) dcos.inv(rotmat[1,1])

rotarc = function(r=1,old=c(1,0),new=c(0,1),...)
{
   arc = function(r,start,end,...)
   {
      rs = seq(start,end,len=3600)
      xc = r*cos(rs)
      yc = r*sin(rs)
      lines(xc,yc,...)
   }
   if (old[2] >= 0) Half1 = 1
   if (old[2] <  0) Half1 = 2
   if (new[2] >= 0) Half2 = 1
   if (new[2] <  0) Half2 = 2


   if ((Half1 == 1) & (Half2 == 1))
   {
     start = acos(old[1]/r)
     end   = acos(new[1]/r)
     arc(r,start,end,...)
   }
   if ((Half1 == 1) & (Half2 == 2))
   {
     start = acos(old[1]/r)
     end   = asin(new[1]/r)+3/2*pi
     arc(r,start,end,...)
   }
   if ((Half1 == 2) & (Half2 == 2))
   {
     start = asin(old[1]/r)+3/2*pi
     end   = asin(new[1]/r)+3/2*pi
     arc(r,start,end,...)
   }
   if ((Half1 == 2) & (Half2 == 1))
   {
     start = asin(old[1]/r)+3/2*pi
     end   = acos(new[1]/r)+2*pi
     arc(r,start,end,...)
   }
}

A = erg1$unrotated[,1:2]
vm = varimax(A)
degree.rot(vm$rotmat)
Old = erg1$varimax.loadings %*% vm$rotmat
h = sqrt(rowSums(erg1$varimax.loadings**2))

plot(erg1$varimax.loadings,xlim=c(-1,1),ylim=c(-1,1),pch=19,col='blue',asp=1,
    main=paste('Varimax rotation (',round(degree.rot(vm$rotmat),1),'°)',sep=''))
abline(h=0)
abline(v=0)
lines(c(vm$rotmat[1,2],vm$rotmat[2,1]),c(-vm$rotmat[1,1],vm$rotmat[1,1]),
      lty=2,col='gray70')
lines(c(-vm$rotmat[1,1],vm$rotmat[1,1]),c(vm$rotmat[2,1],vm$rotmat[1,2]),
      lty=2,col='gray70')
points(Old,col='red')
for (i in 1:6)
{
  rotarc(r=h[i],old=Old[i,],new=erg1$varimax.loadings[i,],lty=2,col='blue')
}
text(erg1$varimax.loadings,labels=rownames(erg1$varimax.loading),pos=3,cex=0.7)
text(Old,labels=rownames(Old),col='red',pos=3,cex=0.7)

lines(c(Old["SOZSI05",1],Old["SOZSI05",1]),c(Old["SOZSI05",2],0),
      lty=2,col='red')
text(Old["SOZSI05",1],0,labels=round(Old["SOZSI05",1],2),pos=1,cex=0.7)
lines(c(Old["SOZSI05",1],0),c(Old["SOZSI05",2],Old["SOZSI05",2]),
      lty=2,col='red')
text(0,Old["SOZSI05",2],labels=round(Old["SOZSI05",2],2),pos=2,cex=0.7)

lines(c(erg1$varimax.loadings["SOZSI05",1],erg1$varimax.loadings["SOZSI05",1]),
      c(erg1$varimax.loadings["SOZSI05",2],0),col='blue')
text(erg1$varimax.loadings["SOZSI05",1],0,
     labels=round(erg1$varimax.loadings["SOZSI05",1],2),pos=1,cex=0.7)
lines(c(erg1$varimax.loadings["SOZSI05",1],0),
      c(erg1$varimax.loadings["SOZSI05",2],erg1$varimax.loadings["SOZSI05",2]),
      col='blue')
text(0,erg1$varimax.loadings["SOZSI05",2],
     labels=round(erg1$varimax.loadings["SOZSI05",2],2),pos=2,cex=0.7)

legend("topleft",legend=c('unrotated','varimax rotated'),
       pch=c(1,19),col=c('red','blue'),cex=0.8)

# Die Abbildung zeigt die Position der Variablen im Faktorraum vor und nach
# einer Varimax-Rotation (Items SozSi02 und FachV06 sind nicht rekodiert). Man
# beachte, dass die Konfiguration der unrotierten Variablen vor und nach Rota-
# tion gleich ist, die Rotation verändert zwar die Ladungsmuster, das rotierte
# Objekt verändert sich dabei jedoch nicht! Nach einer Rotation von 32.4° (gegen
# den Uhrzeigersinn) haben die Ladungen eine deutlich bessere Einfachstruktur.
# So beträgt z.B. die Ladung von SozSi05 in der ursprünglichen Ladungsmatrix .49
# (Faktor 1) und .48 (Faktor 2), während sie nach der Rotation .16 (Faktor 1)
# und .67 (Faktor 2) beträgt.

# ------------------------------------------------------------------------------

# Es fällt auf, dass die Variablengruppen FACH... (Vertrauen fachliche Kompe-
# tenzen) und SOZSI... (soziale Sicherheit) ihre höchsten Ladungen auf zwei
# unterschiedlichen Faktoren haben. Allerdings korreliert Item FACHV06 (Zweifel)
# negativ mit "seinem" Faktor, das gleiche gilt für Item "SOZSI02" (Lampenfie-
# ber). Deshalb sollten zur Berechnung von Scores für beide Dimensionen diese
# Items rekodiert werden (es handelt sich um eine 6-stufige Likertskala von
# 1-6):

ErwStud$SOZSI02 = recode(ErwStud$SOZSI02,c(1:6),c(6:1))
ErwStud$FACHV06 = recode(ErwStud$FACHV06,c(1:6),c(6:1))

# Zwar ist eine Rekodierung von Items für die Faktorenanalyse nicht nötig, sie
# kann aber die Interpretation der Ergebnisse deutlich erleichtern. Deshalb wird
# die Faktorenanalyse noch einmal mit den rekodierten Items wiederholt. Diesmal
# wird auch die Berechnung von Faktorwerten angefordert:

erg2 = pca(ErwStud,fscores=T,sort=T)

# Hier noch einmal die Varimax rotierte Ladungsmatrix, jetzt mit rekodierten
# Items:

erg2$varimax.loadings

# ------------------------------------------------------------------------------

# Zur Veranschaulichung der Bedeutung des Ladungsmusters die Korrelationsmatrix
# der Items, die dabei anhand des Ergebnisses der Faktorenanalyse (vgl. die
# Faktorladungen) sortiert werden:

round(cor(ErwStud[,c(6,4,1,3,2,5)]),3)

# Es zeigt sich, dass die Items innerhalb der Dimensionen deutlich höher als
# zwischen den Dimensionen miteinander korrelieren. Relativ gering sind jeweils
# die Korrelationen des Items FACHV01 und SOZSI05, die auch die jeweils niedrig-
# sten Ladungen auf "ihren" Dimensionen aufweisen.

# ------------------------------------------------------------------------------

# Berechnung von Scores (Itemmittelwerten) für die beiden Dimensionen:
SOZSI = meanscore(ErwStud[,c(2,3,5)])
FACHV = meanscore(ErwStud[,c(1,4,6)])

# Anhand der Varimax rotierten Faktorenlösung lassen sich für jede Person Fak-
# torwerte berechnen. Sie werden im Folgenden den z-standardisierten Item- und
# Skalenwerten (jeweils für beide Gruppen von Variablen getrennt) gegenüber-
# gestellt. Hinweis: Zum Vergleich werden hier Stichproben-z-Werte berechnet:

n = nrow(ErwStud)
k = sqrt(n/(n-1))

round(cbind(k*scale(ErwStud)[,c(6,4,1)],Score.FACHV=k*scale(FACHV)[,],
            FactorScore1=erg2$varimax.scores[,1]),2)
round(cbind(k*scale(ErwStud)[,c(3,2,5)],Score.SOZSI=k*scale(SOZSI)[,],
            FactorScore2=erg2$varimax.scores[,2]),2)

# ------------------------------------------------------------------------------

# Kurzer Exkurs für besonders Interessierte:
#
# Zur Veranschaulichung der Bedeutung der Ladungen und der Faktorwerte lassen
# sich anhand der Gleichung
#
# z(geschätzt) = Ladung(1)*Faktorwert(1) + Ladung(2)*Faktorwert(2)
#
# die geschätzten z-Werte für jedes Item für jede Person berechnen. Matrix-
# algebra ermöglicht dies durch Multiplikation von zwei Matrizen:
#
# Z = F * L'
#
# (L' bedeutet L transponiert, in R wird hierfür die Funktion t() benutzt; die
# Multiplikation von Matrizen wird in R mit dem Operator %*% durchgeführt):

z.geschätzt = erg2$varimax.scores %*% t(erg2$varimax.loadings)

# Am Beispiel des Items FACH04 werden die wahren und die geschätzten z-Werte
# miteinander verglichen:

n = nrow(ErwStud)
Z = sqrt(n/(n-1))*scale(ErwStud[,"FACHV04"])
round(cbind(FACHV04=Z[,],geschätzt=z.geschätzt[,"FACHV04"]),3)

# ------------------------------------------------------------------------------

# Als Kriterium für die Anzahl der Faktoren wurde ein minimaler Eigenwert von 1
# gewählt. Die Eigenwerte der Korrelationsmatrix zeigen, dass dies Kriterium in
# diesem speziellen Fall keine "grenzwertige" Entscheidung darstellt:

erg2$eigenvalues

# Dividiert man die Eigenwerte durch die Anzahl der Items der Korrelations-
# matrix, erhält man den Anteil der durch die jeweiligen Faktoren "aufgeklärte"
# oder repräsentierte Varianz der Items. Ziel ist es, möglichst wenig Faktoren
# mit maximaler Varianz zu extrahieren.

# Es wurden nur zwei Faktoren extrahiert. Die ursprüngliche Ladungsmatrix war:

erg2$unrotated.loadings

# Die Summe der quadrierten Ladungen pro Item (= Kommunalität oder Communality)
# gibt an, wieviel Varianz des Items durch die Faktoren "erklärt" wird.

# Damit die Faktoren besser interpretierbar sind, wurde eine Varimaxrotation
# durchgeführt, wodurch sich die Ladungen besser auf hohe Haupt- und niedrige
# Nebenladungen verteilen. Dabei ändert sich die Summe der quadrierten Ladungen
# (Kommunalität) nicht. Allerdings verteilen sich die Summen der quadrierten
# Ladungen pro Faktor anders, auch wenn die Gesamtvarianz (kumuliert) gleich
# bleibt:

erg2$varimax.SS

# ------------------------------------------------------------------------------

# Eine Varimaxrotation dreht nur das Koordinatenkreuz durch die Ladungen, d.h.
# die Faktoren bleiben dabei senkrecht aufeinander stehen bzw. sind nicht
# korreliert. Die Annahme unkorrelierter Faktoren ist häufig unangemessen,
# weshalb eine schiefwinklige Rotation vorzuziehen ist. Von den unterschiedli-
# chen schiefwinkligen Rotationsmethoden, die korrelierte Faktoren zulassen,
# sind Oblimin und Promax im Allgemeinen am brauchbarsten. Eine schiefwinklige
# Promaxrotation der Items ergibt folgende Pattern- bzw. Ladungsmatrix mit
# korrelierten Faktoren:

erg2$promax.loadings
erg2$corr.factors

# Die Faktoren sind nur schwach korreliert, deshalb sind die "Hauptladungen",
# der Items, die die jeweiligen Dimensionen definieren, nur geringfügig höher
# und die Nebenladungen auf den "fremden" Dimension nur geringfügig niedriger.

# Anhand dieser Lösung lassen sich ebenfalls Faktorwerte berechnen:

erg2$promax.scores

# Die Faktorwerte korrelieren genauso wie die Faktoren:

round(cor(erg2$promax.scores),4)

# Varimaxrotierte Faktorenwerte sind demgegenüber unkorreliert:

round(cor(erg2$varimax.scores),4)

# ------------------------------------------------------------------------------

# Besitzen die via Itemsummen oder -mittelwerten gebildeten Skalenwerte eine
# ähnliche Bedeutung wie die Faktoren, sollte die Korrelation der geschätzten
# Faktorwerte mit den Skalenwerten der jeweiligen Faktoren hoch sein:

round(cor(data.frame(FACHV,SOZSI,erg2$varimax.scores)),4)

# Tatsächlich ist die "Faktorkongruenz" (d.h. die Korrelation zwischen Skalen-
# werten und Faktorwerten) sehr gut. Auch ist die Annahme unkorrelierter
# Faktoren anhand der geringen Korrelation der Skalenwerte untereinander
# begründbar. Das Ergebnis der schiefwinklig rotierten Lösung zeigt, dass das
# Ergebnis in diesem Beispiel sehr ähnlich ausfällt:

round(cor(data.frame(FACHV,SOZSI,erg2$promax.scores)),4)

# ------------------------------------------------------------------------------

# Bisher wurde die Anzahl der Faktoren anhand des Kaiser-Guttman Kriteriums
# (Eigenwert >= 1) bestimmt. Dies häufig benutzte Kriterium hat viele Nachteile
# (häufig werden damit zu viele Faktoren extrahiert). Neben dem Scree-Test
# existiert ein Kriterium, dass sich auf den Vergleich der empirischen Eigen-
# werte mit Eigenwerten stützt, die sich aus Korrelationsmatrizen von Zufalls-
# zahlen ergeben. Es werden dabei genauso viele Fälle und Variablen pro Korrela-
# tionsmatrix zufällig erzeugt, wie im aktuellen Datensatz vorhanden sind. Da
# die Korrelationen nur zufällig von Null abweichen, sollte es keinerlei bedeut-
# same Faktoren zu extrahieren geben. Legt man diese Zufallseigenwerte parallel
# neben die empirischen Eigenwerte (Parallel-Analyse), gelten bei diesem Verfah-
# ren nur so viele Faktoren zu extrahieren, wie es empirische Eigenwerte gibt,
# die größer als die Zufallseigenwerte sind.

# Im Folgenden wird eine Parallelanalyse zur Bestimmung der optimalen Faktoren-
# zahl spezifiziert. Per Voreinstellungen werden Korrelationsmatrizen von 1000
# Stichproben zur Berechnung der Zufallseigenwerte benutzt; mehr Stichproben
# liefern präzisere Ergebnisse, sind aber sehr rechenaufwendig:

erg3=pca(ErwStud,parallel=T)
erg3$eigenvalues

# Auch nach dem Paralleltest-Kriterium werden in diesem Fall zwei Faktoren
# extrahiert, d.h. es ergibt sich die gleiche Lösung wie zuvor.

# Manchmal ist es sinnvoll, den Verlauf der Eigenwerte graphisch darzustellen,
# um besser entscheiden zu können, wie viele Faktoren extrahiert werden sollten.
#
# Hinweis: Da pca() mit dem Argument parallel=T aufgerufen wurden, stehen die
# empirischen Eigenwerte hier nicht in der 1. Spalte von erg3$eigenvalues,
# sondern in der 2. Ohne eine Parallelanalyse stünden sie in der 1. Spalte und
# beim Aufruf von plot() (siehe 1. plot()-Anweisung) müsste ebenfalls Spalte 1
# (via erg$eigenvalues[,1]) spezifiziert werden.

# Beachte, dass über den Menüpunkt "History" des Grafikfensters auch die vor-
# herige Grafik angezeigt werden kann.

plot(erg3$eigenvalues[,2],type='b',col='blue',
     ylim = range(0,erg3$eigenvalues[,2]),ylab='Eigenwert',
     main = 'Eigenwertverlauf')
abline(h=1,lty=2)
par(new=T)
plot(erg3$eigenvalues[,1],type='b',lty=2,col='red',
     ylim = range(0,erg3$eigenvalues[,2]),ylab='Eigenwert')
legend('topright',legend=c('empirische Eigenwerte','Zufallseigenwerte'),
       col=c('blue','red'),lty=c(1,2))

# Die Grafik der Eigenwerte (Scree-Test) zeigt, dass die Eigenwerte ab dem 3.
# Eigenwert tendenziell in eine Linie übergehen bzw. dass sie im "Geröll" des
# Abhangs untergehen. Die Zufallseigenwerte (rot) schneiden die empirischen
# Eigenwerte zwischen dem 2. und 3. Eigenwert, so dass sowohl nach dem Kaiser-
# Guttman Kriterum (minimaler Eigenwert >= 1; waagerechte Linie), dem Scree-
# Test, als auch nach dem Kriterium der Parallelanalyse nur zwei Faktoren zu
# extrahieren sind.

# ------------------------------------------------------------------------------

# Man kann auch aus theoretischen Gründen (a priori) eine bestimmte Faktorenzahl
# erwarten. Auch wenn sich letztlich eine andere Anzahl von Faktoren als besser
# interpretierbar erweist, sollte man dennoch theoretisch erwartete Faktor-
# modelle prüfen, da die aktuellen Ergebnisse immer auch von Besonderheiten
# der jeweiligen Stichprobe beeinflusst sein können. Beispiel:

pca(ErwStud,factors=3,sort=T)

# In dieser Lösung stellt das Item SOZSI05 (zu großer Respekt vor Dozenten)
# einen eigenen Faktor dar, der allerdings hoch mit Faktor 2 (SOZSI) korreliert
# ist.

# ==============================================================================

# Bisher wurde als "Faktorenanalyse" eine Hauptkomponentenanalyse (PCA =
# principal components analysis) durchgeführt, bei der die gesamte Varianz der
# Variablen in weniger Komponenten zerlegt wird. Es gibt gute Gründe, andere
# Verfahren zu benutzen, bei denen nur die den Variablen gemeinsame Varianz in
# Faktoren zerlegt wird. Die Haupkomponentenanalyse dient in erster Linie der
# Datenreduktion (von Variablen mit wenig Messfehlern). Besteht das Ziel darin,
# latente Dimensionen oder Variablen zu identifizieren, sind andere Verfahren
# wie die Hauptachsen-Faktorenanalyse (PAF = principal axes factoring) oder eine
# Maximum-Likelihood-Faktorenanalyse (ML) angemessener. Im Allgemeinen sind die
# Ladungen bei derartigen Verfahren niedriger, da nur die den Variablen gemein-
# same Varianz analysiert wird (der Rest setzt sich aus itemspezifischer Varianz
# und Messfehlern zusammen).

# Maximum-Likelihood Faktorenanalysen haben den Vorteil, dass sich die Anzahl
# der zu extrahieren Faktoren (bzw. die Passung des jeweiligen Faktormodells
# zu den Daten) empirisch prüfen lässt. Darüber hinaus sind sie enger verwandt
# mit so genannten konfirmatorischen Faktorenanalysen, mit denen ein theoretisch
# erwartetes Mess- bzw. Faktormodell geprüft werden kann.

# Beispiel: Bestimmung der optimalen Faktorenzahl mit Hilfe einer Maximum-
# Likelihood Faktorenanalyse. Im ersten Schritt wird eine einfaktorielle Lösung
# (Voreinstellung) spezifiziert:

fa.promax(ErwStud)

# Hinweis: Unter "Uniqueness" werden die Residualvarianzen, d.h. die durch die
# nicht extrahierten Faktoren nicht erklärten Anteile der Varianzen der Varia-
# blen dargestellt. Dies entspricht 1-Kommunalität bei einer Hauptkomponenten-
# analyse.

# Das Ergebnis zeigt, dass die Residualmatrix der Korrelationen signifikant von
# Null verschieden ist (Chi²(9) = 20.8, p = .014: der Chi-Quadrat-Test prüft,
# wie gut die Daten zum Modell passen; ist er signifikant, ist die Hypothese,
# dass das Modell mit den Daten übereinstimmt, zu verwerfen), d.h. das Faktor-
# modell erklärt die Korrelationen der Variablen nicht hinreichend. Ein zwei-
# faktorielles Modell passt demgegenüber sehr gut zu den Daten (Chi²(4) = 0.1,
# p = .998):

fa.promax(ErwStud, factors=2, sort=T)

# Allerdings zeigt sich auch, dass die Lösung keine gute "Einfachstruktur" hat:
# Die Items der Dimension FACHV haben teilweise relativ hohe Ladungen auf dem
# Faktor SOZSI, was nicht durch eine Korrelation der Faktoren erklärbar ist
# (die Promax rotierte Lösung zeigt, dass die Faktoren nur gering miteinander
# korrelieren). Eine gute Faktorenlösung zeichnet sich immer dadurch aus, dass
# die Hauptladungen sehr hoch (idealerweise |1|) und die Nebenladungen sehr
# niedrig (idealerweise 0) sind.

# ==============================================================================

# Schließlich kann die Reliabilität der Skalen untersucht werden, indem Cron-
# bachs Alpha für die Items der jeweiligen (Sub)Skalen berechnet wird:

reliability(ErwStud)

# Cronbachs Alpha der Gesamtskala ist mit .67 unbefriedigend (es sollte minde-
# stens .70, besser .80 und mehr) betragen. Die Analyse zeigt, dass die Skala
# ohne das Item FACHV01 homogener wäre. Auf diese Weise würde aber ein Item
# entfernt, dass in der Faktorenanalyse eine ziemlich hohe Ladung aufwies.
#
# Werden statt dessen die Items der beiden Faktoren getrennt analysisiert, ist
# das Ergebnis für beide Skalen wesentliche besser, obwohl eine größere Anzahl
# von Items in der Regel die Reliabilität einer Skala erhöht. Offenbar sind nur
# die Subskalen in sich homogen:

reliability(ErwStud[,c(6,4,1)])
reliability(ErwStud[,c(3,2,5)])

# Die Subskalen haben gute bzw. befriedigende Reliabilitätskoeffizienten: Alpha
# von FACHV = .84, Alpha von SOZSI = .75. Zwar würde in beiden Fällen je ein
# Item entfernt werden können, was das Alpha noch weiter erhöhen würde. Die
# Reliabilität der Subskalen ist jedoch ausreichend - Items zu entfernen
# bedeutet immer auch ein Veränderung dessen, was die Skala misst. FACHV01 ist
# sicher konzeptuell ein wesentlicher Aspekt von Fachvertrauen, weshalb es in
# der Skala belassen werden sollte. Diskutabel wäre eher Item SOZSI05, weil
# damit nicht eindeutig eine Aussage über die befragte Person selbst sondern
# auch über die übrigen Studenten gemacht wird. Auch die dreifaktorielle
# Faktorenanalyse hat gezeigt, dass dieses Item am schlechtesten zur Skala
# passt.

# Ein häufiger Fehler bei Reliabilitätsanalysen besteht darin, dass nicht zuvor
# alle Items in die gleiche Richtung gepolt wurden. Wäre z.B. das Item SOZSI02
# nicht rekodiert worden, hätte sich folgendes Ergebnis für die Gesamtskala
# ergeben:

ErwStud$SOZSI02 = recode(ErwStud$SOZSI02,c(1:6),c(6:1))
reliability(ErwStud)

# Demnach hätten scheinbar weder SOZSI02 noch SOZSI03 zur Skala gepasst. Diese
# Items daraufhin zu löschen, wäre angesichts der obigen Ergebnisse eine klare
# Fehlentscheidung:

reliability(ErwStud[,-c(2,3)])

# Negative korrigierte Item-Total-Korrlationen wie im obigen Beispiel sind fast
# immer ein Hinweis darauf, dass bestimmte Items hätten rekodiert werden sollen.
# Vor Reliabilitätsanalysen Faktorenanalysen durchzuführen und sorgfältig zu
# interpretieren, kann vor derartigen Fehlentscheidungen bewahren.