32  Schließende Statistik

32.1 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.

Die Formel der Modellgeraden lautet also \(y = 1,208 + 1,055\cdot x\).

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

# lade Testdaten
load(url("https://www.produnis.de/R/data/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)

Denken Sie aber daran, anschließend die Plotausgabe wieder zurückzusetzen:

par(mfrow = c(1, 1))

32.1.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)

par(mfrow = c(1, 1))

32.1.2 nicht-lineare Regression

Häufig folgen die vorhandenen Daten keinem linearen Zusammenhang.

# erzeuge dummy-Werte
x <- c(6, 9, 12, 14, 30, 35, 40, 47, 51, 55, 60)
y <- c(14, 28, 50, 70, 89, 94, 90, 75, 59, 44, 27)

# Plotte die Daten mit Regressionsgeraden
plot(x, y, col="blue")
abline(lm(y~x), col="red", lty=3)

Das Aussehen des Plots lässt schon erahnen, dass kein linearer Zusammenhang vorliegt. Wenn wir ein lineares Modell auf die Daten fitten, ist R2 sehr niedrig.

# Modell fitten
fit <- lm(y~x)
# Modelldaten ausgeben
summary(fit)
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -39.34 -21.99  -2.03  23.50  35.11 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  48.4531    17.3288   2.796   0.0208 *
## x             0.2981     0.4599   0.648   0.5331  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28.72 on 9 degrees of freedom
## Multiple R-squared:  0.0446, Adjusted R-squared:  -0.06156 
## F-statistic: 0.4201 on 1 and 9 DF,  p-value: 0.5331

R2 ist hier 0,0446, was bedeutet, dass nur 4% der Daten durch das Modell erklärt werden können.

Das Plot legt nahe, dass ein quadratischer Zusammenhang besteht, denn die Punkte scheinen einer Parabel zu folgen. Die Formel einer quadratischen Funktion lautet vereinfacht \(y = x^2 + x\). Den rechten Teil der Formel müssen wir innerhalb von lm() verwenden. Leider können wir die Formel nicht direkt in lm() schreiben, weil innerhalb von lm() nach dem ~-Zeichen keine weiteren Rechenoperationen möglich sind. Das liegt daran, dass die Prädiktoren mittels + und - übergeben werden. Unser Trick besteht nun darin, alle Komponeten der Formel als eigenständige Prädiktoren des Modells anzugeben. Dazu müssen wir zunächst ein neues Objekt erzeugen, welches die x2 -Werte enthält:

# quadrierte Werte als eigenes Objekt
x2 <- x^2
# Modell erstellen mit quadratischer Formel
fit <- lm(y ~ x2 + x)
# Modellwerte ausgeben
summary(fit)
## 
## Call:
## lm(formula = y ~ x2 + x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.2484 -3.7429 -0.1812  1.1464 13.6678 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -18.25364    6.18507  -2.951   0.0184 *  
## x2           -0.10120    0.00746 -13.565 8.38e-07 ***
## x             6.74436    0.48551  13.891 6.98e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.218 on 8 degrees of freedom
## Multiple R-squared:  0.9602, Adjusted R-squared:  0.9502 
## F-statistic: 96.49 on 2 and 8 DF,  p-value: 2.51e-06

R2 liegt nun bei 0,9602, was bedeutet, dass 96% der Daten durch das Modell erklärt werden können. Die Funktionsformel lautet \(y = -0,10120\cdot x^2 + 6,74436\cdot x - 18,25364\).

Mit der Funktion predict() lassen sich nun für jedes x die entsprechenden y-Werte des Modells bestimmen (vorhersagen). Der Funktion wird das gefittete Objekt übergeben, sowie eine Liste der gewünschten Prädiktorenwerte. Die so erzeugten Modellwerte können dann ins Plot gezeichnet werden.

vorhersage <- predict(fit, list(x, x2))


plot(x, y, col="blue")
lines(x, vorhersage, col='darkgreen')

Um eine schöne Parabel zu zeichnen, haben wir zu wenige x-Werte. Daher erstellen wir einen neuen Vektor, der kontinuierliche Werte für x enthält. Mit diesem können wir dann die entsprechenden (vielen) y-Werte mittels predict() bestimmen und in das Plot einzeichnen. So entsteht - bei genügend Werten - eine schöne Modellparabel.

# Erstelle kontinuierliche Werte
# für x von 0 bis 60
xref <- seq(0, 60, 0.1)
# berechne dazugehörige y-Modellwerte
yref <- predict(fit, list(x = xref, 
                          x2 = xref^2))

# Plotten
plot(x, y, col="blue")
# Modell-Parabel einzeichnen
lines(xref, yref, col='darkgreen')
abline(lm(y~x), col="red", lty=3)

32.1.2.1 Funktionen höherer Ordnung

Wird ein Zusammenhang mit einer Funktion höherer Ordnung vermutet, gehen wir analog zur quadratischen Funktion vor. In diesem Beispiel nehmen wir an, der Zusammenhang ließe sich vereinfacht durch die Formel \(y = x^3 + x^2 +x\) beschreiben. Die Daten verhalten sich wie folgt:

# erzeuge Testwerte
x <- c(0, 1, 2, 3, 4, 5, 6, 7, 8, 0.5, 1.5, 2.5, 3.5, 4.5, 5.5, 6.5, 7.5)
y <- c(1, 46, 71, 82,  82,  86,  91, 106, 137,  17,  51,  68, 77,  76,  82.1,  87, 108)

# plotten
plot(x, y, col="blue")

# falsche Regressionsgerade
abline(lm(y~x), col="red", lty=3)

summary(lm(y~x))
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -26.6471  -7.6728  -0.3309   8.2743  19.8794 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   27.647      6.387   4.329 0.000596 ***
## x             11.737      1.362   8.619  3.4e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.75 on 15 degrees of freedom
## Multiple R-squared:  0.832,  Adjusted R-squared:  0.8208 
## F-statistic: 74.29 on 1 and 15 DF,  p-value: 3.398e-07

Die Regressionsgerade beschreibt die Daten nicht gut, aber das R2 des linearen Modells ist mit 0,832 recht hoch. Das Modell beschreibt also 83% der vorhandenen Daten.

Versuchen wir nun, die Funktion dritten Grades ins Modell zu übernehmen. Wie bei den quadratischen Funktionen müssen wir die Werte für x3 und x2 in eigene Vektoren schreiben, die wir dann als Prädiktoren innerhalb von lm() übergeben.

# erzeuge Werte für x^3 und x^2
x3 <- x*x*x
x2 <- x^2

# Modell fitten
fit <- lm(y ~ x3 + x2 +x)
summary(fit)
## 
## Call:
## lm(formula = y ~ x3 + x2 + x)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.038 -4.905  1.872  4.254  5.005 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -2.0945     4.0612  -0.516    0.615    
## x3            0.9771     0.1101   8.876 7.04e-07 ***
## x2          -12.6470     1.3414  -9.428 3.55e-07 ***
## x            55.5108     4.5344  12.242 1.64e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.14 on 13 degrees of freedom
## Multiple R-squared:  0.9797, Adjusted R-squared:  0.975 
## F-statistic: 208.8 on 3 and 13 DF,  p-value: 3.038e-11

Das R2 liegt bei 0,9797, was bedeutet, dass knapp 98% der Daten durch unser Modell erklärt werden können. Dieses Modell ist also wesentlich besser als das zuvor berechnete lineare Modell. Die Modellformel lautet \(y = 0,9772\cdot x^3 - 12,647\cdot x^2 + 55,5108\cdot x - 2,0945\).

Um die Funktionskurve des Modells in das Plot zu zeichnen, erstellen wir zunächst einen neuen Vektor, der kontinuierliche Werte für x enthält. Über die Funktion predict() können dann die Modellwerte für y berechnet werden.

# Erstelle kontinuierliche Werte
# für x von 0 bis 8
xref <- seq(0, 8, 0.1)
# berechne dazugehörige y-Modellwerte
yref <- predict(fit, list(x = xref, 
                          x2 = xref^2,
                          x3 = xref*xref*xref))

# Plotten
plot(x, y, col="blue")
# Modell-Parabel einzeichnen
lines(xref, yref, col='darkgreen')
abline(lm(y~x), col="red", lty=3 )

32.1.2.2 Exponentialfunktion

Wird ein exponentieller Zusammenhang vermutet, gehen wir analog zur quadratischen Regression vor.

In diesem Beispiel nehmen wir an, der Zusammenhang ließe sich vereinfacht durch die Formel \(y = a^x\) beschreiben. Die Daten verhalten sich wie folgt:

# Testdaten
x <- c(2, 4, 5, 6, 7, 9, 10, 11, 13, 14, 16, 17, 18, 20, 21, 22, 22, 23, 24, 24, 25)
y <- c(1, 2, 2, 1, 3, 3, 2, 3, 2, 5, 4, 7, 10, 7, 15, 12, 18, 20, 17, 19, 20)

# plotten
plot(x, y, col="blue")
# falsche Regressionsgerade
abline(lm(y~x), col="red", lty=3)

# lineares Modell
summary(lm(y~x))
## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.5601 -2.2566  0.3153  3.0118  4.8952 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.40475    1.60736  -2.740    0.013 *  
## x            0.84824    0.09688   8.756 4.27e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.236 on 19 degrees of freedom
## Multiple R-squared:  0.8014, Adjusted R-squared:  0.7909 
## F-statistic: 76.67 on 1 and 19 DF,  p-value: 4.275e-08

Die Regressionsgerade beschreibt die Daten nicht gut, aber das R2 des linearen Modells ist mit 0,8014 recht hoch. Das Modell beschreibt also 80% der vorhandenen Daten.

Versuchen wir nun, die Exponentialfunktion ins Modell zu übernehmen. Hierfür fitten wir das Modell nicht auf y, sondern auf den natürlichen Logarithmus yon y.

# Modell fitten
fit <- lm(log(y) ~ x)
summary(fit)
## 
## Call:
## lm(formula = log(y) ~ x)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.73254 -0.10440  0.01856  0.24804  0.44867 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.25510    0.16442  -1.551    0.137    
## x            0.12929    0.00991  13.047 6.23e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.331 on 19 degrees of freedom
## Multiple R-squared:  0.8996, Adjusted R-squared:  0.8943 
## F-statistic: 170.2 on 1 and 19 DF,  p-value: 6.232e-11

Mit R2 = 0,8996 beschreibt das Modell die Daten besser, als das zuvor erstellte lineare Modell. Die Funktionsformel lautet \(ln(y) = 0,12929\cdot x -0,2551\).

Um die Funktionskurve des Modells in das Plot zu zeichnen, erstellen wir zunächst einen neuen Vektor, der kontinuierliche Werte für x enthält. Über die Funktion predict() können dann die Modellwerte für y berechnet werden. Da wir das Modell auf log(y) gefittet haben, müssen die Werte von predict() mittels exp() umgewandelt werden.

# Erstelle kontinuierliche Werte
# für x von 0 bis 8
xref <- seq(0, 25, 0.1)
# berechne dazugehörige y-Modellwerte
yref <- exp(predict(fit, data.frame(x = xref)))

# Plotten
plot(x, y, col="blue")
# Modell-Parabel einzeichnen
lines(xref, yref, col='darkgreen')
abline(lm(y~x), col="red", lty=3 )

32.2 Generalisierte lineare Modelle

Mit der Funktion glm() können generalisierte lineare Modelle gebildet werden. Strenggenommen können alle Regressionsmodelle mit glm() erzeugt werden. In R wird sie vor allem verwendet für

  • Poisson-Regressionen
  • logistische Regressionen
  • binomiale Regressionen

Das Standardvorgehen ist in allen Fällen gleich, das gefittete Modell wird in ein Objekt geschrieben

fit <- glm(ZIEL ~ Prädiktor1 + Prädiktor2, family= gaussian)

wobei über den Parameter family die Art des Modells festgelgt wird (z.B. binomial, poisson oder gaussian).

Anschließend werden typischerweise die folgenden Befehle angewendet:

# Modellübersicht
summary(fit)

# Koeffizientenwerte
coef(fit)
#oder
fit$coefficients

# Konfidenzintervalle der Koeffizienten
confint(fit)

# Residuen / Devianz
residuals(fit, type="deviance")

# Vorhersagewerte
predict(fit, type="response")

32.2.1 lineare Regression

Für lineare Regressionen haben wir die Funktion lm() verwendet. Tatsächlich könnten wir die Modelle auch per gml() erstellen. Kommen wir zurück zu dem Beispiel des Lese- und Rechtschreibetests.

# 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"
summary(lm(Rechtschreibung~Lesetest))
## 
## Call:
## lm(formula = Rechtschreibung ~ Lesetest)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.42907 -0.26211  0.04498  0.36332  1.29412 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.20761    0.67774   1.782    0.113    
## Lesetest     1.05536    0.05718  18.458 7.65e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9221 on 8 degrees of freedom
## Multiple R-squared:  0.9771, Adjusted R-squared:  0.9742 
## F-statistic: 340.7 on 1 and 8 DF,  p-value: 7.648e-08

Der Aufruf in glm() erfolgt analog, wobei wir den Parameter family auf gaussian stellen müssen.

# fitte mit glm()
fit <- glm(Rechtschreibung~Lesetest, family = gaussian)
summary(fit)
## 
## Call:
## glm(formula = Rechtschreibung ~ Lesetest, family = gaussian)
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.20761    0.67774   1.782    0.113    
## Lesetest     1.05536    0.05718  18.458 7.65e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for gaussian family taken to be 0.850346)
## 
##     Null deviance: 296.5000  on 9  degrees of freedom
## Residual deviance:   6.8028  on 8  degrees of freedom
## AIC: 30.526
## 
## Number of Fisher Scoring iterations: 2

Die errechneten Parameter stimmen mit lm() überein. Es fällt auf, dass glm() keinen Wert für R2 zeigt, dafür aber einen AIC-Wert liefert.

Eine alternative Modellübersicht bietet die Funktion summ() (mit 2 m) aus dem Paket jtools.

jtools::summ(fit, digits=4)
Observations 10
Dependent variable Rechtschreibung
Type Linear regression
𝛘²(1) 289.6972
Pseudo-R² (Cragg-Uhler) 0.9790
Pseudo-R² (McFadden) 0.6062
AIC 30.5262
BIC 31.4340
Est. S.E. t val. p
(Intercept) 1.2076 0.6777 1.7818 0.1126
Lesetest 1.0554 0.0572 18.4576 0.0000
Standard errors: MLE

32.3 Poisson-Regression

Die Poisson-Regression wird beispielsweise verwendet, um die Anzahl an Ereignissen (Kunden in der Schlange, Patienten in der Notaufnahme) auf ihre Abhängikeit zu bestimmten Prädiktoren zu untersuchen.

Wir verwenden den Testdatensatz “Notaufnahme.txt”, der als Texttabelle im Internet verfügbar ist.

df <- read.table(url("https://www.produnis.de/R/data/Notaufnahme.txt"), 
                 header=TRUE, sep=",")

Er besteht aus 350 Beobachtungen der Patientenanzahl pro Tag in einer Notaufnahme. Zusätzlich wurden in der Spalte Wochentag die Wochentage erhoben. In Spalte Event wurde dichtotom erfasst, ob eine “größere” Veranstaltung in der Stadt stattgefunden hat.

# anzeigen
head(df)
##   Patienten Wochentag Event
## 1         8        Mo  nein
## 2         9        Di  nein
## 3        14        Mi  nein
## 4        10        Do  nein
## 5        10        Fr  nein
## 6        15        Mo  nein

Untersuchen wir nun, inwieweit der Wochentag und das Stattfinden einer größeren Veranstaltung die Anzahl der Notfälle beeinflusst. Das Modell kann wie gewohnt mit glm() erstellt werden, wobei wir den Parameter family=poisson setzen müssen.

# Poissonregressionsmodell erzeugen
fit <- glm(Patienten ~ Wochentag + Event, data = df, family = poisson)

# Zusammenfassung des Modells
summary(fit)
## 
## Call:
## glm(formula = Patienten ~ Wochentag + Event, family = poisson, 
##     data = df)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  2.528011   0.050672  49.889  < 2e-16 ***
## WochentagDo  0.071575   0.057063   1.254   0.2097    
## WochentagFr  0.050956   0.057349   0.889   0.3743    
## WochentagMi  0.119094   0.056421   2.111   0.0348 *  
## WochentagMo  0.003367   0.058026   0.058   0.9537    
## WochentagSa  0.429869   0.065146   6.599 4.15e-11 ***
## WochentagSo  0.365318   0.066489   5.494 3.92e-08 ***
## Eventnein   -0.292176   0.037640  -7.762 8.33e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 512.02  on 349  degrees of freedom
## Residual deviance: 345.57  on 342  degrees of freedom
## AIC: 1830.7
## 
## Number of Fisher Scoring iterations: 4

Das Modell zeigt an, dass die Patientenzahl signifikant an den Wochentagen “Sa” und “So” erhöht ist. Auch der Prädiktor Event ist signifikant. Der negative Prädiktorenwert bei “Eventnein” gibt an, dass sich die Patientenzahl verringert, wenn kein Event stattfindet. In der Modellzusammenfassung ist “Event=ja” der Bezugswert, so dass der Schätzwerte für “Event=nein” ausgegeben wird. Schöner anzusehen wäre es genau andersherum. Um “Event=nein” als Basiswert zu setzen, und den Schätzwert für “Event=ja” anzuzeigen, können wir die Events in einen Factor umwandeln, und dann mittels relevel() den Wert “nein” als Bezugspunkt setzen.

# Wir wollen Event="nein" als Basis
df$Event <- factor(df$Event)
df$Event <- relevel(df$Event, "nein")

# Poissonregressionsmodell neu erzeugen
fit <- glm(Patienten ~ Wochentag + Event, data = df, family = poisson)

# Zusammenfassung des Modells
summary(fit)
## 
## Call:
## glm(formula = Patienten ~ Wochentag + Event, family = poisson, 
##     data = df)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept) 2.235835   0.041828  53.453  < 2e-16 ***
## WochentagDo 0.071575   0.057063   1.254   0.2097    
## WochentagFr 0.050956   0.057349   0.889   0.3743    
## WochentagMi 0.119094   0.056421   2.111   0.0348 *  
## WochentagMo 0.003367   0.058026   0.058   0.9537    
## WochentagSa 0.429869   0.065146   6.599 4.15e-11 ***
## WochentagSo 0.365318   0.066489   5.494 3.92e-08 ***
## Eventja     0.292176   0.037640   7.762 8.33e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 512.02  on 349  degrees of freedom
## Residual deviance: 345.57  on 342  degrees of freedom
## AIC: 1830.7
## 
## Number of Fisher Scoring iterations: 4

Eine alternative Übersicht bietet die Funktion summ() (mit 2 m) aus dem Paket jtools.

jtools::summ(fit, digits=4)
Observations 350
Dependent variable Patienten
Type Generalized linear model
Family poisson
Link log
𝛘²(7) 166.4510
Pseudo-R² (Cragg-Uhler) 0.3798
Pseudo-R² (McFadden) 0.0840
AIC 1830.7077
BIC 1861.5712
Est. S.E. z val. p
(Intercept) 2.2358 0.0418 53.4531 0.0000
WochentagDo 0.0716 0.0571 1.2543 0.2097
WochentagFr 0.0510 0.0573 0.8885 0.3743
WochentagMi 0.1191 0.0564 2.1108 0.0348
WochentagMo 0.0034 0.0580 0.0580 0.9537
WochentagSa 0.4299 0.0651 6.5986 0.0000
WochentagSo 0.3653 0.0665 5.4944 0.0000
Eventja 0.2922 0.0376 7.7624 0.0000
Standard errors: MLE

Die Werte der Koeffizienten sind als Log-Odds dargestellt. Eine bessere Lese- und Interpretierbarkeit bietet die Darstellung als Odds-Ratio. Hierfür muss der Logarithmus mittels exp() aufgehoben werden.

# Zeige die Koeffizienten als Odds-Ratio
exp(fit$coefficients)
## (Intercept) WochentagDo WochentagFr WochentagMi WochentagMo WochentagSa 
##    9.354288    1.074199    1.052277    1.126476    1.003373    1.537056 
## WochentagSo     Eventja 
##    1.440972    1.339339

Die Konfidenzintervalle der Prädiktorenwerte erhalten wir mit der Funktion confint().

confint(fit)
##                    2.5 %    97.5 %
## (Intercept)  2.152783594 2.3167746
## WochentagDo -0.040220978 0.1835209
## WochentagFr -0.061421948 0.1634414
## WochentagMi  0.008603452 0.2298270
## WochentagMo -0.110388958 0.1171302
## WochentagSa  0.301707989 0.5571528
## WochentagSo  0.234302393 0.4950263
## Eventja      0.217989976 0.3655505

Bzw. mit vorgeschalteter exp()-Funktion, um die Werte als Odds-Ratio auszugeben.

# Konfidenzintervalle als Odds Ratio
exp(confint(fit))
##                 2.5 %    97.5 %
## (Intercept) 8.6087884 10.142906
## WochentagDo 0.9605771  1.201440
## WochentagFr 0.9404263  1.177556
## WochentagMi 1.0086406  1.258382
## WochentagMo 0.8954858  1.124266
## WochentagSa 1.3521663  1.745695
## WochentagSo 1.2640267  1.640541
## Eventja     1.2435746  1.441307

Die Funktion tab_model() aus dem sjPlot-Paket erstellt ihrerseits eine schöne Übersichtstabelle, wobei die Werte in Odds-Ratios angegeben und die Konfidenzintervalle enthalten sind.

sjPlot::tab_model(fit)
  Patienten
Predictors Incidence Rate Ratios CI p
(Intercept) 9.35 8.61 – 10.14 <0.001
Wochentag [Do] 1.07 0.96 – 1.20 0.210
Wochentag [Fr] 1.05 0.94 – 1.18 0.374
Wochentag [Mi] 1.13 1.01 – 1.26 0.035
Wochentag [Mo] 1.00 0.90 – 1.12 0.954
Wochentag [Sa] 1.54 1.35 – 1.75 <0.001
Wochentag [So] 1.44 1.26 – 1.64 <0.001
Event [ja] 1.34 1.24 – 1.44 <0.001
Observations 350
R2 Nagelkerke 0.493

Um die Güte des Modells zu bestimmen, können nun weitere Berechnungen vorgenommen werden.

32.3.1 Dispersion

Die Equidistation (Dispersion) des Modells sollte zwischen \(0,9\) und \(1,1\) liegen. Um sie zu bestimmen teilen wir die Residualdeviance durch die Anzahl der Feiheitsgrade.

# berechne Dispersions-Ratio
fit$deviance / fit$df.residual
## [1] 1.010433

Unser Wert liegt innerhalb des Bereichs, d.h. das Modell ist nicht überdispersioniert. Mit der Funktion dispersiontest() aus dem AER-Paket lässt sich zusätzlich ein Signifikanztest auf Überdispersion fahren. Die Nullhypothese besagt, das keine Überdispersion vorliegt.

AER::dispersiontest(fit, trafo=1)
## 
##  Overdispersion test
## 
## data:  fit
## z = -0.64112, p-value = 0.7393
## alternative hypothesis: true alpha is greater than 0
## sample estimates:
##       alpha 
## -0.04479655

Der Test ist nicht signifikant, es liegt keine Überdispersion vor.

Sollte in dem Modell Überdispersion vorliegen, erfolgt der Aufruf von gml() mit dem Parameter family=quasipoisson.

glm(Patienten ~ Wochentag + Event, data = df, family = quasipoisson)

32.3.2 Nullmodell

Mit einem Nullmodell können wir testen, ob unser gefittetes Modell die Daten besser erklären kann. Hierfür ersetzen wir alle Prädiktoren durch die Ziffer \(1\).

# erstelle Null-Modell
fit0 <- glm(Patienten ~ 1, data = df, family = poisson)

Vergleichen wir die AIC-Werte beider Modelle.

# Ziehe von unserem Modell-AIC den Nullmodell-AIC ab
fit$aic - fit0$aic
## [1] -152.451

Das Ergebnis ist negativ. Das beudetet, dass der AIC-Wert des Nullmodells größer ist. Somit ist unser AIC-Modellwert kleiner, und das bedeutet, dass es besser ist als das Nullmodell.

Zusätzlich können wir die beiden Modelle mittels Varianzanalyse (siehe Abschnitt 32.9) untersuchen, wobei wir als Testeinstellung den Chi2-Test wählen.

# ANOVA auf beide Modelle mittels Chi^2-Test
anova(fit0, fit, test="Chisq")
## Analysis of Deviance Table
## 
## Model 1: Patienten ~ 1
## Model 2: Patienten ~ Wochentag + Event
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1       349     512.02                          
## 2       342     345.57  7   166.45 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Unser gefittetes Modell erklärt die Daten signifikant besser als das Nullmodell.

32.3.3 Vorhersagen

Mit der Funktion predict() können wir nun die wahrscheinliche Patientenzahl für Kombinationen aus Wochentag und Event vorhersagen.

# Mit wievielen Patienten ist an einem
# Dienstag mit Event in der Stadt zu rechnen?
predict(fit, data.frame(Wochentag="Di", Event="ja"), type = "response")
##        1 
## 12.52856
# ohne Label "1" ausgeben
as.numeric(predict(fit, data.frame(Wochentag="Di", Event="ja"), 
                   type = "response"))
## [1] 12.52856

Um gut vorbereitet zu sein, runden wir die Patientenzahl mittels ceiling() auf die nächste ganze Zahl auf.

as.numeric(ceiling(predict(fit, data.frame(Wochentag="Di", Event="ja"), 
                           type = "response")))
## [1] 13

Wir können davon ausgehen, dass an einem Dienstag mit Event in der Stadt etwa 13 Patienten in die Notaufnahme kommen werden.

Die Wahrscheinlichkeiten für verschiedene Patientenanzahlen an beliebigen Wochentagen mit oder ohne Stadtevent lassen sich mit der Funktion ppois() bestimmen (siehe Abschnitt 18.4).

# Lambda-Wert für "Dienstag" und Event="ja"
Modelllambda <- predict(fit, data.frame(Wochentag="Di", Event="ja"), 
                        type = "response")

## Wahrscheinlichkeit, dass an einem Dienstag mit Event...
## - bis zu 16 Patienten eingeliefert werden (nicht mehr als 16)
ppois(16, Modelllambda)
## [1] 0.8674935
## - mindestens 16 Patienten eingeliefert werden (16 oder mehr)
1 - ppois(16, Modelllambda)
## [1] 0.1325065
## - exakt 16 Patienten eingeliefert werden
ppois(16, Modelllambda) - ppois(15, Modelllambda)
## [1] 0.0637844

So können wir einen Plot erstellen, auf dem die Wahrscheinlichkeiten angegeben sind, dass an einem Dienstag mit oder ohne Event in der Stadt mindestens x Patienten aufgenommen werden.

# Lambda für Event="nein"
Modelllambda_ohne <- predict(fit, data.frame(Wochentag="Di", Event="nein"), 
                             type = "response")

# Patienten N 1 bis 25
x <- 1:25

# plotten für Event="ja"
plot(x, 1-ppois(x, Modelllambda), 
     type="b", col="blue",
     main = "Wahrscheinlichkeiten für Mindestpatientenanzahl",
     sub="an einem Dienstag mit und ohne Event ",
     ylab="Wahrscheinlichkeit",
     xlab="Mindestanzahl Patienten",  xaxt = "n")
# X-Achsenticks in 2er-Schritten
axis(1, at=seq(0, 25, 2))
# füge Punktlinie für Event="nein" hinzu
lines(x, 1-ppois(x, Modelllambda_ohne),
      type="b", col="red")
# füge Hilfslinie hinzu
abline(h=.5, lty=5, col="grey65")
# füge Legendenbox hinzu
legend(x="topright", inset=0.005, 
       legend=c("ohne Event", "mit Event"),
       col=c("red", "blue"), lwd=6, bg="grey95")

Das geht ebenso für die Wahrscheinlichkeiten der exakten Patientenanzahl.

# Wahrscheinlichkeit für exakte Patienten-N mit Event="ja"
y1 <- ppois(x, Modelllambda) - ppois(x-1, Modelllambda)
# Wahrscheinlichkeit für exakte Patienten-N mit Event="nein"
y2 <- ppois(x, Modelllambda_ohne) - ppois(x-1, Modelllambda_ohne)

# plotte für Event="nein"
plot(x, y2, type="b",  col="red",
     main = "Wahrscheinlichkeiten für exakte Patientenanzahl",
     sub="an einem Dienstag mit und ohne Event ",
     ylab="Wahrscheinlichkeit",
     xlab="Anzahl Patienten",  xaxt = "n")
# X-Achsenticks in 2er-Schritten
axis(1, at=seq(0, 25, 2))
# Hilfslinie für Event="nein"
abline(v=9, lty=5, col="pink", lwd=2)
# Punktlinie für Event="ja" einzeichnen
lines(x, y1, type="b", col="blue")
# Hilfslinie für Event="ja"
abline(v=12, lty=5, col="paleturquoise", lwd=2)
# Legendenbox hinzufügen
legend(x="topright", inset=0.005, 
       legend=c("ohne Event", "mit Event"),
       col=c("red", "blue"), lwd=6, bg="grey95")

32.4 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)
## 
## 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"

