30  Deskriptive Statistik

Zur Beschreibung von Verteilungen eigenen sich Tabellen, Lage- und Streuungskenngrößen.

30.1 Häufigkeiten

Eine Häufigkeitstabelle wird mit der Funktion table() erstellt

# Daten erzeugen
x <- c(1, 3, 5, 7, 7, 7, 6, 34)

# Häufigkeitstabelle
table(x)
## x
##  1  3  5  6  7 34 
##  1  1  1  1  3  1

Für relative Häufigkeiten teilt man einfach durch die Länge der Variable x. So ändert sich der Befehl in:

# relative Häufigkeitstabelle
table(x)/length(x)
## x
##     1     3     5     6     7    34 
## 0.125 0.125 0.125 0.125 0.375 0.125

Kumulierte Häufigkeiten erhält man mit der Funktion cumsum()

# kumulierte Häufigkeitstabelle
cumsum(table(x))
##  1  3  5  6  7 34 
##  1  2  3  4  7  8

Und entsprechend auch kumulierte relative Häufigkeiten

# kumulierte relative Häufigkeitstabelle
cumsum(table(x)/length(x))
##     1     3     5     6     7    34 
## 0.125 0.250 0.375 0.500 0.875 1.000

Eine weitere Option bietet die Funktion freq() aus dem rettyR-Paket. Sie erstellt eine Tabelle mit absoluten und relativen Häufigkeiten sowie relativen Häufikeiten ohne fehlende Wert (NA).

# Paket "prettyR" installieren
install.packages("prettyR", dependencies=TRUE)
# Paket laden
library(prettyR)

# Häufigkeitstabelle mit absoluten und relative Werten erstellen
freq(x)
## 
## Frequencies for x 
##         7    1    3    5    6   34   NA
##         3    1    1    1    1    1    0
## %    37.5 12.5 12.5 12.5 12.5 12.5    0 
## %!NA 37.5 12.5 12.5 12.5 12.5 12.5

30.2 Lagekenngrößen

Mit dem Befehl summary() erhält man eine erste Zusammenfassung der Werteverteilung mit den wichtigsten Lagekenngrößen.

# Daten erzeugen
x <- c(1, 3, 5, 7, 7, 7, 6, 34)

# Zusammenfassung
summary(x)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    1.00    4.50    6.50    8.75    7.00   34.00

Die Funktion describe() aus dem psych-Zusatzpaket liefert noch detailiertere Übersichten:

psych::describe(x)
##    vars n mean    sd median trimmed  mad min max range skew kurtosis   se
## X1    1 8 8.75 10.43    6.5    8.75 1.48   1  34    33 1.69     1.35 3.69

Das Funktioniert auch bei kompletten Datensätzen.

# Zusammenfassung von pf8
psych::describe(pf8)
##                vars   n   mean    sd median trimmed   mad min max range  skew
## Standort*         1 731   2.85  1.61      2    2.81  1.48   1   5     4  0.28
## Alter             2 730  37.00 18.00     28   34.98 12.60  15  86    71  0.74
## Geschlecht*       3 725   1.61  0.49      2    1.63  0.00   1   3     2 -0.40
## Größe             4 730 174.04  9.62    173  173.84 10.38 151 206    55  0.26
## Gewicht           5 720  75.53 16.88     73   73.93 16.31  45 175   130  1.12
## Bildung*          6 506   4.88  1.77      5    4.98  2.97   1   7     6 -0.39
## Beruf*            7 731  40.42 37.04     34   38.58 48.93   1 104   103  0.26
## Familienstand*    8 729   2.11  1.12      2    1.99  1.48   1   6     5  0.83
## Kinder            9 727   0.38  0.79      0    0.18  0.00   0   5     5  2.19
## Wohnort*         10 729   1.41  0.49      1    1.39  0.00   1   2     1  0.36
## Rauchen*         11 729   1.21  0.40      1    1.13  0.00   1   2     1  1.45
## SportHäufig      12 677   2.56  1.88      2    2.41  1.48   0  21    21  1.91
## SportMinuten     13 667  62.42 40.14     60   60.47 44.48   0 400   400  1.53
## SportWie*        14 656   1.64  0.69      2    1.55  1.48   1   3     2  0.62
## SportWarum*      15 454   3.75  1.64      4    3.49  1.48   1   8     7  1.06
## LebenZufrieden   16 724   7.91  1.59      8    8.07  1.48   1  11    10 -1.10
##                kurtosis   se
## Standort*         -1.53 0.06
## Alter             -0.70 0.67
## Geschlecht*       -1.75 0.02
## Größe             -0.26 0.36
## Gewicht            2.40 0.63
## Bildung*          -1.26 0.08
## Beruf*            -1.55 1.37
## Familienstand*     0.48 0.04
## Kinder             4.45 0.03
## Wohnort*          -1.87 0.02
## Rauchen*           0.11 0.01
## SportHäufig       12.75 0.07
## SportMinuten       8.64 1.55
## SportWie*         -0.77 0.03
## SportWarum*        0.79 0.08
## LebenZufrieden     1.86 0.06

Die Lagewerte lassen sich auch einzeln bestimmen.

Mit der Funktion min() wird das Minimum der Verteilung ausgegeben.

# Minimum bestimmen
min(x)
## [1] 1

Dementsprechend liefert max() das Maximum.

# Maximum bestimmen
max(x)
## [1] 34

Die Spannweite wird mit range() abgefragt.

# Spannweite bestimmen
range(x)
## [1]  1 34

Das arithmetische Mittel wird mit der Funktion mean() bestimmt

# x-quer bestimmen
mean(x)
## [1] 8.75

Sind fehlende Werte in der Reihe enthalten, gibt R nur NA zurück.

#  mean für "Gewicht" im Datensatz "pf8"
mean(pf8$Gewicht)
## [1] NA

Beheben kann man dies, indem der Parameter na.rm auf TRUE gesetzt wird.

#  mean für "Gewicht" im Datensatz "pf8"
# entferne NAs
mean(pf8$Gewicht, na.rm=TRUE)
## [1] 75.5325

Der Median wird mit der Funktion median() ermittelt.

# Median bestimmen
median(x)
## [1] 6.5

Auch hier müssen enthaltene NAs ausgeblendet werden:

# Median bestimmen, OHNE NAs
median(pf8$Gewicht, na.rm=TRUE)
## [1] 73

In den R-Hausmitteln existiert keine eigene Funktion zur Bestimmung des Modalwertes. Wir können ihn daher “auf Umwegen” bestimmen. Wir lassen mittels table() eine Häufigkeitstabelle der Variable ausgeben. Diese sortieren wir mittels sort() absteigend. Der Modus ist dann der erste Wert der sortierten Tabelle.

# Daten als Häufigkeitstabelle
table(x)
## x
##  1  3  5  6  7 34 
##  1  1  1  1  3  1
# absteigend sortieren
sort(table(x), decreasing=TRUE)
## x
##  7  1  3  5  6 34 
##  3  1  1  1  1  1
# Modalwert
sort(table(x), decreasing=TRUE)[1]
## 7 
## 3

Es ist der Wert 7 und er kommt 3 mal vor.

Eine weitere Möglichkeit bietet das Zusatzpaket statip. Es enthält die Funktion mfv(), mit welcher der Modalwert bestimmt werden kann.

# Daten als Häufigkeitstabelle
library(statip)

# Modalwert bestimmen
mfv(x)
## [1] 7

Die vier Quartile werden mit der Funktion quantile() ermittelt (ja, der Befehl für Quartile heisst in R quantile()).

# Quartile bestimmen, ohne NAs
quantile(x, na.rm=T)
##   0%  25%  50%  75% 100% 
##  1.0  4.5  6.5  7.0 34.0

