Kapitel 11 Korrelationen

Dies ist eine unvollständige und unkorrigierte Version 2022-12-24.

Wir lernen in diesem Kapitel folgendes:

• Streudiagramm (Scatter Plot)
• Korrelationskoeffizienten (Pearson, Spearman, etc.)
• Regressionslinie

Wir benötigen folgende Pakete:

library(dplyr) # Paket für Datenmanipulation
library(psych) # Paket mit Hilfsmittel für Analysen zu Messeigenschaften (Psychometrics) und Datendarstellungen
library(ggplot2) # Grammar of Graphics 
library(ggiraph) # Hilfmittel zu GGplot für interaktive Visualisierungen
library(jtools)  # Hilft Resultate von Regressionen einfacher darzustellen
library(MASS, exclude="select") # Verschiedene Hilfsfunktionen, dass Paket wurde als Begleitung zum Buch *Modern Applied Statistics with S* geschrieben (S ist die Vorgängerin von R)
library(AER) # für Ökonomische Analysen
library(pscl) # Verschiedene Hilfsmittel für Item Response Analysen und Regressionen (z.B. Count Analysen)
library(rms) # Hilfmisttel für Regressionen von Frank Harrell, wir benutzen es u.a. für Splines
library(faux) # Hilfmittel zur Datensimulation
library(kableExtra) # Tabellen darstellen

Korrelationen können im Prinzip zwischen allen quantitativen Daten (d.h. binären, ordinalen, intervallskalierten, und Proportionalskalen) berechnet werden. Wir benutzen Sie jedoch am meisten für ordinale und intervall- oder proportional skalierte Daten.

Je nach Skalenniveau, Verteilung der Daten und der Art des Zusammenhanges, benutzen wir unterschiedliche Formen der Korrelationen oder andere Zusammenhangsmasse.

• Lineare Zusammenhänge von mindestens Intervallskalierten Daten: Pearson Korrelation
• Lineare Zusammenhänge ordinaler Daten: Spearman Korrelation
• Monotone Zusammenhänge intervallskalierter Daten: Spearman Korrelation
• Zusammenhänge zwischen zwei Variablen, wenn eine kontinuierliche latente Variable als ordinale Variable operationalisiert wurde: Polychorische Korrelation
• Zusammenhänge zwischen zwei Variablen, wenn eine kontinuierliche latente Variable als dichotome Variable operationalisiert wurde: Tetrachorische Korrelation
• Zusammenhänge zwischen intervallskalierter Variable und einer dichotomen Variable: Punktbiseriale Korrelation (Spezialfall Pearson Korrelation)

11.1 Pearson’s r

Korrelationen mit dem Pearson’s Korrelationskoeffizietenten (auch Produkt-Moment-Korrelation genannt) sind geeignet für lineare Zusammenhänge zweier intervallskalierter Variablen. Die Pearson Korrelation wird meistens mit r abgekürzt (Beispiel: r = 0.8).

11.1.1 Voraussetzungen für die Pearson’s Korrelation (gilt immer für beide Variablen, da ja eine Korrelation zwischen zwei Variablen berechnet wird):

• mindestens intervallskalierte Daten
• gepaarte Daten, d.h. die Beobachtung in einer Variable sollte einer Beobachtung der zweiten Variable zugeordnet werden können. Wir können also keine Korrelation zwischen zwei unabhängigen Gruppen berechnen. Wenn wir die Korrelation zwischen Berg Balance Scale Score und der Angst vor dem Stürzen (FES-I) berechnen wollen, muss jeder Patient einen Wert für den BBS haben und einen Wert für den FES-I. Wir können nicht eine Gruppe von Patienten mit dem BBS messen und eine andere Gruppe mit dem FES-I. Praktischer Tipp: Können wir ein Streudiagramm zeichnen, so können wir auch eine Korrelation berechnen.
• Linearer Zusammenhang: Ist der Zusammenhang nicht linear, sollte keine Pearson Korrelation berechnet werden. (Falls Zusammenhang monoton: Spearman’s Korrelation).
  
• Keine Ausreisser: Diese beeinflussen die Korrelation. Wir sollten immer ein Streudiagramm erstellen um Ausreisser zu erkennen.
• Achtung bei mehreren Gruppen (siehe Abbildung 13.5 weiter unten)!
• Homogene Varianzen: Die Streuung sollte in etwa gleich bleiben über die Ausprägungen beider Variablen, doch am besten hilft hier eine Graphik, siehe Abbildung 13.1. Man nennt dies auch Homoskedastizität (Homoscedasticity). Das Gegenteil wäre Heteroskedastizität (Heteroscedasticity)

set.seed(1234)
par(mar = c(4, 4, .1, .1))

id<-1:1000
variable_A<-rnorm(1000, 100,30)

data<-data.frame(id, variable_A)

data<-data %>% 
  mutate(variable_B=variable_A+rnorm(id, 40, 0.299*variable_A))

plot(data$variable_A, data$variable_B, col="#B40CF5")
round(cor(data$variable_A, data$variable_B, use="complete.obs"),2)

## [1] 0.71

data<-data %>% 
  mutate(variable_C=variable_A+rnorm(id, 2, 30))

plot(data$variable_A, data$variable_C, col="#B40CF5")
cor(data$variable_A, data$variable_C)

## [1] 0.7105243

Abbildung 2.8: Linke Graphik zeigt eine nicht-homogene Varianz, die Streuung der Variable B nimmt mit höheren Werten der Variable A zu. Auf der rechten Seite sehen wir eine homogene Varianz. In beiden Varianten ist die wahre Korrelation 0.7.

Visualisieren können wir das auch mit folgender Graphik, wenn wir eine lineare Regression durchführen (Wir werden später sehen, wie der Zusammenhang zwischen der Korrelation und einem linearen Regressionsmodell ist):

par(mar = c(4, 4, .1, .1))
lm.fit_heteroscedasticity<-lm(variable_A~variable_B, data=data)
lmtest::bptest(lm.fit_heteroscedasticity)

## 
##  studentized Breusch-Pagan test
## 
## data:  lm.fit_heteroscedasticity
## BP = 11.205, df = 1, p-value = 0.0008157

plot(lm.fit_heteroscedasticity, which=1)

lm.fit_homoscedasticity<-lm(variable_A~variable_C, data=data)
lmtest::bptest(lm.fit_homoscedasticity)

## 
##  studentized Breusch-Pagan test
## 
## data:  lm.fit_homoscedasticity
## BP = 0.00032745, df = 1, p-value = 0.9856

plot(lm.fit_homoscedasticity, which=1)

Abbildung 3.1: Streudiagramm der Residuen über die geschätzten Werte (fitted values). Linke Graphik zeigt eine nicht-homogene Varianz, die Streuung der Residuen ist nicht gleich über die Bandbreite der gefitteten Werte. Auf der rechten Seite sehen wir eine homogene Varianz.

Zuerst das Beispiel mit Heteroskedastizität: Es gibt auch statistische Tests um zu testen, ob die Varianzen gleich sind. Zum Beispiel der Breush Pagan Test. Hier sehen wir, dass der p-Wert kleiner als 0.05 ist, wir verwerfen deswegen die Nullhypothese, dass die Varianzen homogen sind.

lm.fit_heteroscedasticity<-lm(variable_A~variable_B, data=data)
lmtest::bptest(lm.fit_heteroscedasticity)

## 
##  studentized Breusch-Pagan test
## 
## data:  lm.fit_heteroscedasticity
## BP = 11.205, df = 1, p-value = 0.0008157

Zuerst das Beispiel mit Homoskedastizität: Hier unten sehen wir, dass der P-Wert über 0.05 ist, deswegen verwerfen wir die Nullhypothese (dass wir Homoskedastizität haben) nicht.

lm.fit_homoscedasticity<-lm(variable_A~variable_C, data=data)
lmtest::bptest(lm.fit_homoscedasticity)

## 
##  studentized Breusch-Pagan test
## 
## data:  lm.fit_homoscedasticity
## BP = 0.00032745, df = 1, p-value = 0.9856