32.5 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 & 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/data/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.

32.5.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

32.5.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
## 
## 
## Exponentiated coefficients:
##   Konflikt Zufriedenh 
##  0.5848225  3.1902597

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
## 
## 
## Exponentiated coefficients:
##   Konflikt Zufriedenh 
##  1.5719790  0.2847038

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
## 
## 
## Exponentiated coefficients:
##   Konflikt:1   Konflikt:2   Konflikt:3 Zufriedenh:1 Zufriedenh:2 Zufriedenh:3 
##    1.1577337    1.6264285    1.8791766    0.1651675    0.2836103    0.4338519
# 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

32.6 Konfidenzintervalle

Konfidenzintervalle folgen der Logik, dass zum Punktschätzwert ein gewisser Fehlerbereich addiert und subtrahiert wird. Hierdurch wird ein Intervall mit Unter- und Obergrenze aufgezogen, das mit einer bestimmten Irrtumswahrscheinlichkeit den “wahren Wert” der Grundgesamtheit enthält.

32.6.1 Normalverteilung

Um ein Konfidenzintervall für einen Mittelwert aus einer normalverteilten Wertemenge zu berechnen wird folgende Formel angewendet:

\[\begin{equation} \label{formel:konfidenznormal} \begin{array}{l l} \mu_{u}\\\mu_{o}\\ \end{array} \Big\}= \bar{x} \pm t \cdot \frac{s}{\sqrt{n}} \end{equation}\]

  • \(\bar{x}\) = arithmetisches Mittel

  • \(n\) = Stichprobenumfang

  • \(s\) = Standardabweichung

  • \(t\) = Wert der \(t\)-Verteilung, (mit Freiheitsgraden \(df = n-1\))

