invisible header

link zum pdf

# Statistik 2: R Tutorat
# Übungsskript zum Thema Linearität und Ausreisser
# Datum: 01.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(dplyr)
# install.packages("ggplot2")
library(ggplot2)
# install.packages("olsrr")
library(olsrr)
library(labelled)

 

1. Einleitung

In der Regressionsanalyse untersuchen wir den Zusammenhang zweier metrischer Variablen. Während die Regressionsgerade ein lineares Analyseinstrument ist, sind die damit beschriebenen Zusammenhänge nicht zwingend linear. Das stellt nicht zwingend ein Problem dar, da die Interpretation des Regressionskoeffizienten als “durchschnittlicher Effekt” in den meisten Fällen trotzdem Sinn ergibt. Insbesondere bei starken Abweichungen von der Linearität sind allerdings analytische Modifikationen der Regressionsanalyse geboten, die auf die Modellierung der Nicht-Linearität ausgerichtet sind. Und auch geringgradige Linearitätsabweichungen motivieren Forschende gelegentlich dazu, nach analytischen Alternativen zur Standardregression zu suchen. Als Teil der Dateninspektion bzw. -exploration sollte bei analytischem Interesse an zwei metrischen Variablen immer auch die Form des Zusammenhangs - bzw. der Grad an Linearitätsabweichung - überprüft werden.

Dabei spielen zunächst theoretische Überlegungen eine Rolle. Ist es plausibel, von einem zumindest näherungsweise linearen Zusammenhang zwischen UV und AV auszugehen? Oder gibt es gute theoretische Gründe für die Annahme einer Linearitätsabweichung? In einem zweiten Schritt können diese theoretischen Überlegungen anhand einer visueller Inspektion der bivariaten Verteilung (auf Basis des Streudiagramms) eingeordnet, plausibilisiert oder entkräftet werden. Schliesslich beinhaltet die Exploration der Zusammenhangsform auch einen analytischen Teil, in dem Linearitätsabweichungen auf Basis statistischer Parameter empirisch identifiziert werden.

Im Folgenden beschäftigen wir uns mit dem Zusammenhang zwischen dem Lebensalter und der Klimaverantwortung in der schweizer Bevölkerung. Wir fragen uns, welche Form dieser Zusammenhang hat.

 

2. Data Management

2.1 Selektieren von Merkmalsträgern

Wir führen die Analyse auf Basis der Schweizer Teilstichprobe durch.

ess8_ch <- filter(ess8, cntry == "CH")
dim(ess8_ch)
## [1] 1525  534

Unser Datensatz reduziert sich daher auf 1525 Merkmalsträger.

 

2.2 Inspektion von Variablen

In unseren Analysen arbeiten wir mit den Variablen, agea, ccrdprs, eduyrsund netustm. Wir können entweder im Codebook nachschauen, was diese Variablen messen, oder uns über den attributes()- oder look_for()-Befehl einen raschen Überblick verschaffen.

attributes(ess8_ch$agea)
attributes(ess8_ch$ccrdprs)
attributes(ess8_ch$eduyrs)
attributes(ess8_ch$netustm)

Wir stellen fest:

  • agea misst das Alter in Jahren.

  • ccrdprs misst auf einer 10er-Skala, wie stark sich eine Person dafür verantwortlich fühlt, zur Reduktion des Klimawandels beizutragen. Dabei ist 0 = ‘Not at all’ und 10 = ‘A great deal’. Folglich entsprechen zunehmende Werte grösserer Zustimmung, sprich die Variable ist positiv gepolt.

  • eduyrs misst die Anzahl Bildungsjahre einer Person.

  • netustm misst die Zeit (in Minuten), die eine Person täglich im Internet verbringt.

  • Wir sehen zusätzlich, dass bei allen Merkmalen alle ungültigen Werte als NAs kodiert sind.

     

2.3 Selektieren von Merkmalen

Im Folgenden erstellen wir zwei separate Teildatensätze. Im ersten Teil werden wir uns mit dem Teildatensatz beschäftigen, welcher den Identifyer idno sowie die Variablen agea und ccrdprs beinhaltet. Für den zweiten Teil (Ausreisseranalyse) erstellen wir einen Teildatensatz, welcher den Identifyer idno sowie die Variablen eduyrs und netustm umfasst.

ess8_ch_ss_1 <- select(ess8_ch, identifyer = idno,
                     alter = agea, 
                     klima_ver = ccrdprs)
ess8_ch_ss_2 <- select(ess8_ch, identifyer = idno,
                     eduyrs, internet = netustm)

Da wir jetzt wissen worum es sich bei den einzelnen Merkmalen handelt, können wir ihnen prägnantere Namen zuweisen. Statt in einem zusätzlichen Schritt mit dem rename()-Befehl, können wir dies direkt über den select()-Befehl machen. Die Struktur ist dabei innerhalb des Befehls gleich wie im rename-Befehl: select (daten, neuer name = alter name)

 

3. Univariate Statistik

Die Univariate Statistik ist wichtiger vorbereitender Schritt der analytischen Inspektion und Aufbereitung eines Zusammenhangs. Sie prüft die korrekte Codierung und plausible Verteilung der beteiligten Variablen, gibt Hinweise auf mögliche Ausreisserproblematiken (siehe Punkt 7.) und hilft bei der Einschätzung von Effektgrössen sowie der Organisation der Skalen in den graphischen Darstellungen.

 