• Normalverteilung: Man liest oft, dass für eine Pearson Korrelation die Daten normalverteilt sein sollten. Eigentlich benötigen wir eine bivariate Normalverteilung
• Die bivariate Normalverteilung kann auch getestet werden hier klicken

Bei grösseren Stichproben hilft auch der zentrale Grenzwertsatz wieder (central limit theorem). Was heisst hier gross? Das kommt natürlich auf die Verteilungen der Daten in der Population an, es sollten jedoch mindestens 30 sein. Doch auch mit grossen Stichproben, kann es zu Problemen kommen, wenn wir schiefe Verteilungen haben - und insbesondere, wenn wir Ausreisser haben. Siehe auch Kapitel Problems With Pearson’s Correlation in Wilcox, Rand R.. Introduction to Robust Estimation and Hypothesis Testing (Statistical Modeling and Decision Science) (p. 541). Elsevier Science. Kindle Edition..

Siehe auch diesen Artikel

11.1.2 Formel Pearson’s r

Die Formel ist eigentlich relativ einfach:

• Es ist das Verhältnis zwischen der Kovarianz von X und Y und dem Produkt ihrer Standardabweichungen. An den griechischen Buchstaben sehen wir, dass es hier die Korrelation und die Standardabweichungen der Population sind:

\[\rho_{xy} = \frac{cov(x,y)}{\sigma_x\sigma_y}\] Was ist denn die Kovarianz? Wir benutzen die Kovarianz als Kennzahl eigentlich nie direkt, sie wird jedoch oft im Hintergrund benutzt. Sie ist eine Kennzahl darüber, wie stark ein hoher Wert einer Variable mit einem hohen oder einem tiefen Wert einer anderen Variablen einhergeht. Sie finden eine graphische Darstellung hier (hier klicken).

Die Formel für die Kovarianz

\[ Cov(x,y) = \frac{\sum_{i=1}^{n}(x_i - \overline{x})(y_i - \overline{y})}{n - 1}\]

Wer mag, kann sich dieses kurze Beispiel anschauen. Wir haben 15 Patient:innen mit Daten der Berg Balance Scale und ihrer normalen Gehgeschwindigkeit in Metern pro Sekunde.

Um die Kovarianz zu berechnen benögiten wir (wie wir in der Formel sehen) die beiden Mittelwert, die Differenz von jedem Wert zum Mittelwert, und das Produkt dieser Differerenzen beider Variablen.

data<-rio::import("BBS_normal_Walking_Speed.dta")

data<-data %>% 
  mutate(Mean_BBS=mean(BBS_Total), 
         Mean_Speed=mean(speed)) %>% 
  mutate(Diff_BBS=BBS_Total-Mean_BBS, 
         Diff_Speed=speed-Mean_Speed) %>% 
  mutate(Produkt=Diff_BBS*Diff_Speed)


kbl(data) %>% kable_styling()

BBS_Total	speed	Mean_BBS	Mean_Speed	Diff_BBS	Diff_Speed	Produkt
51	0.6666667	36.46667	0.5866821	14.533333	0.0799846	1.1624428
48	0.7142857	36.46667	0.5866821	11.533333	0.1276036	1.4716954
50	0.9090909	36.46667	0.5866821	13.533333	0.3224088	4.3632664
51	0.7692308	36.46667	0.5866821	14.533333	0.1825487	2.6530410
8	0.0793651	36.46667	0.5866821	-28.466667	-0.5073170	14.4416241
6	0.0571429	36.46667	0.5866821	-30.466667	-0.5295392	16.1332951
16	0.2702703	36.46667	0.5866821	-20.466667	-0.3164118	6.4758955
33	0.3846154	36.46667	0.5866821	-3.466667	-0.2020667	0.7004979
45	1.2500000	36.46667	0.5866821	8.533333	0.6633179	5.6603128
44	0.5882353	36.46667	0.5866821	7.533333	0.0015532	0.0117010
18	0.1388889	36.46667	0.5866821	-18.466667	-0.4477932	8.2692476
50	1.0000000	36.46667	0.5866821	13.533333	0.4133179	5.5935691
47	0.7692308	36.46667	0.5866821	10.533333	0.1825487	1.9228463
31	0.2941177	36.46667	0.5866821	-5.466667	-0.2925644	1.5993522
49	0.9090909	36.46667	0.5866821	12.533333	0.3224088	4.0408576

Jetzt müssen wir nur noch die Summe der Produkte berechnen (für unseren Zähler der Formel).

Summe_der_Produkte<-sum(data$Produkt)
Summe_der_Produkte

## [1] 74.49964

Die Formel sieht jetzt so aus:

\[ cov(BBS, Walking~Speed) = \frac{\sum_{i=1}^{n}(BBS_i - Mean_BBS)\sum_{i=1}^{n}(Speed_i-MeanSpeed))}{15-1} \] Und mit eingesetzten Zahlen:

\[ cov(BBS,Walking~Speed) = \frac{74.49964}{14} = 5.321403 \]

Summe_der_Produkte/14

## [1] 5.321403

Für die Korrelation brauchen wir jetzt noch die Standardabweichungen des BBS und der Gehgeschwindigkeit.

(sd_BBS=sd(data$BBS_Total))

## [1] 16.60407

(sd_Walking_Speed = sd(data$speed))

## [1] 0.3657566

Die Formel für die Korrelation sieht mit den eingesetzten Zahlen so aus:

\[ r = \frac{5.321403}{16.60407*0.3657566} = 0.8762326 \]

Wenn wir das Ergebnis präsentieren, werden wir natürlich runden also zum Beispiel r = 0.88.

In R sieht dies so aus:

5.321403/(16.60407*0.3657566)

## [1] 0.8762326

Oder gerunded so:

round(5.321403/(16.60407*0.3657566),2)

## [1] 0.88

Zur Kontrolle können wir ja auch den R-Befehl cor benutzen:

cor(data$BBS_Total, data$speed)

## [1] 0.8762325

Oder gerundet so:

round(cor(data$BBS_Total, data$speed),2)

## [1] 0.88

Gratuliere! Sie haben soeben eine Kovarianz und eine Korrelation (fast) von Hand gerechnet.

Ein weiteres Beispiel zum Mitrechnen finden Sie hier.

Die Formel oben ist für die Korrelation in der Population, in einer Studie kennen wir natürlich die Populationsparameter nicht, deswegen benutzen wir die Stichprobenwerte.

Meistens findet man diese Formel:

\[ r_{xy} = \frac{n \sum_{i = 1}^{n} (x_i - \overline{x}) (y_i - \overline{y})} {\sqrt{\sum_{i=1}^n(x_i- \overline{x})^{2} \sum_{i=1}^{n}(y_i - \overline{y})^{2}}}\]

Die Formel kann natürlich auch in verschiedenen Transformationen geschrieben werden. Zum Beispiel so:

\[ r = \frac{S_{xy}}{S_x*S_y}\] Hier ist \(S_x\) und \(S_y\) die Standardabweichung der jeweiligen Variablen. \(S_{xy}\) ist: \[ \frac{1}{n-1} \sum_{i=1}^{n} (x_i - \overline{x}) (x_i - \overline{y}) \]

Eine gute Erklärung zur Formel finden Sie hier (hier klicken), oder hier auf Englisch ein gerechnetes Beispiel hier klicken

11.1.3 Zusammenhang zwischen Pearson’s r und Regressionskoeffizienten

Den Regressiononskoeffizienten erhält man durch Multiplikation des Verhältnisses von SD von Y und SD von X mit dem Korrelationskoeffizienten “r”.

Beispiel:

set.seed(123)
x = rnorm(100, 50, 12)
y = x+5+rnorm(100,0, 8.2)

sd_x<-sd(x)
sd_y<-sd(y)
(r<-cor(x,y))

## [1] 0.8000443

plot(x,y)

lm.fit<-lm(y~x)
(beta_coefficient<-lm.fit$coefficients[2])

##         x 
## 0.9641444

Wir können jetzt den Regressionskoeffizienten wie folgt berechnen:

regressionskoeffizient = (sd_y / sd_x) * r
regressionskoeffizient

## [1] 0.9641444

beta_coefficient

##         x 
## 0.9641444

Wir sehen jetzt, dass der berechnete Regressionskoeffizient gleich ist, wie der koeffizient aus der Regressionsgleichung. Wir könnten das jetzt auch wieder zurückrechnen, in dem wir den Regressionskoeffizienten mit dem Quotionten \(\frac{sd(x)}{sd(y)}\) multiplizieren.

pearsons_r = (sd_x/sd_y)*beta_coefficient
pearsons_r

##         x 
## 0.8000443

Pearson’s ist übrigens das gleiche wie der standardisierte Regressionskoeffizient.

11.2 Spearman’s rho

Die Spearman Korrelation ist geeignet, wenn eine Variable ordinalskaliert ist und die zweite ordinal- oder intervallskaliert.

11.2.1 Voraussetzungen Spearman Korrelation

Die Spearman Korrelation wird auch Spearman Rangkorrelation genannt, weil sie nicht mit den Standardabweichungen (respektive Varianzen) arbeitet, sondern mit den Rängen.

• Beispiele Korrelation

Schauen wir uns ein paar Beispiele an.

Zuerst simulieren wir drei korrelierte Variablen:

data <- rnorm_multi(n = 100, 
                  mu = c(0, 20, 20),
                  sd = c(1, 5, 5),
                  r = c(0.9, 0.5, 0.1), 
                  varnames = c("A", "B", "C"),
                  empirical = FALSE)

Eine einfache Variante die Korrelationen aller Variablen eines Datensatzes darzustellen, ist eine Korrelationsmatrix:

psych::pairs.panels(data)

Abbildung 11.1: Matrix mit Korrelationskoeffizienten und Streudiagrammen. Wir sehen im unteren linken Dreieck die Streudiagramme und im oberen rechten Dreieck die Korrelationskoeffizienten (hier Pearson’s r). Wir sagen, dass die Korrelation zwischen Variable A und B 0.88 ist. Das ist eine hohe Korrelation. Zwischen der Variable A und der Variable C ist die Korrelation 0.49, wir könnten sagen, das sei eine mittlere Korrelation. Je näher die Korrelation bei 0 ist, desto kleiner ist die Korrelation. Null bedeutet: Keine Korrelation vorhanden. Die Korrelation kann auch negativ sein. Die rote Linie ist eine Anpassungslinie, sie hilft uns zu sehen, ob der Zusammenhang linear ist (wenn die rote Linie eine gerade Linie ist) oder nicht-linear (wenn die rote Linie gekurvt ist).

Nun fügen wir noch ein paar Variablen hinzu.

data<-data %>% 
  mutate(D=A*-1, 
         E= B^2, 
         F=20*A^2+3*A+20, 
         G=3)

Im nächsten Bild sehen wir, dass zwischen der Variable A und der Variable B eine perfekte negative Korrelation besteht (Korrelation = - 1). Wenn der Zusammenhang nicht linear ist, zum Beispiel zwischen der Variable E und F, kann der Korrelationskoeffizient den Zusammenhang nicht gut genug quantifizieren. Wir sehen auch, dass R eine Warnung ausgibt. Die Variable G hat keine Variabilität, da jeder Beobachtung in der Variable G 3 ist. Deswegen kann R die Korrelation nicht berechnen - und wir sehen in der Graphik fehlende Werte (NA) anstelle der Korrelationskoeffizienten.

psych::pairs.panels(data)

## Warning in cor(x, y, use = "pairwise", method = method): the standard deviation
## is zero

## Warning in cor(x, y, use = "pairwise", method = method): the standard deviation
## is zero

## Warning in cor(x, y, use = "pairwise", method = method): the standard deviation
## is zero

## Warning in cor(x, y, use = "pairwise", method = method): the standard deviation
## is zero

## Warning in cor(x, y, use = "pairwise", method = method): the standard deviation
## is zero

## Warning in cor(x, y, use = "pairwise", method = method): the standard deviation
## is zero

## Warning in cor(x, y, use = "pairwise", method = method): the standard deviation
## is zero

## Warning in cor(x, y, use = "pairwise", method = method): the standard deviation
## is zero

## Warning in cor(x, y, use = "pairwise", method = method): the standard deviation
## is zero

## Warning in cor(x, y, use = "pairwise", method = method): the standard deviation
## is zero

## Warning in cor(x, y, use = "pairwise", method = method): the standard deviation
## is zero

## Warning in cor(x, y, use = "pairwise", method = method): the standard deviation
## is zero

## Warning in min(diff(breaks)): no non-missing arguments to min; returning Inf

Abbildung 11.2: Matrix mit Korrelationskoeffizienten und Streudiagrammen. Nicht alle Zusammenhänge sind linear. Die Korrelationskoeffizienten im oberenn rechten Dreieck sind nur für lineare Zusammenhänge wirklich aussagekräftig. Die Variable G hat keine Variabilität, d.h. die Standardabweichung ist 0, deswegen kann die Korrelation nicht gerechnet werden. Wir könnten auch sagen, dass die Korrelation 0 ist.

11.3 Wann ist eine Korrelation 0?

Die Korrelation ist 0, falls die Regressionslinie keine Steigung aufweist.

par(mar = c(4, 4, .1, .1))
set.seed(1234)
id<-1:10000
variable_a=rnorm(10000,100,40)
variable_b=rnorm(10000, 50, 0.000000000000002)
plot(variable_a, variable_b, ylim=c(0,100))
abline(lm(variable_b ~ variable_a), col="red", lwd=3)
round(cor(variable_a, variable_b),2)

## [1] 0.01

variable_c=rnorm(10000, 50, 10)

plot(variable_a, variable_c, ylim=c(0,100))
abline(lm(variable_c ~ variable_a), col="red", lwd=4)

round(cor(variable_a, variable_c),2)

## [1] -0.01

Abbildung 11.3: Zwei Stichproben aus einer Population in der die Korrelationen zwischen Variable A und B, sowie zwischen A und C 0 ist (d.h. es gibt keine Korrelation in der Population). In beiden Graphiken ist die Korrelation praktisch 0, da die Regressionslinie praktisch keine Steigung aufweist. Wie immer bei zufälligen Daten: in der Stichprobe ist hier die Korrelation nicht exakt 0 und die Regressionslinie hat eine ganz kleine Steigung, da bei einer Stichprobe die Werte ja einen Stichprobenfehler aufweisenn (“sampling error”).

11.4 Wann ist eine Korrelation 1 oder -1?

Eine Korrelation ist 1 oder -1, wenn alle Punkte auf einer Linie sind und die Linie eine Steigung oder ein Gefälle hat. Der Grad der Steigung spielt keine Rolle, so lange es eine Steigung oder ein Gefälle hat.

par(mfrow=c(3,2))
set.seed(1234)

id<-1:1000
Variable_A <-rnorm(1000, 100, 20)
Variable_B <-Variable_A
plot(Variable_A, Variable_B)
abline(lm(Variable_B~Variable_A), col="red", lwd=2)
text(40,100,labels=paste0("r = ",cor(Variable_A, Variable_B)))

Variable_C <-200-Variable_A
plot(Variable_A, Variable_C)
abline(lm(Variable_C~Variable_A), col="red", lwd=2)
text(40,100,labels=paste0("r = ",cor(Variable_A, Variable_C)))

Variable_D<-0.1*Variable_A+0.0001
plot(Variable_A, Variable_D, ylim=c(0,200))
abline(lm(Variable_D~Variable_A), col="red", lwd=2)
text(40,100,labels=paste0("r = ",cor(Variable_A, Variable_D)))

Variable_E<- 100+ (-0.1*Variable_A+0.0001)
plot(Variable_A, Variable_E, ylim=c(0,200))
abline(lm(Variable_E~Variable_A), col="red", lwd=2)
text(40,140,labels=paste0("r = ",cor(Variable_A, Variable_E)))