Nehmen wir an, eine Zufalls-Erhebung über die jährlichen Ausgaben für Pflegehilfsmittel ergab bei 100 Befragten (n = 100) einen Mittelwert \(\bar{x}\)=400,-€ bei einer Standardabweichung von s=20,-€. Das Konfidenzintervall für \(\alpha\)=0.05 wird nun wie folgt bestimmt.

# schreibe gegebene Werte in Variablen
n     <- 100
xquer <- 400
s     <- 20

# 95%KI-Fehler berechnen mit (alpha/2) 
fehler95 <- qt(0.975, df=n-1)*s/sqrt(n)

# Fehlerwert ausgeben
fehler95
## [1] 3.968434

Beachten Sie, dass in der Funktion qt() zur Bestimmung des kritischen t-Wertes die Irrtumswahrscheinlichkeit \(\alpha\) halbiert werden muss, da es sich um eine zweiseitige Fragestellung handelt. So muss qt() nicht mit dem Wert 0.95 (\(\alpha=0.05\)) aufgerufen werden, sondern mit 0.975 (\(\frac{\alpha}{2}=0.025\)).

# untere Grenze des 95%KI
xquer - fehler95
## [1] 396.0316

Die untere Grenze des 95%-Konfidenzintervalls liegt bei 396.0316 Euro.

# obere Grenze des 95%KI
xquer + fehler95
## [1] 403.9684

Die obere Grenze des 95%-Konfidenzintervalls liegt bei 403.9684 Euro.

# Länge des 95%KI
(xquer + fehler95) - (xquer - fehler95)
## [1] 7.936868
# oder schlauer
2*fehler95
## [1] 7.936868

Das 95%-Konfidenzintervall hat eine Länge von 7.936868 Euro.

Für das 99% Konfidenzintervall ändern wir entsprechend um. Auch hier muss, da es sich um eine zweiseitige Fragestellung handelt, mit \(\frac{\alpha}{2}\), also 0.005, gerechnet werden. In der Funktion qt() ändert sich also der Wert auf 0,995.

# 99%KI-Fehler berechnen mit (alpha/2) 
fehler99 <- qt(0.995, df=n-1)*s/sqrt(n)

# zeige fehler99 an
fehler99
## [1] 5.252811
# untere Grenze des 99%KI
xquer - fehler99
## [1] 394.7472

Die untere Grenze des 99%-Konfidenzintervalls liegt bei 394.7472 Euro.

# obere Grenze des 99%KI
xquer + fehler99
## [1] 405.2528

Die obere Grenze des 99%-Konfidenzintervalls liegt bei 405.2528 Euro.

# Länge des 99%KI
(xquer + fehler99) - (xquer - fehler99)
## [1] 10.50562
# oder schlauer
2*fehler99
## [1] 10.50562

Das 99%-Konfidenzintervall hat eine Länge von 10.50562 Euro.

Die Funktion KInormal_a() aus dem jgsbook-Paket nimmt uns die vorgestellten Rechenschritte ab.

# Berechne das 95%-KI für
# n=100, xquer=400 und s=20
jgsbook::KInormal_a(xquer=400, n=100, s=20, alpha=.05)
## Registered S3 method overwritten by 'statip':
##   method         from      
##   predict.kmeans parameters
##                       x          y
## 1 0.95 KI untere Grenze 396.031566
## 2  0.95 KI obere Grenze 403.968434
## 3        0.95 KI Laenge   7.936868
# Berechne das 99%-KI für
# n=100, xquer=400 und s=20
jgsbook::KInormal_a(xquer=400, n=100, s=20, alpha=.01)
##                       x         y
## 1 0.99 KI untere Grenze 394.74719
## 2  0.99 KI obere Grenze 405.25281
## 3        0.99 KI Laenge  10.50562

32.6.1.1 Unterschied von Mittelwerten

Das Konfidenzintervall für den Unterschied von Mittelwerten normalverteilter Daten berechnet sich mit der Formel

\[\begin{equation} \begin{array}{l l} \mu_{u}\\\mu_{o}\\ \end{array} \Big\}= \hat{d} \pm t \cdot \sqrt{\frac{s_{pool}^2}{n_x}+\frac{s_{pool}^2}{n_y}} \end{equation}\]

  • \(\hat{d}\) = Differenz der Mittelwerte \(\bar{x}\),\(\bar{y}\)

  • \(n_{x}, n_{y}\) = die Umfänge der beiden Stichproben

  • \(s_{pool}\) = gepoolte Standardabweichungen = \(\frac{(n_x-1)\cdot s_x^2\ +\ (n_y-1)\cdot s_y^2}{n_x + n_y -2}\)

  • \(t\) = Wert der \(t\)-Verteilung (mit Freiheitsgraden \(df = n_x+n_y -1\))