3.1 Stichprobenstatistik

Wir ermitteln als erstes die univariaten Paramter, welche die zentralen Verteilungseigenschaften unserer Variablen aufgreifen:

summary(ess8_ch_ss_1)
##    identifyer       alter         klima_ver     
##  Min.   :   1   Min.   :15.00   Min.   : 0.000  
##  1st Qu.:1017   1st Qu.:32.00   1st Qu.: 6.000  
##  Median :1827   Median :48.00   Median : 7.000  
##  Mean   :1807   Mean   :47.83   Mean   : 6.861  
##  3rd Qu.:2645   3rd Qu.:62.00   3rd Qu.: 8.000  
##  Max.   :3460   Max.   :94.00   Max.   :10.000  
##                 NA's   :6       NA's   :39
sd(ess8_ch_ss_1$alter, na.rm = TRUE)
## [1] 18.78213
sd(ess8_ch_ss_1$klima_ver, na.rm = TRUE)
## [1] 2.235231
summary(ess8_ch_ss_2)
##    identifyer       eduyrs        internet     
##  Min.   :   1   Min.   : 0.0   Min.   :   3.0  
##  1st Qu.:1017   1st Qu.: 9.0   1st Qu.:  60.0  
##  Median :1827   Median :10.0   Median : 120.0  
##  Mean   :1807   Mean   :11.3   Mean   : 163.6  
##  3rd Qu.:2645   3rd Qu.:13.0   3rd Qu.: 240.0  
##  Max.   :3460   Max.   :26.0   Max.   :1200.0  
##                 NA's   :3      NA's   :338
sd(ess8_ch_ss_2$eduyrs, na.rm = TRUE)
## [1] 3.496909
sd(ess8_ch_ss_2$internet, na.rm = TRUE)
## [1] 163.1974
Im Durchschnitt haben SchweizerInnen in der Stichprobe ein Alter von 47.83 Jahren mit einer Standardabweichung von 18.78 Jahren, eine Klimaverantwortung von 6.86 mit einer Standardabweichung von 2.24, eine Ausbildung von 11.30 Jahren mit einer Standardabweichung von 3.50 Jahren und verbringen an einem Tag durchschnittlich 163.60 Minuten (ca. 2.7 Stunden) im Internet mit einer nahezu gleich grossen Standardabweichung von 163.20 Minuten.

 

3.2 Graphische Darstellung

Nun wollen wir die univariaten Verteilungen graphisch darstellen. Die Verteilungen der metrischen Variablen eduyrs, alter und internet stellen wir im Histogramm dar. Für die ordinalskalierte, “quasi-metrische” Variable klima_ver nehmen wir ein Säulendiagramm. Diese univariaten Darstellungen sind ebenfalls Bestandteil der Dateninspektion. Sie helfen, die analyserelevanten Variablen besser kennenzulernen und ggf. Probleme der Variablen zu identifzieren:

ggplot(ess8_ch_ss_1, aes(alter))+
  geom_histogram(bins = 15, fill = "lightblue", color = "blue")+
  theme_bw()+
  scale_x_continuous(breaks = seq(15,100,5))+
  labs(title = "Verteilung des Alters", 
       caption = "Data: ESS(2016), Teilstichprobe CH, N = 1525", 
       y = "Häufigkeit", 
       x ="Alter")

ggplot(ess8_ch_ss_1, aes(klima_ver))+
  geom_bar(aes(y = ..count../sum(..count..)), fill = "lightblue", color = "blue")+
  scale_y_continuous(labels = scales::percent)+
  scale_x_continuous(breaks = seq(0,10), labels = c("Gar nicht","1", "2", "3", "4", "5", "6", "7", "8", "9", "Sehr"))+
  labs(title ="Ich fühle mich verantwortlich dafür den Klimawandel zu reduzieren",
       caption = "Data: ESS(2016), Teilstichprobe CH, N = 1525", 
       y = "Häufigkeit in Prozenten", 
       x ="Klimaverantwortung")+
  theme_bw()

ggplot(ess8_ch_ss_2, aes(eduyrs))+
  geom_histogram(bins = 15, fill = "lightblue", color = "blue")+
  theme_bw()+
  scale_x_continuous(breaks = seq(0,30,3))+
  labs(title ="Verteilung der Bidlungsjahre", 
       caption = "Data: ESS(2016), Teilstichprobe CH, N = 1525", 
       y = "Häufigkeit", 
       x ="Bildungsjahre")

ggplot(ess8_ch_ss_2, aes(internet))+
  geom_histogram(bins = 20, fill = "lightblue", color = "blue")+
  theme_bw()+
  scale_x_continuous(breaks = seq(0,1200,120), minor_breaks = seq(0,1200,60))+
  labs(title ="Verteilung der täglichen Internetnutzung", 
       caption = "Data: ESS(2016), Teilstichprobe CH, N = 1525", 
       y = "Häufigkeit", 
       x ="tägliche Internetnutzung in Minuten")

 

4. Linearitätsdiagnose

4.0 Vorüberlegungen

