invisible header

link zum pdf

# Statistik 2: R Tutorat
# Übungsskript zum Thema Inferenzstatistik
# Datum: 12.04.2025
# AutorIn: XXX

 

Vorbereitung: Einlesen der Daten und Aktivieren der nötigen Packages

# Working directory setzen (z.B. "c:\daten")
setwd("mein_laufwerk/mein_datenverzeichnis")
# Daten einlesen
library(haven)
datensatz <- read_dta("daten.dta")
# install.packages("dplyr")
library(tidyverse)
library(ggplot2) 
library(labelled)
library(ggeffects)

 

1. Hypothesen formulieren

Uns beschäftigt der Zusammenhang zwischen Bildung und dem zeitlichen Umfang der Internetnutzung.

Unterschiedliche Begründungen und Vermutungen sind hier zwar denkbar, wir finden es aber plausibel, dass Personen mit mehr Bildung tendenziell mehr professionelle Anwendungsbezüge haben und das Internet somit umfangreicher nutzen. Aus diesen Überlegungen können wir eine empirisch prüfbare Hypothese (auch Forschungs- oder Alternativhypothese genannt) ableiten:

Hypothese 1: Bildung hat einen positiven Einfluss auf den zeitlichen Umfang der Internetnutzung.

Der H1 (auch Forschungshypothese oder Alternativhypothese genannt) steht immer eine Nullhypothese gegenüber, selbst wenn diese oft nicht explizit formuliert wird:

Nullhypothese: Bildung hat keinen positiven Einfluss auf den zeitlichen Umfang der Internetnutzung.

Damit wir die tatsächliche forschungspraktische Relevanz eines Koeffizienten bewerten können, reicht es nicht aus, diesen inhaltlich über die Interpretation und Visualisierung des Koeffizienten anzuschauen. Wir müssen zusätzlich seine statistische Bedeutsamkeit ermitteln, indem wir den Grad an Überzufälligkeit des Zusammenhangs in der Stichprobe anhand des p-Wertes ermitteln. Diese Identifikation von Überzufälligkeit ist das zentrale Element der Inferenzstatistik. Hypothesen bilden dabei die Grundlage. Unterschreitet der p-Wert einen hinreichend tiefen Schwellenwert, so drücken wir die statistische Bedeutsamkeit durch Ablehnung der Nullhypothese aus.

 

2. Analyse

2.1 Datenmanagement – Inspektion und Selektion

Wir messen die Konstrukte Bildung und Internetnutzung in unserer Hypothese mit den Variablen eduyrs und netustm. Zudem beschränken wir unsere Stichprobe auf den Teildatensatz der Schweiz. Wir verschaffen uns einen Überblick über die beiden Variablen und reduzieren unseren Datensatz entsprechend.

ess8_CH <- filter(ess8, cntry == "CH") 
look_for(ess8_CH, "eduyrs")
##  pos variable label                          col_type missing
##  336 eduyrs   Years of full-time education ~ dbl+lbl  3      
##                                                              
##                                                              
##  values            
##  [NA(b)] Refusal   
##  [NA(c)] Don't know
##  [NA(d)] No answer
look_for(ess8_CH, "netustm")
##  pos variable label                     col_type missing values                
##  9   netustm  Internet use, how much t~ dbl+lbl  338     [NA(a)] Not applicable
##                                                          [NA(b)] Refusal       
##                                                          [NA(c)] Don't know    
##                                                          [NA(d)] No answer
ess8_CH <- select(ess8_CH, internet = netustm, eduyrs, idno)
summary(ess8_CH)
##     internet          eduyrs          idno     
##  Min.   :   3.0   Min.   : 0.0   Min.   :   1  
##  1st Qu.:  60.0   1st Qu.: 9.0   1st Qu.:1017  
##  Median : 120.0   Median :10.0   Median :1827  
##  Mean   : 163.6   Mean   :11.3   Mean   :1807  
##  3rd Qu.: 240.0   3rd Qu.:13.0   3rd Qu.:2645  
##  Max.   :1200.0   Max.   :26.0   Max.   :3460  
##  NA's   :338      NA's   :3
sd(ess8_CH$eduyrs, na.rm = TRUE)
## [1] 3.496909
sd(ess8_CH$internet, na.rm = TRUE)
## [1] 163.1974