Nehmen wir an, in einer Studie wurde die Wirkung von zwei Medikamenten auf den Magnesiumgehalt des Blutserums (in mg/dl) untersucht. 13 Personen erhielten Medikament A und 10 Personen Medikament B. Nach 6 Wochen wurden folgende Daten bestimmt:

-\(\bar{x}_{A}\) = 2,22 mg/dl

  • \(\bar{x}_{B}\) = 2,70 mg/dl

  • \(s_{A}\) = 0,255 mg/dl

  • \(s_{B}\) = 0,306 mg/dl

  • \(n_{A}\) = 13

  • \(n_{B}\) = 10

Das Konfidenzintervall für den Unterschied der Mittelwerte bei \(\alpha=0,05\) wird nun wie folgt bestimmt (Achtung, da es sich um eine zweiseitige Fragestellung handelt, müssen wir mit \(\frac{\alpha}{2} =0,025\) rechnen).

# erzeuge die Daten
xa <- 2.22
xb <- 2.7
sa <- 0.255
sb <- 0.306
na <- 13
nb <- 10

# Differenz der Mittelwerte
d   <- xb - xa

# ausgeben
d
## [1] 0.48
# berechne die gepoolten Standardabweichungen
s_pool = ((na-1)*sa^2 + (nb-1)*sb^2) / (na+nb-2)

# ausgeben
s_pool
## [1] 0.07728686
#Fehlerquote berechnen
fehler95 <- qt(0.975, df=na+nb-1) * sqrt(s_pool^2/na + s_pool^2/nb)

# ausgeben
fehler95
## [1] 0.06741865
# untere Grenze des 95%KI
d - fehler95
## [1] 0.4125813

Die untere Grenze des 95%-Konfidenzintervalls liegt bei 0.4125813 mg/dl.

# obere Grenze des 95%KI
d + fehler95
## [1] 0.5474187

Die obere Grenze des 95%-Konfidenzintervalls liegt bei 0.5474187 mg/dl.

# Länge des Intervalls
2*fehler95
## [1] 0.1348373

Das 95%-Konfidenzintervall hat eine Länge von 0.1348373 mg/dl.

Für das 99%-Konfidenzintervall wird der Wert für die Funktion qt() auf 0.995 gesetzt (da es sich um eine zweiseitige Fragestellung handelt, müssen wir mit \(\frac{\alpha}{2} = 0,005\) rechnen).

#Fehlerquote berechnen
fehler99 <- qt(0.995, df=na+nb-1) * sqrt(s_pool^2/na + s_pool^2/nb)

# ausgeben
fehler99
## [1] 0.09163373
# untere Grenze des 99%KI
d - fehler99
## [1] 0.3883663

Die untere Grenze des 99%-Konfidenzintervalls liegt bei 0.3883663 mg/dl.

# obere Grenze des 99%KI
d + fehler99
## [1] 0.5716337

Die obere Grenze des 99%-Konfidenzintervalls liegt bei 0.5716337 mg/dl.

# Länge des Intervalls
2*fehler99
## [1] 0.1832675

Das 99%-Konfidenzintervall hat eine Länge von 0.1832675 mg/dl.

Die Funktion KInormal_u() aus dem jgsbook-Paket nimmt uns auch hier die vorgestellten Rechenschritte ab.

# Berechne das 95%-KI
jgsbook::KInormal_u(xa, sa, na, xb, sb, nb, alpha=0.05)
##                           x         y
## 1 Differenz der Mittelwerte 0.4800000
## 2     0.95 KI untere Grenze 0.4125813
## 3      0.95 KI obere Grenze 0.5474187
## 4            0.95 KI Laenge 0.1348373
# Berechne das 99%-KI
jgsbook::KInormal_u(xa, sa, na, xb, sb, nb, alpha=0.01)
##                           x         y
## 1 Differenz der Mittelwerte 0.4800000
## 2     0.99 KI untere Grenze 0.3883663
## 3      0.99 KI obere Grenze 0.5716337
## 4            0.99 KI Laenge 0.1832675

32.6.2 Binomialverteilung

Zur Berechnung von Konfidenzintervallen für Anteilswerte wird folgende Formel verwendet.

\[\begin{equation} \label{formel:konfidenzbinomial} \begin{array}{l l} p_{o}\\p_{u}\\ \end{array} \Big\}= k \pm c \cdot \sqrt{\frac{k \cdot (1-k)}{n}} \end{equation}\]

  • \(n\) = Umfang der Stichprobe

  • \(k\) = geschätzter Anteil

  • \(c\) = Quartilswert für Irrtumswahrscheinlichkeit \(\alpha\)

Nehmen wir an, in einer klinischen Studie wurde an 150 Patienten die Dekubitusprävalenz ermittelt. Dabei wurde eine Prävalenz von 35% festgestellt. Die Konfidenzintervalle für den “wahren Wert” der Grundgesamtheit wird wie folgt bestimmt.

# Daten erzeugen
n <- 150
k <- 0.35

# Fehler berechnen mit (alpha/2)!
fehler95 <- qnorm(0.975)*sqrt(k*(1-k)/n)

# anzeigen
fehler95
## [1] 0.07632963

Beachten Sie, dass in der Funktion qnorm() die Irrtumswahrscheinlichkeit \(\alpha\) halbiert werden muss, da es sich um eine zweiseitige Fragestellung handelt. So muss qnorm() nicht mit dem Wert 0.95 (\(\alpha=0.05\)) aufgerufen werden, sondern mit 0.975 (\(\frac{\alpha}{2}=0.025\)).

# untere Grenze des KI
k - fehler95
## [1] 0.2736704

Die untere Grenze des 95%-Konfidenzintervalls liegt bei 0.2736704 = 27,36704%.

# obere Grenze des KI
k + fehler95
## [1] 0.4263296

Die obere Grenze des 95%-Konfidenzintervalls liegt bei 0.4263296 = 42,63296%.

# Länge des KI
2*fehler95
## [1] 0.1526593

Das 95%-Konfidenzintervall hat eine Länge von 0.1526593 = 15,26593%.

Für das 99%-Konfidenzintervall muss der Wert der Funktion qnorm() auf 0.995 (\(\frac{\alpha}{2}=0,005\)) angepasst werden.

# Fehler berechnen mit (alpha/2)!
fehler99 <- qnorm(0.995)*sqrt(k*(1-k)/n)

# anzeigen
fehler99
## [1] 0.1003141
# untere Grenze des KI
k - fehler99
## [1] 0.2496859

Die untere Grenze des 99%-Konfidenzintervalls liegt bei 0.2496859 = 24,96859%.

# obere Grenze des KI
k + fehler99
## [1] 0.4503141

Die obere Grenze des 99%-Konfidenzintervalls liegt bei 0.4503141 = 45,03141%.

# Länge des KI
2*fehler99
## [1] 0.2006283

Das 99%-Konfidenzintervall hat eine Länge von 0.2006283 = 20,06283%.

Die Funktion KIbinomial_a() aus dem jgsbook-Paket nimmt uns auch hier die vorgestellten Rechenschritte ab.

# Berechne das 95%-KI
jgsbook::KIbinomial_a(p=0.35, n=150, alpha=0.05)
##                       x         y
## 1 0.95 KI untere Grenze 0.2736704
## 2  0.95 KI obere Grenze 0.4263296
## 3        0.95 KI Laenge 0.1526593
# Berechne das 99%-KI
jgsbook::KIbinomial_a(p=0.35, n=150, alpha=0.01)
##                       x         y
## 1 0.99 KI untere Grenze 0.2496859
## 2  0.99 KI obere Grenze 0.4503141
## 3        0.99 KI Laenge 0.2006283

32.6.2.1 Unterschied von Anteilswerten

Zur Berechnung der Konfidenzintervalle bei Unterschieden von Anteils- werten wird die folgende Formel verwendet.

\[\begin{equation} \label{formel:konfidenzunterschiedbinomial} \begin{array}{l l} d_{o}\\d_{u}\\ \end{array} \Big\}= \hat{d} \pm c \cdot \sqrt{\frac{\hat{p}_{1} \cdot (1 - \hat{p}_{1})}{n_{1}} + \frac{\hat{p}_{2} \cdot (1 - \hat{p}_{2})}{n_{2}}} \end{equation}\]

  • \(\hat{d}\) = die Differenz der Anteile

  • \(n_{1}, n_{2}\) = die Umfänge der beiden Stichproben

  • \(\hat{p}_{1}\), \(\hat{p}_{2}\) = die geschätzten Anteile in beiden Stichproben

  • \(c\) = Quartilswert für Irrtumswahrscheinlichkeit \(\alpha\)

Nehmen wir an, in Bundesland A gaben 25 von 100 Befragten an, Zuzahlungen bei Pflegehilfsmitteln leisten zu müssen. In Bundesland B wurden 150 Personen befragt, von denen 60 von Zuzahlungen betroffen waren. Bundesland A hat

dementsprechend einen Anteil von 0,25 (25%) und Bundesland B einen Anteil von 0,40 (40%) Zuzahlern. Der Punktschätzwert für die Differenz beträgt folglich 0,40-0,25 = 0,15.

Die Konfidenzintervalle für die wahre Differenz d kann nun wie folgt bestimmt werden.

# Daten erzeugen
n1 <- 100
p1 <- 0.25
  
n2 <- 150
p2 <- 0.4

# Differenz der Anteile
d <- p2-p1

Beachten Sie, dass in der Funktion qnorm() die Irrtumswahrscheinlichkeit \(\alpha\) halbiert werden muss, da es sich um eine zweiseitige Fragestellung handelt. So muss qnorm() nicht mit dem Wert 0.95 (\(\alpha=0.05\)) aufgerufen werden, sondern mit 0.975 (\(\frac{\alpha}{2}=0.025\)).

# 95%KI-Fehler berechnen mit (alpha/2) 
fehler95 <- qnorm(0.975)*sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2)

#anzeigen
fehler95
## [1] 0.1155382
# untere Grenze des 95%KI
d - fehler95
## [1] 0.03446183

Die untere Grenze des 95%-Konfidenzintervalls liegt bei 0.03446183 = 3,446183%.

# obere Grenze des 95%KI
d + fehler95
## [1] 0.2655382

Die obere Grenze des 95%-Konfidenzintervalls liegt bei 0.2655382 = 26,55382%.

# Länge des KI
2*fehler95
## [1] 0.2310763

Die Länge des 95%-Konfidenzintervalls liegt bei 0.2310763 = 23,10763%.

Für das 99%-Konfidenzintervall muss der Wert der Funktion qnorm() auf 0.995 (\(\frac{\alpha}{2}=0,005\)) angepasst werden.

# 99%KI-Fehler berechnen mit (alpha/2) 
fehler99 <- qnorm(0.995)*sqrt(p1*(1-p1)/n1 + p2*(1-p2)/n2)

#anzeigen
fehler99
## [1] 0.1518429
# untere Grenze des 99%KI
d - fehler99
## [1] -0.001842898

Die untere Grenze des 99%-Konfidenzintervalls liegt bei -0.0018429 = -0,18429%.

# obere Grenze des 99%KI
d + fehler99
## [1] 0.3018429

Die obere Grenze des 99%-Konfidenzintervalls liegt bei 0.3018429 = 30,18429%.

# Länge des KI
2*fehler99
## [1] 0.3036858

Die Länge des 99%-Konfidenzintervalls liegt bei 0.3036858 = 30,36858%.

Wie Sie sehen, liegt die untere Grenze des 99%-Konfidenzintervalls im negativen Bereich. Errechnet sich ein Intervall, welches die 0 einschließt, bedeutet dies, dass der wahre Unterschied auch 0 sein könnte.