Welche verschiedenen Zusammenhangsformen sind auf Basis der beiden Variablen alter und klima_ver denkbar?:

  • Zunehmendes Alter könnte möglicherweise zu einer grösseren Klimaverantwortung führen. Begründen könnten wir diese Vermutung damit, dass erst mit zunehmendem Wissen zur Situation und der Auswirkungen der eignenen Aktionen sich das Gefühl der persönlichen Verantwortung entwickeln kann. Mit zunehmendem Alter würde zudem auch die Verantwortung für die nächsten Generationen (Kinder und Enkelkinder) hinzukommen.

  • Eine weitere Vermutung könnte sein, dass zunehmendes Alter zu geringerer Klimaverantwortung führt. Dies wäre begründbar damit, dass ältere Menschen weniger mit den Konsequenzen des Klimawandels zu kämpfen haben werden und daher nicht so direkt damit konfrontiert sind als jüngere. Somit könnte das Verantwortlichkeitsgefühl mit zunehmendem Alter schwächer werden.

  • Es wäre allerdings auch ein nicht-linearer Zusammenhang denkbar: Die Lebenszufriedenheit könnte beispielsweise bis zu einem gewissen Alter -z.B. bis typischerweise der Auszug der eigenen Kinder erfolgt - zunehmen und dann wieder abnehmen.

 

4.1 Visuelle Inspektion

(a) Streudiagramm

Wir möchten zur Inspektion den Zusammenhang aussagekräftig visualisieren und starten mit einem einfachen Streudiagramm.

ggplot(ess8_ch_ss_1, 
         aes(x = alter, y = klima_ver))+
  geom_point()+
    scale_x_continuous(breaks = seq(15,100,5))+
    scale_y_continuous(breaks = seq(0,10,1))+
    theme_classic()+
    labs(title = "Klimaverantwortung nach Alter", 
         y = "Klimaverantwortung",
         x = "Alter",
         subtitle = "Ich fühle mich gar nicht (0) - stark (10) verantwortlich dafür den Klimawandel zu reduzieren",
         caption = "ESS(2016), Teilstichprobe CH, N = 1525")

Das Streudiagramm hat aufgrund von Overplotting keine Aussagekraft, wir modifizieren daher mit den üblichen Tricks ‘jitter’ und ‘alpha’:

ggplot(ess8_ch_ss_1, 
         aes(x = alter, y = klima_ver))+
    geom_jitter(alpha = 0.2, size =2)+
    scale_x_continuous(breaks = seq(15,100,5))+
    scale_y_continuous(breaks = seq(0,10,1))+
    geom_smooth(method = "lm", se = F, color = "red")+
    theme_classic()+
    labs(title = "Klimaverantwortung nach Alter", 
         y = "Klimaverantwortung",
         x = "Alter",
         subtitle = "Ich fühle mich gar nicht (0) - stark (10) verantwortlich dafür den Klimawandel zu reduzieren",
         caption = "ESS(2016), Teilstichprobe CH, N = 1525")