Wir haben nun einen reduzierten Datensatz der ausschliesslich die Variablen netustm und eduyrs, sowie den Identifyer idno enthält. Im Rahmen des select()-Befehls haben wir gleichzeitig auch die Variable zur Internetnutzung umbenannt. Der Output des look_for-Befehls legt zudem die Messung der verwendeten Variablen offen: eduyrs misst die Anzahl an Bildungsjahren, internet misst die tägliche Internetnutzung in Minuten. Beide Variablen sind metrisch skaliert.

Alle Variablen scheinen korrekt eingelesen. Es gibt zudem keine Kategorien für fehlende Werte, die wir noch zu NAs rekodieren müssten. internet enthält Fälle mit einem Maximum von 1200 Minuten, was unplausibel hoch erscheint. Hier wäre daher eigentlich eine Ausreisseranalyse und ggf. Datenbereinigung sinnvoll, wir beschäftigen uns jedoch hier nicht weiter damit.

Für die Analyse beschränken wir uns auf die Schweiz, oben umgesetzt über den filter-Befehl. Der Datensatz enthält damit noch 1525 Merkmalsträger. D

2.2 Regressionsanalyse

Nun zur Regressionsanalyse mithilfe des lm()-Befehls.

fit <- lm(internet ~ eduyrs, data = ess8_CH)
summary(fit)
## 
## Call:
## lm(formula = internet ~ eduyrs, data = ess8_CH)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -204.19 -105.20  -52.06   44.16 1011.66 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   72.646     16.171   4.492 7.74e-06 ***
## eduyrs         7.713      1.313   5.875 5.48e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 161.1 on 1182 degrees of freedom
##   (341 Beobachtungen als fehlend gelöscht)
## Multiple R-squared:  0.02838,    Adjusted R-squared:  0.02755 
## F-statistic: 34.52 on 1 and 1182 DF,  p-value: 5.48e-09

Dem Summary-Output können wir u.a. folgende Werte entnehmen:

  • Regressionskoeffizient: Der Koeffizient (7.713) zeigt uns einen positiven Zusammenhang zwischen den Bildungsjahren und der Nutzungszeit an: Mit jedem zusätzlichen Bildungsjahr steigt also die Nutzungszeit im Schnitt um 7 Minuten und 43 Sekunden. Ordnen wir den Koeffizienten zusätzlich vor dem Hintergrund der Standardabweichung ein: Mit jedem zusätzlichen Bildungsjahr steigt die Internetnutzung um 0.047 Standardabweichungen. Ein MA nutzt also das Internet im Schnitt etwa 30 Minuten länger als ein BA - ein Unterschied, der immerhin etwa ein Fünftel einer typischen Mittelwertabweichung ausmacht. Es handelt sich also eindeutig um einen inhaltlich bedeutsamen Effekt. (Gleichwohl wissen wir aus der letzten Einheit, dass wir einen Disclaimer in die Auswertung einbringen müssen, da der Zusammenhang relativ stark durch Ausreisserwertepaare getrieben ist).

  • Standardfehler: Hinter dem Koeffizienten können wir den Standardfehler ablesen. Er beträgt 1.313. Achtung: Standardfehler und Standardabweichung dürfen nicht verwechselt werden. Beide messen Variation, beziehen sich allerdings auf unterschiedliche Ebenen: Die Standardabweichung misst die Variation innerhalb einer Stichprobe, der Standardfehler misst die Variation zwischen verschiedenen Stichproben bzw. ihren Kennwerten und gibt somit den erwartbaren Stichprobenfehler an.

  • t-Wert: Der t-Wert hinter dem Standardfehler beträgt 5.875.

  • p-Wert: Ganz rechts befindet sich der p-Wert, er beträgt 5.48e-09 (0.00000000548). Das Signifikanzniveau ist in Sternchen hinter dem p-Wert ablesbar.