Die Funktion KIbinomail_u() aus dem jgsbook-Paket nimmt uns auch hier die vorgestellten Rechenschritte ab.

# Berechne das 95%-KI
jgsbook::KIbinomial_u(p1, n1, p2, n2, alpha=0.05)
##                       x          y
## 1 Differenz der Anteile 0.15000000
## 2 0.95 KI untere Grenze 0.03446183
## 3  0.95 KI obere Grenze 0.26553817
## 4        0.95 KI Laenge 0.23107635
# Berechne das 99%-KI
jgsbook::KIbinomial_u(p1, n1, p2, n2, alpha=0.01)
##                       x            y
## 1 Differenz der Anteile  0.150000000
## 2 0.99 KI untere Grenze -0.001842898
## 3  0.99 KI obere Grenze  0.301842898
## 4        0.99 KI Laenge  0.303685796

32.7 Signifikanztests

32.7.1 \(\chi^2\)-Test

Der \(\chi^2\)-Test wird mit der Funktion chisq.test() durchgeführt. Hierfür müssen die Daten als Kreuztabelle vorliegen.

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

# erstelle Kreuztabelle aus Alter und Geschlecht
xtabelle <- xtabs(~ nw$age + nw$sex)

# Führe Chiquadrat-Test durch
chisq.test(xtabelle)
## Warning in chisq.test(xtabelle): Chi-Quadrat-Approximation kann inkorrekt sein
## 
##  Pearson's Chi-squared test
## 
## data:  xtabelle
## X-squared = 43.886, df = 44, p-value = 0.4765

Der p-Wert ist nicht signifikant, das heisst, das kein Zusammenhang gefunden werden konnte.

Die Warnung bedeutet, dass der p-Wert aus sehr kleinen Zellenhäufigkeiten (kleiner 5) geschätzt wurde, und daher fehlerhaft sein kann (wir sollten daher einen exakten Test durchführen, zum Beispiel den Fisher-Test).

Ein gutes Anleitungsvideo zum \(\chi^2\)-Test findet sich bei Youtube: https://www.youtube.com/watch?v=YJUuyaC0x48

32.7.2 Fisher’s Exakttest

Der exakte Fisher-Test wird mit der Funktion fisher.test() durchgeführt.

# Führe Fisher-Test durch
fisher.test(xtabelle)
Fisher's Exact Test for Count Data with simulated p-value (based on 2000 replicates)
data: nw$sex and nw$age
p-value = 0.3848
alternative hypothesis: two.sided


Der exakte Fisher-Test bestätigt, dass es keinen Unterschied gibt.

32.7.3 t-Test

Der t-Test ist wohl einer der bekanntesten Signifikanztests. In R wird er mit der Funktion t.test() ausgeführt.

32.7.3.1 Einstichproben t-Test

Beim Einstichproben t-Test wird der Mittelwert der Grundgesamtheit über den Parameter mu mitgegeben. Schauen wir, ob im Nachtwachendatensatz das Alter von einem gedachten Mittelwert 40 in der Grundgesamtheit abweicht.

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

# Führe t-Test für "Alter" auf Mittelwert 40 durch
t.test(nw$age, mu=40)
## 
##  One Sample t-test
## 
## data:  nw$age
## t = 6.8274, df = 275, p-value = 5.524e-11
## alternative hypothesis: true mean is not equal to 40
## 95 percent confidence interval:
##  43.00908 45.44744
## sample estimates:
## mean of x 
##  44.22826

Der Test ist signifikant, das bedeutet, das Alter im Datensatz nw weicht von einem Mittelwert in der Grundgesamtheit von 40 ab. Es liegt ein Unterschied vor.

Zur Beschreibung der Effektstärke des Unterschieds können wir Cohen’s d berechnen. Dies kann über das Zusatzpaket lsr und der Funktion cohensD() erfolgen.

# aktiviere Zusatzpaket "lsr"
library(lsr)

# berechne Effektsärke Cohen's D
cohensD(nw$age, mu=40)
## [1] 0.4109627

Dieser Wert kann nun nach Cohen’s Angaben (1992) interpretiert werden (0.2=kleiner Effekt, 0.5=mittel, 0.8=stark).

Unsere Effektgröße ist mit 0.4109627 eher “mittel”.

Ein gutes Anleitungsvideo zum Einstichproben t-Test findet sich bei Youtube: https://www.youtube.com/watch?v=7vcyjuSwfzI

32.7.4 unabhängiger t-Test

Für den Vergleich der Mittelwerte bei unabhängige Stichproben werden die Werte beider Gruppen als Variable übergeben. So können wir z.B. schauen, ob sich die Datensätze nw und epa hinsichtlich des Alters unterscheiden

# t-test auf Alter in Datensatz "epa" und "nw"
t.test(nw$age, epa$age)
## 
##  Welch Two Sample t-test
## 
## data:  nw$age and epa$age
## t = -19.316, df = 854.37, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -20.65610 -16.84545
## sample estimates:
## mean of x mean of y 
##  44.22826  62.97903

Das Ergebnis ist signifikant, die Gruppen unterscheiden sich also.

Standardmäßig wird ein zweiseitiger Test durchgeführt. Über den Parameter alternative kann dies geändert werden.

# t-test einseitig, kleiner
t.test(nw$age, epa$age, alternative = "less")
# t-test einseitig, größer
t.test(nw$age, epa$age, alternative = "greater")
# t-test zweiseitig
t.test(nw$age, epa$age, alternative = "two.sided")

Die Funktion nimmt auch die Syntax erklärt durch entgegen. Wir schauen, ob sich das Alter im Datensatz epa in den Geschlechtergruppen unterscheidet, also Alter erklärt durch Geschlecht.

# t-test Alter zwischen Geschlechtern im Datensatz "epa"
t.test(age ~ sex, data=epa)
## 
##  Welch Two Sample t-test
## 
## data:  age by sex
## t = -4.7535, df = 595.03, p-value = 2.51e-06
## alternative hypothesis: true difference in means between group m and group w is not equal to 0
## 95 percent confidence interval:
##  -9.923268 -4.120822
## sample estimates:
## mean in group m mean in group w 
##        59.79646        66.81851

Das Ergebnis ist signifikant. Auch hier lässt sich wieder die Effektstärke Cohen’s D berechnen.

# Effektsärke 
lsr::cohensD(age ~ sex, data=epa)
## [1] 0.3837919

Der Effekt ist mit 0.3837919 eher klein.

Ein gutes Anleitungsvideo zum unabhängigen t-Test findet sich bei Youtube: https://www.youtube.com/watch?v=3Yz-XNRasGM

32.7.4.1 abhängiger t-Test

Bei gepaarten Werten muss ein abhängiger t-Test berechnet werden. Dies erfolgt über den Parameter paired, der auf TRUE gesetzt wird. Auch hier können wir wieder die Richtung (größer, kleiner, zweiseitig) über den Parameter alternative verändern.

# erzeuge Testwerte
t0 <- rnorm(100, mean=12, sd=1)
t1 <- rnorm(100, mean=11, sd=3)

# gepaarter t-Test
t.test(t0, t1, paired=TRUE)
## 
##  Paired t-test
## 
## data:  t0 and t1
## t = 4.5818, df = 99, p-value = 1.345e-05
## alternative hypothesis: true mean difference is not equal to 0
## 95 percent confidence interval:
##  0.8011115 2.0249811
## sample estimates:
## mean difference 
##        1.413046
# gepaarter  t-test einseitig, kleiner
t.test(t0, t1, paired=TRUE, alternative = "less")
# gepaarter t-test einseitig, größer
t.test(t0, t1, paired=TRUE, alternative = "greater")
# gepaarter t-test zweiseitig
t.test(t0, t1, paired=TRUE, alternative = "two.sided")

Die Effektstärke Cohen’s D lässt sich entsprechend berechnen, indem der Parameter method auf paired gesetzt wird:

# Effektsärke 
lsr::cohensD(t0, t1, method="paired")
## [1] 0.4581845

Der Effekt ist mit 0.4581845 eher gering.

Ein gutes Anleitungsvideo zum verbundenen t-Test findet sich bei Youtube: https://www.youtube.com/watch?v=H15M_8gh1Ok

32.7.5 F-Test

Der F-Test wird mit der Funktion var.test() berechnet. So können wir prüfen, ob die Varianz des Alters im Datensatz epa in den Geschlechtergruppen unterschiedlich ist.

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

# F-Test 
var.test(age ~ sex , data=epa)
## 
##  F test to compare two variances
## 
## data:  age by sex
## F = 0.98352, num df = 338, denom df = 280, p-value = 0.8816
## alternative hypothesis: true ratio of variances is not equal to 1
## 95 percent confidence interval:
##  0.7848842 1.2294811
## sample estimates:
## ratio of variances 
##           0.983521

Der Test ist nicht signifikant, es liegt kein Unterschied in den Varianzen vor.

32.7.6 Levene-Test

Ebenso kann mit dem Levene-Test auf Varianzunterschied geprüft werden. Hierfür rufen wir die Funktion leveneTest() aus dem Zusatzpaket car auf.

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

# Levene-Test 
car::leveneTest(epa$age, epa$sex)
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1  0.0232 0.8789
##       618

Der Test ist nicht signifikant, es liegt kein Unterschied in den Varianzen vor. Ein gutes Anleitungsvideo zum verbundenen Levene-Test findet sich bei Youtube: https://www.youtube.com/watch?v=717UWGnyQGA

32.7.7 Mann-Whitney-U-Test

Dies ist die nicht-parametrische alternative zum t-Test für unabhängige Stichproben. Er wird mit der Funktion wilcox.test() durchgeführt, wobei der Parameter paired auf FALSE gesetzt wird (ist die Standardeinstellung).

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

# Man-Whitney-U-Test
wilcox.test(epa$age ~ epa$sex)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  epa$age by epa$sex
## W = 35862, p-value = 1.151e-07
## alternative hypothesis: true location shift is not equal to 0

Das Ergebnis ist signifikant, das Alter im Datensatz epa unterscheidet sich zwischen den Geschlechtern.

Mit dem Parameter conf.int kann das Konfidenzintervall des Schätzwertes mit ausgegeben werden.

# Man-Whitney-U-Test mit Konfidenzintervall
wilcox.test(epa$age ~ epa$sex, conf.int=TRUE)
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  epa$age by epa$sex
## W = 35862, p-value = 1.151e-07
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -9.999986 -4.999952
## sample estimates:
## difference in location 
##              -7.000029

Ein gutes Anleitungsvideo zum Mann-Whitney-U-Test finden Sei bei Youtube: https://www.youtube.com/watch?v=4lFyRcXoJB8

32.7.8 Wilcoxon-Test

Dies ist die nicht-parametrische alternative zum abhängigen t-Test. Er wird mit der Funktion wilcox.test() durchfegührt, wobei der Parameter paired auf TRUE gesetzt werden muss. Andernfalls rechnet R einen Man-Whitney-U-Test.

# erzeuge Testwerte
t0 <- sample(70:140, 100, replace=T)
t1 <- sample(74:170, 100, replace=T)

wilcox.test(t0, t1, paired=TRUE)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  t0 and t1
## V = 805, p-value = 5.628e-09
## alternative hypothesis: true location shift is not equal to 0

Mit dem Parameter conf.int kann das Konfidenzintervall des Schätzwertes mit ausgegeben werden.

# erzeuge Testwerte
t0 <- sample(70:140, 100, replace=T)
t1 <- sample(74:170, 100, replace=T)

wilcox.test(t0, t1, paired=TRUE, conf.int=TRUE)
## 
##  Wilcoxon signed rank test with continuity correction
## 
## data:  t0 and t1
## V = 1205.5, p-value = 5.748e-06
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -28.00003 -12.50000
## sample estimates:
## (pseudo)median 
##      -20.49995