Die Verteilung der Punkte lässt keinen eindeutigen Zusammenhang erkennen. Die horizontal verlaufende Regressionsgerade zeigt an, dass die durchschnittliche Steigung der Klimaverantwortung mit zunehmendem Alter näherungsweise Null ist (vgl. lm(klima_ver ~ alter, data = ess8_ch_ss_1). Dass sich hinter diesem Koeffizienten ein parabelartiger bzw. nicht-linearer Zusammenhang verbirgt, liegt zwar aufgrund des visuellen Befundes nicht unbedingt nah, kann aber auch nicht ausgeschlossen werden.

(b) Binned Scatterplot

Der Binned Scatterplot ist eine weitere Alternative zum klassischen Streudiagramm, die insb. im Kontext der Linearitätsdiagnose oft hilfreich ist.

library(binsreg)
ess_sample_bins<-as.data.frame(ess8_ch_ss_1)
ess_sample_bins$klima_mean <- as.numeric(ess_sample_bins$klima_ver)
ess_sample_bins$alter_g <- as.numeric(ess_sample_bins$alter)
binsreg(data = ess_sample_bins, x = alter_g, y = klima_mean) 
## Call: binsreg
## 
## Binscatter Plot
## Bin/Degree selection method (binsmethod)  =  IMSE direct plug-in (select # of bins)
## Placement (binspos)                       =  Quantile-spaced
## Derivative (deriv)                        =  0
## 
## Group (by)                         =  Full Sample
## Sample size (n)                    =  1482
## # of distinct values (Ndist)       =  79
## # of clusters (Nclust)             =  NA
## dots, degree (p)                   =  0
## dots, smoothness (s)               =  0
## # of bins (nbins)                  =  12

Die Datenpunkte repräsentieren hier nicht mehr einzelne Personen, sondern altersgruppenspezifische Mittelwerte. Der Befehl binsreg() bildet also zunächst im Hintergrund gleich grosse Gruppen von Einheiten auf Basis der UV. Hier sind dies die Altersgruppen 15-21, 22-26, 27-33, … 75-94. Beachte, dass die Breite der Klassenintervalle varriert, und zwar je nach Belegungsdichte im jeweiligen Altersspektrum. Damit stellt der Befehl sicher, dass eine sinnvolle Anzahl an Gruppen mit identischer Gruppengrösse gebildet wird. Der Default zur Anzahl an Gruppen bzw. “Bins” kann dabei durch zusätzliche Argumente (siehe Dokumentation zum Befehl) modifiziert werden.

Die einzelnen Punkte repräsentieren nun die Mittelwerte der AV innerhalb der jeweiligen Altersgruppen. Unter den 15- bis 21-jährigen liegt der Mittelwert der Klimaverantwortung demnach bei etwa 6,5, unter den 38- bis 43-jährigen bei etwa 7,3 usw.. Durch die starke Kompression der Dateninformation und der damit Verbundenen Einebnung der Streuung der AV kristallisiert sich auf Basis dieses geplotteten Mittelwertvergleichs nun doch relativ eindeutig ein umgekehrt U-förmiger Verlauf heraus.

Die Schwächen des binsreg- Befehl sind einerseits die unflexible Variablenverarbeitung (beide Variablen mussten vorab explizit als “numeric” definiert werden, obgleich sie schon entsprechend angelegt waren), die etwas zu starke Verengung der y-Achse, die den Verlauf überdramatisiert darstellt, sowie - damit verbunden - die beschränkten Möglichkeiten zur graphischen Aufbereitung und Annotation. So gibt es z.B. keine Möglichkeit, Titel etc. hinzuzufügen. Dazu muss relativ aufwendig zunächst ein ggplot-Objekt auf Basis des binsreg-Objektes kreiert werden:

# Convert variables to numeric
alter_num <- as.numeric(as.character(ess_sample_bins$alter))
klima_ver_num <- as.numeric(as.character(ess_sample_bins$klima_ver))

# Run the binned scatterplot regression analysis and store the output
binsreg_output <- binsreg(data = ess_sample_bins, x = alter_num, y = klima_ver_num)
# Create the ggplot object from the binsreg output
bins_plot <- ggplot(data = ggplot_build(binsreg_output$bins_plot)$data[[1]], aes(x = x, y = y)) +
  geom_point(size=3) +
  ggtitle ("Binned Scatterplot: Mittlere Klimaverantwortung nach Altersgruppen") +
  ylim(4, 9) +
    labs(y = "Mittlere Klimaverantwortung",
         x = "Alter (Gruppiert)",
         caption = "ESS(2016), Teilstichprobe CH, N = 1525")
bins_plot

 

4.3 Multigruppenanalyse

Im Rahmen unserer Vorüberlegungen (vgl. 4.0) haben wir festgehalten, dass es sich beim Alter und der Klimaverantwortung auch um einen nicht-linearen Zusammenhang handeln könnte, bei dem der Alterseffekt in verschiedenen Bereichen der Altersverteilung unterschiedlich gerichtet ist. Hierfür bietet zwar das Streudiagram keinen unmittelbaren Anhaltspunkt, trotzdem ist es im Rahmen der Linearitätsexploration immer sinnvoll, auch analytische Techniken anzuwenden, da auch die beste Visualisierung bestimmte Nuancen der gemeinsamen Verteilung oft verschleiert - zumal der Binned Scatterplot relativ deutlich einen umgekehrt U-förmigen Zusammenhang aufzeigt.

Zur analytischen Inspektion nutzen wir das Diagnoseintstument der Multigruppenanalyse. Wir vergleichen dafür nun in getrennten Analysen den Alterseffekt in verschiedenen Bereichen der Altersverteilung. Als Erstes müssen wir dazu festlegen, bei welchem Wert der Altersvariable wir den cut setzen bzw. wo wir unseren Datensatz aufteilen. Dazu haben wir verschiedene Möglichkeiten:

  • Theoretische Begründung (Gibt es einen theoretisch begründen Altersscheitelpunkt bei der Klimaverantwortung?)

  • Datengetriebene Identifikation: Scheitelpunkt (Eignet sich natürlich nur, wenn dieser auch visuell erkennbar ist - was hier nicht der Fall ist)

  • Konvention: Mean

  • Konvention: Median (Hat den Vorteil, dass dadurch (etwa) zwei gleich grosse Gruppen entstehen.)

  • Konvention: Quartile (…also mehr als zwei Gruppen: erlaubt es feinere Facetten der Zusammenhangsform zu erkennen.)

Für unsere jetzige Analyse nehmen wir den Median. Das heisst: Wir erstellen zwei neue Teildatensätze ess8_split1 und ess8_split2, die nur Merkmalsträger enthalten, welche entweder höchstens das Medianalter haben oder oberhalb des Medians liegen. Das machen wir mit dem bereits bekannten filter() Befehl - Achtung: Bei UV die wenige Skalenpunkten haben (10-er Skala) oder sehr schwerpunktlastig sind (z.B. Bildungsjahre) ist es oft sinnvoll, den Median in beide Splitts mit einzubeziehen!).

summary(ess8_ch_ss_1$alter)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max.    NA's 
##   15.00   32.00   48.00   47.83   62.00   94.00       6
  ess8_split1 <- filter(ess8_ch_ss_1, alter <= 48)
  ess8_split2 <- filter(ess8_ch_ss_1, alter > 48)

Jetzt können wir zwei separate, altersgruppenspezifische Regressionsmodelle berechnen, vergleichen und auf dem Scatterplot abbilden.

model_split1 <- lm(klima_ver ~ alter, data = ess8_split1) 
model_split2 <- lm(klima_ver ~ alter, data = ess8_split2) 
  
