Kapitel 22 Beispiele zum Kapitel “Mit Statistik Fragen beantworten und Hypothesen testen”
library(tidyverse)## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.2 ──
## ✔ ggplot2 3.4.0 ✔ purrr 0.3.4
## ✔ tibble 3.1.8 ✔ dplyr 1.0.10
## ✔ tidyr 1.2.1 ✔ stringr 1.5.0
## ✔ readr 2.1.3 ✔ forcats 0.5.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()library(kableExtra)##
## Attaching package: 'kableExtra'
##
## The following object is masked from 'package:dplyr':
##
## group_rowslibrary(describedata)Wir wollen uns einen neuen Computer kaufen. Sollen wir einen Mac oder einen PC kaufen? Wenn wir uns in der Klasse umschauen, benutzen etwa die Hälfte der Studierenden einen Mac, die anderen einen PC.
Die beste Studie wäre wohl eine randomisierte kontrollierte Studie, in der Studierenden zufällig ein PC oder ein Mac zugeteilt würde. Doch diese Studie wird es wohl aus kostengründen nicht geben.
Wir können aber eine einfache Befragung durchführen. Wie zufrieden seid ihr mit eurem Computer?
Wir haben das gemacht und schauen uns einmal die Analyse der Daten an.
data<-rio::import("mac_versus_pc_2021.csv")
data<-data %>%
dplyr::select(record_id, gender_1.factor, computer_work_1.factor, satisfaction_work_vas) %>%
filter(!is.na(satisfaction_work_vas))Zusammenfassung pro Computerart
# length2: New version of length which can handle NA's: if na.rm==T, don't count them;
# Dank an: http://www.cookbook-r.com/Manipulating_data/Summarizing_data/
length2 <- function (x, na.rm=FALSE) {
if (na.rm) sum(!is.na(x))
else length(x)
}
summary_table<-data %>%
group_by(computer_work_1.factor) %>%
summarize(n_satisfaction=length2(., na.rm=TRUE),
mean_satisfaction=round(mean(satisfaction_work_vas, na.rm=TRUE),1),
sd_satisfaction=round(sd(satisfaction_work_vas, na.rm=TRUE),1))
kbl(summary_table, caption='Mittelwert und Standardabweichung der Zufriedenheit mit dem Computer, aufgeteilt nach PC oder Mac.') %>% kable_styling()| computer_work_1.factor | n_satisfaction | mean_satisfaction | sd_satisfaction |
|---|---|---|---|
| Mac OS | 208 | 91.3 | 10.4 |
| Windows | 208 | 80.7 | 21.9 |
22.1 Hypothesen testen
Wir können jetzt die Nullhypothese (H0) testen, dass in der Population der Mittelwert der Zufriedenheit bei den PCs und den Macs gleich ist. \(H0: \mu PC = \mu Mac\) . Die Alternativhypothese wäre, das der Populationsmittelwert nicht gleich ist.
Die Statistiker:innen müssen nun überlegen, was der geeigneteste Test ist, d.h. bei welchem Test es nicht zu einer zu hohen Typ-I Fehlerhäufigkeit oder Typ-II Fehlerhäufigkeit kommt. Ganz grob gesagt gibt es zwei Gruppen von Tests: die parametrischen Tests und die nichtparametrischen Tests. Wenn wir mindestens intervallskalierte Daten haben, die in der Population normalverteilt sind oder wir einge genügend hohe Stichprobenzahl haben, können wir parametrische Tests auswählen. Diese Tests benutzen die Parameter \(\mu\) und \(\sigma\). Falls wir ordinale Daten haben oder die intervall- oder Ratioskalierten Daten in der Population nicht normalverteilt sind und wir eine zu kleine Stichproben haben, so benutzen wir besser nichtparametrische Tests, die zum Beispiel mit Rängen arbeiten. Man nennt diese Tests auch Verteilungsfreie Tests.
Nun ist es jedoch nicht so einfach herauszufinden, ob die Daten in der Population normalverteilt sind.
Es gibt statistische Tests dafür, doch diese sind nicht ganz unproblematisch (darauf gehen wir hier nicht ein). Wir sollten aber zumindest die Daten visuell überprüfen.
22.2 Visuelle Überprüfung der Verteilung
Schauen wir die Verteilung der Zufriedenheit an. Was häufig gemacht wird, ist die abhängige Variable zu überprüfen.
hist(data$satisfaction_work_vas)
Abbildung 3.2: Histogramm der Variable Zufriedenheit. Wenn wir die Daten hier sehen, glauben wir eher nicht, dass die Daten in der Population normalverteilt sind.
Besser wäre es, die Variable geschichtet nach Computer anzuschauen.
hist(data$satisfaction_work_vas[data$computer_work_1.factor=="Windows"],col=rgb(1,0,0,1/4), breaks=20)
hist(data$satisfaction_work_vas[data$computer_work_1.factor=="Mac OS"], col=rgb(0,0,1,1/4), breaks=20, add=TRUE)
Abbildung 5.5: Histogramm der Variable Zufriedenheit, geschichtet nach Computer. Wenn wir die Daten hier sehen, glauben wir eher nicht, dass die Daten in der Population normalverteilt sind.
Wir können eine noch umfassendere Analyse machen, ob unsere Analyse für diese Daten geeignet sind. Die nächsten Graphiken geben uns Hinweise über mögliche Probleme:
- In der Graphik Resdiuals vs Fitted sollte die rote Linie auf 0 und ganz horizontal liegen.
- Im Q-Q Plot sollten die Punkte auf der diagonalen Linie liegen - eine Abweichung am Ende der Linie is jedoch akzeptierbar - siehe auch weiter unten im Text. Eine Abweichung deutet auf eine Abweichung von der Normalverteilung hin.
- Die mit Scale-Location betitelte Graphik zeigt, ob die Residuen über alle Gruppen ähnliche Variabilität haben (“Homoskedastizität”, homoscedasticity). Die rote Linie sollte horizontal sein.
- Die Graphik Residuals vs Leverage zeigt uns den Einfluss der Ausreisser auf die Analyes. Wir schauen, ob es Datenpunkte ausserhalb der rot gestrichelten Linie hat.
Eine ausführlichere Erklärunge finden Sie hier (hier klicken)
fit.lm<-lm(satisfaction_work_vas~computer_work_1.factor, data=data)
plot(fit.lm)
Wir können auch noch statistische Tests durchführen für die Nullhypothese, dass die Residuen normalverteilt sind.
Wir müssen leider hier die Nullhypothese verwerfen.
library(olsrr)##
## Attaching package: 'olsrr'## The following object is masked from 'package:datasets':
##
## riversols_test_normality(fit.lm)## Warning in ks.test.default(y, "pnorm", mean(y), sd(y)): ties should not be
## present for the Kolmogorov-Smirnov test## -----------------------------------------------
## Test Statistic pvalue
## -----------------------------------------------
## Shapiro-Wilk 0.8108 0.0000
## Kolmogorov-Smirnov 0.1888 0.0491
## Cramer-von Mises 4.2537 0.0000
## Anderson-Darling 2.6946 0.0000
## -----------------------------------------------Bei grösseren Stichproben wäre dies eigentlich nicht so ein Problem. Wann ist nun eine Stichprobe gross genug? Da gibt es keine klare Grenze. Viele nennen 30 als eine Grenze. Wie gross ist unsere Stichprobe?
psych::describeBy(data$satisfaction_work_vas, data$computer_work_1.factor)##
## Descriptive statistics by group
## group: Mac OS
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 17 91.29 10.4 96 92.27 5.93 68 100 32 -1.09 -0.19 2.52
## ------------------------------------------------------------
## group: Windows
## vars n mean sd median trimmed mad min max range skew kurtosis se
## X1 1 35 80.74 21.87 86 84.79 14.83 21 100 79 -1.52 1.65 3.7Wir haben 17 und 35, also könnten wir uns auf den zentralen Grenzwertsatz und das Gesetz der grossen Zahlen verlassen.
Der zentrale Grenzwertsatz: Die Stichprobenverteilung ist annähernd normalverteilt, wenn der Stichprobenumfang groß genug ist (beispielsweise > 30).
Die Stichprobenverteilung ist: Wenn wir unendlich viele Stichproben nehmen und jeweils für jede Stichproben den Mittelwert berechnen. Die Stichprobenverteilung ist dann die Verteilung dieser Stichprobenmittelwerte.
Das Gesetz der grossen Zahlen sagt, dass der Mittelwert einer Stichprobe tendenziell näher beim Populationsmittelwert ist, je grösser die Stichprobe ist.
Wir könnten uns jetzt also auf den zentralen Grenzwertsatz beruhen und sagen, dass wir mehr als n = 30 haben und dass somit die Nicht-Normalverteilung kein grosses Problem ist.
Eine möglichkeit wäre auch noch, die Daten zu transformieren, damit sie eher normalverteilt sind.
Die nächste Graphik zeigt - mit der Hilfe des Befehls gladder aus dem describedata Paket, eine Serie von Transformationen.
- Achtung: Funktioniert nur mit positiven Zahlen.
gladder(data$satisfaction_work_vas)## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## ℹ The deprecated feature was likely used in the describedata package.
## Please report the issue to the authors.
Abbildung 22.1: Histogramme und Density-Kurve nach unterschiedlcihen Transformationen der Daten. Wir sehen hier, dass eigentlich keine der Transformationen eine optimale Lösung bietet.
Wir können auch statistische Tests zu allen Transformationen anzeigen lassen. Ein hoher P-Wert wäre hier besser, d.h. ein hoher P-Wert bedeutet weniger Evidenz gegen eine Normalverteilung in der Population.
ladder(data$satisfaction_work_vas)## # A tibble: 9 × 3
## Transformation statistic p.value
## <fct> <dbl> <dbl>
## 1 cubic 0.899 3.35e- 4
## 2 square 0.861 2.26e- 5
## 3 identity 0.762 8.58e- 8
## 4 sqrt 0.685 2.78e- 9
## 5 log 0.596 1.01e-10
## 6 1/sqrt 0.508 5.98e-12
## 7 inverse 0.432 6.74e-13
## 8 1/square 0.329 4.88e-14
## 9 1/cubic 0.279 1.52e-1422.2.1 Nichtparametrische Tests
Wir können auch nichtparametrische Tests durchführen, in unserem Falle haben wir ja zwei unabhängige Gruppen (Zufriedenheit bei Studierenden mit Mac und bei Studierenden mit PC), so können wir den Mann-Whitney-U-Test durchführen. Eine sehr gute deutsche Übersicht zu diesem Test finden Sie hier (hier klicken).
In R geht dies wie folgt (in R wird der Test Wilcoxon rank sum test genannt):
wilcox.test(satisfaction_work_vas~computer_work_1.factor, data = data, exact = FALSE, correct = FALSE, conf.int = TRUE)##
## Wilcoxon rank sum test
##
## data: satisfaction_work_vas by computer_work_1.factor
## W = 383.5, p-value = 0.09006
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
## -5.399237e-05 1.499996e+01
## sample estimates:
## difference in location
## 6Der Test sagt uns, dass es keine statistisch signifikanten Unterschied in der Zufriedenheit zwischen den beiden Computern gibt. Die passende Graphik zu diesem Test ist die Boxplot Graphik.
plot<-ggplot(data=data, aes(x=computer_work_1.factor, y=satisfaction_work_vas))+
geom_boxplot(aes(fill=computer_work_1.factor))
print(plot)
Abbildung 21.6: Boxplot Zufriedenheit mit dem Computer auf der Arbeit, separat für PC und Mac.
Es gibt auch ein Paket, dass uns hilft, den P-Wert direkt in die Graphik zu plotten (was wir natürlich auch manuell tun könnten).
library(ggpubr)
plot+ stat_compare_means()
Abbildung 22.2: Boxplot Zufriedenheit mit dem Computer auf der Arbeit, separat für PC und Mac.
Wie schon gesagt, könnten wir uns auf den zentralen Grenzsatz berufen und einen parametrischen Tests durchführen. Da wir zwei unabhängige Gruppen haben, wäre dies hier ein T-Test für unabhängige Gruppen (auch ungepaarter T-Test genannt).
t.test(data$satisfaction_work_vas~data$computer_work_1.factor, paired=FALSE)##
## Welch Two Sample t-test
##
## data: data$satisfaction_work_vas by data$computer_work_1.factor
## t = 2.3574, df = 49.999, p-value = 0.02236
## alternative hypothesis: true difference in means between group Mac OS and group Windows is not equal to 0
## 95 percent confidence interval:
## 1.561514 19.541007
## sample estimates:
## mean in group Mac OS mean in group Windows
## 91.29412 80.74286Wir können wieder das ggpubr Paket benutzen, um den p-Wert direkt in die Graphik zu plotten.
plot+stat_compare_means(method = "t.test")
Abbildung 22.3: Boxplot Zufriedenheit mit dem Computer auf der Arbeit, separat für PC und Mac.
In unserem Beispiel ist jedoch der Unterschied zwischen dem nichtparametrischen Test und dem parametrischen Test so gross, weil wir zwei ungleichgrosse Gruppen haben, die nicht die gleichegrosse Varianz haben. Hier wählt r als Standardtest beim T-Test den Welch Test, der nicht voraussetzt, dass die Varianzen in beiden Gruppen gleich gross sind. Dadurch hat dieser Tests eine höhere statistische Power und der P-Wert wird in unserem Beispiel kleiner.
Doch oft gibt eine andere Problematik eine Antwort auf diese Frage: Da wir keine randomiserte Studie durchgeführt haben, könnte es sein, dass unser Resultat confounded ist. Vielleicht ist der Anteil von Macs bei den Frauen grösser als bei den Männern und vielleicht beurteilen die Frauen die Zufriedenheit eher höher als die Männer. Wäre dies der Fall, hätten wir ein Confounding und dann müssten wir für die Variable Geschlecht kontrollieren. Dies können wir mit dem nichtparametrischen Test nicht (*das ist jetzt gelogen, aber sagen wir mal: nicht so einfach, siehe Harrell, Abschnitt 7.6, (hier klicken)). Wenn wir jedoch einen parametrischen Ansatz wählen, können wir anstelle des T-Testes eine lineare Regression durchführen und dort für die Variable geschlecht kontrollieren. Zuerst ohne kontrollieren für das Geschlecht. Hier haben wir nicht den gleichen P-Wert wie beim T-Test, da der T-Test nicht voraussetzt, dass die Varianzen in beiden Gruppen gleich sind, die lineare Regression jedoch schon. Siehe hier
lm.fit<-lm(satisfaction_work_vas~computer_work_1.factor, data=data)
jtools::summ(lm.fit, confint=TRUE)| Observations | 52 |
| Dependent variable | satisfaction_work_vas |
| Type | OLS linear regression |
| F(1,50) | 3.54 |
| R² | 0.07 |
| Adj. R² | 0.05 |
| Est. | 2.5% | 97.5% | t val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 91.29 | 82.05 | 100.54 | 19.84 | 0.00 |
| computer_work_1.factorWindows | -10.55 | -21.82 | 0.71 | -1.88 | 0.07 |
| Standard errors: OLS |
Und jetzt mit dem Geschlecht im Modell:
lm.fit<-lm(satisfaction_work_vas~computer_work_1.factor+gender_1.factor, data=data)
jtools::summ(lm.fit, confint=TRUE)| Observations | 52 |
| Dependent variable | satisfaction_work_vas |
| Type | OLS linear regression |
| F(3,48) | 1.51 |
| R² | 0.09 |
| Adj. R² | 0.03 |
| Est. | 2.5% | 97.5% | t val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 93.11 | 83.11 | 103.11 | 18.73 | 0.00 |
| computer_work_1.factorWindows | -10.79 | -22.29 | 0.71 | -1.89 | 0.07 |
| gender_1.factorMale | -6.17 | -18.30 | 5.96 | -1.02 | 0.31 |
| gender_1.factorRather not say | 0.18 | -28.05 | 28.41 | 0.01 | 0.99 |
| Standard errors: OLS |
In unserem Beispiel ändert sich das Resultat nicht - das Geschlecht ist hier kein Confounder.
Wir sehen jedoch auch, dass der P-Wert hier unter der Voraussetzung der gleichen Varianz nicht viel anders ist als der P-Wert vom nichtparametrischen Test. Der Unterschied zwischen dem P-Wert des nichtparametrischen Tests und des T-Tests mit der “ungleichen Varianz” Variante (d.h. Welch Test) wegen der Tatsache zustande kam, dass wir zwei ungleich grosse Gruppen haben, die stark unterschiedliche Varianzen aufweisen.
Take Home Message: Es ist nicht ganz einfach, einen klaren Ablauf vorzuschlagen. Der Autor würde in diesem Fall ohne schlechtes Gewissen den T-Test präsentieren, aber sicher gäbe es da auch Kritiker:innen, die sagen würden, dass ein nichtparametrischer Test hier angebrachter wäre.
22.3 Vergleich T-Test - Lineare Regression
Der nächste Abschnitt ist nur für sehr Interessierte, er kann aber hilfreich sein, um nicht übereinstimmende Resultate zu erklären.
- Wenn wir zwei Gruppen so behandeln, als hätten sie die gleiche Varianz in der abhängigen Variable, so wird der T-Test und die lineare Regression das gleiche Resultat ergeben (dazu müssen wir die option var.equal = TRUE angeben):
t.test(data$satisfaction_work_vas~data$computer_work_1.factor, paired=FALSE, var.equal=TRUE)##
## Two Sample t-test
##
## data: data$satisfaction_work_vas by data$computer_work_1.factor
## t = 1.8814, df = 50, p-value = 0.06575
## alternative hypothesis: true difference in means between group Mac OS and group Windows is not equal to 0
## 95 percent confidence interval:
## -0.7131705 21.8156915
## sample estimates:
## mean in group Mac OS mean in group Windows
## 91.29412 80.74286Der P-Wert ist hier 0.06575.
Hier die lineare Regression:
lm.fit<-lm(satisfaction_work_vas~computer_work_1.factor, data=data)
summary(lm.fit)##
## Call:
## lm(formula = satisfaction_work_vas ~ computer_work_1.factor,
## data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -59.743 -3.993 4.706 10.507 19.257
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 91.294 4.601 19.842 <2e-16 ***
## computer_work_1.factorWindows -10.551 5.608 -1.881 0.0657 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 18.97 on 50 degrees of freedom
## Multiple R-squared: 0.06611, Adjusted R-squared: 0.04743
## F-statistic: 3.54 on 1 and 50 DF, p-value: 0.06575Wir sehen, dass wir den gleichen P-Wert für den Computer erhalten (p = 0.06575).
Wenn wir jedoch die Standardeinstellung beim T-Test in R lassen, führt R automatisch einen T-Test durch, der nicht annimmt, dass die Varianzen der beiden Gruppen gleich sind (Welch-Test).
t.test(data$satisfaction_work_vas~data$computer_work_1.factor, paired=FALSE, var.equal=FALSE)##
## Welch Two Sample t-test
##
## data: data$satisfaction_work_vas by data$computer_work_1.factor
## t = 2.3574, df = 49.999, p-value = 0.02236
## alternative hypothesis: true difference in means between group Mac OS and group Windows is not equal to 0
## 95 percent confidence interval:
## 1.561514 19.541007
## sample estimates:
## mean in group Mac OS mean in group Windows
## 91.29412 80.74286Hier ist der P-Wert 0.02236. Diese Variante des T-Tests können wir nicht ganz so einfach mit einer linearen Regression reproduzieren. Wir müssen hier ein gemischtes Modell benutzen (das müssen Sie nicht erklären können und ist erst Inhalt in Master oder noch späteren Ausbildungen). Sie finden hier eine Erklärung oder eine andere Variante hier.
library(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:dplyr':
##
## collapsesummary(gls(satisfaction_work_vas ~ computer_work_1.factor, data = data, weights=varIdent(form = ~ 1 | computer_work_1.factor)))## Generalized least squares fit by REML
## Model: satisfaction_work_vas ~ computer_work_1.factor
## Data: data
## AIC BIC logLik
## 441.0183 448.6664 -216.5091
##
## Variance function:
## Structure: Different standard deviations per stratum
## Formula: ~1 | computer_work_1.factor
## Parameter estimates:
## Windows Mac OS
## 1.0000000 0.4756618
##
## Coefficients:
## Value Std.Error t-value p-value
## (Intercept) 91.29412 2.523077 36.18364 0.0000
## computer_work_1.factorWindows -10.55126 4.475717 -2.35745 0.0224
##
## Correlation:
## (Intr)
## computer_work_1.factorWindows -0.564
##
## Standardized residuals:
## Min Q1 Med Q3 Max
## -2.7316768 -0.1825690 0.2861010 0.8368697 0.8805118
##
## Residual standard error: 21.8704
## Degrees of freedom: 52 total; 50 residualOder diese Variante mit einem anderen Paket (hier müssen wir einen Umweg machen, siehe Erklärung hier):
library(lmerTest)## Loading required package: lme4## Loading required package: Matrix##
## Attaching package: 'Matrix'## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack##
## Attaching package: 'lme4'## The following object is masked from 'package:nlme':
##
## lmList##
## Attaching package: 'lmerTest'## The following object is masked from 'package:lme4':
##
## lmer## The following object is masked from 'package:stats':
##
## stepdata <- transform(data,
obs=factor(1:nrow(data)),
dummy=as.numeric(computer_work_1.factor=="Windows"))
summary(lmer(satisfaction_work_vas ~ computer_work_1.factor + (dummy-1|obs), data=data,
control=lmerControl(check.nobs.vs.nlev = "ignore",
check.nobs.vs.nRE = "ignore")))## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: satisfaction_work_vas ~ computer_work_1.factor + (dummy - 1 |
## obs)
## Data: data
## Control:
## lmerControl(check.nobs.vs.nlev = "ignore", check.nobs.vs.nRE = "ignore")
##
## REML criterion at convergence: 433
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.23919 -0.08684 0.13609 0.41883 0.83687
##
## Random effects:
## Groups Name Variance Std.Dev.
## obs dummy 370.1 19.24
## Residual 108.2 10.40
## Number of obs: 52, groups: obs, 52
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 91.294 2.523 16.000 36.184 <2e-16 ***
## computer_work_1.factorWindows -10.551 4.476 49.999 -2.357 0.0224 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## cmptr_w_1.W -0.56422.4 Wie sieht ein QQ-Plot bei einer Stichrpobe aus, die aus normalverteilten Daten gezogen wurde?
data<-rnorm(20, mean=100, sd=12)
hist(data)
Abbildung 22.4: Hisogramm einer Stichprobe aus einer normalverteilten Population
qqnorm(data)
qqline(data)
Abbildung 11.5: Quantil-Quantil-Plot (QQPlot) einer Stichprobe aus einer normalverteilten Population. Wir sehen, dass auch hier an den Enden eine Abweichung von der diagonalen Linie vorhanden ist.
library(ggpubr)
ggqqplot(data)## Warning: The following aesthetics were dropped during statistical transformation: sample
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: sample
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
Abbildung 1.11: Quantil-Quantil-Plot (QQPlot) mit Konfidenzintervall einer Stichprobe aus einer normalverteilten Population. Wir sehen, dass auch hier an den Enden eine Abweichung von der diagonalen Linie vorhanden ist, die Punkte verlassen jedoch das Konfidenzband nicht.
Wir können den Shapiro-Wilk Test benutzen, um die Nullhypothese zu testen, dass die Daten aus einer normalverteilten Population stammen.
shapiro.test(data)##
## Shapiro-Wilk normality test
##
## data: data
## W = 0.97117, p-value = 0.7794Der p-Wert ist nicht kleiner als 0.2 (hier nehmen wir oft 0.2 als Grenze und nicht 0.05), deshalb müssen wir die Hypothese, dass die Daten aus einer Normalverteilung stammen, nicht verwerfen.
Ziehen wir noch einmal eine Stichprobe:
set.seed(123)
data<-rnorm(80, mean=100, sd=12)library(ggpubr)
ggqqplot(data)## Warning: The following aesthetics were dropped during statistical transformation: sample
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: sample
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
Abbildung 4.8: Quantil-Quantil-Plot (QQPlot) mit Konfidenzintervall einer Stichprobe aus einer normalverteilten Population. Wir sehen, dass auch hier an den Enden eine Abweichung von der diagonalen Linie vorhanden ist, die Punkte verlassen jedoch das Konfidenzband nicht.
shapiro.test(data)##
## Shapiro-Wilk normality test
##
## data: data
## W = 0.99471, p-value = 0.9867Jetzt ziehen wir eine zufällige Stichprobe aus einer nichtnormalverteilten Population.
mean=4585
sd=7470
data<-rgamma(20,shape=(mean/sd)^2, rate=mean/sd^2 )hist(data)
Abbildung 11.6: Hisogramm einer Stichprobe aus einer nichtnormalverteilten Population
library(ggpubr)
ggqqplot(data)## Warning: The following aesthetics were dropped during statistical transformation: sample
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
## The following aesthetics were dropped during statistical transformation: sample
## ℹ This can happen when ggplot fails to infer the correct grouping structure in
## the data.
## ℹ Did you forget to specify a `group` aesthetic or to convert a numerical
## variable into a factor?
Abbildung 11.7: Quantil-Quantil-Plot (QQPlot) mit Konfidenzintervall einer Stichprobe aus einer nichtnormalverteilten Population. Wir sehen, dass hier an den Enden eine Abweichung von der diagonalen Linie vorhanden ist, die Punkte verlassen das Konfidenzband
shapiro.test(data)##
## Shapiro-Wilk normality test
##
## data: data
## W = 0.78142, p-value = 0.0004624Beim Shapiro-Wilk Test sehen wir, dass wir die Nullhypothese (die Daten stammen aus einer normalverteilten Population) verwerfen sollten. Wir glauben also nicht, dass die Daten normalverteilt sind.
22.5 Drei Gruppen vergleichen
Wenn wir zwei Gruppen haben, können wir eine lineare Regression oder einen T-Test (parametrisch) durchführen, oder einen nichtparametrischen Test (Wilcoxon). Wenn wir jedoch mehr als zwei Gruppen haben, können wir den T-Test und den Wilcoxon Test nicht mehr durchführen. Falls die Daten aus einer normalverteilung stammten, wurde früher meistens eine Varianzanalyse (Analysis of Variance, kurz ANOVA) durchgeführt. Wir führen jedoch lieber ein lineares Modell mit einer Faktor-Variable durch. Die Faktor-Variable gibt uns die Gruppenzugehörigkeit an.
Schauen wir uns ein Beipsiel an.
data<-rio::import("sit_to_stand_center.csv")Aus einer Studie von Anne-Gabrielle Mittaz Hager und ihrem Team haben wir Daten zum Five-Time-Sit-to-Stand Test aus drei Regionen: Oberwallis (VG), Mittell- und Unterwallis (VR für Valais Romande), sowie Waadt (VD).
summarytools::freq(data$center)## Frequencies
## data$center
## Type: Character
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## VD 81 20.30 20.30 20.30 20.30
## VG 47 11.78 32.08 11.78 32.08
## VR 271 67.92 100.00 67.92 100.00
## <NA> 0 0.00 100.00
## Total 399 100.00 100.00 100.00 100.00Hier eine Zusammenfassung der Daten Danke an den Author der Funktion length2
# New version of length which can handle NA's: if na.rm==T, don't count them
length2 <- function (x, na.rm=FALSE) {
if (na.rm) sum(!is.na(x))
else length(x)
}
table_ftst<-data %>%
group_by(center) %>%
summarize(n = length2(five_time_sts_t0m0, na.rm=TRUE),
mean=mean(five_time_sts_t0m0, na.rm=TRUE),
sd=sd(five_time_sts_t0m0, na.rm=TRUE))
kbl(table_ftst, caption='Zusammenfassung des 5xAufstehen Tests pro Gebiet') %>% kable_styling()| center | n | mean | sd |
|---|---|---|---|
| VD | 74 | 16.33108 | 6.823223 |
| VG | 44 | 15.25341 | 5.561068 |
| VR | 233 | 17.09410 | 6.363446 |
Da wir eine grosse Stichprobe haben, müssen wir nicht zwingend testen, ob die Daten aus einer Normalverteilung stammen und wir können einen parametrischen Test durchführen. Wir werden drei Varianten zeigen, um die Vergleichbarkeit mit anderen Statistikprogrammen wie SPSS zu gewährleisten.
- Variante 1: ANOVA (Standardeinstellung R)
- Variante 2: ANOVA (Einstellung, die zu den gleichen Resultate führen wie SPSS)
- Variante 3: Lineare Regression
Zuerst transformieren wir die Variable center noch zu einer Faktor Variable
data$center<-factor(data$center)22.5.1 Variante 1
one.way.anova <- aov(five_time_sts_t0m0 ~ center, data = data)
summary(one.way.anova)## Df Sum Sq Mean Sq F value Pr(>F)
## center 2 138 69.17 1.704 0.183
## Residuals 348 14123 40.58
## 48 observations deleted due to missingnessDas Resultat zeigt uns, dass es keinen statistisch signifikanten Unterschied zwischen den drei Gruppen gibt.
par(mfrow=c(2,2))
plot(one.way.anova)
par(mfrow=c(1,1))Wir können auch noch die Unterschiede zwischen den Gruppen anzeigen. In unserem Beispiel würden wir dies jedoch normalereise nicht tun, da es keinen signifikanten Unterschied zwischen den Gruppen gab.
tukey.plot.test<-TukeyHSD(one.way.anova)
plot(tukey.plot.test, las = 1)
22.5.2 Variante 2
Damit wir mit der Variante 2 die gleichen Resultate wie das Statistikprogramm SPSS erhalten, müssen wir eine Standardeinstellung ändern. Die Kontraste müssen wir zu Helmert Kontraste ändern.
options(contrasts = c("contr.helmert", "contr.poly"))Nun können wir die Anova mit einer linearen Regression durchführen. Damit wir Type-III Sum of Squares erhalten (wie in SPSS), müssen wir das Paket car benutzen.
Siehe dazu diese Seite.
lm_anova <- lm(five_time_sts_t0m0~center, data=data)
# anova(lm_anova) # Das würde die Standard Sum of Squares Einstellung ergeben, d.h. das gleiche Resultat wie Variante 1.
car::Anova(lm_anova,type=3)## Anova Table (Type III tests)
##
## Response: five_time_sts_t0m0
## Sum Sq Df F value Pr(>F)
## (Intercept) 58462 1 1440.5450 <2e-16 ***
## center 138 2 1.7043 0.1834
## Residuals 14123 348
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Viele Statistiker ziehen heutzutage die Analyse mit einem linearen Modell vor, da diese Analyse direkt die Mittelwert der einzelnen Gruppen anzeigt. Wir müssen jedoch zuerst wieder die Standardeinstellung für die Kontraste ändern (weil wir die weiter oben abgeändert haben).
- Standardmässig ist contr.treatment, bei dieser Einstellung ist der Intercept der Mittelwert der Referenzgruppe (bei uns VD) und mit den Koeffizienten der beiden anderen Levels können dann die Mittelwerte der beiden Gruppen berechnet werden.
- Beim Beispiel zuvor haben wir die Kontrastmethode auf helmert gesetzt (contr.helmert), bei dieser Variante ist der Intercept nicht der Mittelwert der Referenzgruppe, sondern der Durchschnitt der drei Levels. Mit den Koeffizienten kann danach wieder der Mittelwert zwischen zwei Gruppen berechnet werden siehe hier für Details. Deshalb wollen wir nicht diese Variante und müsssen wieder zur Standardmethode wechseln:
options(contrasts = c("contr.treatment", "contr.poly"))Der Mittelwert der Gruppe VD (ist die Referenzgruppe) sehen wir im Intercept (also 16.3311 Sekunden). Der Mittelwert für die Gruppe1 (VG) ist 16.3311 -1.0777 = 15.2534, und für die Gruppe2 (VR) = 16.3311 + 0.7630 = 17.0941.
summarytools::freq(data$center)## Frequencies
## data$center
## Type: Factor
##
## Freq % Valid % Valid Cum. % Total % Total Cum.
## ----------- ------ --------- -------------- --------- --------------
## VD 81 20.30 20.30 20.30 20.30
## VG 47 11.78 32.08 11.78 32.08
## VR 271 67.92 100.00 67.92 100.00
## <NA> 0 0.00 100.00
## Total 399 100.00 100.00 100.00 100.00lm_fit <- lm(five_time_sts_t0m0~center, data=data)
summary(lm_fit)##
## Call:
## lm(formula = five_time_sts_t0m0 ~ center, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.224 -4.075 -1.409 2.213 36.906
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.3311 0.7406 22.053 <2e-16 ***
## centerVG -1.0777 1.2127 -0.889 0.375
## centerVR 0.7630 0.8501 0.898 0.370
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.37 on 348 degrees of freedom
## (48 observations deleted due to missingness)
## Multiple R-squared: 0.0097, Adjusted R-squared: 0.004009
## F-statistic: 1.704 on 2 and 348 DF, p-value: 0.1834jtools::summ(lm_fit, confint=TRUE)| Observations | 351 (48 missing obs. deleted) |
| Dependent variable | five_time_sts_t0m0 |
| Type | OLS linear regression |
| F(2,348) | 1.70 |
| R² | 0.01 |
| Adj. R² | 0.00 |
| Est. | 2.5% | 97.5% | t val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 16.33 | 14.87 | 17.79 | 22.05 | 0.00 |
| centerVG | -1.08 | -3.46 | 1.31 | -0.89 | 0.37 |
| centerVR | 0.76 | -0.91 | 2.43 | 0.90 | 0.37 |
| Standard errors: OLS |
Manuell können wir die Mittelwert der drei Gruppen so berechnen
(VD = 16.3311)## [1] 16.3311(VG = 16.3311-1.0777)## [1] 15.2534(VR = 16.3311+ 0.7630)## [1] 17.0941Wir können uns die Mittelwerte der Gruppen auch so ausgeben lassen:
lm_fit <- lm(five_time_sts_t0m0~0+center, data=data)
levels(data$center)## [1] "VD" "VG" "VR"summary(lm_fit)##
## Call:
## lm(formula = five_time_sts_t0m0 ~ 0 + center, data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -9.224 -4.075 -1.409 2.213 36.906
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## centerVD 16.3311 0.7406 22.05 <2e-16 ***
## centerVG 15.2534 0.9604 15.88 <2e-16 ***
## centerVR 17.0941 0.4173 40.96 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 6.37 on 348 degrees of freedom
## (48 observations deleted due to missingness)
## Multiple R-squared: 0.8741, Adjusted R-squared: 0.873
## F-statistic: 805.4 on 3 and 348 DF, p-value: < 2.2e-16Oder so:
emmeans::emmeans(lm_fit, "center")## center emmean SE df lower.CL upper.CL
## VD 16.3 0.741 348 14.9 17.8
## VG 15.3 0.960 348 13.4 17.1
## VR 17.1 0.417 348 16.3 17.9
##
## Confidence level used: 0.95- Take Home Message: Alle drei Varianten sind korrekt - sie beantworten einfach nicht exakt die gleiche Fragestellung. Wir wählen in R eigentlich immer die dritte Variante.
22.6 Beispiel Risikofgruppen
Wir schauen uns jetzt noch die Variante 3 mit einem Beispiel an, wo es einen statistisch signifikanten Unterschied zwischen den Gruppen gibt. Wir teilen die Personen in drei Gruppen ein mit der vorletzten Version des STEADI-Tools, mit dem die älteren Personen in niedriges Sturzrisiko, mittleres Sturzrisiko und höheres Sturzrisiko eingeteilt werden. Die abhängige Variable ist der Five-Time Sit-to-Stand test. Personen die den Test nicht durchführen konnten, bekamen in unserer Analyse den Wert 60 Sekunden (damit sie trotzdem analysiert werden konnten).
library(dplyr)
data<-rio::import("sit_to_stand_risk_category.csv")
data<-data %>%
mutate(FTSTS_notpossible60_t0m0=ifelse(FTSTS_notpossible121_t0m0>120, 60, FTSTS_notpossible121_t0m0),
STEADI_CDC_0LR_1MR_2HR_t0m0=factor(STEADI_CDC_0LR_1MR_2HR_t0m0, levels=c(0,1,2), labels=c("Low Risk", "Moderate Risk" ,"High Risk")))
boxplot(data$FTSTS_notpossible60_t0m0~data$STEADI_CDC_0LR_1MR_2HR_t0m0)
lm.fit<-lm(FTSTS_notpossible60_t0m0~factor(STEADI_CDC_0LR_1MR_2HR_t0m0), data=data)
jtools::summ(lm.fit, confint=TRUE)| Observations | 377 (28 missing obs. deleted) |
| Dependent variable | FTSTS_notpossible60_t0m0 |
| Type | OLS linear regression |
| F(2,374) | 21.97 |
| R² | 0.11 |
| Adj. R² | 0.10 |
| Est. | 2.5% | 97.5% | t val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 11.01 | 7.44 | 14.57 | 6.07 | 0.00 |
| factor(STEADI_CDC_0LR_1MR_2HR_t0m0)Moderate Risk | 10.78 | 6.42 | 15.14 | 4.86 | 0.00 |
| factor(STEADI_CDC_0LR_1MR_2HR_t0m0)High Risk | 13.78 | 9.69 | 17.87 | 6.63 | 0.00 |
| Standard errors: OLS |
means<-emmeans::emmeans(lm.fit, "STEADI_CDC_0LR_1MR_2HR_t0m0")
means## STEADI_CDC_0LR_1MR_2HR_t0m0 emmean SE df lower.CL upper.CL
## Low Risk 11.0 1.81 374 7.44 14.6
## Moderate Risk 21.8 1.28 374 19.27 24.3
## High Risk 24.8 1.02 374 22.79 26.8
##
## Confidence level used: 0.95plot(means)
Abbildung 11.16: Mittelwerte des Five-Time Sit-to-Stand Tests in den drei Sturzrisikogruppen.