invisible header

link zum pdf

# Statistik 2: R Tutorat
# Übungsskript zum Thema Inferenzstatistik
# Datum: 26.04.2024
# 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-Werts 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 Ablehung 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 complet~ 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
##  9   netustm  Internet use, how much time on ~ dbl+lbl  338    
##                                                                
##                                                                
##                                                                
##  values                
##  [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 vorhergesagte 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 eines typischen Unterschiedes zwischen zwei zufälligen Personen 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, color = "red", fill = "orange", 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.lvl = 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]", interval = "confidence", ci.lvl = 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]

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 vorrhergesagt.
Die dritte Spalte informiert zum Grenzverlauf bzw. zu den Grenzwerte der 95%-Bandes: Bei 2 Bildungsjahren verläuft das 95% KonfidenzbandBand zwischen 61 und 115 Minuten.

3.3 Vorhersageband

Das Vorhersageband zeigt den Bereich an, in dem 95% aller Werte liegen bzw. in dem ein einzelener Wert mit 95% Sicherheit liegt. Es berücksichtigt zusätzlich zur Unsicherheit ob der korrekten Regressionsgerade auch die Streuung der tatsächlichen Werte um die Regressionsgerade (daher die 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 hier niedriger (z.B. 50% oder 75%) anzusetzen. Wir können hierfür den ggpredict()-Befehl nutzen.

library (ggeffects)
ggpredict(fit, 
          terms = "eduyrs[2, 8, 14, 20, 26]", 
          interval = "prediction", 
          ci.lvl = 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]

Der tatsächliche tägliche Nutzungswert für eine Person mit 13 Bildungsjahren liegt mit 50% Sicherheit zwischen 64.23 und 281.60 Minuten.

Wie das Konfidenzband können wir das Vorhersageband um unsere Regressionsgerade mit ggplot anzeigen lassen, es braucht dafür jedoch noch einen zusätzlichen Schritt.

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. Wenn wir im Rahmen des terms-Argumentes keine konkreten x-Werte definieren, sucht der Befehl automatisch eine sinnvolle Anzahl an Wertepaaren aus:

predictions <- ggpredict(fit, 
                         terms = "eduyrs", 
                         interval = "prediction", 
                         ci.lvl = 0.50)

Der hier erstellte Datensatz predictions enthält nun zu mehreren X-Werten entsprechende Y-Vorhersagen sowie die unteren und oberen Y-Grenzwerte des Vorhersagebandes. Durch Verweis auf diese Datenpunkte im geom_smooth lassen sich nun oberer und unterer Grenzverlauf (bzw. optimierte Verlaufskurven über die Grenzwerte) im ursprünglichen Plot plot1 ergänzen. Per method=“loess” - gleichzeitig auch der Default des geom_smooth - lassen sich die optimierten Verlaufskurven einbinden:

plot2 <- plot1 +
  geom_smooth(method = "loess", 
              data = predictions, 
              aes(x = x, y = conf.high),
              size = 0.5, 
              color = "brown",
              linetype = "longdash")+
  geom_smooth(method = "loess", 
              data = predictions, 
              size = 0.5, 
              color = "brown",
              aes(x = x, y = conf.low), 
              linetype = "longdash")+
  labs(title = "Bildung und Internetnutzung", 
       y = "Internetnutzung in Minuten/Tag",
       x = "Bildungsjahre",
       caption = "ESS(2016), Teilstichprobe CH, N=1184.
       \n Regressionsgerade mit 95-Prozent-Konfidenzband
       und \n 50-Prozent Vorhersageband.")
plot2

Hinweis: Das Vorhersageband wird in wissenschaftlichen Publikation – anders als das Konfidenzband – selten dargestellt (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. 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”