Es lassen sich beliebige Quantile berechnen.

# Quartile bestimmen, ohne NAs
quantile(x, probs = c(.333, .666))
## 33.3% 66.6% 
## 5.331 7.000

Und auch die Art der Berechnung lässt sich über den Parameter type festlegen. Beispielsweise verwendet type=6 die Methode von SPSS.

# Quartile bestimmen, ohne NAs
# so wie SPSS
quantile(x, probs = c(.333, .666), type=6)
## 33.3% 66.6% 
## 4.994 7.000

30.3 Streuungskenngrößen

Die wichtigste Streuungskenngröße ist die Standardabweichung. Sie wird bestimmt mittels der Funktion sd().

# Standardabweichung bestimmen
sd(x, na.rm=T)
## [1] 10.43004

Die Varianz der Stichprobe wird mit var() berechnet.

# Varianz von x bestimmen
var(x, na.rm=T)
## [1] 108.7857

Der Interquartilsabstand (oder einfach nur Quartilsabstand) wird mit der Funktion IQR() bestimmt.

# Quartilsabstand bestimmen
IQR(x, na.rm=T)
## [1] 2.5

Ebenso wie bei quantile() lässt sich über den Parameter type die Berechnungsart ändern. Mit type=6 nutzt R die Methode wie in SPSS:

# Quartilsabstand bestimmen
# Methode wie in SPSS
IQR(x, na.rm=T, type=6)
## [1] 3.5

30.4 Kreuztabellen

Mit der Funktion table() lassen sich auch Kreuztabellen erstellen, indem einfach beide Variablen übergeben werden.

# erzeuge Daten
Punktwert <- c(2 , 12 , 7 , 15 , 10 , 4 , 13 , 9 , 16 , 19)
Geschlecht <- rep(c("maennlich" , "weiblich") ,5)

# Kreuztabelle
table(Punktwert, Geschlecht)
##          Geschlecht
## Punktwert maennlich weiblich
##        2          1        0
##        4          0        1
##        7          1        0
##        9          0        1
##        10         1        0
##        12         0        1
##        13         1        0
##        15         0        1
##        16         1        0
##        19         0        1

Ebenso steht die Funktion xtabs() zur Verfügung.

# Kreuztabelle mit "xtabs()"
xtabs(~ Punktwert + Geschlecht)
##          Geschlecht
## Punktwert maennlich weiblich
##        2          1        0
##        4          0        1
##        7          1        0
##        9          0        1
##        10         1        0
##        12         0        1
##        13         1        0
##        15         0        1
##        16         1        0
##        19         0        1

Siehe Abschnitt 34.1 für weitere Beispiele zur Kreuztabelle.

30.5 z-Transformation

Die z-Transformation erfolgt mit der Funktion scale().

# Wertereihe
Punktwerte <- c(2 , 12 , 7 , 15 , 10 , 4 , 13 , 9 , 16 , 19)

# z-Transformation
scale(Punktwerte)
##             [,1]
##  [1,] -1.6183421
##  [2,]  0.2418212
##  [3,] -0.6882604
##  [4,]  0.7998702
##  [5,] -0.1302114
##  [6,] -1.2463094
##  [7,]  0.4278376
##  [8,] -0.3162278
##  [9,]  0.9858866
## [10,]  1.5439356
## attr(,"scaled:center")
## [1] 10.7
## attr(,"scaled:scale")
## [1] 5.375872

30.6 Korrelation

Die allgemeine Funktion zur Berechnung von Korrelationen heisst cor(). Die Korrelation nach Pearson berechnet sich mit der Funktion cor() und deren Parameter method="pearson".

# erzeuge Daten
Lesetest <- c(2 , 12 , 7 , 15 , 10 , 4 , 13 , 9 , 16 , 19)
Rechtschreibung <- c(3 , 14 , 9 , 17 , 12 , 4 , 16 , 12 , 18 , 20)

# Pearson Maßkorrelationskoeffizient
cor(Lesetest, Rechtschreibung, method="pearson")
## [1] 0.9884616

Mit der Funktion cor.test() kann direkt ein Siginifkanztest mitgefahren werden

# Pearson Maßkorrelationskoeffizient
# mit Signifikanztest
cor.test(Lesetest, Rechtschreibung, method="pearson")
## 
##  Pearson's product-moment correlation
## 
## data:  Lesetest and Rechtschreibung
## t = 18.458, df = 8, p-value = 7.648e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9502085 0.9973659
## sample estimates:
##       cor 
## 0.9884616

Sollten fehlende Wert (NA) enthalten sein, gibt der Befehl nur ein NA zurück.

# Pearson Maßkorrelationskoeffizient
# mit Datensatz "pf8", da sind NAs enthalten!!!
with(pf8, cor(Alter,Größe, method="pearson"))
## [1] NA

Mit dem Parameter use="complete.obs" (nutze nur vollständige Beobachtungen) können wir dies beheben.

# Pearson Maßkorrelationskoeffizient
# mit Datensatz "pf8", verwende nur "komplette" Daten (ohne NA)
with(pf8, cor(Alter,Größe, method="pearson", use="complete.obs"))
## [1] -0.1067484

Für die Korrelationsberechnung nach Spearmen (Spearman’s \(\rho\)) wird der Parameter method auf spearman gesetzt.

# Spearman's Korrelationsberechnung
cor(Lesetest, Rechtschreibung, method="spearman")
## [1] 0.9969651

Auch hier müssen NAs mit dem Parameter use="complete.obs" unterdrückt werden.

# Spearman's Rho
# mit Datensatz "pf8", verwende nur "komplette" Daten (ohne NA)
with(pf8, cor(Alter,Größe, method="spearman", use="complete.obs"))
## [1] -0.05920716

Auch für Spearman’s Rho kann ein Signifikanztest mitgefahren werden, indem cor.test() aufgerufen wird.

# Spearman's Korrelationsberechnung
# mit Signifikanztest
cor.test(Lesetest, Rechtschreibung, method="spearman")
## Warning in cor.test.default(Lesetest, Rechtschreibung, method = "spearman"):
## Kann exakten p-Wert bei Bindungen nicht berechnen
## 
##  Spearman's rank correlation rho
## 
## data:  Lesetest and Rechtschreibung
## S = 0.50076, p-value = 3.698e-10
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## 0.9969651

Die Warnmeldung kann ignoriert werden, der p-Wert wird in der Ausgabe angezeigt.

Da cor.test() ein Objekt der Klasse Liste zurückliefert, kann auf die einzelnen Komponenten des Tests zugegriffen werden.

# schreibe Testergebnis in Variable "test"
test <- cor.test(Lesetest, Rechtschreibung, method="spearman")

# zeige nur p-Wert an
test$p.value
## [1] 3.698093e-10
# zeige nur Teststatistik an
test$statistic
##         S 
## 0.5007599
# zeige nur Spearmans's Rho an
test$estimate
##       rho 
## 0.9969651

30.7 Regressionen

Die Funktion zur Berechnung linearer Regressionen ist lm() (für “linear model”). Sie folgt der Syntax x erklärt durch y, was durch eine Tilde ~ ausgedrückt wird.

# erzeuge Daten
Lesetest <- c(2 , 12 , 7 , 15 , 10 , 4 , 13 , 9 , 16 , 19)
Rechtschreibung <- c(3 , 14 , 9 , 17 , 12 , 4 , 16 , 12 , 18 , 20)

# lineare Regression "Rechtschreibung erklärt durch Lesetest"
lm(Rechtschreibung~Lesetest)
## 
## Call:
## lm(formula = Rechtschreibung ~ Lesetest)
## 
## Coefficients:
## (Intercept)     Lesetest  
##       1.208        1.055