2.3 Standardfehler

Der Standardfehler (SE) gibt die durchschnittliche Abweichung eines Parameters der Stichprobe vom wahren Parameterwert in der Grundgesamtheit an. In unserem Fall: Wir können erwarten, dass der wahre Anstieg der Nutzungsdauer pro Bildungsjahr um 1.3 Minuten vom Regressionskoeffizienten in der Stichprobe abweicht. Wir müssen also damit rechnen, dass die wahre Steigung der Nutzungsdauer pro Bildungsjahr um etwa 1.3 Minuten grösser oder kleiner ausfällt als 7.7.

2.4 t-Wert

Das Verhältnis von Koeffizient zum SE wird durch den t-Wert angegeben.

\(t=\frac{b}{SE}=\frac{7.713}{1.313}\)

Er misst die Grösse des Koeffizienten in Einheiten des Standardfehlers bzw. den «Sicherheitsabstand» des Koeffizienten von der 0. Der t-Wert wird in euren Auswertungen zumeist nicht berichtet.

2.5 p-Wert

Der p-Wert zeigt uns, wie wahrscheinlich der vorgefundene Stichprobenkoeffizient ist, wenn es in Wirklichkeit keinen Zusammenhang zwischen UV und AV gäbe, also demnach die (beidseitige Variante) der Nullhypothese richtig wäre.

Als Konvention gelten Werte \(p<0.05\) oder \(p<0.01\) als Schwellen für die Feststellung statistischen Signifikanz. (Achtung: Diese Konventionen gelten u.U. nicht disziplin- und kontextübergreifend). Im Regressionsoutput wird zudem noch ein Wert von \(p<0.001\) aufgelistet. Dem Signifikanzniveau entsprechend werden Sterne verteilt. Achtet darauf, dass ihr keine normative Bewertung des p-Wertes vornehmt, hohe p-Werte sind nicht “schlecht” und tiefe p-Werte nicht “gut”. Sie sind zuallererst statistische Kennzahlen bzw. neutrale statistische Befunde.

2.6 Hypothesenevaluation

Wir können die Nullhypothese, dass es keinen positiven Zusammenhang zwischen den beiden Variablen gibt, ablehnen. Die Forschungshypothese, dass Bildung den zeitlichen Umfang der Internetnutzung positiv beeinflusst, wird inferenzstatistisch gestützt. Hinweis: Auch wenn die von R postulierte Forschungshypothese zweiseitig ausgerichtet ist, können wir den abgeleiteten p-Wert einer Konvention folgend für die Prüfung unserer einseitigen Hypothese verwenden. Letztlich wird hierdurch die Schwelle für die Verwerfung der H0 höher gelegt (siehe Folien letzte Vorlesung).

Der Summary-Output zeigt uns jedoch nur einen Teil der wichtigen inferenzstatistischen Parameter an, einige weitere müssen wir zusätzlich berechnen.

 

3. Weitere inferenzstatistische Parameter

Weitere wichtige inferenzstatistische Parameter sind:

3.1 Konfidenzintervall des Koeffizienten

Das Konfidenzintervall des Koeffizienten zeigt an, zwischen welchen Werten der wahre Koeffizient in der Grundgesamtheit mit einer bestimmten Sicherheit (z.B. 95%) liegt. Wir können die Grenzen mit dem Befehl confint ermitteln, das Argument level können wir dabei beliebig variieren.

confint(fit, level= 0.95)
##                 2.5 %    97.5 %
## (Intercept) 40.918306 104.37282
## eduyrs       5.137214  10.28828