Ein gutes Anleitungsvideo zum Wilcoxon-Test finden Sie bei Youtube: https://www.youtube.com/watch?v=4lFyRcXoJB8

32.7.9 Shapiro-Wilk-Normalverteilungstest

Zum Test auf Normalverteilung kann der Shapiro-Wilk-Test über die Funktion shapiro.test() verwendet werden. Die Nullhypothese lautet, dass eine Normalverteilung vorliegt. Testen wir, ob im Datensatz nw das Alter (age) normalverteilt ist.

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

# Teste "age" auf Normalverteilung
shapiro.test(nw$age)
## 
##  Shapiro-Wilk normality test
## 
## data:  nw$age
## W = 0.96636, p-value = 4.677e-06

Das Ergebnis ist signifikant, die Daten sind nicht normalverteilt.

Erzeugen wir normalverteilte Werte zur Gegenprobe.

# Test mit 100 zufälligen und normalverteilten Werten
shapiro.test(rnorm(100))
## 
##  Shapiro-Wilk normality test
## 
## data:  rnorm(100)
## W = 0.99307, p-value = 0.8921

Das Ergebnis ist (wie erwartet) nicht signifikant, das heisst, unsere Zufallszahlen sind wirklich normalverteilt.

Der Shapiro-Wilk-Test funktioniert nur bis zu einer Stichprobengröße von \(n=5000\). Besteht der Datensatz aus mehr Fällen, kann auf den Kolmogorov-Smirnov-Test ausgewichen werden.

32.7.10 Kolmogorov-Smirnov-Test

Der Kolmogorov-Smirnov-Test prüft die Nullhypothese, dass die Daten einer bestimmten Verteilung (z.B. der Normalverteilung) folgen. In R ist er über die Funktion ks.test() implementiert.

Prüfen wir nun also mittels Kolmogorov-Smirnov-Test, ob im Datensatz nw das Alter (age) normalverteilt ist.

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

# Teste, ob Alter normalverteilt ist
ks.test(nw$age, "pnorm")
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  nw$age
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided

Der Test ist signifikant, das heisst, es liegt keine Normalverteilung vor.

Die Testparameter können zudem auf die Lage- und Streuungskenngrößen der Wertereihe angepasst werden.

# passe auf Kenngrößen an
ks.test(nw$age, "pnorm", 
        mean = mean(nw$age, na.rm = TRUE), 
        sd = sd(nw$age, na.rm = TRUE))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  nw$age
## D = 0.11477, p-value = 0.001391
## alternative hypothesis: two-sided

Auch dieser Test ist signifikant, das heisst, es liegt keine Normalverteilung vor.

Machen wir nun den Gegentest mit normalverteilten Werten:

# erzeuge normalverteilte Dummywerte
dummy <- rnorm(10000)

# Test auf Normalverteilung
ks.test(dummy, "pnorm", 
        mean = mean(dummy), 
        sd = sd(dummy))
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dummy
## D = 0.0059544, p-value = 0.8703
## alternative hypothesis: two-sided

Der Test ist nicht signifikant, das heisst, es liegen normalverteilte Werte vor.

32.8 Post-Hoc-Tests

Post-hoc-Tests werden verwendet, um paarweise Vergleiche nach einem initialen Test (wie ANOVA (siehe Abschnitt 32.9)) durchzuführen, der signifikante Unterschiede zwischen Gruppenmittelwerten anzeigt. Durch die Mehrfachvergleiche (multiples Testen) entsteht das erhöhte Risiko von falsch-positiven Ergebnissen (Typ-I-Fehler).

Post-Hoc-Tests helfen dabei, spezifische Gruppen zu identifizieren, die sich voneinander unterscheiden, während sie das Problem der Mehrfachtestung berücksichtigen.

R bietet mehrere Funktionen zur Durchführung von Post-hoc-Tests. In diesem Abschnitt bschränken wir uns auf die weit verbreiteten Funktionen reporttools::pairwise.fisher.test(), pairwise.wilcox.test() und pairwise.t.test().

32.8.1 Paarweise Chiquadrat

In R-base ist keine Funktion implementiert, die einen paarweisen Chiquadrat-Test durchführt. Das mag daran liegen, dass der Chiquadrat-Test als (weniger genaue) Alternative zu dem rechenaufwändigen exakten Fisher-Test entwickelt wurde. Mit der heutigen Computertechnologie ist es jedoch in den meisten Fällen kein Problem, den exakten Test nach Fisher auch auf älterer Hardware anzuwenden.

32.8.2 Paarweise Fisher-Test

Der pairwise.fisher.test aus dem Paket reporttools wird verwendet, um paarweise Vergleiche zwischen Gruppen durchzuführen, wenn nonimale Daten vorliegen.

Als Beispiel untersuchen wir im Datensatz pf8, ob sich die Geschlechter je nach Standort unterscheiden.

# verwende pf8 Datensatz
pf8 <- jgsbook::pf8

# führe exakten Test durch
reporttools::pairwise.fisher.test(pf8$Geschlecht, pf8$Standort, p.adjust.method = "holm")
## 
##  Pairwise comparisons using 
## 
## data:  pf8$Geschlecht and pf8$Standort 
## 
##           Rheine  Münster Bahn    Ladbergen
## Münster   1.00000 -       -       -        
## Bahn      1.00000 1.00000 -       -        
## Ladbergen 0.56615 0.24097 0.32261 -        
## Internet  0.00262 0.02282 0.09245 0.00065  
## 
## P value adjustment method: holm

Als Ergebnis erhalten wir eine Tabelle mit p-Werten. Es ist zu erkennen, dass sich die Geschlechterverteilungen zwischen den Standortvergleichen Internet-Rheine, Internet-Münster und Internet-Landbergen unterscheiden.

32.8.3 Paarweise Wilcoxon-Test (Mann-Whitney-U)

Die Funktion pairwise.wilcox.test() wird verwendet, um paarweise Vergleiche zwischen Gruppen durchzuführen, wenn die Daten nicht normalverteilt sind. Dieser Test basiert auf dem Wilcoxon-Rangsummentest, auch bekannt als Mann-Whitney-U-Test, und ist eine nichtparametrische Methode, die keine Annahmen über die Verteilung der Daten erfordert.

Als Beispiel untersuchen wir im Datensatz pf8, ob sich das Alter der Probanden je nach Standort unterscheidet.

pairwise.wilcox.test(pf8$Alter, pf8$Standort, p.adjust.method = "bonferroni")
## 
##  Pairwise comparisons using Wilcoxon rank sum test with continuity correction 
## 
## data:  pf8$Alter and pf8$Standort 
## 
##           Rheine  Münster Bahn    Ladbergen
## Münster   1.00000 -       -       -        
## Bahn      0.55704 1.00000 -       -        
## Ladbergen 0.56286 1.00000 1.00000 -        
## Internet  0.02780 3.5e-05 0.00027 0.00062  
## 
## P value adjustment method: bonferroni

Als Ergebnis erhalten wir eine Tabelle mit p-Werten. Es ist zu erkennen, dass sich die Altersverteilungen zwischen den Standortvergleichen Internet-Rheine, Internet-Münster und Internet-Landbergen unterscheiden.

32.8.4 Paarweise t-Test

Die Funktion pairwise.t.test() wird verwendet, um paarweise Vergleiche zwischen Gruppen durchzuführen, wenn die Daten normalverteilt sind.

# erzeuge Testdaten
Gruppe <- factor(rep(c("Handball", "Fussball", "Volleyball"), each = 10))
Punkte <- c(rnorm(10, mean = 5), rnorm(10, mean = 6), rnorm(10, mean = 7))

# Führe paarweise t-Tests durch
pairwise.t.test(Punkte, Gruppe, p.adjust.method = "holm")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  Punkte and Gruppe 
## 
##            Fussball Handball
## Handball   0.116    -       
## Volleyball 0.391    0.025   
## 
## P value adjustment method: holm

Als Ergebnis erhalten wir eine Tabelle mit p-Werten. Es ist zu erkennen, dass sich die Punkteverteilungen zwischen Handball-Fussball, und Handball-Volleyball unterscheiden.

32.9 Varianzanalysen

32.9.1 ANOVA mit Messwiederholung

Eine ANOVA mit Messwiederholung kann wie folgt in R durchgeführt werden. Wir nutzen den Testdatensatz Messwiederholung.rda.

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

Der Datensatz besteht aus 200 Probanden, bei denen zu drei Zeitpunkten ein Wert erhoben wurde.

head(Messwiederholung)
##      Name        t1       t2       t3
## 1  Brycen 11.728110 8.190290 6.202859
## 2   Waail  9.781605 8.646142 6.311455
## 3 Abigail 11.937967 9.386623 4.577478
## 4   Holly  4.597950 8.204762 7.250849
## 5  Sheila  8.277670 8.322449 6.952457
## 6 Rebekah  7.891258 9.104856 7.191639

Zunächst müssen wir die Daten ins long table-Format überführen.

# Überführe die Daten ins "long table"-Format
df <- Messwiederholung %>%
  pivot_longer(cols=t1:t3,
               names_to = "Zeitpunkt",
               values_to = "Wert") %>%
  mutate(Zeitpunkt = factor(Zeitpunkt),
         Name = factor(Name))

# anschauen
head(df)
## # A tibble: 6 × 3
##   Name   Zeitpunkt  Wert
##   <fct>  <fct>     <dbl>
## 1 Brycen t1        11.7 
## 2 Brycen t2         8.19
## 3 Brycen t3         6.20
## 4 Waail  t1         9.78
## 5 Waail  t2         8.65
## 6 Waail  t3         6.31

Zunächst verschaffen wir uns einen deskriptiven Überblick über die Daten.

# mal deskriptiv anschauen
df %>%
  group_by(Zeitpunkt) %>%
  summarise(M = mean(Wert),
            SD = sd(Wert),
            Q1 = quantile(Wert, probs = 0.25),
            Q2 = quantile(Wert, probs = 0.5),
            Q3 = quantile(Wert, probs = 0.75),
            IQR = IQR(Wert),
            N = n())
## # A tibble: 3 × 8
##   Zeitpunkt     M    SD    Q1    Q2    Q3   IQR     N
##   <fct>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 t1         9.88 2.14   8.41  9.77 11.3   2.91   200
## 2 t2         7.99 0.935  7.24  8.09  8.60  1.36   200
## 3 t3         6.92 0.868  6.31  6.89  7.51  1.21   200

Das geht auch graphisch:

# mal graphisch anschauen
boxplot(df$Wert ~ df$Zeitpunkt)

ez::ezPlot(df, Wert, Name, within = Zeitpunkt, x = Zeitpunkt)

# mit Regressionsgerade
ez::ezPlot(df, Wert, Name, within = Zeitpunkt, x = Zeitpunkt) +
  geom_smooth(aes(x = as.numeric(Zeitpunkt)), method = 'lm', se = F)
## `geom_smooth()` using formula = 'y ~ x'

Die Plots legen nahe, dass es einen Unterschied in den Zeitpunktgruppen gibt.

Wir überprüfen mittels Shapiro-Test und QQ-Plot, ob die Daten normalverteilt sind.

## qq-Plot (Normalverteilung)
ggpubr::ggqqplot(df, x="Wert", 
                 facet.by = "Zeitpunkt", 
                 color = "Zeitpunkt",
                 scales = "free")

# Teste "Wert" auf Normalverteilung
shapiro.test(df$Wert[df$Zeitpunkt=="t1"])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Wert[df$Zeitpunkt == "t1"]
## W = 0.99364, p-value = 0.5481
shapiro.test(df$Wert[df$Zeitpunkt=="t2"])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Wert[df$Zeitpunkt == "t2"]
## W = 0.98766, p-value = 0.08025
shapiro.test(df$Wert[df$Zeitpunkt=="t3"])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Wert[df$Zeitpunkt == "t3"]
## W = 0.99116, p-value = 0.2625

Die Testergebnisse sind nicht signifikant, also liegt Normalverteilung vor.