Hierbei ist Intercept der Schnitt durch die Y-Achse.

Wenn Sie über den Parameter data das Datenframe angeben, können Sie auf die Variablen referenzieren, ohne $ wenden zu müssen.

# lade Testdaten
load(url("https://www.produnis.de/R/nw.RData"))

# entweder Variablen per "$" referenzieren
lm(nw$workyear ~ nw$worknight)
## 
## Call:
## lm(formula = nw$workyear ~ nw$worknight)
## 
## Coefficients:
##  (Intercept)  nw$worknight  
##       9.9143        0.8291
# oder Datenframe mit angeben
lm(workyear ~ worknight, data=nw)
## 
## Call:
## lm(formula = workyear ~ worknight, data = nw)
## 
## Coefficients:
## (Intercept)    worknight  
##      9.9143       0.8291

Es ist üblich, Modelle in ein Objekt zu speichern und über die summary()- Funktion aufzurufen. Die Informationen der Modellzusammenfassung sind so detailierter.

# speichere Modell in "fit"
fit <- lm(workyear~worknight, data=nw)
# schaue mit "summary()" an
summary(fit)
## 
## Call:
## lm(formula = workyear ~ worknight, data = nw)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -12.889  -5.743  -2.351   3.660  26.427 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.91425    0.74088   13.38   <2e-16 ***
## worknight    0.82912    0.05958   13.91   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.463 on 274 degrees of freedom
## Multiple R-squared:  0.4141, Adjusted R-squared:  0.4119 
## F-statistic: 193.6 on 1 and 274 DF,  p-value: < 2.2e-16

Außerdem lässt sich das Modell so leicht plotten.

plot(fit)

Dies erzeugt insgesamt 4 Plots (Residualplot, Q-Q, Scale-Location und Einfluss). Durch drücken der [Enter]-Taste schalten Sie zum nächsten Plot weiter.

Sie können aber auch alle 4 Plots in einem ausgeben lassen.

# plotte alle 4 auf einer Seite
par(mfrow = c(2, 2), oma = c(0, 0, 2, 0))
plot(fit)

30.7.1 multiple lineare Regressionen

Multiple lineare Regressionen werden ebenfalls mit der Funktion lm() berechnet. Hierbei werden die erklärenden Variablen per + Zeichen aneinandergereiht.

# multiple lineare Regression 
# "workyear" erklärt durch "worknight" UND "age"
fit <- lm(workyear ~ worknight+age, data=nw)
summary(fit)#
## 
## Call:
## lm(formula = workyear ~ worknight + age, data = nw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -14.7077  -3.9748  -0.3066   3.5936  19.9569 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.99856    1.72557  -3.476 0.000591 ***
## worknight    0.57707    0.05714  10.100  < 2e-16 ***
## age          0.41614    0.04195   9.921  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.41 on 273 degrees of freedom
## Multiple R-squared:  0.5693, Adjusted R-squared:  0.5662 
## F-statistic: 180.5 on 2 and 273 DF,  p-value: < 2.2e-16

Ebenfalls können hier die 4 Plots erzeugt werden.

# plotte alle 4 auf einer Seite
par(mfrow = c(2, 2), oma = c(0, 0, 2, 0))
plot(fit)

30.7.2 Logistische Regression

Als Beispiel dienen Daten der Challenger-Katastrophe (am 28. Januar 1986 explodierte die Raumfähre Challenger 73 Sekunden nach dem Start).

Der Ausfall eines oder mehrerer Dichtungsringe in einer der seitlichen Feststoffrakten wurde als Grund der Katastrophe ermittelt. Probleme mit den Dichtungsringen sind auch vor- her am Boden aufgetreten. So liegen Erhebungsdaten vor, die den Defekt ei- nes Rings in Zusammenhang mit der Außentemperatur darstellen können: Je kälter es draußen ist, desto höher die Ausfallwahrscheinlichkeit eines Rings. Bedauerlicher Weise wurden diese Daten erst rückwirkend ausgewertet. Diese Daten übertragen wir wie folgt in R:

Temp <- c(66, 67, 68, 70, 72, 75, 76, 79, 53, 58, 70, 75, 67, 67, 69, 70, 73, 76, 78, 81, 57, 63, 70)
Ausfall <- factor(c(0, 0,0 ,0 ,0 ,0 ,0 ,0 ,1 ,1 ,1 ,1 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,0 ,1 ,1, 1))
shuttle <- data.frame(Temp, Ausfall)
colnames(shuttle) <- c("Temp", "Ausfall")

Die Kodierung des factor Ausfall lautet: 0=kein Ausfall | 1=Ausfall Die Einheit für Temp lautet Fahrenheit.

Der Aufruf zur logistischen Regression (Ausfall erklärt durch Temp) lautet:

fit <- glm(Ausfall~Temp,data=shuttle, family=binomial)
summary(fit)
## 
## Call:
## glm(formula = Ausfall ~ Temp, family = binomial, data = shuttle)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.0611  -0.7613  -0.3783   0.4524   2.2175  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  15.0429     7.3786   2.039   0.0415 *
## Temp         -0.2322     0.1082  -2.145   0.0320 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 28.267  on 22  degrees of freedom
## Residual deviance: 20.315  on 21  degrees of freedom
## AIC: 24.315
## 
## Number of Fisher Scoring iterations: 5

Die Odds Ratio (OR) für \(\beta\)Temp= -0,2322 errechnet sich mit der exp()-Funktion:

exp(-0.2322)
## [1] 0.7927875

Dies liefert das Ergebnis OR = 0.7927875. Die Interpretation lautet: Mit jedem höhren Grad an Temparatur ändert sich die Chance auf einem Ausfall um den Faktor 0.7927. Das bedeutet, dass die Chancen eines Ausfalls niedriger werden, je höher die Außentemperatur steigt.

In der Nacht vor dem Challengerstart konnte eine Temperatur von 37° Fahrenheit gemessen werden. Auf Grundlage der logistischen Regression lässt sich nun die Ausfallwahrscheinlichkeit bei 37° Außentemperatur berechnen:

predict(fit, data.frame(Temp=37),type="response")
##         1 
## 0.9984265

Das Modell sagt bei einer Außentemperatur von 37° eine Ausfallwahrscheinlichkeit von 99,84% voraus.

Jetzt kann man dies für jede mögliche Temperatur durchmachen, und schon erhält man dieses schöne Plot:

y <- predict(fit, data.frame(Temp=30:90), type="response")
plot(30:90, y, type="l", ylab="Ausfallswahrscheinlichkeit", xlab="Temperatur in F")
abline(v=seq(20,100,1), h=seq(-0.2,1.2,0.05),col="gray") #Mathegitter
abline(v=37, col="red") # 37-Grad
abline(h=0.5, col="black") # Mittellinie
points(30:90, y, type="l",lwd=2, col="blue") # Schätzwerte in "blau"

30.7.3 Ordinale Regression

Zur Bestimmung des Einflusses von Prädiktoren auf ein ordinales Ziel stehen verschiedene Modelle der ordinalen Regression zur Verfügung. Namentlich werden hier vor allem Proportional Odds Modelle und Continuation Ratio Modelle verwendet.

Eine verständliche Einführung in ordinale Modelle findet sich im Open Access Artikel von große Schlarmann und Galatsch (große Schlarmann and Galatsch 2014).

Für die ordinalen Regressionen verwenden wir das “VGAM”, welches zunächst in R installiert werden muss

install.packages("VGAM", dependencies=T)

Als Beispiel dient der folgende Datensatz:

load(url("http://www.produnis.de/R/OrdinalSample.RData"))
mydata <- ordinalSample
head(mydata)
##    Konflikt Zufriedenh Geschlecht Stimmung
## 1       4.2       2.50          1  maessig
## 6       1.8       3.00          1 sehr gut
## 15      3.4       1.75          1  maessig
## 16      4.0       2.50          2  maessig
## 22      2.2       2.50          1      gut
## 23      3.4       2.25          1  maessig
  • Die Variable “Stimmung” dient hier als ordinales Ziel mit den Ausprägungen “schlecht” < “maessig” < “gut” < “sehr gut”.

  • Die Variable “Konflikt” dient als Prädiktor. Sie gibt an, wieviele Konflikte derzeit im Arbeitsleben vorliegen.

  • Die Variable “Zufriedenh” beschreibt, wie zufrieden Probanden mit ihrem Job sind.

  • Die Variable “Stimmung” soll nun durch “Konflikt” und “Zufriedenh” beschrieben werden.

Wir wollen wissen, wie stark Stimmung durch Konflikt und Zufriedenh beinflusst wird. Das heisst, in unserem Modell soll die Variable Stimmung durch Konflikt und Zufriedenh beschrieben werden.

30.7.3.1 Proportional Odds Modell

Ein Proportional Odds Modell errechnet sich leicht mit dem VGAM-Paket und der Funktion vglm():

library(VGAM)
pom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=T))

# oder abgekürzt
pom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=propodds)
summary(pom)
## 
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = propodds, 
##     data = mydata)
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1   0.6210     0.6482   0.958   0.3381    
## (Intercept):2  -1.7703     0.6540  -2.707   0.0068 ** 
## (Intercept):3  -4.0578     0.6759  -6.004 1.93e-09 ***
## Konflikt       -0.5762     0.0970  -5.940 2.86e-09 ***
## Zufriedenh      1.3643     0.2021   6.751 1.47e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y>=2]), logitlink(P[Y>=3]), 
## logitlink(P[Y>=4])
## 
## Residual deviance: 920.5493 on 1240 degrees of freedom
## 
## Log-likelihood: -460.2746 on 1240 degrees of freedom
## 
## Number of Fisher scoring iterations: 5 
## 
## No Hauck-Donner effect found in any of the estimates
## 
## 
## Exponentiated coefficients:
##   Konflikt Zufriedenh 
##  0.5620549  3.9129275

Mit dem Modell können nun weitere Parameter errechnet werden:

# Koeffizienten
pom.coef <- (coef(summary(pom)));pom.coef
##                 Estimate Std. Error    z value     Pr(>|z|)
## (Intercept):1  0.6209962 0.64821816  0.9580049 3.380603e-01
## (Intercept):2 -1.7702692 0.65404383 -2.7066523 6.796539e-03
## (Intercept):3 -4.0577657 0.67586741 -6.0037895 1.927645e-09
## Konflikt      -0.5761558 0.09700358 -5.9395317 2.858375e-09
## Zufriedenh     1.3642858 0.20208974  6.7508909 1.469400e-11

Die Odds Ratio errechnet sich mit der exp()-Funktion:

# Odds Ratio
pom.odds <- exp(coef(pom));pom.odds
## (Intercept):1 (Intercept):2 (Intercept):3      Konflikt    Zufriedenh 
##     1.8607808     0.1702871     0.0172876     0.5620549     3.9129275
# Devianz und AIC
pom.devi <- deviance(pom);pom.devi
## [1] 920.5493
pom.aic  <- AIC(pom);pom.aic
## [1] 930.5493
# logLikelihood
pom.ll <- logLik(pom);pom.ll
## [1] -460.2746
# 0-Modell (fuer pseudo R^2)
p0 <- vglm(Stimmung ~ 1, data=mydata, family=propodds)
p0.ll <- logLik(p0)
# R^2 McFadden
pom.mcfad <- as.vector(1- (pom.ll/p0.ll));pom.mcfad
## [1] 0.1240613
# R^2 Cox&Snell
N <- length(mydata[,1]) # Anzahl der Fälle
pom.cox <- as.vector(1 - exp((2/N) * (p0.ll - pom.ll)));pom.cox
## [1] 0.2696034
# R^2 Nagelkerke
pom.nagel <- as.vector((1 - exp((2/N) * (p0.ll - pom.ll))) / (1 - exp(p0.ll)^(2/N))) ;pom.nagel
## [1] 0.2928789

Das Proportional Odds Modell geht von der “equal slopes assumption” (auch: “proportional odds assumption”) aus. Diese Annahme muss geprüft werden, bevor das Modell als gültig angesehen werden kann. Dies geht mit dem VGAM-Paket recht einfach. Der Befehl vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=propodds) ist eine Abkürzung für vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=T)). Der Parameter parallel=TRUE/FALSE stellt ein, ob das Modell mit equal slopes (=TRUE), oder ohne equal slopes assumption (=FALSE) erstellt werden soll. Zur Überprüfung der “Equal Slopes Assumption” erstellt man jeweils ein TRUE und ein FALSE-Modell, und führt dann einen Likelihood-Test durch. Die Nullhypothese lautet, dass equal slopes vorliegen.

# Modell OHNE equal slopes
npom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=F))

# log-likelihood-Test auf Equal Slopes Assumption
lrtest(pom,npom)# Test
## Likelihood ratio test
## 
## Model 1: Stimmung ~ Konflikt + Zufriedenh
## Model 2: Stimmung ~ Konflikt + Zufriedenh
##    #Df  LogLik Df  Chisq Pr(>Chisq)
## 1 1240 -460.27                     
## 2 1236 -456.46 -4 7.6202     0.1065

Der Test ist mit p=0,1065 nicht signifikant. Das heisst, es liegen equal slopes vor. Das Modell darf also verwendet werden.

Alternativ kann dieser Test auch auf die Devianz der Modelle bestimmt werden:

# Devianz-Test auf Equal Slopes Assumption
pom.pdevi = (1 - pchisq(deviance(pom) - deviance(npom), df=df.residual(pom)-df.residual(npom)));pom.pdevi
## [1] 0.106526

Ist der Test signifikant, liegen keine “equal slopes” vor. Hier kann ein Partial Proportional Odds Modell erstellt werden, welches die “equal slopes” nur für bestimmte Prädiktoren annimmt. Mit dem VGAM-Paket kann recht einfach bestimmt werden, wie ein solches Modell erstellt werden soll. Hierfür erweitertet man den “parallel”-Switch wie folgt:

# Equal Slopes NUR für "Konflikt"
ppom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=T~Konflikt-1));ppom
## 
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cumulative(parallel = T ~ 
##     Konflikt - 1), data = mydata)
## 
## 
## Coefficients:
## (Intercept):1 (Intercept):2 (Intercept):3      Konflikt  Zufriedenh:1 
##     0.6614816     1.5373699     2.7337389     0.5705897    -1.9684047 
##  Zufriedenh:2  Zufriedenh:3 
##    -1.2805318    -0.8865225 
## 
## Degrees of Freedom: 1245 Total; 1238 Residual
## Residual deviance: 914.4929 
## Log-likelihood: -457.2464
# Nochmals EqualSlopes NUR für "Konflikt"
ppom <- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cumulative(parallel=F~Zufriedenh));ppom
## 
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cumulative(parallel = F ~ 
##     Zufriedenh), data = mydata)
## 
## 
## Coefficients:
## (Intercept):1 (Intercept):2 (Intercept):3      Konflikt  Zufriedenh:1 
##     0.6614816     1.5373699     2.7337389     0.5705897    -1.9684047 
##  Zufriedenh:2  Zufriedenh:3 
##    -1.2805318    -0.8865225 
## 
## Degrees of Freedom: 1245 Total; 1238 Residual
## Residual deviance: 914.4929 
## Log-likelihood: -457.2464
  • parallel=T~Konflikt-1 bedeutet, dass equal slopes nur für Konflikt angenommen wird.

  • parallel=F~Zufriedenh bedeutet, dass equal slopes nur für Zufriedenh nicht angenommen wird.

Beide Befehle bedeuten also das selbe. Daher ist die R-Ausgabe bei beiden Befehlen gleich.

Eine Koeffizientenübersicht erhält man per

ppom.ce <- coef(summary(ppom));ppom.ce
##                 Estimate Std. Error    z value     Pr(>|z|)
## (Intercept):1  0.6614816 0.85866412  0.7703613 4.410856e-01
## (Intercept):2  1.5373699 0.71997824  2.1353005 3.273647e-02
## (Intercept):3  2.7337389 0.98084390  2.7871294 5.317723e-03
## Konflikt       0.5705897 0.09715456  5.8730094 4.279541e-09
## Zufriedenh:1  -1.9684047 0.34535219 -5.6997024 1.200167e-08
## Zufriedenh:2  -1.2805318 0.23585229 -5.4293805 5.654999e-08
## Zufriedenh:3  -0.8865225 0.32542569 -2.7241933 6.445877e-03

Mit einem kleinen Script können Konfidenzintervalle und andere Werte ausgegeben werden:

get.ci <- function(x){
  back <-cbind( x[1], # estimate
                (x[1] - (1.96 * x[2])), #lCI
                (x[1] + (1.96 * x[2])), #uCI
                x[2], x[3], # SD und z
                (2*(1 -pnorm(abs(x[3]))) ),# p-wert
                exp(x[1])) 
  colnames(back) <- c("Estimate", "lCI", "uCI","SD","z","p-value", "OR")
  return(back)
}

Der Aufruf erfolgt z.B. so:

get.ci(as.numeric(ppom.ce[4,]))
##       Estimate       lCI       uCI         SD        z      p-value      OR
## [1,] 0.5705897 0.3801667 0.7610126 0.09715456 5.873009 4.279541e-09 1.76931

30.7.3.2 Continuation Ratio Model

Um ein Continuation Ratio Modell zu rechnen, muss der Parameter family auf cratio gesetzt werden.

# CR-Modell MIT PO-Assumption  (Vorwärts)
crm<- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cratio(parallel=T,reverse=F))
summary(crm)
## 
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cratio(parallel = T, 
##     reverse = F), data = mydata)
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1  0.92610    0.57405   1.613   0.1067    
## (Intercept):2 -1.14132    0.57924  -1.970   0.0488 *  
## (Intercept):3 -3.01660    0.60706  -4.969 6.72e-07 ***
## Konflikt      -0.53645    0.08639  -6.209 5.32e-10 ***
## Zufriedenh     1.16010    0.17944   6.465 1.01e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y>1|Y>=1]), logitlink(P[Y>2|Y>=2]), 
## logitlink(P[Y>3|Y>=3])
## 
## Residual deviance: 924.3824 on 1240 degrees of freedom
## 
## Log-likelihood: -462.1912 on 1240 degrees of freedom
## 
## Number of Fisher scoring iterations: 6 
## 
## No Hauck-Donner effect found in any of the estimates

Dies berechnet standardmäßig die Vorwärts-Methode. Möchte man die Rückwärts-Methode rechnen, setzt man den Parameter reverse auf TRUE.

# CR-Modell MIT PO-Assumption (Rückwärts)
crm<- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cratio(parallel=T,reverse=TRUE))
summary(crm)
## 
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cratio(parallel = T, 
##     reverse = TRUE), data = mydata)
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1  0.12822    0.57259   0.224 0.822809    
## (Intercept):2  2.08751    0.58094   3.593 0.000326 ***
## (Intercept):3  4.04245    0.60636   6.667 2.62e-11 ***
## Konflikt       0.45234    0.08488   5.329 9.88e-08 ***
## Zufriedenh    -1.25631    0.18108  -6.938 3.98e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y<2|Y<=2]), logitlink(P[Y<3|Y<=3]), 
## logitlink(P[Y<4|Y<=4])
## 
## Residual deviance: 920.4627 on 1240 degrees of freedom
## 
## Log-likelihood: -460.2314 on 1240 degrees of freedom
## 
## Number of Fisher scoring iterations: 5 
## 
## No Hauck-Donner effect found in any of the estimates

Hat man die passende Methode gewählt, folgen die selben Befehle wie beim Proportional Odds Modell. Zunächst muss überprüft werden, ob “equal slopes” vorliegen:

# CR-Modell OHNE PO-Assumption (Rückwärts)
ncrm<- vglm(Stimmung ~ Konflikt+Zufriedenh, data=mydata, family=cratio(parallel=F,reverse=T));summary(ncrm)
## 
## Call:
## vglm(formula = Stimmung ~ Konflikt + Zufriedenh, family = cratio(parallel = F, 
##     reverse = T), data = mydata)
## 
## Coefficients: 
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept):1   2.3649     1.1740   2.014 0.043969 *  
## (Intercept):2   1.9998     0.8117   2.464 0.013753 *  
## (Intercept):3   2.4670     1.1190   2.205 0.027474 *  
## Konflikt:1      0.1465     0.1811   0.809 0.418633    
## Konflikt:2      0.4864     0.1194   4.073 4.65e-05 ***
## Konflikt:3      0.6308     0.1689   3.736 0.000187 ***
## Zufriedenh:1   -1.8008     0.3919  -4.596 4.32e-06 ***
## Zufriedenh:2   -1.2602     0.2578  -4.889 1.01e-06 ***
## Zufriedenh:3   -0.8351     0.3382  -2.469 0.013543 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Names of linear predictors: logitlink(P[Y<2|Y<=2]), logitlink(P[Y<3|Y<=3]), 
## logitlink(P[Y<4|Y<=4])
## 
## Residual deviance: 914.7924 on 1236 degrees of freedom
## 
## Log-likelihood: -457.3962 on 1236 degrees of freedom
## 
## Number of Fisher scoring iterations: 6 
## 
## No Hauck-Donner effect found in any of the estimates
# loglikelihood test auf equal slopes assumption
lrtest(ncrm,crm)# Test
## Likelihood ratio test
## 
## Model 1: Stimmung ~ Konflikt + Zufriedenh
## Model 2: Stimmung ~ Konflikt + Zufriedenh
##    #Df  LogLik Df  Chisq Pr(>Chisq)
## 1 1236 -457.40                     
## 2 1240 -460.23  4 5.6703     0.2252

Der Test ist nicht signifikant (p=0,23). Das bedeutet, dass die Annahme der equal slopes für die hier getesteten Rückwärts-Modelle beibehalten werden kann.

Nun können weitere Modellparameter bestimmt werden:

# 0-Modell (fuer pseudo R^2)
c0 <- vglm(Stimmung ~ 1, data=mydata, family=cratio(parallel=T))
c0.ll <- logLik(c0)
crm.ll <- logLik(crm)

# R^2 Nagelkerke
crm.nagel <- as.vector((1 - exp((2/N) * (c0.ll - crm.ll))) / (1 - exp(c0.ll)^(2/N))) ;crm.nagel
## [1] 0.2930443
# Devianz
crm.devi <- deviance(crm); crm.devi
## [1] 920.4627
# AIC
crm.aic  <- AIC(crm); crm.aic
## [1] 930.4627

Mit unserer Funktion von oben können wir uns die Modellparameter der Koeffizienten ausgeben lassen.

