invisible header

link zum pdf

# Statistik 2: R Tutorat
# Übungsskript zum Thema Inferenzstatistik
# Datum: 26.04.2022
# 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 aufstellen

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_ss <- select(ess8_CH, internet = netustm, eduyrs, idno)
summary(ess8_CH_ss)
##     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_ss$eduyrs, na.rm = TRUE)
## [1] 3.496909
sd(ess8_CH_ss$internet, na.rm = TRUE)
## [1] 163.1974
ess8_noNA <- na.omit(ess8_CH_ss)

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. Durch Entfernung von Personen mit NAs beziehungsweise durch die Erstellung eines neuen Datensatzes ohne NAs mithilfe des na.omit-Befehls reduziert sich der Datensatz von 1525 auf 1184 Merkmalsträger. Wichtig ist, dass der Ausschluss von NAs erst ganz zum Schluss nach Erstellung des Teildatensatzes und Inspektion aller relevanten Variablen vorgenommen wird. Ansonsten besteht die Gefahr, dass auch analyserelevante Merkmalsträger aus dem Datensatz und somit aus der Regressionsanalyse ausgeschlossen werden.

2.2 Regressionsanalyse

Nun zur Regressionsanalyse mithilfe des lm()-Befehls.

fit <- lm(internet ~ eduyrs, data = ess8_noNA)
summary(fit)
## 
## Call:
## lm(formula = internet ~ eduyrs, data = ess8_noNA)
## 
## 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
## 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 (üblichweise 95%) liegt. Wir können die Grenzen mit dem Befehl confint ermitteln:

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. Genauer: Mit 95% Sicherheit steigt mit jedem zusätzlichen Bildungsjahr die Dauer der durchschnittlichen Internetnutzung in der GG zwischen 5 min 08 und 10 min 17. 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 95% 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_noNA, aes(x = eduyrs, y = internet))+
  geom_jitter(alpha = 0.2, height = 5) +
  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 = "Bildung und Internetnutzung in der Schweiz",
       subtitle = "Regressionsgerade mit 95%-Konfidenzband",
  y = "Internetnutzung in Minuten/Tag",
  x = "Bildungsjahre", 
  caption = "ESS(2016), Teilstichprobe CH, N = 1184")
plot1

Das Argument se = TRUE erzeugt das Konfidenzband, mit dem Argument fill können wir zudem seine Farbe wählen. Per Default wird ein 95-Prozent-Konfidenzband dargestellt, der Bereich kann aber über das level-Argument variiert werden (das ist allerdings unüblich).
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 y=0 und y=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).

3.3 Vorhersageband

Das 95% Vorhersageband zeigt den Bereich an, in dem 95% aller Werte der GG liegen bzw. in dem ein einzelener Wert der GG 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.

Wie das Konfidenzband können wir das Vorhersageband um unsere Regressionsgerade im ggplot darstellen, 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 alle Merkmalsträger im Datensatz berechnen und erstellen dazu ein neues Objekt, was wir hier predictions nennen:

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

Schliesslich können wir unseren ursprünglichen Plot plot1 mit den beiden Intervallgrenzen aus der hier generieren “predictions”-Tabelle ergänzen:

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

Querschnitte des Vorhersagebandes lassen sich auf konkrete Ausprägungen von X beziehen. Zur Ermittlung solcher Vorhersageintervalle nutzen wir ebenfalls den ggpredict()-Befehl:

ggpredict(fit, terms = "eduyrs[13]", interval = "prediction", ci.lvl = 0.5)
## # Predicted values of Internet use, how much time on typical day, in minutes
## 
## eduyrs | Predicted |          50% CI
## ------------------------------------
##     13 |    172.91 | [64.19, 281.63]

Der tatsächliche tägliche Nutzungswert für eine Person mit 13 Bildungsjahren liegt mit 50% Sicherheit zwischen 64.23 und 281.60 Minuten.
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”