Variable_F<- rnorm(1000, 50, 2)
plot(Variable_A, Variable_F, ylim=c(0,200))
abline(lm(Variable_F~Variable_A), col="red", lwd=2)
text(40,140,labels=paste0("r = ",round(cor(Variable_A, Variable_F),2)))


Variable_G<- rnorm(1000, 50, 30)
plot(Variable_A, Variable_G, ylim=c(0,200))
abline(lm(Variable_G~Variable_A), col="red", lwd=2)
text(40,140,labels=paste0("r = ",round(cor(Variable_A, Variable_G),2)))

Abbildung 11.4: Vier Streudiagramme mit perfekten Korrelationen von 1 oder -1. r = Pearson’s r, d.h. Pearson Korrelationskoeffizient. Die untersten beiden Streudiagramme zeigen zwei Regressionslinien, die praktisch keine Steigung aufweisen, deswegen ist die Korrelation beinahe 0, obschon die Streuung nicht sehr gross ist.

11.5 Nicht-lineare Zusammenhänge

11.5.1 Motone Zusammenhänge

Nicht immer können wir in einem Streudiagramm eine gerade Linie durch die Punkte ziehen. Machnmal ist der Zusammenhang nicht linear, aber er die eine Variable nimmt stetig zu, wenn die andere grösser wird. Dies nennt man positiver monotoner Zusammenhang. Ein negativ monotoner Zusammenhang bedeutet, dass, wenn eine Variable grösser wird, die andere kleiner wird. Beispiel für positiv monoton: mehr Kraft in den Beinen, höhere Gehgeschwindigkeit.

Für einen negativen monotonen Zusammenhang schauen wir uns ein echtes Beispiel an:

Für dieses Beispiel nutzen wir Daten vom London Marathon von dieser Webseite (hier klicken). Mehr Informationen finden Sie hier dazu.

Beispiel für negativ monotoner Zusammenhang: mehr Trainingskilometer, weniger Laufzeit beim Marathon.

Bei diesem Beispiel sehen wir, dass der Zusammenhang nicht linear ist. Der Zusammenhang ist jedoch negativ monoton. Was heist das: Je mehr Kilometer man trainiert pro Woche, desto weniger lange hat man um den Marathon zu laufen. Die Zeit nimmt mit zunehmenden Trainingskilometer stetig ab.

ggplot(data=data2, aes(x=training_kilometer_per_week, y=marathon_time_hours))+
  geom_point(aes(colour=marathon_time_hours*training_kilometer_per_week), show.legend = FALSE)+
  scale_y_continuous(breaks=seq(0,10,by=0.5))+
  theme_classic()+
  geom_smooth(method="loess", se=FALSE, colour="red")+
  labs(x="Training per week in kilometers", y="Marathon time in hours")

## `geom_smooth()` using formula = 'y ~ x'

Abbildung 11.5: Scatterplot and smoothed loess line (red) for the training kilometers peer week and the time in hours at the London Marathon. We see non-linear, but negative monotonic association. Data from: https://blog.scottlogic.com/2017/02/28/london-marathon-training-visualisation.html.

Schauen wir uns einmal die Zusammenhänge zwischen den Variablen A bis F an.

psych::pairs.panels(data[,c("A", "B", "E", "F")])

Abbildung 1.11: Korrelationsmatrix zwischen Variablen A, B, E und F

Wir sehen, dass alle Streudiagramme mit der Variable F nicht-lineare Zusammenhänge aufzeigen.

Ausserdem sehen wir einen nicht-linearen Zusammenhang zwischen der Variable B und E, was uns ja nicht erstaunt, da E ja das Quadrat von B ist.

11.5.2 Quantifizierung der Zusammenhänge

Für den Zusammenhang A und F sollten wir keinen Korrelationskoeffizienten berechnen, da die Korrelationskoeffizienten diesen stark nicht-linearen und nicht monotonen Zusammenhang nicht quantifizieren können. Was heisst monotoner Zusammenhang: In einem positiven monotonen Zusammenhang bewegen sich die Werte der beiden Variablen in die gleiche Richtung, d.h. wenn eine Variable einen höheren Wert zeigt, hat die andere Variable auch einen höheren Wert. Der Zusammenhang muss jedoch nicht linear sein. Ein quadratischer Zusammenhang wie zwischen B und E ist monoton, aber nicht linear. Die Variable F hat mit keiner anderen Variable einen monotonen Zusammenhang.

Wenn wir intervallskalierte Daten haben, und die beiden Variablen linear miteinander verbunden sind (wenn man eine gerade Linie durch die Punkte verbinden kann), dann können wir die Pearson’s Produkt-Moment-Korrelation berechnen (Pearson product-moment correlation coefficient). Wenn eine der beiden Variablen oder beide nur ordinal skaliert sind, oder der Zusammenhang nicht linear, aber immer noch monoton (d.h. stetig steigend oder sinkend), dann benutzen wir die Spearman Korrelation. In unserem Beispiel können wir die Pearson Korrelation berechnen.

Für die Korrelation zwischen den Variablen B und E würden wir also eher den Spearmanschen Korrelationskoeffizienten berechnen. Diese Korrelation wird auch Rangkorrelation nach Spearman genannt.

Berechnen wir einmal die Pearson Korrelation und die Rangkorrelation nach Spearman:

Zuerst Pearson’s r:

cor.test(data$B, data$E, method="pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  data$B and data$E
## t = 57.981, df = 98, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.9788359 0.9903966
## sample estimates:
##       cor 
## 0.9857354

Wir sehen, dass die Korrelation (Pearson’s r) 0.983 ist.

Spearman’s rho:

cor.test(data$B, data$E, method="spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  data$B and data$E
## S = 0, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho 
##   1

Das Spearman’s rho ist 1.

Das Pearson’s r ist also etwas kleiner als das Spearman’s rho. Welches ist nun korrekter? Da wir ja den Zusammenhang deterministisch erstellt haben, also ohne einen zufälligen Fehler hinzuzufügen (\(E = B^{2}\)), haben wir ja in Wirklichkeit einen perfekten Zusammenhang: Kennen wir E, kennen wir auch B, ohne Fehler. Deswegen ist hier die Spearman Rangkorrelation die bessere Wahl.

Wie sieht es für die beiden Korrelationen zwischen A und B aus?

Pearson’s r:

cor.test(data$A, data$B, method="pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  data$A and data$B
## t = 18.692, df = 98, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.8316511 0.9203753
## sample estimates:
##       cor 
## 0.8837134

und Spearman’s rho:

cor.test(data$A, data$B, method="spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  data$A and data$B
## S = 20076, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##      rho 
## 0.879532

Da wir hier einen linearen Zusammenhang haben und die Variablen A und B intervallskaliert sind, würden wir hier eher den Pearson Korrelationskoeffizienten wählen. Wir sollten aber diese Wahl treffen, bevor wir die Resultate anschauen.

Wie würde es bei einem exponentiellen Zusammenhang aussehen?

data<-data %>% 
  mutate(H=2^B)

Hier das Streudiagramm mit der logarithmierten Y-Achse.

ggplot(data, aes(x=B, y=H))+
  geom_point()+
  theme_classic()+
  scale_y_continuous(trans='log10')

Abbildung 11.6: Streudiagramm einer exponentiellen Steigung. Der Zusammenhang erscheint hier nur linear, da die Y-Achse logarithmiert wurde.

Hier mit der nicht transformierten Achse:

ggplot(data, aes(x=B, y=H))+
  geom_point()+
  theme_classic()

Abbildung 11.7: Streudiagramm einer exponentiellen Steigung.

Hier die beiden Korrelationskoeffizienten.

Pearson’s r:

cor.test(data$B, data$H, method="pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  data$B and data$H
## t = 4.2774, df = 98, p-value = 4.4e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.2171384 0.5501939
## sample estimates:
##       cor 
## 0.3966399

Spearman’s rho:

cor.test(data$B, data$H, method="spearman")

## 
##  Spearman's rank correlation rho
## 
## data:  data$B and data$H
## S = 0, p-value < 2.2e-16
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
## rho 
##   1

Hier sehen wir einen extremen Unterschied. Welcher Wert ist hier korrekter? Wir haben wieder einen deterministischen Zusammenhang ohne zufälligen Fehler, das heisst, wir haben einen perfekten Zusammenhang. Deswegen ist hier der Rangkorrelationskoeffizient nach Spearman die bessere Wahl.

11.6 Vor dem Rechnen von Korrelationskoeffizienten immer die Daten visualisieren

Schauen wir uns einmal folgendes reales Beispiel an.

Die Daten sind mit WebplotDigitizer aus dem folgenden Artikel Fantozzi, S., Cortesi, M., Giovanardi, A., Borra, D., Di Michele, R., & Gatta, G. (2020). Effect of walking speed during gait in water of healthy elderly. Gait & Posture, 82, 6-13. ausgelesen worden (Aus Figure 4C).

Wir haben einerseits die dimensionslose Gehgeschwindigkeit, die wie folgt berechnet wurde: \(\frac{walking speed}{\sqrt{g*leg_length}})\)., anderseits die Schrittdistanz in Zentimetern.

Es wurden junge und ältere Personen gemessen. Das Beispiel soll illustrieren, dass man nicht immer Korrelationen über zwei Gruppen berechnen sollte. In diesem Fall ist die Korrelation über beide Gruppen viel höher als die Korrelation in den beiden Gruppen. Diese hohe Korrelation über alle Personen ist nicht korrekt, dass wir nicht eine homogene Gruppe haben, sondern zwei, die sich deutlich in der Schrittlänge unterscheiden.

data<-rio::import("Fantozzi_Figure_4_c_walking_speed_stride_distance.csv")
names(data)

## [1] "dimensionless_walking_speed" "Stride_distance_cm"         
## [3] "age_group"

ggplot(data=data, aes(x=dimensionless_walking_speed, y=Stride_distance_cm))+
  geom_point(aes(colour=age_group))+
  geom_smooth(aes(colour=age_group), method="lm")+
  geom_smooth(method="lm")+
  ggpubr::stat_cor(label.y=120, color="blue")+
  ggpubr::stat_cor(aes(colour=age_group))+
  theme_classic()

## Registered S3 methods overwritten by 'broom':
##   method            from  
##   tidy.glht         jtools
##   tidy.summary.glht jtools

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Abbildung 11.8: Streudiagramm zwischen der dimensionslosen Gehgeschwindigkeit und der Schrittdistanz bei jungen und älteren Personen. Wir sehen, dass die Korrelation über alle Personen deutlich höher ist, als die Korrelationen in den Subgruppen. In solchen Situationen sollten wir nicht Korrelationen über alle berechnen.

Im nächsten Beispiel sollten wir auch keine Korrelation über beide Gruppen rechnen.

data<-rio::import("Fantozzi_Figure_2_D.csv")
names(data)

## [1] "dimensionless_walking_speed"                  
## [2] "knee_flexion_extension_at_foot_strike_degrees"
## [3] "group"

Wir haben hier die dimensionslose Gehgeschindigkeit und die Kniewindkel beim Aufsetzen des Fusses bei zwei Bedingungen, d.h. Gehen im Wasser und Gehen an Land. Bei unserem Beispiel vernachlässigen wir die Tatsache, dass wir hier nicht unabhängige Gruppen haben (nicht unabhängig weil die gleichen Personen an Land und im Wasser gehen).

ggplot(data=data, aes(x=dimensionless_walking_speed, y=knee_flexion_extension_at_foot_strike_degrees))+
  geom_point(aes(colour=group))+
    geom_smooth(aes(colour=group), method="lm")+
  geom_smooth(method="lm")+
  ggpubr::stat_cor(label.y=-10, color="blue")+
  ggpubr::stat_cor(aes(colour=group))+
  theme_classic()

## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'

Abbildung 11.9: Streudiagramm mit dimensionslosen Gehgeschwindigkeit und Kniewinkel beim Aufsetzen des Fusses. Es wurde das Gehen im Wasser und das Gehen an Land untersucht. Wir sehen, dass die Korrelation über beide Konditionen negativ ist, innerhalb der beiden Konditionen jedoch eine positive Korrelation besteht.

11.7 Reduzierte Streuung führt zu einer kleineren Korrelation

Wenn wir einen reduzierten Bereich auswählen, kann dies die Streuung reduzieren, was zu einer reduzierten Korrelation führt.

Für das erste Beispiel benutzen wir einen Datensatz, der öffentlich zugänglich ist. Der veröffentlichte Artikel zu diesen Daten finden Sie mit diesem Link (hier klicken).

Informationen zu diesen Daten finden Sie mit diesem Link.

data<-rio::import("https://research-data.westernsydney.edu.au/default/rdmp/pubrecord/a28a55a0519311ecb15399911543e199/pubattach/be4b3a8e7e4a418e8da61f03775c9146?pubId=160811c0519411ecb15399911543e199", format="xlsx")
#skimr::skim(data)

Wir schauen uns die Korrelation zwischen dem Alter und dem Score in der Short Physiocal Performance Battery an.

Eine Beschreibung des Short Physical Performance Battery finden Sie hier (hier klicken).

Zuerst die Korrelation über den ganzen Altersbereich:

age_sppb<-data %>% 
  select(age, t0_SPPB_score)

ggplot(data=age_sppb, aes(x=age, y=t0_SPPB_score))+
  geom_point()+
  geom_smooth(method="lm")+
  theme_classic()+
  ggpubr::stat_cor()

## `geom_smooth()` using formula = 'y ~ x'

Abbildung 11.10: Streudiagramm mit Fehler: Wir sehen, dass die Variable Alter noch fehlerhafte Werte hatten. Wir müssen die Graphik neu machen ohne diese Werte.

• Wir sehen, dass wir Fehler in der Variable age haben, Alter 0 kann ja nicht sein. Wir schliessen diese Werte aus der Analyse aus.

ggplot(data=age_sppb[age_sppb$age>0,], aes(x=age, y=t0_SPPB_score))+
  geom_point()+
  geom_smooth(method="lm")+
  theme_classic()+
  ggpubr::stat_cor()

## `geom_smooth()` using formula = 'y ~ x'

Abbildung 11.11: Streudiagramm mit ganzem Altersbereich. Korrelation ist -0.4.

Nun wählen wir nur Personen im Alter zwischen 70 und 80 aus und wir sehen in der unteren Graphik, dass die Korrelation kleiner wird (d.h. näher bei Null ist). Dies ist der Fall, weil wir die Variabilität im Alter reduziert haben. Der Korrelationskoeffizient ist hier nicht falsch, wir müssen einfach in der Interpretation aufpassen, dass wir nicht falsche Schlüsse ziehen.

ggplot(data=age_sppb[age_sppb$age>70 & age_sppb$age<80,], aes(x=age, y=t0_SPPB_score))+
  geom_point()+
  geom_smooth(method="lm")+
  theme_classic()+
  ggpubr::stat_cor()

## `geom_smooth()` using formula = 'y ~ x'

Abbildung 11.12: Streudiagramm nur mit Personen zwischen 70 und 80 Jahren. Da wir jetzt weniger Streuung in der Variable Alter haben, ist auch die Korrelation geringer.

11.8 Weitere Beispiele aus der Literatur

Wir werden für dieses Beispiel die Variablen für den Trail Making Test und für die Short Physical Performance Battery benutzen. Informationen zum Trail Making Test finden sie hier (hier klicken) oder hier (hier klicken).

Eine Beschreibung des Short Physical Performance Battery finden Sie hier (hier klicken).

Zusätzlich wählen wir noch die Anzahl Stürze aus.

df<-data %>% 
  select(Participant_ID,contains("TMT"), contains("SPPB"), number_falls, age)

Hier ein grober Überblick über alle Variablen, die wir ausgewählt haben.

psych::pairs.panels(df)

Abbildung 11.13: Korrelationsmatrix zwischen allen Variablen

Zuerst möchten wir etwas genauer den Zusammenhang zwischen den Resultaten der beiden Subtests des Trail Making Tests untersuchen (A und B).

Im Teil A muss die Versuchsperson eine Spur durch eine Reihe von Kreisen zu ziehen, die von 1-25 nummeriert sind. Im Teil B muss die Versuchsperson eine Spur durch Zahlen und Buchstaben ziehen, und zwar so, dass sie zuerst zur 1 zeigt, dann zu A, danach zu 2, und weiter zur B, und so weiter.

Es wird die Zeit gemessen, die die Versuchsperson benötigt, um den Test zu beenden. Tiefere Werte sind bessere Werte.

Eine detaillierte Beschreibung der Instruktionen finden Sie hier.

Wir erstellen ein Streudiagramm und sehen, wie die beiden Variablen zusammen variieren. Jeder Punkt entspricht einer Person. Sie können mit der Maus auf irgendeinen Punkt zeigen, dann sehen Sie welche Person es ist, und was die Werte in den beiden Tests sind. Wir sehen, dass bei höheren Werten im Test A die Werte im Test B auch tendenziell höher sind.

Wir sehen also, dass es eine positive Korrelation zwischen den beiden Variablen gibt.

names(df)

## [1] "Participant_ID" "t0_TMT_A"       "t0_TMT_B"       "t12_TMT_A"     
## [5] "t12_TMT_B"      "t0_SPPB_score"  "t12_SPPB_score" "number_falls"  
## [9] "age"

df<-df %>% 
  mutate(show=paste0("ID= ",Participant_ID,";","Test A=",t0_TMT_A, "; Test B=",t0_TMT_B ))


gg_point=ggplot(df, aes(x=t0_TMT_A, y=t0_TMT_B))+
  geom_point_interactive(aes(tooltip = show ), colour="#AD50B5", alpha=0.3)+
  theme_classic()

girafe(ggobj = gg_point)

Abbildung 11.14: Interaktives Streudiagramm. Sie können einen Punkt anwählen, um die Werte zu sehen

Wir können jetzt die Korrelation auch noch quantifizieren. Wenn wir intervallskalierte Daten haben, und die beiden Variablen linear miteinander verbunden sind (wenn man eine gerade Linie durch die Punkte verbinden kann), dann können wir die Pearson’s Produkt-Moment-Korrelation berechnen (Pearson product-moment correlation coefficient). Wenn eine der beiden Variablen oder beide nur ordinal skaliert sind, oder der Zusammenhang nicht linear, aber immer noch monoton (d.h. stetig steigend oder sinkend), dann benutzen wir die Spearman Korrelation. In unserem Beispiel können wir die Pearson Korrelation berechnen.

cor.test(df$t0_TMT_A, df$t0_TMT_B)

## 
##  Pearson's product-moment correlation
## 
## data:  df$t0_TMT_A and df$t0_TMT_B
## t = 17.788, df = 528, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.5559430 0.6627472
## sample estimates:
##       cor 
## 0.6121292

Wir können auch durch das Streudiagramm eine Regressionslinie zeichnen.

ggplot(df, aes(x=t0_TMT_A, y=t0_TMT_B))+
  geom_point()+
  geom_smooth(method="lm")+
  theme_classic()

## `geom_smooth()` using formula = 'y ~ x'

Abbildung 11.15: Streudiagramm mit Regressionslinie

Übrigens ist die Korrelation auch in einer linearen Regression ersichtlich. Der standardisierte Koeffizient einer Regression entspricht dem Pearson Korrelationskoeffizient. Der normale Output einer linearen Regrssion zeigt jedoch den unstandardisierten Koeffizienten. Der Koeffizient ist in der Einheit der Variablen, d.h. für eine Einheit mehr in der unabhängigen Variable, ändert sich die abhängige Varialbe um den Faktor des Koeffizienten. Der Koeffizient entspricht der Steigung der Regressionslinie (siehe Graphik oben).

lm_fit<-lm(t0_TMT_A~t0_TMT_B, data=df)
summary(lm_fit)

## 
## Call:
## lm(formula = t0_TMT_A ~ t0_TMT_B, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -82.334  -7.763  -1.605   6.196  79.843 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 24.149099   1.183444   20.41   <2e-16 ***
## t0_TMT_B     0.149271   0.008392   17.79   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 13.25 on 528 degrees of freedom
## Multiple R-squared:  0.3747, Adjusted R-squared:  0.3735 
## F-statistic: 316.4 on 1 and 528 DF,  p-value: < 2.2e-16

Es gibt verschiedene Möglichkeiten, sich die standardisierten Koeffizienten anzeigen zu lassen. Wir benutzen hier das Paket jtools und den Befehl summ. Zuerst die Ausgabe unstandardisiert.

jtools::summ(lm_fit, scale=FALSE,transform.response = FALSE, confint=TRUE)

Observations	530
Dependent variable	t0_TMT_A
Type	OLS linear regression

F(1,528)	316.40
R²	0.37
Adj. R²	0.37

	Est.	2.5%	97.5%	t val.	p
(Intercept)	24.15	21.82	26.47	20.41	0.00
t0_TMT_B	0.15	0.13	0.17	17.79	0.00
Standard errors: OLS

Und hier standardisiert

jtools::summ(lm_fit, scale=TRUE,transform.response = TRUE, confint=TRUE)

Observations	530
Dependent variable	t0_TMT_A
Type	OLS linear regression

F(1,528)	316.40
R²	0.37
Adj. R²	0.37

	Est.	2.5%	97.5%	t val.	p
(Intercept)	-0.00	-0.07	0.07	-0.00	1.00
t0_TMT_B	0.61	0.54	0.68	17.79	0.00
Standard errors: OLS; Continuous variables are mean-centered and scaled by 1 s.d.

Wir sehen jetzt, dass der standardisierte Koeffizient der Korrelation entspricht.

• Achtung: Damit der standardisierte Koffizient dem Korrelationskoeffizienten entspricht, muss nicht nur die unabhängige Variable standardisiert werden (was wir mit der option scale=TRUE machen), sondern auch die abhängige Variable (was wir mit der Option transform.response = TRUE erreichen.). Siehe dazu (hier klicken)

Wenn wir nur eine Variable im Modell haben, können wir auch die Wurzel des \(R^{2}\) nehmen, dies entspricht auch der Korrelation zwischen der abhängigen und unabhängigen Variable. In unserem Beispiel ist \(R^{2}\) 0.3747, die Wurzel davon ist 0.6121274

.

11.9 Korrelation zwischen Trail Making Test und SPPB Score

Wir sehen, dass wir hier eine sinkende Regressionslinie haben. Das bedeutet, dass mit steigendem Trail Making B Testresultat der Score in der Short Physical Performance Battery sinkt. Zur Erinnerung: Höhere Werte im Trail Making Test deuten auf eine schlechtere kognitive Funktionsfähigkeit hin.

ggplot(data=df, aes(x=t0_TMT_B, y=t0_SPPB_score))+
  geom_point()+
  geom_smooth(method="lm")+
  theme_classic()

## `geom_smooth()` using formula = 'y ~ x'

Abbildung 11.16: Streudiagramm mit Regressionslinie zwischen Trail Making Test B und dem Score in der Short Physical Performance Battery

Für die Berechnung des Korrelationskoeffizienten müssen wir zwischen dem Pearson und dem Spearman Korrelationskoeffizient auswählen. Beide wären hier möglich, wir könnten jedoch argumentieren, dass der SPPB Score eher ordinalskaliert ist und nicht intervallskaliert. Berechnen wir doch beide Koeffizienten:

cor.test(df$t0_TMT_A, df$t0_SPPB_score)

## 
##  Pearson's product-moment correlation
## 
## data:  df$t0_TMT_A and df$t0_SPPB_score
## t = -9.6948, df = 528, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.4587113 -0.3139517
## sample estimates:
##        cor 
## -0.3887279

cor.test(df$t0_TMT_A, df$t0_SPPB_score, method='spearman')

## Warning in cor.test.default(df$t0_TMT_A, df$t0_SPPB_score, method = "spearman"):
## Cannot compute exact p-value with ties

## 
##  Spearman's rank correlation rho
## 
## data:  df$t0_TMT_A and df$t0_SPPB_score
## S = 33152235, p-value = 1.839e-15
## alternative hypothesis: true rho is not equal to 0
## sample estimates:
##       rho 
## -0.336097

Wir sehen, dass die beiden Korrelationskoeffizienten in diesem Beispiel sehr ähnlich sind.

Jetzt interessieren wir uns für den Zusammenhang zwischen dem Alter und der Short Physical Performance Battery.

ggplot(data=df, aes(x=age, y=t0_SPPB_score))+
  geom_point()+
  geom_smooth(method="loess")+
  theme_classic()

## `geom_smooth()` using formula = 'y ~ x'

Abbildung 11.17: Streudiagramm SPPB Score über das Alter. Wir sehen, dass es fehlerhafte Altersdaten gibt, die wir korrigieren müssen.

Wir sehen, dass es Fehler in den Daten haben muss, da wir zwei Personen mit Alter 0 haben. Wir werden diese zuerst ausschliessen.

df<-df %>% 
  mutate(age=ifelse(age<18, NA, age))

Jetzt plotten wir die Daten noch einmal und zeichnen eine Regressionslinie in rot und eine angepasste Linie (blau) ein. Wir sehen, dass wahrscheinlich der Zusammenhang zwischen dem Alter und dem SPPB-Score nicht linear ist.

ggplot(data=df, aes(x=age, y=t0_SPPB_score))+
  geom_point()+
  geom_smooth(method="loess")+
  geom_smooth(method="lm", color="red")+
  theme_classic()

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 2 rows containing non-finite values (`stat_smooth()`).

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 2 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 2 rows containing missing values (`geom_point()`).

Abbildung 11.18: Streudiagramm mit Regressionslinie (rot) und Anpassungslinie (blau)

Wenn wir das Streudiagramm anschauen, sind wir nicht ganz sicher, ob der Zusammenhang linear ist. Wie können wir das testen?

Wir können unterschiedliche Modelle vergleichen.

• Normales Modell
• Modell mit Alter plus Alter im Quadrad (dafür erstellen wir noch die Variable age2 indem wir das Alter quadrieren)
• Modell mit Polynom-Funktion zweiten Grades (was im Prinzip das gleiche ist wie im Modell 2) Siehe hier.
• Modell mit Splines mit drei Knoten Für eine illustration von Splines und Knoten siehe hier und für eine Diskussion über die Anzahl Knoten siehe hier.

df<-df %>% 
  filter(!is.na(age))
df<-df %>% 
  mutate(age2=age^2)

Einfaches lineares Modell.

m1<-lm(t0_SPPB_score~age,data=df)
summary(m1)

## 
## Call:
## lm(formula = t0_SPPB_score ~ age, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.0011 -0.8014  0.3923  1.1956  2.8021 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 18.234820   0.785487   23.21   <2e-16 ***
## age         -0.100411   0.009987  -10.05   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.625 on 526 degrees of freedom
## Multiple R-squared:  0.1612, Adjusted R-squared:  0.1596 
## F-statistic: 101.1 on 1 and 526 DF,  p-value: < 2.2e-16

df<-cbind(df, residm1 = resid(m1), fittedm1 = fitted(m1))

m2<-lm(t0_SPPB_score~age+age2, data=df)
summary(m2)

## 
## Call:
## lm(formula = t0_SPPB_score ~ age + age2, data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.1557 -0.8535  0.3553  1.0198  3.2832 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -10.515365   7.030563  -1.496 0.135342    
## age           0.645637   0.181597   3.555 0.000412 ***
## age2         -0.004799   0.001167  -4.114 4.51e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.6 on 525 degrees of freedom
## Multiple R-squared:  0.1874, Adjusted R-squared:  0.1843 
## F-statistic: 60.54 on 2 and 525 DF,  p-value: < 2.2e-16

anova(m1,m2)

## Analysis of Variance Table
## 
## Model 1: t0_SPPB_score ~ age
## Model 2: t0_SPPB_score ~ age + age2
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    526 1388.2                                  
## 2    525 1344.8  1    43.361 16.928 4.508e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

df<-cbind(df, residm2 = resid(m2), fittedm2 = fitted(m2))

Polinomiale Regerssion mit Alter und Polynom zweiten Grades, entspricht im Prinzip dem Model m2

m3<-lm(t0_SPPB_score~poly(age,2), data=df)
summary(m3)

## 
## Call:
## lm(formula = t0_SPPB_score ~ poly(age, 2), data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.1557 -0.8535  0.3553  1.0198  3.2832 
## 
## Coefficients:
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    10.36932    0.06965 148.872  < 2e-16 ***
## poly(age, 2)1 -16.33374    1.60049 -10.205  < 2e-16 ***
## poly(age, 2)2  -6.58494    1.60049  -4.114 4.51e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.6 on 525 degrees of freedom
## Multiple R-squared:  0.1874, Adjusted R-squared:  0.1843 
## F-statistic: 60.54 on 2 and 525 DF,  p-value: < 2.2e-16

anova(m1,m2, m3)

## Analysis of Variance Table
## 
## Model 1: t0_SPPB_score ~ age
## Model 2: t0_SPPB_score ~ age + age2
## Model 3: t0_SPPB_score ~ poly(age, 2)
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1    526 1388.2                                  
## 2    525 1344.8  1    43.361 16.928 4.508e-05 ***
## 3    525 1344.8  0     0.000                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

df<-cbind(df, residm3 = resid(m3), fittedm3 = fitted(m3))

Modell mit Restricted Cubic Spline, dazu benutzen wir das rms Paket und die Funktion ols.

m4<-rms::ols(t0_SPPB_score~rcs(age,4),data=df)
df<-cbind(df, residm4 = resid(m4), fittedm4 = fitted(m4))
summary.lm(m4)

## 
## Call:
## rms::ols(formula = t0_SPPB_score ~ rcs(age, 4), data = df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.9818 -0.8315  0.2824  0.9893  3.1177 
## 
## Coefficients:
##           Estimate Std. Error t value Pr(>|t|)    
## Intercept  7.39749    2.53414   2.919 0.003661 ** 
## age        0.05511    0.03626   1.520 0.129119    
## age'      -0.31242    0.08671  -3.603 0.000345 ***
## age''      1.06088    0.39795   2.666 0.007916 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.593 on 524 degrees of freedom
## Multiple R-squared:  0.1965, Adjusted R-squared:  0.1919 
## F-statistic: 42.72 on 3 and 524 DF,  p-value: < 2.2e-16

Wir können uns alle Modelle in einer Graphik anschauen.

ggplot(data=df,aes(x=age))+
  geom_jitter(aes(y=t0_SPPB_score), alpha=0.3)+
  geom_line(aes(y=fittedm1), colour="red", alpha=1, linetype="dotted", size=1)+
  geom_line(aes(y=fittedm2), colour="pink", alpha=1, linetype="twodash", size=2)+
  geom_line(aes(y=fittedm3), colour="blue",linetype="dashed", alpha=0.5, size=1 )+
  geom_line(aes(y=fittedm4), colour="purple",linetype="twodash", alpha=1, size=1.4 )+
  theme_classic()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.

Abbildung 11.19: Rot = normale lineare Regression, Pink = Polinomiale Regression (Polynom 2. Grades), Blau = Modell mit Alter und Alter im Quadrad, Violett = Spline mit 3 Knoten.

Hier ein Video über Splines:

11.10 Zusammenhang zwischen Trail Making Test, Short Physical Performance Battery und der Anzahl Stürze (Count Data).

Als Übersicht, können wir uns die Korrelations- und Streudiagramm Matrix zwischen allen Variablen anschauen. Diese ist jedoch mit Vorsicht zu geniessen.

psych::pairs.panels(df[,2:9])

Abbildung 11.20: Korrelationsmatrix

Betrachten wir das Streudiagramm zwischen der Short Physical Performance Battery und der Anzahl Stürze genauer an.

ggplot(data=df, aes(x=t0_SPPB_score, y=number_falls))+
  geom_point()+
  geom_smooth(method="loess")+
  theme_classic()

## `geom_smooth()` using formula = 'y ~ x'

## Warning: Removed 7 rows containing non-finite values (`stat_smooth()`).

## Warning: Removed 7 rows containing missing values (`geom_point()`).

Abbildung 11.21: Streudiagramm zwischen SPPB und Anzahl der Stürze.

Diese Graphik macht nicht wahnsinnig viel Sinn. Am besten analysieren wir den Zusammenhang zwischen diesen beiden Variablen mit einer Count-Regression. Siehe hier für ein Beispiel einer Count-Regression.

pr_fit<-glm(number_falls~t0_SPPB_score, family=poisson, data=df)
summary(pr_fit)

## 
## Call:
## glm(formula = number_falls ~ t0_SPPB_score, family = poisson, 
##     data = df)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.8265  -1.2401  -1.1815   0.2518   7.9640  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    0.80206    0.25259   3.175   0.0015 ** 
## t0_SPPB_score -0.09681    0.02468  -3.923 8.74e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 1015.9  on 520  degrees of freedom
## Residual deviance: 1001.5  on 519  degrees of freedom
##   (7 observations deleted due to missingness)
## AIC: 1534.8
## 
## Number of Fisher Scoring iterations: 6

Wir erstellen hier ein paar Diagnosen, um zu evaluieren, ob die Poisson Regression sinnvoll ist: Falls Sie mehr über Diagnosen bei Poisson Regression wissen möchten, klicken Sie hier.

Da wir sehr viele Patienten mit Null Stürzen haben, sollten wir testen, ob ein Zero inflation Modell besser wäre.

summarytools::freq(df$number_falls)

## Frequencies  
## df$number_falls  
## Type: Numeric  
## 
##               Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
##           0    297     57.01          57.01     56.25          56.25
##           1    142     27.26          84.26     26.89          83.14
##           2     41      7.87          92.13      7.77          90.91
##           3     17      3.26          95.39      3.22          94.13
##           4      8      1.54          96.93      1.52          95.64
##           5      6      1.15          98.08      1.14          96.78
##           6      4      0.77          98.85      0.76          97.54
##           7      1      0.19          99.04      0.19          97.73
##           8      1      0.19          99.23      0.19          97.92
##          11      1      0.19          99.42      0.19          98.11
##          15      3      0.58         100.00      0.57          98.67
##        <NA>      7                               1.33         100.00
##       Total    528    100.00         100.00    100.00         100.00

Wir vergleichen den AIC Wert der beiden Modelle. Tiefe AIC sind besser. Wir sehen hier, dass das Zero Inflation Modell besser wäre.

library(pscl)
zir_fit <- zeroinfl(number_falls~t0_SPPB_score|t0_SPPB_score, data=df, dist="poisson")
AIC(pr_fit, zir_fit)

##         df      AIC
## pr_fit   2 1534.844
## zir_fit  4 1438.830

Wir betrachten, ob wir overdispersion haben, siehe dazu dieses Video:

library(AER)
deviance(pr_fit)/pr_fit$df.residual

## [1] 1.929699

dispersiontest(pr_fit)

## 
##  Overdispersion test
## 
## data:  pr_fit
## z = 2.7185, p-value = 0.003279
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion 
##   3.379106

Der Wert bei dispersion sollte nahe bei 0 sein. Hier ist er über 3, also sollten wir eher nicht eine Poisson Regression berechnen. Eigentlich wussten wir das schon vorher, da bei Stürzen meistens die negative binomiale Regression empfohlen wird.

nbr_fit <- glm.nb(number_falls ~ t0_SPPB_score, data = df)
summary(nbr_fit)

## 
## Call:
## glm.nb(formula = number_falls ~ t0_SPPB_score, data = df, init.theta = 0.6467214557, 
##     link = log)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -1.2874  -1.0063  -0.9723   0.1666   3.9067  
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)  
## (Intercept)    0.81450    0.41600   1.958   0.0502 .
## t0_SPPB_score -0.09803    0.03996  -2.453   0.0142 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Negative Binomial(0.6467) family taken to be 1)
## 
##     Null deviance: 460.99  on 520  degrees of freedom
## Residual deviance: 454.90  on 519  degrees of freedom
##   (7 observations deleted due to missingness)
## AIC: 1302.1
## 
## Number of Fisher Scoring iterations: 1
## 
## 
##               Theta:  0.6467 
##           Std. Err.:  0.0891 
## 
##  2 x log-likelihood:  -1296.0920

Wir könnten noch dieses Modell mit dem Zero-Inflation Modell vergleichen und wir sehen, dass das negative binomiale Modell besser ist.

AIC(nbr_fit, zir_fit)

##         df      AIC
## nbr_fit  3 1302.092
## zir_fit  4 1438.830

Wir könnten auch noch einen Schritt weiter gehen und ein zero inflated negative binomial Modell anschauen siehe hier{target=“_blank”. Oder auch hier und auch das nächste Kapitel hier.

df<-df %>% 
  filter(!is.na(number_falls))
zirnb_fit <- zeroinfl(number_falls~t0_SPPB_score |t0_SPPB_score, data=df, dist="negbin",start = list(count = c(0,0), zero = c(-0, -0))) # the number of starting values in the list does correspond to the number of variables in the model - same for the zero values. 
summary(zirnb_fit)

## 
## Call:
## zeroinfl(formula = number_falls ~ t0_SPPB_score | t0_SPPB_score, data = df, 
##     dist = "negbin", start = list(count = c(0, 0), zero = c(-0, -0)))
## 
## Pearson residuals:
##     Min      1Q  Median      3Q     Max 
## -0.6835 -0.5925 -0.5791  0.1789 11.8935 
## 
## Count model coefficients (negbin with log link):
##               Estimate Std. Error z value Pr(>|z|)   
## (Intercept)    0.81450    0.41843   1.947  0.05159 . 
## t0_SPPB_score -0.09803    0.04019  -2.439  0.01473 * 
## Log(theta)    -0.43584    0.13781  -3.163  0.00156 **
## 
## Zero-inflation model coefficients (binomial with logit link):
##               Estimate Std. Error z value Pr(>|z|)
## (Intercept)     -9.220   1008.638  -0.009    0.993
## t0_SPPB_score   -0.732    185.170  -0.004    0.997
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Theta = 0.6467 
## Number of iterations in BFGS optimization: 34 
## Log-likelihood:  -648 on 5 Df

AIC(nbr_fit, zirnb_fit)

##           df      AIC
## nbr_fit    3 1302.092
## zirnb_fit  5 1306.092

vuong(nbr_fit,zirnb_fit)

## Vuong Non-Nested Hypothesis Test-Statistic: 
## (test-statistic is asymptotically distributed N(0,1) under the
##  null that the models are indistinguishible)
## -------------------------------------------------------------
##               Vuong z-statistic             H_A p-value
## Raw                4.056070e-01 model1 > model2 0.34252
## AIC-corrected      5.228992e+04 model1 > model2 < 2e-16
## BIC-corrected      1.635555e+05 model1 > model2 < 2e-16

Wir sehen, dass das zeroinflated negative binomiale Modell das höhere AIC hat, wir wählen deshalb eher das normale negative binomiale Modell aus. Wir sehen das auch im Vuong Test.

Hier noch kurz zur Erklärung des zero inflated negative binomial models:

Der Teil, der als Count model coefficients bezeichnet wird, modelliert die Anzahl Stürze (wenn die Anzahl nicht 0 ist). Wir sehen, dass der SPPB Score einen signifikanten Einfluss auf die Anzahl Stürze hat.

Der zweite Teil, der als Zero-inflation model coefficients bezeichnet ist, zeigt uns, dass der SPPB Score nicht signifikant hilft vorherzusagen, wer zu den Personen mit 0 Stürzen gehört.

Ein Beispiel mit dieser Analyse kann hier gefunden werden (hier klicken).