# Konfidenzinztervalle
crm.ce <- coef(summary(crm));crm.ce
##                 Estimate Std. Error    z value     Pr(>|z|)
## (Intercept):1  0.1282227 0.57259420  0.2239329 8.228095e-01
## (Intercept):2  2.0875126 0.58094099  3.5933299 3.264788e-04
## (Intercept):3  4.0424508 0.60636000  6.6667504 2.615293e-11
## Konflikt       0.4523353 0.08488357  5.3288915 9.881398e-08
## Zufriedenh    -1.2563060 0.18108066 -6.9378255 3.981811e-12
get.ci(as.numeric(crm.ce[4,]))
##       Estimate       lCI       uCI         SD        z      p-value       OR
## [1,] 0.4523353 0.2859635 0.6187071 0.08488357 5.328892 9.881398e-08 1.571979
get.ci(as.numeric(crm.ce[5,]))
##       Estimate       lCI        uCI        SD         z      p-value        OR
## [1,] -1.256306 -1.611224 -0.9013879 0.1810807 -6.937825 3.981704e-12 0.2847038

30.8 Überlebenszeitanalysen

30.8.1 Kaplan-Meier-Analysen

Zur Überlebenszeitanalyse nach dem Kaplan-Meier-Verfahren steht das Zusatzpaket survival zur Verfügung. Zum Ausprobieren der Funktionen liegt dem Paket der Datensatz ovarian bei.

# Beschreibung des Datensatzes lesen
?survival::ovarian
# speichere Testdatensatz "ovarian" in Objekt "ov"
ov <- survival::ovarian
head(ov)
##   futime fustat     age resid.ds rx ecog.ps
## 1     59      1 72.3315        2  1       1
## 2    115      1 74.4932        2  1       1
## 3    156      1 66.4658        2  1       2
## 4    421      0 53.3644        2  2       1
## 5    431      1 50.3397        2  1       1
## 6    448      0 56.4301        1  1       2

Bevor wir mit der Analyse beginnen, sollten wir den Datensatz noch aufbereiten. Wir überführen alle kategorialen Variablen in einen factor.

### Datensatz aufbereiten
## alle kategorialen Variablen
# rx ist die Behandlungsgruppe
ov$rx <- factor(ov$rx,
                levels = c("1", "2"),
                labels = c("A", "B"))

# resid.ds sind Nebenerkrankungen
ov$resid.ds <- factor(ov$resid.ds,
                      levels = c("1", "2"),
                      labels = c("nein", "ja"))

# ecog.ps ist der ECOG-Status
ov$ecog.ps <- factor(ov$ecog.ps,
                     levels = c("1", "2"),
                     labels = c("gut", "schlecht"))

Als ersten Analyseschritt erzeugen wir mit der Funktion Surv() ein “Survivalobjekt”. Der Funktion muss eine Zeitvariable mit der Beobachtungszeit sowie eine dichotome Statusvariable mit dem Ereignis (0 = nicht eingetreten; 1 = eingetreten) übergeben werden. Im ovarian-Datensatz ist die Zeit in futime und der Status in fustat gespeichert. Der Funktionsaufruf lautet dementsprechend

# Survivalobjet erzeugen
survObjekt <- survival::Surv(time=ov$futime, event=ov$fustat)

Mit der Funktion survfit wird nun das Kaplan-Meier-Modell gefittet. Der Funktion muss das eben erstellte Survivalobjekt übergeben werden.

### Kaplan-Meier
#  ~1 bedeutet "Gesamtmodell", ohne Gruppierung
km_fit <- survival::survfit(survObjekt ~ 1, data=ov)

Das Modell kann als Tabelle oder Plot ausgegeben werden.

## Tabelle anschauen
#  "survival"-Spalte ist die Überlebenswahrscheinlichkeit
summary(km_fit)
## Call: survfit(formula = survObjekt ~ 1, data = ov)
## 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    59     26       1    0.962  0.0377        0.890        1.000
##   115     25       1    0.923  0.0523        0.826        1.000
##   156     24       1    0.885  0.0627        0.770        1.000
##   268     23       1    0.846  0.0708        0.718        0.997
##   329     22       1    0.808  0.0773        0.670        0.974
##   353     21       1    0.769  0.0826        0.623        0.949
##   365     20       1    0.731  0.0870        0.579        0.923
##   431     17       1    0.688  0.0919        0.529        0.894
##   464     15       1    0.642  0.0965        0.478        0.862
##   475     14       1    0.596  0.0999        0.429        0.828
##   563     12       1    0.546  0.1032        0.377        0.791
##   638     11       1    0.497  0.1051        0.328        0.752

Die Spalte survival gibt die entsprechenden Überlebenswahrscheinlichkeiten an.

## plot()  (sehr hässlich)
plot(km_fit)

Etwas hübschere Plots lassen sich mit der Funktion autoplot() erzeugen. Damit die Funktion die Modelldaten versteht, muss das Zusatzpaket ggfortify geladen werden.

# etwas hübscher
library(ggfortify)
autoplot(km_fit)

Darüber hinaus bietet das Paket survminer die Funktion ggsurvplot():

## survminer ggsurvplot
survminer::ggsurvplot(km_fit)

Das Modell lässt sich auch für Gruppenvergleiche erstellen. Die Variable ov$rx unterteilt den Datensatz in zwei Behandlungsgruppen. Das Überlebensmodell pro Gruppe lässt sich wie folgt erzeugen:

### Gruppierungen
# nach Behandlungsgruppe "rx" gruppieren
km_fit <- survival::survfit(survObjekt ~ rx, data=ov)

Mittels summary() werden die Überlebenstabellen pro Gruppe ausgegeben.

summary(km_fit)
## Call: survfit(formula = survObjekt ~ rx, data = ov)
## 
##                 rx=A 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##    59     13       1    0.923  0.0739        0.789        1.000
##   115     12       1    0.846  0.1001        0.671        1.000
##   156     11       1    0.769  0.1169        0.571        1.000
##   268     10       1    0.692  0.1280        0.482        0.995
##   329      9       1    0.615  0.1349        0.400        0.946
##   431      8       1    0.538  0.1383        0.326        0.891
##   638      5       1    0.431  0.1467        0.221        0.840
## 
##                 rx=B 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##   353     13       1    0.923  0.0739        0.789        1.000
##   365     12       1    0.846  0.1001        0.671        1.000
##   464      9       1    0.752  0.1256        0.542        1.000
##   475      8       1    0.658  0.1407        0.433        1.000
##   563      7       1    0.564  0.1488        0.336        0.946

Im Plot werden die Kaplan-Meier-Kurven beider Gruppen angezeigt.

# plot (hässlich)
plot(km_fit, conf.int=T) # mit Konfidenzintervallen

# plot (hübscher)
plot(km_fit, conf.int=TRUE, col=c("red", "blue"), yscale=100,
     xlab = "Überlebenszeit in Tagen", ylab = "Prozent Überleben")
legend("bottomright", title="Gruppe", c("1", "2"),
       fill=c("red", "blue"))

Die Funktion ggsurvplot() liefert auf Wunsch den p-Wert des Gruppenunterschieds mit:

## survminer ggsurvplot mit p-Wert
survminer::ggsurvplot(km_fit, pval=TRUE, conf.int = TRUE)

Der p-Wert lässt sich mit dem Lograng-Test auch selbst berechnen. Hierfür steht die Funktion survdiff() zur Verfügung.

# p-Wert selbst berechnen
## Lograng-Test
survival::survdiff(survObjekt ~rx, data=ov)
## Call:
## survival::survdiff(formula = survObjekt ~ rx, data = ov)
## 
##       N Observed Expected (O-E)^2/E (O-E)^2/V
## rx=A 13        7     5.23     0.596      1.06
## rx=B 13        5     6.77     0.461      1.06
## 
##  Chisq= 1.1  on 1 degrees of freedom, p= 0.3

30.8.2 Cox-Regression