Wie können wir nun diese Werte interpretieren? Wir erkennen, dass der wahre Koeffizient der Grundgesamtheit mit 95%-iger Sicherheit zwischen 5.14 und 10.29 liegt. Diese Grenzen für das 95%-Konfidenzintervall können wir genau wie in der univariaten Statistik näherungsweise mit einer Daumenregel berechnen:

\(KI=b\pm 2*SE\)

In diesem Fall: \(KI = 7.7\pm 2.6\)

Hinweis: Bei kleineren Stichproben kann diese Berechnung ungenau sein.

3.2 Konfidenzband

Das Konfidenzband zeigt den Bereich an, in dem die wahre Regressionsgerade mit 95%-Sicherheit verläuft. Das Konfidenzband wird von ggplot() automatisch generiert.

plot1 <- ggplot(ess8_CH, aes(x = eduyrs, y = internet))+
  geom_jitter(alpha = 0.2, height = 10) +
  scale_x_continuous(breaks = seq(0,26,2))+
  coord_cartesian(ylim = c(0,450)) +
  geom_smooth(method = "lm", se = TRUE, level = 0.95)+
  theme_bw()+
  labs(title = "Streudiagramm: Bildung und Internetnutzung", 
  y = "Internetnutzung in Minuten/Tag",
  x = "Bildungsjahre", 
  caption = "ESS(2016), Teilstichprobe CH, N = 1184.\n Regressionsgerade mit 95-Prozent-Konfidenzband.")
plot1

Das Argument se = TRUE erzeugt das Konfidenzband, mit dem Argument fill können wir zudem noch eine Farbe wählen, mit der das Band gefüllt wird. Per Default wird ein 95-Prozent-Konfidenzband dargestellt, der Bereich kann aber über das level-Argument variiert werden (das ist allerdings unüblich). Hierbei gilt: Je grösser der Konfidenzwert, desto höher der angelegte Sicherheitsstandard, desto breiter das Konfidenzband.

Hier wurde zudem der Anzeigebereich des Plots mittels coord_cartesian(ylim = c(0,450)) eingeschränkt. Visualisiert wird damit nur der Ausschnitt auf der vertikalen Achse zwischen 0 und 450. Mit diesem Argument werden die Werte ausserhalb des ausgeschnittenen Bereichs ausgeblendet, aber nicht ausgeschlossen (was bei ähnlichen Befehlen wie xlim() der Fall wäre).

Die Grenzen bzw. der Grenzverlauf des Konfidenzbandes kann leicht mit ggpredict nachvollzogen werden.

library (ggeffects)
ggpredict(fit, terms = "eduyrs[2, 8, 14, 20, 26]", interval = "confidence", ci_level = 0.95)
## # Predicted values of Internet use, how much time on typical day, in minutes
## 
## eduyrs | Predicted |         95% CI
## -----------------------------------
##      2 |     88.07 |  61.23, 114.91
##      8 |    134.35 | 120.94, 147.75
##     14 |    180.62 | 169.82, 191.43
##     20 |    226.90 | 203.85, 249.95
##     26 |    273.18 | 235.45, 310.91

Auf die Optionen kann hier auch verzichtet werden, da sie default sind.

library (ggeffects)
ggpredict(fit, terms = "eduyrs[2, 8, 14, 20, 26]")
## # Predicted values of Internet use, how much time on typical day, in minutes
## 
## eduyrs | Predicted |         95% CI
## -----------------------------------
##      2 |     88.07 |  61.23, 114.91
##      8 |    134.35 | 120.94, 147.75
##     14 |    180.62 | 169.82, 191.43
##     20 |    226.90 | 203.85, 249.95
##     26 |    273.18 | 235.45, 310.91

