Modul 6 - Inferenzstatistik

Sebastian Schnettler und Svenja Heinecke 2024

 

 

Einleitung

Da dieses Modul zur Inferenzstatistik hauptsächlich darauf abzielt, Ihnen den Hintergrund für das Verständnis der inferenzstatistischen Implikationen der Regressionsanalyse zu vermitteln, ist diese R-Sitzung sehr kurz gehalten und enthält neben dieser Übungsdatei keine weiteren Aufgaben mit Lösungen. Weiterhin verfügen Sie mittlerweile über die grundlegenden R-Kenntnisse, die es Ihnen ermöglichen sollten, dem hier vorgestellten Material problemlos zu folgen. So zeigen wir hier nur überblicksartig, wie Sie einige inferenzstatistische Aufgaben erledigen können.

Beispieldaten laden

Auch hier laden wir den Datensatz GSSvocab:

library(car) # für den GSSvocab-Datensatz
data(GSSvocab)
head(GSSvocab)

 

##        year gender nativeBorn ageGroup educGroup vocab age educ
## 1978.1 1978 female        yes    50-59    12 yrs    10  52   12
## 1978.2 1978 female        yes      60+   <12 yrs     6  74    9
## 1978.3 1978   male        yes    30-39   <12 yrs     4  35   10
## 1978.4 1978 female        yes    50-59    12 yrs     9  50   12
## 1978.5 1978 female        yes    40-49    12 yrs     6  41   12
## 1978.6 1978   male        yes    18-29    12 yrs     6  19   12

Konfidenzintervalle

 Punktschätzer und Standardfehler

Um den Punktschätzer zu erhalten, müssen Sie die jeweilige Statistik von Interesse kennen. Dies kann die aus einer Tabelle gewonnene Häufigkeit oder ein Maß für die zentrale Tendenz sein, das z. B. mit der Funktion median() oder mean() ermittelt wird.

Zur Berechnung des Standardfehlers muss man wissen, wie man eine Standardabweichung (sd()) berechnet und wie man den Stichprobenumfang und die Quadratwurzel (sqrt()) des Stichprobenumfangs erhält.

Die Funktion length() gibt Ihnen die Anzahl der Fälle in einem Vektor an. Die Funktion dim() liefert Ihnen Informationen über die Anzahl der Fälle und Variablen in einem Datensatz. Wenn Sie nur die Anzahl der Fälle (der erste Wert, der von dim() zurückgegeben wird) erhalten möchten, würden Sie schreiben: dim(meineDaten)[[1]].

Ein Beispiel folgt nach dem nächsten Unterabschnitt.

Kritischer Wert

Traditionell war es üblich, kritische Werte für gewünschte statistische Signifikanzniveaus aus Verteilungstabellen in Statistikbüchern nachzuschlagen. R bietet dafür eine Reihe von Funktionen, mit denen Sie Werte aus bekannten statistischen Verteilungen ziehen können. So gibt es z.B. eine Funktion qnorm der Standardnormalverteilung, mit der Sie einen Quantilwert angeben können und die den entsprechenden kritischen Wert für dieses Quantil zurückgibt.

Beispielsweise entsprechen in der Standardnormalverteilung das 2,5 %-Quantil und das 97,5 %-Quantil dem kritischen Wert von |1,96| bei einem statistischen Signifikanzniveau von 0,05:

qnorm(p = c(0.025, 0.975), mean = 0, sd = 1)

## [1] -1.959964  1.959964

Ebenso gibt es eine ganze Reihe weiterer Funktionen für andere Verteilungen, z.B. qt() für die t-Verteilung, qf() für die f-Verteilung und qpois() für die Poisson-Verteilung. Um eine Liste der in R verfügbaren Verteilungen zu erhalten, geben Sie ?distributions ein.

Der folgende Befehl liefert zum Beispiel den kritischen t-Wert bei den 10%- und 90%-Quantilen der t-Verteilung mit 10 Freiheitsgraden:

qt(p = c(0.1,0.9), df = 10)

## [1] -1.372184  1.372184

Für jede Verteilung stehen weitere nützliche Funktionen zur Verfügung: Der führende Buchstabe “r” steht für eine Funktion, die zufällig n-Werte aus der jeweiligen Verteilung zieht. Der folgende Code zieht zum Beispiel zufällig 1000 Fälle aus einer Standardnormalverteilung und stellt sie als Histogramm dar:

x <- rnorm(n = 1000, mean = 10, sd = 2)
hist(x) # Basis-Funktion in R für ein Histogramm

 

Der führende Buchstabe “p” leitet Funktionen ein, die es ermöglichen, einen beliebigen Wert der jeweiligen Verteilung einzugeben und zu erfahren, welches Quantil diesem Wert entspricht. Zum Beispiel:

pnorm(1.96)

## [1] 0.9750021

Zusammenfassung mit Beispiel

Verwenden wir also all diese Informationen, um das 95%-Konfidenzintervall für das arithmetische Mittel des Abschneidens im Vokabeltest im Datensatz GSSvocab zu erhalten:

x <- GSSvocab$vocab
mean.vocab <- mean(x, na.rm=T)
crit.low <- qnorm(p=0.025)
crit.hi <- qnorm(p=0.975)
st.error <- sd(x, na.rm=T) / sqrt(length(x))
CI.low <- mean.vocab + crit.low * st.error
CI.hi <- mean.vocab + crit.hi * st.error
CI <- c(CI.low, mean.vocab, CI.hi)
names(CI) <- c("95%CI.lo", "Mean", "95%CI.hi")
CI

## 95%CI.lo     Mean 95%CI.hi 
## 5.973929 5.998219 6.022510

Wenn Ihnen das zu umständlich erscheint, gibt es Pakete, die Konfidenzintervalle für Sie berechnen. Zum Beispiel gibt es die Funktion MeanCI()im Zusatzpaket DescTools, die ein Konfidenzintervall um den jeweiligen arithmetischen Mittelwert berechnet:

library(DescTools)
MeanCI(x, na.rm=T)

##     mean   lwr.ci   upr.ci 
## 5.998219 5.973340 6.023099

Die Werte weichen geringfügig von unseren eigenen Werten ab. Dies lässt sich auf Rundungsfehler zurückführen.

Hypothesentests

t.test in R

R bietet eine Reihe von Funktionen zur Durchführung von Null-Hypothesentests. Der t-Test kann mit der Funktion t.test() aus dem Zusatzpaket stats durchgeführt werden, welches bereits mit der Standardinstallation von R geliefert wird. Durch Anpassung der Funktionsargumente können Sie zwischen einem t-Test mit einer Stichprobe, zwei Stichproben und gepaarten Stichproben wählen.

Wenn Sie zum Beispiel im GSSvocab-Datensatz testen möchten, ob sich die mittleren Testergebnisse im Vokabeltest zwischen Männern und Frauen unterscheiden, können Sie einen t-Test mit zwei Stichproben durchführen:

t.test(vocab ~ gender, alternative = "two.sided",
       data = GSSvocab)

## 
##  Welch Two Sample t-test
## 
## data:  vocab by gender
## t = 3.0642, df = 25344, p-value = 0.002185
## alternative hypothesis: true difference in means between group female and group male is not equal to 0
## 95 percent confidence interval:
##  0.02836647 0.12907480
## sample estimates:
## mean in group female   mean in group male 
##             6.032289             5.953568

Im Output dieses t-Tests wird Ihnen dazu direkt die Information für die Alternativ-Hypothese mitgeliefert: “alternative hypothesis: true difference in means between group female and group male is not equal to 0”.

Das Ergebnis zeigt dann, dass Sie die Nullhypothese, dass es keinen Unterschied zwischen den beiden Gruppen gibt, verwerfen können, da der p-Wert unter 0,05 liegt. Damit wird die Alternativ-Hypothese angenommen und wir gehen von Gruppenunterschieden zwischen Frauen und Männern aus. Allerdings ist dieser Mittelwertunterschied sehr gering.

Weitere Tests

Andere Tests, die in R zur Verfügung stehen, sind bspw.:

• chisq.test() zum Testen der Nullhypothese, dass es keinen statistischen Zusammenhang zwischen zwei kategorialen Variablen gibt. Zum Beispiel:

chisq.test(x = GSSvocab$ageGroup,
           y = GSSvocab$educGroup)

## 
##  Pearson's Chi-squared test
## 
## data:  GSSvocab$ageGroup and GSSvocab$educGroup
## X-squared = 1474.6, df = 16, p-value < 2.2e-16