Neben der Kaplan-Meier-Analyse können die Daten auch per Cox-Regression analysiert werden. Der Vorteil besteht darin, dass Prädiktoren verschiedensten Datenniveaus in die Analyse integriert werden können. Das Paket survival bietet hierzu die Funktion coxph().

## Cox-Regression
cox <- survival::coxph(survObjekt ~ rx+resid.ds+age+ecog.ps, data=ov)
# Modell ausgeben
cox
## Call:
## survival::coxph(formula = survObjekt ~ rx + resid.ds + age + 
##     ecog.ps, data = ov)
## 
##                     coef exp(coef) se(coef)      z       p
## rxB             -0.91450   0.40072  0.65332 -1.400 0.16158
## resid.dsja       0.82619   2.28459  0.78961  1.046 0.29541
## age              0.12481   1.13294  0.04689  2.662 0.00777
## ecog.psschlecht  0.33621   1.39964  0.64392  0.522 0.60158
## 
## Likelihood ratio test=17.04  on 4 df, p=0.001896
## n= 26, number of events= 12

Die Spalte exp(coef) entspricht der Hazard-Ratio, mit welcher Richtung und Stärke des jeweiligen Einflusses interpretiert werden kann.

Das Modell kann mittels ggforest() geplottet werden:

# Forestplot der Hazard-Ratio
survminer::ggforest(cox, data=ov)

30.9 Faktorenanalyse

Die beiden am Häufigsten genutzen Verfahren sind die Hauptkomponentenanalyse (PCA) und die Hauptachsenanalyse (PAF). Sie unterscheiden sich in den Varianzanteilen, die in den Items enthalten sein können. Die PCA versucht, die gesamte Varianz auf Faktoren zurückzuführen, die PAF die gemeinsame Varianz der Items. In psychologischen Anwendungen (= am Menschen) wird daher eher die PAF verwendet.

30.9.1 explorative Faktorenanalyse

Bei der explorativen Faktorenanalyse (EFA) erfolgt die Auswertung (häufig) nicht theoriegeleitet. Ziel ist es vielmehr zu erfahren, mit wievielen Faktoren die Daten bestmöglich erklärt werden können.

Wir verwenden den Datensatz Faktorenbogen, der wie folgt geladen werden kann.

load(url("https://www.produnis.de/R/Faktorenbogen.rda"))

# anschauen
psych::describe(Faktorenbogen)
##         vars   n  mean    sd median trimmed   mad min max range  skew kurtosis
## age        1 150 53.11 18.82   53.0   54.16 17.05   1 101   100 -0.47     0.29
## gender*    2 150  1.65  0.55    2.0    1.64  0.00   1   3     2  0.02    -0.87
## A          3 150  4.41  2.03    5.0    4.40  2.97   0  10    10  0.05    -0.43
## B          4 150  1.92  0.84    2.0    1.84  1.48   1   5     4  0.76     0.43
## C          5 150  7.91  6.19    7.0    7.22  5.93   0  30    30  1.12     1.31
## D          6 150 12.23  5.05   12.0   11.95  4.45   3  27    24  0.51    -0.04
## E          7 150 14.77  2.39   15.0   14.69  2.97  10  21    11  0.28    -0.13
## F          8 150  5.00  1.32    5.0    5.03  1.48   2   8     6 -0.12    -0.67
## G          9 150 10.35  4.08    9.5   10.06  3.71   3  24    21  0.78     0.76
## H         10 150 15.92  2.93   16.0   15.84  2.97  10  24    14  0.21    -0.39
## I         11 150 85.02 68.19   65.0   75.93 51.15   0 356   356  1.40     2.05
## J         12 150  5.02  2.25    5.0    4.93  1.48   1  11    10  0.33    -0.12
## K         13 150  5.44  2.08    6.0    5.38  2.97   1  12    11  0.31    -0.10
## L         14 150  1.47  0.65    1.0    1.36  0.00   1   4     3  1.20     0.86
##           se
## age     1.54
## gender* 0.04
## A       0.17
## B       0.07
## C       0.51
## D       0.41
## E       0.20
## F       0.11
## G       0.33
## H       0.24
## I       5.57
## J       0.18
## K       0.17
## L       0.05

Die Items für die Faktorenanalyse befinden sich in den Variablen A bis L.

items <- Faktorenbogen %>%
  select(A:L)

Mit dem Bartlett-Test kann überprüft werden, ob die Items miteinander korrelieren (das sollten sie tun).

psych::cortest.bartlett(items)$p.value
## R was not square, finding R from data
## [1] 8.843226e-151

Das Ergebnis ist signifikant.

Mittels Kaiser-Meyer-Olkin-Verfahren (KMO) wird geprüft, wie gut die Daten für eine Faktorenanalyse geeignet sind, wobei Werte zwischen 0 (für “nicht geeignet”) und 1 (für “perfekt geeignet”) angenommen werden. Die Minimum Average Partial (MSA) ist ein Maß der internen Konsistenz, und nimmt ebenfalls Werte zwischen 0 (nicht konsistent) und 1 (perfekt konsistent) an. Beide Werte sollte größer als 0,5 sein.

# berechne KMO und MSA
psych::KMO(items)
## Kaiser-Meyer-Olkin factor adequacy
## Call: psych::KMO(r = items)
## Overall MSA =  0.69
## MSA for each item = 
##    A    B    C    D    E    F    G    H    I    J    K    L 
## 0.36 0.72 0.66 0.68 0.60 0.73 0.24 0.73 0.75 0.62 0.78 0.61

Die “Overall MSA” entspricht dem KMO.

Mittels Screeplot kann die Faktorenzahl graphisch entschieden werden. Dabei sollte der letzte “gültige” Faktor vor dem deutlichen Knick zu einer abflachenden Gerade genutzt werden. Darüber hinaus sollten die Eigenwerte der Faktoren größer 1 sein.

Mit R-Hausmitteln kann ein Screeplot wie folgt erzeugt werden.

plot(eigen(cor(items))$values, type="b")
abline(h=1)

Schöner geht es mit der Funktion fa.parallel().

psych::fa.parallel(items)

## Parallel analysis suggests that the number of factors =  3  and the number of components =  3

Es werden 3 Faktoren/Komponenten vorgeschlagen.

Die Funktion berechnet das Screeplot sowohl für PAC (Sternchen) als auch PAF (Dreieck). Mit dem Parameter fa = "fa" kann die PAC-Methode entfernt werden.

psych::fa.parallel(items, fa ="fa")

## Parallel analysis suggests that the number of factors =  3  and the number of components =  NA

Die Daten lassen sich am besten mit 3 Faktoren erklären.

30.9.2 konfirmatorische Faktorenanalyse

Ist die Anzahl der Faktoren bekannt (z.B. theoriegeleitet), wird die Hauptachsenanaylse wie folgt durchgeführt.

# Hauptachsen-Faktorenanalyse (PAF)
# mit 3 Faktoren, Varimax-Rotation
fit <- psych::fa(items, nfactors=3,
                 rotate = "varimax")