Wir dürfen eine ANOVA mit Messwiederholung rechnen. Mit R-Hausmittel könnten wir so vorgehen:

# ANOVA mit Messwiederholung 
fit <- aov(Wert ~ Zeitpunkt + Error(Name), data=df)
summary(fit)
## 
## Error: Name
##            Df Sum Sq Mean Sq F value Pr(>F)
## Residuals 199  395.8   1.989               
## 
## Error: Within
##            Df Sum Sq Mean Sq F value Pr(>F)    
## Zeitpunkt   2  898.0   449.0   212.3 <2e-16 ***
## Residuals 398  841.6     2.1                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Das Ergebnis ist signifikant. Mit der Funktion pairwise.t.test() können wir nun schauen, welche Gruppenunterschiede es konkret gibt:

# geht auch mit R-Hausmitteln so:
pairwise.t.test(df$Wert, df$Zeitpunkt, 
                p.adjust="bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  df$Wert and df$Zeitpunkt 
## 
##    t1      t2     
## t2 < 2e-16 -      
## t3 < 2e-16 1.1e-12
## 
## P value adjustment method: bonferroni

Die Unterschiede sind zwischen allen Gruppen signifikant.

Das Zusatzpaket rstatix bietet mit der Funktion anova_test() eine alternative Vorgehensweise.

# ANOVA mit Messwiederholung 
fit <- rstatix::anova_test(data = df, dv = Wert, wid = Name,
                           within = Zeitpunkt, effect.size = "pes")

# auf Sphärizität prüfen
fit$`Mauchly's Test for Sphericity`
##      Effect     W        p p<.05
## 1 Zeitpunkt 0.668 4.29e-18     *

Mauchly’s Test ist signifikant, das heisst, es liegt keine Sphärizität vor. R hat in diesem Falle bereits die Freiheitsgrade automatisch korrigiert (mittels Greenhouse-Geisser-Verfahren).

Schauen wir uns die ANOVA an:

rstatix::get_anova_table(fit)
## ANOVA Table (type III tests)
## 
##      Effect DFn    DFd      F        p p<.05   pes
## 1 Zeitpunkt 1.5 298.73 212.34 2.18e-48     * 0.516

Die ANOVA ist signifikant. Jetzt können wir analysieren, welche konkreten Gruppenunterschiede vorliegen.

# post-hoc Analyse: paarweise Gruppenvergleiche
df %>%
  rstatix::pairwise_t_test(Wert ~ Zeitpunkt, paired = TRUE,
                           p.adjust.method = "bonferroni") %>%
  as.data.frame()
##    .y. group1 group2  n1  n2 statistic  df        p    p.adj p.adj.signif
## 1 Wert     t1     t2 200 200  11.24384 199 5.02e-23 1.51e-22         ****
## 2 Wert     t1     t3 200 200  18.28475 199 1.79e-44 5.37e-44         ****
## 3 Wert     t2     t3 200 200  11.26870 199 4.23e-23 1.27e-22         ****

Die Unterschiede sind zwischen allen Gruppen signifikant.

Die Effektstärke kann wie folgt berechnet werden.

# Effektstärke (nur solche, die post-hoc signifikant sind, interpretieren)
df %>%
  rstatix::cohens_d(Wert ~ Zeitpunkt, paired = TRUE) %>%
  as.data.frame()
##    .y. group1 group2   effsize  n1  n2 magnitude
## 1 Wert     t1     t2 0.7950598 200 200  moderate
## 2 Wert     t1     t3 1.2929274 200 200     large
## 3 Wert     t2     t3 0.7968171 200 200  moderate

Achtung, nur solche Werte interpretieren, die oben signifkant waren. In unserem Beispiel sind alle Werte signifikant, darum dürfen wir alle Effekte interpretieren.

32.9.2 Friedman ANOVA

Sind die Werte nicht normalverteilt, kann ein Friedman-ANOVA gerechnet werden. Wir nutzen den Testdatensatz MarioANOVA.rda. Er enthält Informationen über 47 Super Mario Charaktere, deren Fehlerquote beim Spielen eines schweren Levels zu 3 Zeitpunkten erfasst wurde.

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

# anschauen
head(MarioANOVA)
## # A tibble: 6 × 8
##   Name             Alter Kingdom          Geschlecht BadGuy    t1    t2    t3
##   <fct>            <dbl> <fct>            <fct>      <lgl>  <dbl> <dbl> <dbl>
## 1 Mario               35 Mushroom Kingdom männlich   FALSE  0.372 0.472 0.257
## 2 Luigi               35 Mushroom Kingdom männlich   FALSE  1.27  1.34  1.14 
## 3 Prinzessin Peach    25 Mushroom Kingdom weiblich   FALSE  6.15  1.40  1.98 
## 4 Bowser              45 Bowser's Castle  männlich   TRUE   4.85  0.648 1.27 
## 5 Toad                20 Mushroom Kingdom männlich   FALSE  1.34  0.174 1.77 
## 6 Toadette            19 Mushroom Kingdom weiblich   FALSE  0.727 0.413 2.09

In den Variablen t1 bis t3 liegen die Werte der 3 Messzeitpunkte. Auch hier überführen wir die Daten ins long table-Format

df <- MarioANOVA %>%
  select(Name, t1:t3) %>%
  pivot_longer(cols=t1:t3,
               names_to = "Zeitpunkt",
               values_to = "Wert")

Schauen wir uns die Daten zunächst deskriptiv an.

# mal deskriptiv anschauen
df %>%
  group_by(Zeitpunkt) %>%
  summarise(M = mean(Wert),
            SD = sd(Wert),
            Q1 = quantile(Wert, probs = 0.25),
            Q2 = quantile(Wert, probs = 0.5),
            Q3 = quantile(Wert, probs = 0.75),
            IQR = IQR(Wert),
            N = n())
## # A tibble: 3 × 8
##   Zeitpunkt     M    SD    Q1    Q2    Q3   IQR     N
##   <chr>     <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <int>
## 1 t1        2.75  4.61  0.750 1.34   2.46 1.71     47
## 2 t2        0.762 0.585 0.356 0.509  1.15 0.791    47
## 3 t3        0.901 0.701 0.266 0.650  1.34 1.07     47

Das geht auch graphisch wie oben beschrieben:

# mal graphisch anschauen
boxplot(df$Wert ~ df$Zeitpunkt)

ez::ezPlot(df, Wert, Name, within = Zeitpunkt, x = Zeitpunkt)

# mit Regressionsgerade
ez::ezPlot(df, Wert, Name, within = Zeitpunkt, x = Zeitpunkt) +
  geom_smooth(aes(x = as.numeric(Zeitpunkt)), method = 'lm', se = F)

Nun prüfen wir mittels QQ-Plot und Shapiro-Test auf Normalverteilung.

## qq-Plot (Normalverteilung)
ggpubr::ggqqplot(df, x="Wert", 
                 facet.by = "Zeitpunkt", 
                 color = "Zeitpunkt",
                 scales = "free")

# Teste "Wert" auf Normalverteilung
shapiro.test(df$Wert[df$Zeitpunkt=="t1"])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Wert[df$Zeitpunkt == "t1"]
## W = 0.52189, p-value = 3.729e-11
shapiro.test(df$Wert[df$Zeitpunkt=="t2"])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Wert[df$Zeitpunkt == "t2"]
## W = 0.89427, p-value = 0.0004751
shapiro.test(df$Wert[df$Zeitpunkt=="t3"])
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Wert[df$Zeitpunkt == "t3"]
## W = 0.90101, p-value = 0.0007756

Die Werte für t1 bis t3 sind signifikant, also nicht normalverteilt.

Wir führen daher eine Friedman-ANOVA mit Messwiederholung durch.

# Friedman-Test
rstatix::friedman_test(Wert ~ Zeitpunkt | Name, data=df)
## # A tibble: 1 × 6
##   .y.       n statistic    df       p method       
## * <chr> <int>     <dbl> <dbl>   <dbl> <chr>        
## 1 Wert     47      12.8     2 0.00165 Friedman test

Das Ergebnis ist signifikant, es gibt also einen Unterschied.

# post-hoc-test
df %>%
  rstatix::wilcox_test(Wert ~ Zeitpunkt, paired=T,
                       p.adjust.method = "bonferroni")
## # A tibble: 3 × 9
##   .y.   group1 group2    n1    n2 statistic        p    p.adj p.adj.signif
## * <chr> <chr>  <chr>  <int> <int>     <dbl>    <dbl>    <dbl> <chr>       
## 1 Wert  t1     t2        47    47       917 0.000102 0.000306 ***         
## 2 Wert  t1     t3        47    47       882 0.000536 0.002    **          
## 3 Wert  t2     t3        47    47       506 0.546    1        ns

Der Unterschied existiert zwischen t1 und t2 sowie zwischen t1 und t3.

Die Effktstärke wird wie folgt berechnet:

# Effektstärke
rstatix::friedman_effsize(Wert ~ Zeitpunkt | Name, data = df)
## # A tibble: 1 × 5
##   .y.       n effsize method    magnitude
## * <chr> <int>   <dbl> <chr>     <ord>    
## 1 Wert     47   0.136 Kendall W small

Der Effekt des Gesamtmodells ist eher gering. Schauen wir uns die einzelnen Gruppeneffekte an.

Achtung, nur solche Werte interpretieren, die oben signifkant waren. In unserem Beispiel ist das der Unterschied von t1 zu t2 sowie von t1 zu t3..

# nur die signifikanten anschauen (t1 -> t2 und t1 -> t3)
df %>%
  rstatix::wilcox_effsize(Wert ~ Zeitpunkt, paired = TRUE)
## # A tibble: 3 × 7
##   .y.   group1 group2 effsize    n1    n2 magnitude
## * <chr> <chr>  <chr>    <dbl> <int> <int> <ord>    
## 1 Wert  t1     t2      0.545     47    47 large    
## 2 Wert  t1     t3      0.491     47    47 moderate 
## 3 Wert  t2     t3      0.0895    47    47 small

Der Effekt von t1 nach t2 ist “groß”, der von t1 nach t3 “moderat”.

32.9.3 einfaktorielle ANOVA

Eine einfaktorielle Varianzanalyse (ANOVA) wird mit der Funktion aov() durchgeführt. Sie wird nach der Syntax x erklärt durch y aufgerufen. Wie bei anderen Modellen ist es ratsam, die Ergebnisse der Berechnung in ein Objekt zu speichern, und dieses per summary() anzuschauen.

# Lade Testdatensatz
load(url("https://www.produnis.de/R/data/pf8.RData"))

# erstelle ANOVA  "Alter" erklärt durch "Standort"
fit <- aov(Alter ~ Standort, data=pf8)

# schaue Modellzusammenfassung an
summary(fit)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## Standort      4  10557  2639.3    8.48 1.09e-06 ***
## Residuals   725 225633   311.2                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 1 Beobachtung als fehlend gelöscht

Das Ergebnis ist signifikant, wir haben unterschiedliche Varianzen im Alter erklärt durch Standort.

Wir wissen aber noch nicht, innerhalb welche Standorte die Unterschiede bestehen. Dazu führen wir post hoc Tests durch, in denen alle Gruppen paarweise verglichen werden. Wegen dieses multiplen Testens muss eine Fehlerkorrektur für Alpha durchgeführt werden. Dazu steht die Funktion pairwise.t.test() zur Verfügung, die mit den selben Angaben aufgerufen wird.

pairwise.t.test(pf8$Alter, pf8$Standort)
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  pf8$Alter and pf8$Standort 
## 
##           Rheine Münster Bahn   Ladbergen
## Münster   0.5042 -       -      -        
## Bahn      0.4354 1.0000  -      -        
## Ladbergen 0.5432 1.0000  1.0000 -        
## Internet  0.0090 1.8e-05 0.0002 0.0059   
## 
## P value adjustment method: holm

Die gewählte Fehlerkorrektur ist standardmäßig holm. Möchten Sie auf die Bonferroni-Korrekturmethode wechseln, übergeben Sie den Parameter p.adjust="bonferroni".

pairwise.t.test(pf8$Alter, pf8$Standort, p.adjust="bonferroni")
## 
##  Pairwise comparisons using t tests with pooled SD 
## 
## data:  pf8$Alter and pf8$Standort 
## 
##           Rheine  Münster Bahn    Ladbergen
## Münster   1.00000 -       -       -        
## Bahn      0.72572 1.00000 -       -        
## Ladbergen 1.00000 1.00000 1.00000 -        
## Internet  0.01292 1.8e-05 0.00022 0.00733  
## 
## P value adjustment method: bonferroni

Die signifikanten Unterschiede befinden sich zwischen dem Standort Internet und allen anderen Standorten.

Ein anschauliches Video zur ANOVA in R findet sich bei Youtube: https://www.youtube.com/watch?v=GctGncQrJew

32.10 Überlebenszeitanalysen

32.10.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

32.10.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)