Die zweite Spalte (“Predicted”) gibt hier den Verlauf der Regressionsgerade wieder bzw. enthält vorhergesagte Werte: Für eine Person mit 2 Bildungsjahren werden laut Regressionsgerade 88 tägliche Internetminuten vorhergesagt.
Die dritte Spalte informiert zum Grenzverlauf bzw. zu den Grenzwerte der 95%-Bandes: Bei 2 Bildungsjahren verläuft das 95% Konfidenzband zwischen 61 und 115 Minuten.

3.3 Vorhersageband und Vorhersageintervall

Das 50%-Vorhersageband markiert den Bereich der Verteilung, der 50% aller Werte enthält (in dem also ein weiterer Fall mit 50% Sicherheit liegt). Es berücksichtigt zusätzlich zur Unsicherheit ob der korrekten Regressionsgerade auch und insbesondere die Streuung der tatsächlichen Werte um die Regressionsgerade (daher besteht auch eine Verwandtschaft zum R2).

Bei gleicher Sicherheitsstufe ist es immer deutlich breiter als das Konfidenzband und geht sogar oft in den irrealen Wertebereich. Es ist daher häufig sinnvoll, den Sicherheitswert beim Vorhersageintervall niedriger (z.B. 50%, in Anlehnung an das geläufige univariate Konzept des Interquartilabstands). Fünfzig Prozent sind ein leicht verständlicher Grenzmarker (die Hälfte liegt drin, die Hälfte liegt draussen) und korrespondieren zudem gut zur Idee des Kerns einer (bedingten) Verteilung. Fünfzigprozentige Sicherheit wird zudem in vielen Entscheidungsfällen eine risikoabwägende Wahl gut stützen (in eine Aktie aufgrund der vorhergesagten Entwicklung investieren oder nicht? einen Kunden am Self-Check-Out aufgrund seines Scanmusters kontrollieren oder nicht? eine Fussballspielerin ob ihrer prognostizierten Leistung verpflichten oder nicht?).

Auch andere, höhere Sicherheitswerte wären natürlich möglich und kontextabhängig auch sinnvoller: Bei der Vorhersage von Rentenerträgen würden viele von uns sicherlich gerne auch unwahrscheinlichere Szenarien (als die mittleren 50%) mit in das erwartungsleitende Prognoseintervall einbinden. Gleichwohl überlasten Vorhersageintervalle mit Sicherheitswerten von 95% und 99% – die sich also an einschlägige inferenzstatistische Standards anlehnen – oft die Erklärungskraft des Regressionsfit. Sie führen dann zu nichtssagenden oder gar irrealen Prognosen.

Gleichzeitig Grundlage und wichtigste Ableitung des Vorhersagebandes ist das Vorhersageintervall. Dieses gibt den Bereich um einen durch die Regression vorhergesagten y-Wert an, in dem sich die reale y-Ausprägung eines Falls mit bestimmter Sicherheit befindet. Vorhersageintervalle können über die einfache Vorhersagefunktion von ggpredict() abgerufen werden:

library (ggeffects)
ggpredict(fit, 
          terms = "eduyrs[2, 8, 14, 20, 26]", 
          interval = "prediction", 
          ci_level = 0.50)
## # Predicted values of Internet use, how much time on typical day, in minutes
## 
## eduyrs | Predicted |         50% CI
## -----------------------------------
##      2 |     88.07 | -20.99, 197.13
##      8 |    134.35 |  25.58, 243.11
##     14 |    180.62 |  71.89, 289.36
##     20 |    226.90 | 117.94, 335.86
##     26 |    273.18 | 163.74, 382.62
## 
## Intervals are prediction intervals. Use `interval = "confidence"` to
##   return regular confidence intervals.

Der tatsächliche tägliche Nutzungswert für eine Person mit 14 Bildungsjahren liegt mit 50% Sicherheit zwischen 71.9 und 289.4 Minuten. Praktisch: Der Output enthält eine Nachricht, die uns darauf hinweist, dass wir Vorhersageintervalle, keine Konfidenzintervalle, berechnet haben.