• binom.test() für den Binomialtest, mit dem die Nullhypothese, dass es keinen Unterschied in den Prozentsätzen gibt, getestet werden kann.

Die Tests liefern etwas unterschiedliche Ausgaben, aber die Logik ist immer die gleiche. Es wird eine Teststatistik berechnet und ein p-Wert angegeben.

Da es Tausende von benutzergenerierten Paketen in R gibt, stehen die Chancen gut, dass Sie einen in R implementierten Test für fast jeden Zweck finden.

Bootstrapping

Der folgende Text ist optional und kann übersprungen oder zu einem späteren Zeitpunkt behandelt werden.

Wenn Sie Bootstrapping für eine beliebige Statistik anwenden wollen, können Sie das Paket boot verwenden. Die allgemeine Logik, auch wenn die R-Programmierung Kenntnisse einiger fortgeschrittener Funktionen erfordert, die in diesem Kurs nicht behandelt werden, ist einfach und kann für die Berechnung eines Bootstrap-Standardfehlers für das mittlere Abschneiden im Vokabeltest illustriert werden.

Mit der Funktion sample() können Sie zufällig Werte aus einem Vektor ziehen, entweder mit oder ohne Ersetzung. Im Fall von Bootstrapping wollen wir eine Stichprobe mit Ersetzung durchführen und eine Anzahl von Bootstrap-Stichproben aus unserer empirischen Stichprobe ziehen, die den gleichen Umfang haben wie die empirische Stichprobe.

Der folgende Code würde eine solche Bootstrap-Stichprobe aus unserer Variable x, die wir oben erstellt haben und die die Werte der Variablevocab aus dem “GSSvocab”-Datensatz enthält, ziehen.

boot1 <- sample(x = x, size = length(x), replace = TRUE)

Nun müssten wir dies immer und immer wieder tun, sagen wir tausendmal, und jede daraus resultierende Bootstrap-Stichprobe speichern. Glücklicherweise gibt es eine Funktion, mit der wir dies auf sehr einfache Weise und ohne 1000 Zeilen Code schreiben zu müssen tun können. Die Funktion lapply() erlaubt es, eine Funktion, die durch das Funktionsargument FUN = ... spezifiziert wird, auf jedes Element eines Vektors oder einer Liste anzuwenden, und gibt eine Liste mit so vielen Einträgen wie die ursprüngliche Liste zurück. Zu diesem Zweck schreiben wir die ursprüngliche Beispielfunktion leicht um.

Unsere neue Funktion ist einfach eine Kopie der Funktion sample() mit einigen voreingestellten Werten:

myboot <- function(x) sample(
  x = GSSvocab$vocab,
  size = length(GSSvocab$vocab),
  replace = TRUE)

Und nun können wir diese Funktion auf jede Zahl im Vektor 1:1000 anwenden, was eine andere Art zu sagen ist: “Mache dies tausendmal”:

boot1000 <- lapply(1:1000, myboot)

Die Ausgabe ist eine Liste mit 1000 Bootstrap-Stichproben, die jeweils die gleiche Länge wie der Eingabevektor (unsere empirische Stichprobe) haben.

Jetzt können wir einen Mittelwert für jede dieser Bootstrap-Stichproben berechnen, indem wir wieder die Funktion lapply() oder besser die Funktion sapply() verwenden, die eine Kopie der Funktion lapply() ist, mit dem Unterschied, dass sie versucht, die Ausgabe zu vereinfachen, sodass wir einen Vektor von Mittelwerten und nicht eine Liste von Mittelwerten erhalten:

boot1000means <- sapply(boot1000, FUN = mean, na.rm=T)
hist(boot1000means) # Histogramm der Bootstrap-Mittelwerte

 

Schließlich können wir die Standardabweichung der Bootstrap-Verteilung der Mittelwerte berechnen:

sd(boot1000means, na.rm=T)

## [1] 0.0127156

---

Dieses Modul ist Teil des Kurses: Dein SoWi-Statistik-Kurs: Von der Theorie zur Praxis in R - Materialien zum Selbstlernen und Unterrichten. Carl von Ossietzky Universität Oldenburg. Für eine eventuelle Vervielfältigung, Veränderung oder Weitergabe der Kursinhalte beachten Sie bitte die Lizenzbestimmungen nach CC BY 4.0 (siehe dazu auch die Modulbeschreibung dieses Kurses).