32.11 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.

32.11.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/data/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.

32.11.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.
## 
## df null model =  66  with the objective function =  6.38 with Chi Square =  920.48
## df of  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 n.obs is  150 with the empirical chi square  20.2  with prob <  0.96 
## The total n.obs 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] "CFI"           "RMSEA"         "BIC"           "SABIC"        
## [29] "r.scores"      "R2"            "valid"         "score.cor"    
## [33] "weights"       "rotation"      "hyperplane"    "communality"  
## [37] "communalities" "uniquenesses"  "values"        "e.values"     
## [41] "loadings"      "model"         "fm"            "rot.mat"      
## [45] "Structure"     "method"        "scores"        "R2.scores"    
## [49] "r"             "np.obs"        "fn"            "Vaccounted"   
## [53] "ECV"
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

32.11.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"   "valid"        "chi"          "EPVAL"        "R2"          
## [16] "objective"    "residual"     "rms"          "factors"      "dof"         
## [21] "null.dof"     "null.model"   "criteria"     "STATISTIC"    "PVAL"        
## [26] "weights"      "r.scores"     "rot.mat"      "Vaccounted"   "Structure"   
## [31] "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.

32.12 Fallzahlkalkulation

32.12.1 klinische Studien

Zur Fallzahlkalkulation bei klinischen Studien (z.B. RCTs) kann das Zusatzpaket pwr installiert werden. Mit der Funktion pwr.t.test() können alle Werte berechnet werden. Ihr müssen die entsprechenden Werte für das Signifikanzlevel, Power und Cohen’s d (Effektstärke) übergeben werden.

# erechne Fallzahl mit
# Cohens d = 0.1
# alpha = 0.05
# Power = 0.8
pwr::pwr.t.test(d=0.1, sig.level=0.05, power=0.8, type="two.sample", alternative="two.sided")
## 
##      Two-sample t test power calculation 
## 
##               n = 1570.733
##               d = 0.1
##       sig.level = 0.05
##           power = 0.8
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

Über die Parameter type und alternative lassen sich die Testbedingungen festlegen.

# Lade Paket
library(pwr)

# unabhängige Stichprobe, zweiseitig
pwr.t.test(d=0.1, sig.level=0.05, power=0.8, type="two.sample", alternative="two.sided")

# Einstichproben-Test, kleiner
pwr.t.test(d=0.1, sig.level=0.05, power=0.8, type="one.sample", alternative="less")

# gepaarte Stichprobe, größer
pwr.t.test(d=0.1, sig.level=0.05, power=0.8, type="paired", alternative="greater")

Möchten Sie nicht immer die gesamte Testausgabe angezeigt bekommen, sondern nur den Wert für n, können Sie dies mit einem Dollarzeichen $ hinter dem Funktionsaufruf spezifizieren.

# zeige nur "n"
pwr.t.test(d=0.1, sig.level=0.05, power=0.8, type="two.sample")$n
## [1] 1570.733
# runde auf nächste ganze Zahl
ceiling(pwr.t.test(d=0.1, sig.level=0.05, power=0.8, type="two.sample")$n)
## [1] 1571

Sie benötigen 1571 Probanden pro Untersuchungsgruppe.

32.12.2 Umfragen

Soll eine Befragung (ein Survey) durchgeführt werden, kann die benötigte Fallzahl über das Zusatzpaket samplingbook errechnet werden.

Geht es um numerische Werte, wird die Funktion sample.size.mean() verwendet. Ihr übergibt man Fehlerspanne e (ein Präzisionswert, der besagt, wie breit das Konfidenzintervall - nach oben sowie nach unten - sein darf, für gewöhnlich 5% bzw. 0.05), das Konfidenzniveau, die Standardabweichung S innerhalb der Grundgesamtheit (füre gewöhnlich 50% bzw. 0.5) sowie die Populationsgröße. Standardmäßig ist die Populationsgröße auf N=Inf (für unendlich) gesetzt.

# Berechne Fallzahl für einen Survey
# Konfidenzniveau auf 95%, Fehlerbereich 5%.
# In der Grundgesamtheit hat der Wert eine
# Standardabweichung von S=0.5
samplingbook::sample.size.mean(e=0.05, S=0.5, level=0.95 )
## 
## sample.size.mean object: Sample size for mean estimate
## Without finite population correction: N=Inf, precision e=0.05 and standard deviation S=0.5
## 
## Sample size needed: 385

Ist die Grundgesamtheit bekannt, kann dies mit N angegeben werden.

# Die Grundgesamtheit ist 1.500 Personen groß
samplingbook::sample.size.mean(e=0.05, S=0.5, level=0.95, N=1500 )
## 
## sample.size.mean object: Sample size for mean estimate
## With finite population correction: N=1500, precision e=0.05 and standard deviation S=0.5
## 
## Sample size needed: 306

Sollen Prozentwerte kategorialer Variablen erfragt werden (z.B. der Anteil weiblicher Pflegepersonen), wird die Funktion sample.size.prop() verwendet. Ihr werden ebenso e, N und level (Konfidenzniveau) übergeben. Statt der Standardabweichung wird über den Parameter P die prozentuale Verteilung in der Grundgesamtheit angegeben.

# Berechne Fallzahl für einen Survey
# unsere Werte düren maximal 5% abweichen
# Konfidenzniveau auf 95%
# In der Grundgesamtheit hat der Wert eine
# Verteilung von 65%
# Die Grundgesamtheit ist 1.500 Personen groß
samplingbook::sample.size.prop(e=0.05, P=0.65, level=0.95, N=1500)
## 
## sample.size.prop object: Sample size for proportion estimate
## With finite population correction: N=1500, precision e=0.05 and expected proportion P=0.65
## 
## Sample size needed: 284

32.12.2.1 Tabellen

Um ein Gefühl für die benötigten Stichprobengrößen zu erhalten, können wir Tabellen erzeugen, in denen die benötigten Fallzahlen beispielsweise in Abhängigkeit zur Populationsgröße angezeigt werden. Die Tabelle soll Werte für metrische und kategoriale Items anzeigen, wobei wir jeweils das 90%-, 95% und 99%-Konfidenzintervall einsehen wollen. In diesem Beispiel gehen wir für metrische Daten von 30% Streuung aus (S=0.3) und setzen die Fehlerspanne e auf 5%. Für kategoriale Items nehmen wir eine Verteilung von 50% an (P=0.5).

# gewünschte Populationsgrößen
Population <- c(100, 200, 300, 400, 500, 600, 700, 800, 900, 1000,
                1500, 2000, 4000, 6000, 8000, 10000, Inf)

# überführe in Datenframe
df <- data.frame(Population)

Mit der Funktion sapply() können wir sample.size.mean() und sample.size.prop() auf jeden Wert des Vectors Population anwenden. Die Ergebnisse schreiben wir als neue Variable in unser Datenframes df.

## Wende sample.size.mean() auf jeden Wert in "Population" an
# 90%-Konfidenzintervall
df$ConA10 <- sapply(Population, FUN=function(x) samplingbook::sample.size.mean(e=0.05, S=0.3, level=0.90, N=x)$n)
# 95%-Konfidenzintervall
df$ConA05 <- sapply(Population, FUN=function(x) samplingbook::sample.size.mean(e=0.05, S=0.3, level=0.95, N=x)$n)
# 99%-Konfidenzintervall
df$ConA01 <- sapply(Population, FUN=function(x) samplingbook::sample.size.mean(e=0.05, S=0.3, level=0.99, N=x)$n)

## Wende sample.size.prop() auf jeden Wert in "Population" an
# 90%-Konfidenzintervall
df$CatA10 <- sapply(Population, FUN=function(x) samplingbook::sample.size.prop(e=0.05, P=0.5, level=0.9, N=x)$n)
# 95%-Konfidenzintervall
df$CatA05 <- sapply(Population, FUN=function(x) samplingbook::sample.size.prop(e=0.05, P=0.5, level=0.95, N=x)$n)
# 99%-Konfidenzintervall
df$CatA01 <- sapply(Population, FUN=function(x) samplingbook::sample.size.prop(e=0.05, P=0.5, level=0.99, N=x)$n)

Die Variablennamen ConA10 - CatA01 sind etwas kryptisch und für eine “schöne” Tabelle eher ungeeignet. Daher ändern wir die Variablennamen mittels colnames() in leserliche Spaltennamen um. Zuvor speichern wir aber unser Datenframe in ein neues Objekt df2, denn mit den veränderten Variablennamen können wir später viel schlechter weiterarbeiten.

# Mache Sicherungskopie von df
df2 <- df

# ändere die Variablennamen für eine "schöne" Tabelle
Kopfzeile <- c("Population N", "metrisch 90%", "metrisch 95%", "metrisch 99%",
               "kategor 90%", "kategor 95%", "kategor 99%")
colnames(df) <- Kopfzeile

Mit den neuen Variablennamen sieht die Tabelle nun so aus:

Population N metrisch 90% metrisch 95% metrisch 99% kategor 90% kategor 95% kategor 99%
100 50 59 71 74 80 87
200 66 82 109 115 132 154
300 74 95 133 143 169 207
400 79 103 150 162 196 250
500 82 109 162 176 218 286
600 84 113 171 187 235 316
700 86 116 179 196 249 341
800 87 118 184 203 260 363
900 88 120 189 209 270 382
1000 89 122 193 213 278 399
1500 92 127 207 230 306 461
2000 93 130 214 239 323 499
4000 96 134 226 254 351 570
6000 96 136 230 259 362 598
8000 97 136 232 262 367 613
10000 97 137 234 264 370 623
Inf 98 139 239 271 385 664

Und um auf das nächste Kapitel vorzugreifen, können wir die Daten der Tabelle in einem einfachen Plot graphisch darstellen.

# Erzeuge den Grund-Plot mit 1 y-Variable
plot(Population, df2$ConA01, type="b", ylim=c(0,670), 
        ylab="Fallzahl", col="darkgreen", pch=19, lty=2)

# Füge weitere y-Variablen mit anderer Farbe hinzu
lines(Population, df2$ConA05, col="seagreen", type="b", pch=19, lty=1)
lines(Population, df2$ConA10, col="lightgreen", type="b", pch=19, lty=3)
lines(Population, df2$CatA01, col="darkblue", type="b", pch=18, lty=2)
lines(Population, df2$CatA05, col="blue", type="b", pch=18, lty=1)
lines(Population, df2$CatA10, col="lightblue", type="b", pch=18, lty=3)

# Füge eine Legendenbox hinzu
legend(x=0, y=670,
       legend = c("metrisch 90%", "metrisch 95%", "metrisch 99%",
               "kategori 90%", "kategori 95%", "kategori 99%"),
       col = c("lightgreen", "seagreen", "darkgreen",
               "lightblue", "blue", "darkblue"), 
       lty=c(3,1,2), cex=0.8)

Sieht kompliziert aus? Warten Sie es ab, es wird alles im folgenden Kapitel erklärt.