summary(model_split1)
## 
## Call:
## lm(formula = klima_ver ~ alter, data = ess8_split1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.2566 -1.0880  0.3195  1.4460  3.6427 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.90764    0.26513  22.282  < 2e-16 ***
## alter        0.02810    0.00784   3.584 0.000359 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.163 on 756 degrees of freedom
##   (18 Beobachtungen als fehlend gelöscht)
## Multiple R-squared:  0.01671,    Adjusted R-squared:  0.01541 
## F-statistic: 12.85 on 1 and 756 DF,  p-value: 0.0003594
summary(model_split2)
## 
## Call:
## lm(formula = klima_ver ~ alter, data = ess8_split2)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.5461 -1.1720  0.4539  1.5053  3.9089 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  9.583184   0.535634  17.891  < 2e-16 ***
## alter       -0.041573   0.008254  -5.037 5.99e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.254 on 722 degrees of freedom
##   (19 Beobachtungen als fehlend gelöscht)
## Multiple R-squared:  0.03394,    Adjusted R-squared:  0.03261 
## F-statistic: 25.37 on 1 and 722 DF,  p-value: 5.986e-07
ggplot(ess8_ch_ss_1, aes(alter, klima_ver))+
    geom_jitter(alpha = 0.2, size =2)+
    scale_x_continuous(breaks = seq(15,100,5))+
    scale_y_continuous(breaks = seq(0,10,1))+
    geom_smooth(method = "lm", se = F, color = "green", data = ess8_split1)+
    geom_smooth(method = "lm", se = F, color = "red", data = ess8_split2)+ 
    theme_classic()+ labs(title = "Klimaverantwortung nach Alter", 
                  y = "Klimaverantwortung", x = "Alter",
                  subtitle= "Ich fühle mich gar nicht (0) - stark (10) verantwortlich dafür den Klimawandel zu reduzieren",
         caption = "grün: Alter<=43,n=581 / rot: Alter>43, n=574. ESS(2016), 
         Teilstichprobe CH, N = 1525")

Vergleichen wir die Regressionskoeffizienten der beiden Teilgruppen, so bemerken wir: Das Vorzeichen der (statistisch signifikanten) Koeffizienten unterscheidet sich zwischen den beiden Altersintervallen. Der Zusammenhang zwischen Alter und Klimaverantwortung ist positiv in der ersten Hälfte der Altersverteilung und negativ in der zweiten Hälfte der Altersverteilung. Man kann also sagen, dass die Klimaverantwortung bei jungen Personen mit jedem zusätzlichen Lebensjahr steigt und bei älteren Personen mit jedem zusätzlichen Lebensjahr sinkt. Diese Beobachtung deutet darauf hin, dass es zwischen dem Alter und der Klimaverantwortung sehr wohl einen Zusammenhang gibt. Es handelt sich dabei allerdings um einen nicht-linearen Zusammenhang. Dieser scheint zudem substantieller Natur zu sein, daher müssen wir ihn bei der regressionsbasierten Modellierung des Zusammenhangs berücksichtigen.

 

5. Modellierung der Linearitätsabweichung

Wird eine Abweichung von der Linearität diagnostiziert und für relevant befungen, gibt es im Rahmen der Regressionsanalyse mehrere spezifische Analyseoptionen. Neben der Multigruppenanalyse (die zugleich Diagnose- wie auch Modellierungstool ist) gehören die intrinsische Linearisierung durch Quadrierung oder Logarithmierung zu den prominentesten Varianten.

Wir entscheiden uns in diesem Fall für eine Quadrierung, weil wir uns aus theoretischen Überlegungen einen umgekehrt U-förmigen Zusammenhang vorstellen können. Diese Entscheidung ist zudem durch die analytische Inspektion per Multigruppenanalyse gut abgestützt.

 

5.1. Regression mit quadriertem Term

Wir erstellen in unserem Datensatz eine neue Variable mit dem quadrierten Alter.

 ess8_ch_ss_1$alter_sqr <- ess8_ch_ss_1$alter^2

Jetzt berechnen wir das neue Modell, indem wir im lm()-Befehl die quadrierte Altersvariable ergänzen.

# model_sqr <- lm(AV~UV+UV_quadriert, data = data)
model_sqr <- lm(klima_ver~alter_sqr+alter, data = ess8_ch_ss_1)
summary(model_sqr)
## 
## Call:
## lm(formula = klima_ver ~ alter_sqr + alter, data = ess8_ch_ss_1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.2294 -1.1417  0.3478  1.5043  4.0435 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.7966541  0.3530216  13.587  < 2e-16 ***
## alter_sqr   -0.0010240  0.0001584  -6.465 1.37e-10 ***
## alter        0.0998233  0.0156638   6.373 2.47e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.206 on 1479 degrees of freedom
##   (43 Beobachtungen als fehlend gelöscht)
## Multiple R-squared:  0.0275, Adjusted R-squared:  0.02619 
## F-statistic: 20.91 on 2 and 1479 DF,  p-value: 1.106e-09

Anstatt extra eine neue quadrierte Variable zu generieren, können wir auch im lm-Befehl automatisch eine entsprechende Variable anlegen und integrieren lassen.

model_sqr1 <- lm(klima_ver ~ alter + I(alter^2), data = ess8_ch_ss_1)
summary(model_sqr1)
## 
## Call:
## lm(formula = klima_ver ~ alter + I(alter^2), data = ess8_ch_ss_1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.2294 -1.1417  0.3478  1.5043  4.0435 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.7966541  0.3530216  13.587  < 2e-16 ***
## alter        0.0998233  0.0156638   6.373 2.47e-10 ***
## I(alter^2)  -0.0010240  0.0001584  -6.465 1.37e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.206 on 1479 degrees of freedom
##   (43 Beobachtungen als fehlend gelöscht)
## Multiple R-squared:  0.0275, Adjusted R-squared:  0.02619 
## F-statistic: 20.91 on 2 and 1479 DF,  p-value: 1.106e-09

Ausserdem können wir die quadratische Regressionskurve auch grafisch leicht darstellen. Hierfür fügen wir unserem ggplot() folgendes hinzu: geom_smoth(method = "lm", formula = y~poly(x,2))

 ggplot(ess8_ch_ss_1, aes(alter, klima_ver))+
     geom_jitter(alpha = 0.2, size =2)+
     scale_x_continuous(breaks = seq(15,100,5))+
     scale_y_continuous(breaks = seq(0,10,1))+
     geom_smooth(method = "lm", se = F, color = "blue", formula = y ~ poly(x,2))+
   theme_classic()+
     labs(title = "Klimaverantwortung nach Alter (mit quadratischem Fit)", 
         y = "Klimaverantwortung",
         x = "Alter",
         subtitle = "Ich fühle mich gar nicht (0) - stark (10) verantwortlich dafür den Klimawandel zu reduzieren",
          caption = "ESS(2016), Teilstichprobe CH, N=1155")

Auf Grundlage der Koeffizienten und der Grafik können wir Aussagen über den Zusammenhang treffen:

  • Der Befund eines nicht-linearen Zusammenhangs erhärtet sich: Der Koeffizient des quadrierten Terms ist signifikant von 0 verschieden, das R^2 der quadratischen Funktion substanziel gegenüber der linearen erhöht: Das model_sqr erklärt 2.75% der gesamten Varianz der Klimaverantwortung, während dies beim linearen fit deutlich weniger als 1% sind. Insofern hat das Modell mit quadriertem Term eine deutlich bessere Vorhersagekraft. Wir können also davon ausgehen, dass dieses Modell den vorliegenden Zusammenhang besser beschreibt. Merke: Das R^2 kann ebenfalls in vergleichender Perspektive als Linearitätsdiagnostetool eingesetzt werden!
  • Der negative Koeffizient der quadrierten Altersvariable verweist auf eine umgekehrt U-förmige Schätzkurve.

Hinweis: Hätten wir uns stattdessen für die Logarithmierung entschieden, so könnten wir einen logarithmischen Zusammenhang mit formula = y~ln(x) statt formula = y~poly(x,2) visualisieren.

 

5.2 Interpretationsstrategien

Um die Eigenschaften der quadratischen Regressionskurve zu beschreiben, können wir den Scheitelpunkt der Parabel bestimmen. Wir definieren dazu die Funktion und ermitteln dann das Maximum per Ableitung oder mittels des optimize- Befehls in R.

Dazu erstellen wir eine Funktion anhand der Regressionsgleichung:

alter_funktion <- function(x){4.80+0.10*x-0.0010*x^2}

Diese können wir von Hand ableiten, gleich Null setzen und nach x auflösen. Oder automatisiert in R:

optimize(alter_funktion, interval = c(15, 95), maximum = T)
## $maximum
## [1] 50
## 
## $objective
## [1] 7.3

Wichtig: Man muss immer einen Bereich interval= bestimmen, in dem ein Scheitelpunkt gesucht werden soll. Nimm jeweils einen Bereich, von dem du dir sicher bist, dass er den Scheitelpunkt enthält. Wichtig ist auch, dass du richtig angibst ob ein Maximum oder Minimum gesucht werden soll.

Der Scheitelpunkt liegt bei 50. Dies bedeutet, dass 50 jährige Personen im Mittel die grösste Verantwortung gegenüber dem Klima verspüren.

Die Form des Zusammenhangs können wir zudem über geschickt ausgewählte Vorhersagen verfügbar machen. Diese können wiederum auch nochmals helfen, die Bdeutung der quadratischen Modellierung im vergleich zur linearen zu verdeutlichen. Ausgehend von Scheitelpunkt sowie ein paar Aussagekräftigen und relevante Altersstufen lässt sich vorhersagebasiert ein guter Eindruck von der Form des Zusammenhangs sowie der Effektgrösse vermitteln.

  library(ggeffects)
  sd (ess8_ch_ss_1$klima_ver, na.rm=TRUE)
## [1] 2.235231
  ggpredict(model_sqr1, terms = "alter[20, 35, 50, 65, 80]")
## # Predicted values of To what extent feel personal responsibility to reduce climate change
## 
## alter | Predicted |       95% CI
## --------------------------------
##    20 |      6.38 | [6.14, 6.63]
##    35 |      7.04 | [6.89, 7.18]
##    50 |      7.23 | [7.07, 7.39]
##    65 |      6.96 | [6.80, 7.11]
##    80 |      6.23 | [5.93, 6.53]

Der quadratische fit sagt das Maximum der Klimaverantwortung für 50-jährige voraus und eine um fast einen Skalenpunkt (bzw. mehr 0,3 Standardabweichungseinheiten) geringere Klimaverantwortung an den Rändern der Altersverteilung. Obgleich der Verantwortungsanstieg also bis zum 50. Lebensjahr im Mittel stetig positiv ist, ist er im jungen Erwachsenenleben (Anstieg um zwei Drittel Skalenpunkte zwischen 20 und 35) deutlich stärker ausgebildet als im mittleren Alter (Anstieg um weniger als einen Fünftel Skalenpunkt zwischen 35 und 50). Auf der anderen Seite der Verteilung gilt dies Umgekehrt. Ausweislich des quadratischen Regressionsfits zeigt sich also eine relativ hohes und stabiles Plateau der Klimaverantwortung in der Lebensmitte (35-65), gerahmt von einem fulminaten Auf- und Abstieg in den rahmenden, frühen und späten Lebensphasen.

Zuletzt vergleichen wir noch die R-Quadrate von linearer und quadratischer Regression:

model_lin<-lm(klima_ver ~ alter, data = ess8_ch_ss_1) 
summary(model_lin)
## 
## Call:
## lm(formula = klima_ver ~ alter, data = ess8_ch_ss_1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.8866 -0.8726  0.1424  1.1429  3.1496 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 6.8426286  0.1586060  43.142   <2e-16 ***
## alter       0.0004828  0.0030855   0.156    0.876    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.236 on 1480 degrees of freedom
##   (43 Beobachtungen als fehlend gelöscht)
## Multiple R-squared:  1.654e-05,  Adjusted R-squared:  -0.0006591 
## F-statistic: 0.02448 on 1 and 1480 DF,  p-value: 0.8757
summary(model_sqr1)
## 
## Call:
## lm(formula = klima_ver ~ alter + I(alter^2), data = ess8_ch_ss_1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.2294 -1.1417  0.3478  1.5043  4.0435 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.7966541  0.3530216  13.587  < 2e-16 ***
## alter        0.0998233  0.0156638   6.373 2.47e-10 ***
## I(alter^2)  -0.0010240  0.0001584  -6.465 1.37e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.206 on 1479 degrees of freedom
##   (43 Beobachtungen als fehlend gelöscht)
## Multiple R-squared:  0.0275, Adjusted R-squared:  0.02619 
## F-statistic: 20.91 on 2 and 1479 DF,  p-value: 1.106e-09
Das R-Quadrat der quadratischen Regression ist deutlich Grösser als das der linearen und validiert somit zusätzlich die quadratische Regression in diesem Fall als angemessen.

 

7. Ausreisser

Ob ein Modell einen tatsächlich vorliegenden Zusammenhang beschreibt hängt immer auch von dem verwendeten Datensatz ab. Es kann vorkommen, dass vereinzelte Ausprägungen (sogenannte Ausreisser) eine besonders starke Streuung vom Erwartungswert vorweisen und somit das Regressionsmodell auch besonders stark beeinflussen. Im folgenden Abschnitt wird beschrieben, wie du ein Regressionsmodell auf Ausreisser testen kannst. Da die im vorherigen Abschnitt verwendeten Variablen keine grosse Streuung zulassen, verwenden wir für die Ausreisseranalyse im Folgenden die Variablen eduyrs und netustm.
Wir berechnen zuerst das lineare Modell.

net_model <- lm(internet ~ eduyrs, data = ess8_ch_ss_2)
net_model
## 
## Call:
## lm(formula = internet ~ eduyrs, data = ess8_ch_ss_2)
## 
## Coefficients:
## (Intercept)       eduyrs  
##      72.646        7.713

 

7.1 Inspektion des Streudiagramms

Wir inspizieren zunächst visuell, ob im Streudiagramm Ausreisser zu erkennen sind. Wir benutzen hier einen Scatterplot, da es sich bei Ausreissern um einzelne Messungen handelt, die in Heatmap nicht gut repräsentiert werden.

ggplot(ess8_ch_ss_2, 
         aes(x = eduyrs, y = internet))+
    geom_jitter(alpha = 0.2, size = 2)+
    scale_x_continuous(breaks = seq(0,26,2))+
    scale_y_continuous(breaks = seq(0,1200,100))+
    geom_smooth(method = "lm", se =F, color = "red")+
    theme_classic()+
    labs(title = "Zeit im Internet nach Bildung", 
         caption = "ESS(2016), Teilstichprobe CH, N = 1184", 
         y = "Zeit im Internet pro Tag in Minuten", 
         x = "Bildungsjahre")
## Warning: Removed 341 rows containing non-finite values (`stat_smooth()`).
## Warning: Removed 341 rows containing missing values (`geom_point()`).

Zwischen 10 und 20 Jahren sind einige relativ hohe Werte für Internetnutzung zu erkennen. Bei diesen könnte es sich möglicherweise um Ausreisser handeln, welche die Robustheit bzw. Stabilität von Regressionsergebnissen einschränken und über relevante Sachverhalte in den Daten hinwegtäuschen könnten.

 

7.2 Graphische Darstellung der dfbetas

Formal betrachen wir solche Messungen als Ausreisser, die einen unbotmässig grossen Einfluss auf das Regressionsergebnis haben. Den Einfluss einer Messung auf den Regressionskoeffizienten bestimmen wir über das standardisierte Einflussmass dfbetas. Überschreitet das dfbetas einer Messung die kritische Grenze nach Belsley et al. - also \(\frac{2}{\sqrt{n}}\) - verweist dies auf einen übermässigen Einfluss der entsprechenden Messung und bezeugt somit deren Aussreissereigenschaft. Unser Teildatensatz beinhaltet 1184 Merkmalsträger mit gültigen Werten auf beiden Variablen (n = 1184). Der Grenzwert beträgt folglich etwa +-0.06. Mit dem Befehl ols_plot_dfbetas vom Package olsrr können wir uns schnell graphisch einen Überblick von möglichen Ausreissern verschaffen.

library(olsrr)
ols_plot_dfbetas(net_model)

Wir sehen unsere Vermutung nach der visuellen Inspektion im Streudiagramm bestätigt: Es gibt einige Ausreisser (Punkte die den Grenzwert überschreiten). Folglich ist eine Re-Analyse ohne Ausreisser sinnvoll, um die Robustheit der Regressionsanalyse zu überprüfen.

 

7.3 Datensatz mit dfbetas erweitern

Mit folgendem Befehl ergänzen wir unseren Teildatensatz um eine neue Variable, welche allen Merkmalsträgern mit gültigen Werten das standardisierte Einflussmass dfbetas zuweist:

ess8_ch_mitdfb <- cbind(filter(ess8_ch_ss_2, !is.na(eduyrs),!is.na(internet)), data.frame(betas = dfbetas(net_model)))

Hinweis: Indem wir betas= in den Befehl schreiben setzen wir einen Präfix, der vor unseren Variablen kommt, fest.

Anschliessend können wir alle Ausreisser entfernen, indem wir Merkmalsträger, die zu hohe oder zu tiefe dfbetas haben, ausschliessen. Wir erinnern uns: der Grenzwert beträgt 0.06. Das heisst wir wollen nur Merkmalsträger mit dfbetas zwischen -0.06 und 0.06 in den Robustheitstest einschliessen. Wir können nun nach dieser Bedingung die Merkmalsausprägungen herausfilter()n.

ess8_noOut <- filter(ess8_ch_mitdfb, betas.eduyrs > -0.06 & betas.eduyrs < 0.06)

 

7.4 Re-Analyse ohne Ausreisser

Mit unserem erstellten Datensatz ohne Ausreisser ess8_noOut ermitteln wir einen neuen Regressionskoeffizienten und vergleichen diesen mit dem Regressionskoeffizienten der Hauptanalyse.

noOut_model<- lm(internet ~ eduyrs, data = ess8_noOut)
summary(net_model)
## 
## Call:
## lm(formula = internet ~ eduyrs, data = ess8_ch_ss_2)
## 
## 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
summary(noOut_model)
## 
## Call:
## lm(formula = internet ~ eduyrs, data = ess8_noOut)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -147.11  -85.96  -41.09   47.06  937.76 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   97.627     14.065   6.941 6.62e-12 ***
## eduyrs         3.718      1.177   3.159  0.00163 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 125.9 on 1105 degrees of freedom
## Multiple R-squared:  0.008949,   Adjusted R-squared:  0.008052 
## F-statistic: 9.977 on 1 and 1105 DF,  p-value: 0.001628
sd(ess8_ch_ss_2$internet, na.rm = TRUE)
## [1] 163.1974

Zwischen den zwei Modelle unterscheidet sich der Koeffizient um knappe vier Einheiten, das Vorzeichen ändert sich allerdings nicht. Der Koeffizient verringert sich also um etwas mehr als die Hälfte. Aus inhalticher Perspektive ist das Regressionsergebnis somit nicht robust gegenüber einem Ausreisserausschluss.

Schauen wir uns die Robustheit zusätzlich aus statistischer Perspektive an: Der p-Wert indiziert in den Ausgangsergebnissen sowie auch unter Ausschluss der Ausreissern, dass es einen signifikanten Zusammenhang zwischen den Bildungsjahren und der täglichen Internetnutzung gibt. Anders als aus inhaltlicher Perspektive ist das Ergebnis aus statistischer Perspektive also robust gegenüber einem Ausreisserausschluss.

Konklusion: Wir würden in der Darstellung unserer Analyse den Hinweis ergänzen, dass die Ergebnisse aus inhaltlicher Perspektive nur bedingt robust gegenüber Ausreisserausschluss sind, aus statistischer Perspektive allerdings schon.

Zum Schluss können wir die beiden Modelle grafisch darstellen:

ggplot(ess8_ch_ss_2, aes(eduyrs, internet))+
  geom_jitter(alpha = 0.2, size = 2)+
  scale_x_continuous(breaks = seq(0, 30, 5))+
  geom_smooth(method = "lm", se = F, aes(color = "a"))+
  geom_smooth(method = "lm", se = F, aes(color = "b"),  data = ess8_noOut)+
  scale_color_manual(values = c("black", "orange"), 
                     labels = c("mit Ausreisser", "ohne Ausreisser"),
                     name = "")+
  labs(title = "Zeit im Internet nach Bildung", 
       caption = "ESS(2016), Teilstichprobe CH, N=1155", 
       y = "Zeit im Internet pro Tag in Minuten", 
       x = "Bildungsjahre")

 

Hier gehts weiter zur Übung III

 

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”