fit
## Factor Analysis using method =  minres
## Call: psych::fa(r = items, nfactors = 3, rotate = "varimax")
## Standardized loadings (pattern matrix) based upon correlation matrix
##     MR1   MR2   MR3     h2    u2 com
## A  0.00  0.08  0.03 0.0069 0.993 1.2
## B -0.01  0.83 -0.08 0.6997 0.300 1.0
## C -0.05  0.88 -0.01 0.7770 0.223 1.0
## D  0.96 -0.09 -0.07 0.9266 0.073 1.0
## E -0.02  0.10  0.88 0.7820 0.218 1.0
## F -0.04 -0.12  0.59 0.3655 0.635 1.1
## G  0.00  0.00  0.07 0.0048 0.995 1.0
## H  0.94 -0.02 -0.08 0.8855 0.115 1.0
## I -0.03  0.78 -0.05 0.6107 0.389 1.0
## J -0.08  0.08  0.88 0.7866 0.213 1.0
## K  0.89  0.03 -0.01 0.7896 0.210 1.0
## L -0.10  0.09  0.14 0.0393 0.961 2.6
## 
##                        MR1  MR2  MR3
## SS loadings           2.61 2.13 1.94
## Proportion Var        0.22 0.18 0.16
## Cumulative Var        0.22 0.39 0.56
## Proportion Explained  0.39 0.32 0.29
## Cumulative Proportion 0.39 0.71 1.00
## 
## Mean item complexity =  1.2
## Test of the hypothesis that 3 factors are sufficient.
## 
## The degrees of freedom for the null model are  66  and the objective function was  6.38 with Chi Square of  920.48
## The degrees of freedom for the model are 33  and the objective function was  0.32 
## 
## The root mean square of the residuals (RMSR) is  0.03 
## The df corrected root mean square of the residuals is  0.05 
## 
## The harmonic number of observations is  150 with the empirical chi square  20.2  with prob <  0.96 
## The total number of observations was  150  with Likelihood Chi Square =  45.01  with prob <  0.079 
## 
## Tucker Lewis Index of factoring reliability =  0.971
## RMSEA index =  0.049  and the 90 % confidence intervals are  0 0.083
## BIC =  -120.34
## Fit based upon off diagonal values = 0.99
## Measures of factor score adequacy             
##                                                    MR1  MR2  MR3
## Correlation of (regression) scores with factors   0.98 0.94 0.94
## Multiple R square of scores with factors          0.96 0.88 0.88
## Minimum correlation of possible factor scores     0.92 0.76 0.77

Die Ausgabe von fit ist ziemlich lang. Wie gewohnt besteht die Möglichkeit, mittels $ auf die einzelnen Ergebnisse zuzugreifen. Eine Übersicht erhält man mittels names().

names(fit)
##  [1] "residual"      "dof"           "chi"           "nh"           
##  [5] "rms"           "EPVAL"         "crms"          "EBIC"         
##  [9] "ESABIC"        "fit"           "fit.off"       "sd"           
## [13] "factors"       "complexity"    "n.obs"         "objective"    
## [17] "criteria"      "STATISTIC"     "PVAL"          "Call"         
## [21] "null.model"    "null.dof"      "null.chisq"    "TLI"          
## [25] "RMSEA"         "BIC"           "SABIC"         "r.scores"     
## [29] "R2"            "valid"         "score.cor"     "weights"      
## [33] "rotation"      "hyperplane"    "communality"   "communalities"
## [37] "uniquenesses"  "values"        "e.values"      "loadings"     
## [41] "model"         "fm"            "rot.mat"       "Structure"    
## [45] "method"        "scores"        "R2.scores"     "r"            
## [49] "np.obs"        "fn"            "Vaccounted"
fit$communality
##           A           B           C           D           E           F 
## 0.006940720 0.699704328 0.776989353 0.926609287 0.782014553 0.365479522 
##           G           H           I           J           K           L 
## 0.004811051 0.885466837 0.610670890 0.786557566 0.789638724 0.039273937
fit$loadings
## 
## Loadings:
##   MR1    MR2    MR3   
## A                     
## B         0.833       
## C         0.880       
## D  0.956              
## E         0.103  0.878
## F        -0.118  0.592
## G                     
## H  0.938              
## I         0.779       
## J                0.879
## K  0.888              
## L -0.104         0.141
## 
##                  MR1   MR2   MR3
## SS loadings    2.605 2.129 1.940
## Proportion Var 0.217 0.177 0.162
## Cumulative Var 0.217 0.394 0.556

Mit diesen Informationen können die Items ein bisschen umsortiert…

  • Faktor 1: D, H, K
  • Faktor 2: B, C, I
  • Faktor 3: E, F, J

… und ein Korrelationsmatrix erstellt werden.

# ändere Reihenfolgt
items2 <- items %>%
  select(D, H, K, B, C, I, E, F, J)

# Korreltationsmatrix plotten
corrplot::corrplot(corr = cor(items2), # Korrelationsmatrix
         method = "color",   # Ausprägung farblich kodieren
         addCoef.col = "black", # schreibe Werte in schwarz
         number.cex = 0.7)   # Schriftgröße

30.9.3 Hauptkomponentenanalyse

Die Hauptkomponentenanalyse kann mit der Funktion principal() durchgeführt werden.

# Hauptkomponentenanalyse (PCA)
# mit 3 Hauptkomponenten und
# Varimax-Rotation
fit2 <- psych::principal(items, nfactors=3, rotate="varimax")

# anschauen
fit2
## Principal Components Analysis
## Call: psych::principal(r = items, nfactors = 3, rotate = "varimax")
## Standardized loadings (pattern matrix) based upon correlation matrix
##     RC1   RC2   RC3    h2    u2 com
## A  0.01  0.13  0.05 0.019 0.981 1.3
## B -0.01  0.89 -0.08 0.793 0.207 1.0
## C -0.05  0.90  0.00 0.821 0.179 1.0
## D  0.96 -0.08 -0.07 0.928 0.072 1.0
## E -0.03  0.10  0.89 0.803 0.197 1.0
## F -0.03 -0.15  0.76 0.606 0.394 1.1
## G  0.01  0.00  0.12 0.015 0.985 1.0
## H  0.95 -0.02 -0.07 0.916 0.084 1.0
## I -0.04  0.86 -0.05 0.752 0.248 1.0
## J -0.09  0.08  0.89 0.806 0.194 1.0
## K  0.94  0.04  0.00 0.884 0.116 1.0
## L -0.14  0.14  0.24 0.095 0.905 2.3
## 
##                        RC1  RC2  RC3
## SS loadings           2.74 2.44 2.26
## Proportion Var        0.23 0.20 0.19
## Cumulative Var        0.23 0.43 0.62
## Proportion Explained  0.37 0.33 0.30
## Cumulative Proportion 0.37 0.70 1.00
## 
## Mean item complexity =  1.2
## Test of the hypothesis that 3 components are sufficient.
## 
## The root mean square of the residuals (RMSR) is  0.06 
##  with the empirical chi square  61.27  with prob <  0.002 
## 
## Fit based upon off diagonal values = 0.96
# welche $-Werte sind verfügbar?
names(fit2)
##  [1] "values"       "rotation"     "n.obs"        "communality"  "loadings"    
##  [6] "fit"          "fit.off"      "fn"           "Call"         "uniquenesses"
## [11] "complexity"   "chi"          "EPVAL"        "R2"           "objective"   
## [16] "residual"     "rms"          "factors"      "dof"          "null.dof"    
## [21] "null.model"   "criteria"     "STATISTIC"    "PVAL"         "weights"     
## [26] "r.scores"     "rot.mat"      "Vaccounted"   "Structure"    "scores"
# Faktorenladungen
fit2$loadings
## 
## Loadings:
##   RC1    RC2    RC3   
## A         0.128       
## B         0.887       
## C         0.905       
## D  0.958              
## E         0.102  0.890
## F        -0.146  0.764
## G                0.121
## H  0.954              
## I         0.865       
## J                0.890
## K  0.940              
## L -0.135  0.143  0.238
## 
##                  RC1   RC2   RC3
## SS loadings    2.744 2.435 2.260
## Proportion Var 0.229 0.203 0.188
## Cumulative Var 0.229 0.432 0.620

Ein passendes Video von Daniela Keller findet sich bei Youtube unter https://www.youtube.com/watch?v=NFPGQcq1fO8.