Das Vorhersageintervall ist ein Querschnitt aus dem Vorhersageband. Die obere und untere Grenzlinie des Vorhersagebandes ergibt sich also aus den verknüpften unteren und oberen Grenzen der Vorhersageintervalle. Achtung: Gestalt und Breite des Vorhersagebandes wird durch die Erklärungskraft des Regressionsfits und somit den Koeffizienten R2 dominiert. Es scheint daher oft linear - ist es aber nicht: Tatsächlich fliesst aber auch inferenzstatistische Information zur Parametersicherheit (welche in der Mitte der x-Variable grösser und an den Rändern schwächer ist) mit in das Vorhersageband bzw. die Vorhersageintervalle ein, daher sollten die (leider nicht automatisch in den ggplot integrierbaren) Grenzen des Vorhersagebandes nicht als linear modelliert werden (auch wenn sie meist linear erscheinen). Die hier vorgeschlagene Einbindungsroutine orientiert sich an der Integrationslogik der Konfidenzbänder über (geom_smooth) (bzw. den unter dessen Haube ablaufenden Routinen).

Zuerst müssen wir die obere und untere Grenze für das Vorhersageband für mehrere Werte von x berechnen und erstellen dazu ein neues Objekt predictions. Da die Grenzen eines Vorhersageintervalles nicht exakt linear Verlaufen, ist es redlich, zunächst möglichst viele, möglichst nah beianderliegende Vorhersageintervalle zu ermitteln (deren Grenzpunkte dann im zweiten Schritt zum Vorhersageband verknüpft werden).

predictions <- ggpredict(fit, 
                         terms = "eduyrs[0:26 by=0.05]", # Ermittle im Wertebereich von x insgesamt 520 Vorhersagen für y 
                         interval = "prediction", # Ergänze die Vorhersagen um die unteren und oberen Grenzen der Vorhersageintervalle 
                         ci_level = 0.50) # Die Vorhersageintervalle sollen eine reale Ausprägung mit 50% Sicherheit einfangen

Der hier erstellte Datensatz predictions enthält nun zu mehreren x-Werten entsprechende y-Vorhersagen sowie deren unteren und oberen Grenzwerte der jeweiligen Verhersageintervalle. Dieser werden nun per geom_ribbon() im ggplot() miteinander verknüpft und in eine durch sie umschlossene Fläche - das Vorhersageband - verwandelt.

plot2 <- plot1 +
  geom_ribbon(data = predictions, 
              aes(x = x, ymin = conf.low, ymax = conf.high),
              fill = "blue", 
              alpha = 0.1,
              inherit.aes = FALSE) +  
  labs(title = "Bildung und Internetnutzung in der Schweiz", 
       subtitle = "Regressionsgerade mit 
       95-Prozent-Konfidenzband und 50-Prozent Vorhersageband",
       y = "Internetnutzung in Minuten/Tag",
       x = "Bildungsjahre",
       caption = "ESS(2016), Teilstichprobe CH, N=1184.")
plot2

Hinweis: Das Vorhersageband wird in wissenschaftlichen Publikation – anders als das Konfidenzband – selten dargestellt (seine inhaltlichen Implikationen allerdings trotzdem häufig mit dem Konfidenzband verwechselt). Das Vorhersageband ist aber in vielen datenwissenschaftlichen Anwendungen wichtig, bei denen die Vorhersage (und nicht wie bei uns das Vorhersageprinzip bzw. der Hypothesentest) im Vordergrund steht (Bruce et al. 2020, Practical Statistics for Data Scientists, S. 163). Auch würde unter Praxisbedingungen ergänzend zu dem Hypothesentest noch eine Linearitätsprüfung und Ausreisseranalyse durchgeführt werden.

 

Hier gehts weiter zur Übung IV

 

logo.knit

Conforti, E., Siefart, F., De Min, N., Dürr, R., Hofer, L., Rauh, S., Senn, S., Strassmann Rocha, D., Giesselmann, M. (2023): “Regressionsanalysen mit R”