Regression mit kategorialen UV

link zum pdf

# Statistik 2: R Tutorat
# Übungsskript zum Thema ANOVA
# Datum: 31.05.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)
ess8 <- read_dta("ESS8e02_2.dta")

#install.packages("dplyr")
library(dplyr)
#install.packages("table1")
library(table1)
#install.packages("tidyverse")
library(tidyverse)
#install.packages("stargazer")
library(stargazer)
#install.packages("visreg")
library(visreg)
#install.packages("ggplot2")
library(ggplot2)
#install.packages("summarytools")
library(summarytools)
#install.packages("labelled")
library(labelled)

In dieser Einheit lernst du, kategoriale unabhängige Variablen in die Regressionsanalyse einzubinden. Nach dem aus der Vorlesung bekannten Analyseschema ist für diese Merkmalskonstellation (UV kategorial, AV kontinuierlich) eigentlich der Mittelwertvergleich die passende Analysetechnik. Dieser lässt sich im Rahmen der Regressionsanalyse leicht nachbilden. Zentrale präparierende Datenoperation ist dabei die Faktorisierung der unabhängigen kategorialen Variablen.

Der ersten Beispielanalyse legen wir folgende ad-hoc formulierte Hypothese zugrunde: Männer weisen im Schnitt eine höhere Lebenszufriedenheit auf als Frauen (H1). Wir beginnen also mit einer dichotom ausgeprägten unabhängigen Variable (Geschlecht). Im zweiten Abschnitt widmen wir uns dann einer multikategorialen unabhängigen Variable (Erwerbsstatus). In diese zweite Analyse integrieren wir zudem noch weitere Kontrollmerkmale (Alter, Bildung).

Der Titel dieser Einheit spielt übrigens darauf an, dass die bivariate Analyse einer kategorialen unabhängigen Variablen im Rahmen eines linearen Modells oft auch “Varianzanalyse” (oder “Analysis of Variance” - kurz: “ANOVA”) genannt wird.

1. Vorbereitungen

Wir starten mit der Variablenselektion: wie immer beschränken wir den Datensatz auf die für uns relevanten Variablen. Hier sind dies stflife und gndr für die erste Analyse, sowie für die zweite Analyse zusätzlich agea, eduyrs und emplrel.

ess8 <- select(ess8, idno, stflife, gndr, agea, eduyrs, emplrel)

Mit look_for(ess8) oder attributes() können wir uns schnell die Variableneigenschaften erschliessen:

attributes(ess8$gndr)

## $label
## [1] "Gender"
## 
## $format.stata
## [1] "%12.0g"
## 
## $class
## [1] "haven_labelled" "vctrs_vctr"     "double"        
## 
## $labels
##      Male    Female No answer 
##         1         2        NA

attributes(ess8$emplrel)

## $label
## [1] "Employment relation"
## 
## $format.stata
## [1] "%31.0g"
## 
## $class
## [1] "haven_labelled" "vctrs_vctr"     "double"        
## 
## $labels
##                        Employee                   Self-employed 
##                               1                               2 
## Working for own family business                  Not applicable 
##                               3                              NA 
##                         Refusal                      Don't know 
##                              NA                              NA 
##                       No answer 
##                              NA

Die Variablen gndr und emplrel sind kategorial, während die anderen drei Variablen (quasi-)metrisches Skalenniveau aufweisen.
Trotz des kategorialen Niveaus sind gndr und emplrel im double-Format und somit als numerische Variablen angelegt.
In der Variable gndr sind dabei die Ausprägungen Männlich / Weiblich mit den Codes 1 und 2 besetzt.
In der Variable emplrel sind die Ausprägungen Employee / Self-employed / Working for own family business mit den Codes 1, 2 und 3 besetzt.

Rekodierungen werden zweckmässig bereits vor der Faktorisierung der kategorialen Variablen durchgeführt. Sie sind danach zwar auch noch möglich, dann aber meist etwas komplizierter. Allerdings sind in dieser Anwendung keine Rekodierungen notwendig: Fehlende Werte der Analysevariablen sind stimmig mit NA belegt, und Umpolungen oder Klassifizierungen der Variablen sind nicht notwendig bzw. nicht sinnvoll.

Die Faktorisierung der kategorialen Variable lässt sich ebenfalls an verschiedenen Positionen in die Aufbereitungsroutine eintakten, jeweils verbunden mit spezifischen Vor- und Nachteilen. Wir faktorisieren hier direkt nach dem grundlegenden Variablenmanagement:

ess8$emplrel <- as_factor(ess8$emplrel)
ess8$gndr <- as_factor(ess8$gndr)

Nun zur Variableninspektion:

summary(ess8)

##       idno              stflife              gndr            agea       
##  Min.   :        1   Min.   : 0.000   Male     :21027   Min.   : 15.00  
##  1st Qu.:     1208   1st Qu.: 6.000   Female   :23351   1st Qu.: 34.00  
##  Median :     2589   Median : 8.000   No answer:    9   Median : 49.00  
##  Mean   : 31545782   Mean   : 7.155                     Mean   : 49.14  
##  3rd Qu.:    11058   3rd Qu.: 9.000                     3rd Qu.: 64.00  
##  Max.   :551603139   Max.   :10.000                     Max.   :100.00  
##                      NA's   :187                        NA's   :155     
##      eduyrs                                 emplrel     
##  Min.   : 0.00   Employee                       :35074  
##  1st Qu.:11.00   Self-employed                  : 4842  
##  Median :13.00   Working for own family business:  743  
##  Mean   :13.04   Not applicable                 : 3527  
##  3rd Qu.:16.00   Refusal                        :  113  
##  Max.   :54.00   Don't know                     :   71  
##  NA's   :424     No answer                      :   17

Alle Variablen scheinen korrekt eingelesen und haben plausible Verteilungen. Allerdings sind, als Kollateralschaden der Faktorisierung, manche fehlende Werte in den Variablen gndr und emplrel nicht mehr als NA dargestellt, sondern bilden substanzielle Kategorien (“Refusal”, “No Answer”…). Auch wenn die entsprechenden Kategorienbezeichnungen auf fehlende Werte verweisen, werden sie von R nicht als solche erkannt und folglich auch nicht entsprechend verarbeitet. Diese Kategorien sollten daher in “echte” NAs rekodiert werden. Dieses kann über die einfache Standardsyntax zur Variablenrekodierung (ess8$gndr[ess8$gndr=="No answer"] <- NA) erfolgen, zur Rekodierung in NAs bietet sich jedoch zudem das spezifische Kommando na_if an:

ess8$gndr <- na_if(ess8$gndr, "No answer")
ess8$emplrel <- na_if(ess8$emplrel, "Refusal")
ess8$emplrel <- na_if(ess8$emplrel, "Don't know")
ess8$emplrel <- na_if(ess8$emplrel, "No answer")
ess8$emplrel <- na_if(ess8$emplrel, "Not applicable")
summary (ess8)

##       idno              stflife              gndr            agea       
##  Min.   :        1   Min.   : 0.000   Male     :21027   Min.   : 15.00  
##  1st Qu.:     1208   1st Qu.: 6.000   Female   :23351   1st Qu.: 34.00  
##  Median :     2589   Median : 8.000   No answer:    0   Median : 49.00  
##  Mean   : 31545782   Mean   : 7.155   NA's     :    9   Mean   : 49.14  
##  3rd Qu.:    11058   3rd Qu.: 9.000                     3rd Qu.: 64.00  
##  Max.   :551603139   Max.   :10.000                     Max.   :100.00  
##                      NA's   :187                        NA's   :155     
##      eduyrs                                 emplrel     
##  Min.   : 0.00   Employee                       :35074  
##  1st Qu.:11.00   Self-employed                  : 4842  
##  Median :13.00   Working for own family business:  743  
##  Mean   :13.04   Not applicable                 :    0  
##  3rd Qu.:16.00   Refusal                        :    0  
##  Max.   :54.00   (Other)                        :    0  
##  NA's   :424     NA's                           : 3728

Nun passt es: Im Zuge unsere Post-Faktorisierungs Bereinigung sind die “falschen” Kategorien für fehlende Werte zwar noch in den Variablen angelegt, deren Belegungen aber sinnvollerweise nun in die “richtige” NA-Kategorie verschoben.

Nun erstellen wir mit table1() eine gut formatierte Stichprobenstatistik. Diese geht über die einfache Darstellung von summary() hinaus und ist auf Publizierbarkeit ausgelegt. Sie gehört zu jeder publizierten Regressionsanalyse, ganz unabhängig vom modellierten Variablentyp. Super-Praktisch: Das table1()-Kommando verarbeitet faktorisierte Variablen optimal, weisst also für diese automatisch Häufigkeiten (statt Mittelwertsstatistiken) aus. Achtet immer darauf, dass ihr für die Übersichtsstatistik den Teildatensatz mit der vollständigen Stichprobe (also vor der Fallzahlharmonisierung) verwendet!

table1::table1(~ gndr + emplrel + stflife + agea + eduyrs, data = ess8)

	Overall (N=44387)
Gender
Male	21027 (47.4%)
Female	23351 (52.6%)
No answer	0 (0%)
Missing	9 (0.0%)
Employment relation
Employee	35074 (79.0%)
Self-employed	4842 (10.9%)
Working for own family business	743 (1.7%)
Not applicable	0 (0%)
Refusal	0 (0%)
Don't know	0 (0%)
No answer	0 (0%)
Missing	3728 (8.4%)
How satisfied with life as a whole
Mean (SD)	7.15 (2.09)
Median [Min, Max]	8.00 [0, 10.0]
Missing	187 (0.4%)
Age of respondent, calculated
Mean (SD)	49.1 (18.6)
Median [Min, Max]	49.0 [15.0, 100]
Missing	155 (0.3%)
Years of full-time education completed
Mean (SD)	13.0 (3.85)
Median [Min, Max]	13.0 [0, 54.0]
Missing	424 (1.0%)

Etwas lästig: Die vormaligen Kategorien für fehlende Werte sind als Kategorienrelikte weiterhin in den kategorialen Variablen angelegt, auch wenn wir oben die entsprechenden Werte in NA rekodiert haben. Mit fct_drop() können wir diese Phantomkategorien entfernen:

ess8$gndr <- fct_drop(ess8$gndr)
ess8$emplrel <- fct_drop(ess8$emplrel)

table1::table1(~ gndr + emplrel + stflife + agea + eduyrs, data = ess8)

	Overall (N=44387)
Gender
Male	21027 (47.4%)
Female	23351 (52.6%)
Missing	9 (0.0%)
Employment relation
Employee	35074 (79.0%)
Self-employed	4842 (10.9%)
Working for own family business	743 (1.7%)
Missing	3728 (8.4%)
How satisfied with life as a whole
Mean (SD)	7.15 (2.09)
Median [Min, Max]	8.00 [0, 10.0]
Missing	187 (0.4%)
Age of respondent, calculated
Mean (SD)	49.1 (18.6)
Median [Min, Max]	49.0 [15.0, 100]
Missing	155 (0.3%)
Years of full-time education completed
Mean (SD)	13.0 (3.85)
Median [Min, Max]	13.0 [0, 54.0]
Missing	424 (1.0%)

Obwohl wir zwei separate Analysen durchführen, nehmen wir - ausnahmsweise und formal nicht ganz korrekt - die Stichprobenharmonisierung von Einheiten mit fehlenden Werten auf Basis aller noch vorhandenen Variablen (also analyseübergreifend) vor. Durch diesen überzogenen Fallzahlausschluss kommt es somit zu leichten Ergebnisverzerrungen.

ess8_noNA <- na.omit(ess8)

Nun haben wir alle Vorbereitungsschritte zur Durchführung und Dokumentation der Regressionsanalyse abgeschlossen.

2. Regressionsanalyse mit dichotomer UV

Ein Störmerkmal lässt sich für das Geschlecht nicht sinnvoll herleiten, daher schätzen wir lediglich ein einfaches bivariates Regressionsmodell:

bi_model <- lm(stflife ~ gndr, data = ess8_noNA)
summary(bi_model)

## 
## Call:
## lm(formula = stflife ~ gndr, data = ess8_noNA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.1758 -1.1758  0.8242  1.8242  2.8716 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  7.17577    0.01500 478.282   <2e-16 ***
## gndrFemale  -0.04736    0.02084  -2.272   0.0231 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.085 on 40076 degrees of freedom
## Multiple R-squared:  0.0001288,  Adjusted R-squared:  0.0001039 
## F-statistic: 5.164 on 1 and 40076 DF,  p-value: 0.02307

Der Koeffizient der faktorisierten Geschlechtervariable zeigt nun an, wie sich die Lebenszufriedenheit der Frauen im Mittel von der von Männern unterscheidet: Die durchschnittliche Lebenszufriedenheit von Frauen liegt etwa 0.05 Skalenpunkte unterhalb der von Männern. Dasselbe Ergebnis erzeugt auch ein einfacher Mittelwertvergleich, z.B. im Rahmen des t-test: t.test(stflife ~ gndr, data = ess8_noNA, var.equal = T).

Wegen p<0.05 wird die Hypothese eines Geschlechterunterschieds formal-statistisch gestützt. Allerdings zeigt der Abgleich mit der Standardabweichung der AV sd(ess8_noNA$stflife), dass der mittlere Zufriedenheitsunterschied zwischen den Geschlechtern gerade einmal gut 1/40 eines typischen Unterschiedes in der Zufriedenheitsskala ausmacht, also aus inhaltlicher Perspektive eher als marginal eingeordnet werden muss.

Wichtig: Die Kategorie, auf die sich der Unterschiedswert bezieht (hier: Männer), wird auch “Referenzkategorie” genannt. Per Default setzt R die Kategorie mit dem geringer codierten Wert in der unfaktorisierten Ausgangsvariable (oder dem niedrigeren Anfangsbuchstaben) als Referenzkategorie. Sprich, wären statt der Männer die Frauen in der unfaktorisierten Ausgangsvariable mit “1” (und Männer mit “2”) kodiert, wären Frauen die Referenzkategorie der Analyse. Der betragsgleiche Koeffizient würde dann ein positives Vorzeichen - und somit die identische Information zum geschlechtsspezifischen Mittelwertunterschied - tragen.

Für die tabellarische Darstellung des Regressionsfits nutzen wir wieder das schon bekannte stargazer-Package:

stargazer(bi_model, type = "text", 
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("Bivariates Modell"), 
          single.row = TRUE, 
          omit.stat = c("f", "ser", "adj.rsq"), 
          digits = 2, digits.extra = 5,
          star.cutoffs = c(.05, .01, .001),
          notes = "Daten: ESS(2016), Standardfehler in Klammern", 
          title = "Lineares Modell: Geschlecht und Lebenszufriedenheit",
          out = "bi_model.doc")

## 
## Lineares Modell: Geschlecht und Lebenszufriedenheit
## =========================================================
##                                                          
##                           Bivariates Modell              
## ---------------------------------------------------------
## gndrFemale                  -0.05* (0.02)                
## Constant                    7.18*** (0.02)               
## ---------------------------------------------------------
## Observations                    40,078                   
## R2                              0.0001                   
## =========================================================
## Note:                       *p<0.05; **p<0.01; ***p<0.001
##              Daten: ESS(2016), Standardfehler in Klammern

Wichtig: Achtet bei der nun noch notwendigen Weiterverarbeitung und Formatierung der Tabelle darauf, dass die Referenzkategorie in der Tabelle kenntlich wird (siehe unten, multivariates Beispiel).

3. Regressionsanalyse mit multikategorialer UV

Als nächstes analysieren wir den Einfluss des Arbeitsverhältnisses auf die Lebenszufriedenheit. Ad-hoc formulieren wir die Hypothese, dass Arbeit im Familienbetrieb die zufriedenheitsstiftendste aller Erwerbsformen ist (H2) .

Die Inspektion und Faktorisierung haben wir bereits im Rahmen der analyseübergreifenden Vorbereitungen durchgeführt. Wir können uns trotzdem über look_for(ESS8_noNA) oder auch attributes(ESS8_noNA$emplrel) vergewissern, dass es sich bei emprel um eine kategoriale (bzw. nominal skalierte) Variable handelt. Außerdem nehmen wir in diesem Fall eine Referenzmodifikation vor: Wir verändern über relevel() die Referenzkategorie, welche in den folgenden Analysen, entsprechend der Hypothesenformulierung, durch die Ausprägung “Working for own family business” besetzt werden soll:

ess8_noNA$emplrel <- relevel(ess8_noNA$emplrel, ref = "Working for own family business")

Tipp: Grundsätzlich ist es immer sinnvoll, die Zuweisung der Referenzkategorie auf die Formulierung der Hypothese abzustimmen. Wenn, wie in unserem Fall, die Hypothese auf Unterschiede von Familienbetrieben zu allen anderen Erwerbsformen abstellt, bildet ein Regressionsmodell mit “Familienbetrieb” als Referenzkategorie konsistent diese Unterschiede empirisch ab. Falls die unabhängige Variable nicht hypothesenrelevant ist (sondern lediglich als Kontrollvariable integriert wird) ergibt es meistens Sinn, entweder die Modalkategorie als Referenzkategorie zu besetzen oder, bei ordinalen Variablen, diejenige mit dem niedrigsten Rang.

Wir starten in die Analyse mit einem bivariaten Regressionsmodell (Achtung: auch in dieser Anwendung kommt es durch den analyseübergreifenden Fallzahlausschluss wieder zu leichten Ergebnisverzerrungen):

bi_model <- lm(stflife ~ emplrel, data = ess8_noNA)
summary(bi_model)

## 
## Call:
## lm(formula = stflife ~ emplrel, data = ess8_noNA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.4672 -1.1160  0.6412  1.6412  2.8840 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           7.46721    0.07699  96.988  < 2e-16 ***
## emplrelEmployee      -0.35118    0.07780  -4.514 6.39e-06 ***
## emplrelSelf-employed -0.10842    0.08271  -1.311     0.19    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.083 on 40075 degrees of freedom
## Multiple R-squared:  0.001841,   Adjusted R-squared:  0.001791 
## F-statistic: 36.95 on 2 and 40075 DF,  p-value: < 2.2e-16

Das Regressionsergebnis weist wiederum die Mittelwertdifferenzen zur Referenzkategorie (hier: Working for own family business) aus. Tabellarisch aufbereitet mit Stargazer:

stargazer(bi_model, type = "text", 
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("Bivariates Modell"),
          single.row = TRUE, 
          omit.stat = c("f", "ser", "adj.rsq"), 
          digits = 2, digits.extra = 5,
          star.cutoffs = c(.05, .01, .001),
          notes = "Daten: ESS(2016), Standardfehler in Klammern", 
          title = "Der Effekt des Anstellungsverhältnis auf die Lebenszufriedenheit",
          out = "bi_model_2.doc")

## 
## Der Effekt des Anstellungsverhältnis auf die Lebenszufriedenheit
## =================================================================
##                                                                  
##                                   Bivariates Modell              
## -----------------------------------------------------------------
## emplrelEmployee                    -0.35*** (0.08)               
## emplrelSelf-employed                 -0.11 (0.08)                
## Constant                            7.47*** (0.08)               
## -----------------------------------------------------------------
## Observations                            40,078                   
## R2                                      0.002                    
## =================================================================
## Note:                               *p<0.05; **p<0.01; ***p<0.001
##                      Daten: ESS(2016), Standardfehler in Klammern

Entsprechend unserer Vermutung haben Angestellte und Selbstständige im Mittel eine geringere Lebenszufriedenheit als in Familienbetrieben tätige ArbeitnehmerInnen. Die durchschnittliche Lebenszufriedenheit von FamilienbetrieblerInnen liegt gut einen drittel Skalenpunkt über der von Angestellten und gut einen zehntel Skalenpunkt über der von Selbstständigen. Allerdings ist der zweite Unterschied nicht statistisch signifikant.

Auch dieses Ergebnis reproduziert letztlich einen einfachen Mittelwertvergleich der Lebenszufriedenheit zwischen den verschiedenen Berufsgruppen:

ess_gr<-group_by(ess8_noNA, emplrel)
summarise_at(ess_gr, vars(stflife), list(stflife=mean), na.rm=T)

## # A tibble: 3 × 2
##   emplrel                         stflife
##   <fct>                             <dbl>
## 1 Working for own family business    7.47
## 2 Employee                           7.12
## 3 Self-employed                      7.36

Warum aber setzen wir den Mittelwertvergleich im Regressionsrahmen um, obwohl dies komplizierter ist als das einfache Berechnen und Vergleichen kategorienspezifischer Mittelwerte? Der wichtigste Grund ist, dass wir im Rahmen einer multiplen Regressionsanalyse weitere Variablen integrieren und so insbesondere Störmerkmale kontrollieren können.

Mit Blick auf den vorliegenden Zusammenhang gehen wir davon aus, dass das Alter und die Bildung als Störmerkmale auf den Zusammenhang zwischen Anstellungsverhältnis und Lebenszufriedenheit einwirken können. Diese Störmerkmale integrieren wir daher als Drittvariablen (agea und eduyrs) in unsere Regressionsanalyse - und berechnen somit einen bereinigten Mittelwertunterschied der Lebenszufriedenheit zwischen FamilienbetrieblerInnen und Angestellten/Selbstständigen, der, anders als der unbereinigte, keinen Alters- oder Bildungseffekt transportiert.

multi_model <- lm(stflife ~ emplrel + agea + eduyrs, data = ess8_noNA)
summary (multi_model)

## 
## Call:
## lm(formula = stflife ~ emplrel + agea + eduyrs, data = ess8_noNA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.3707 -1.1161  0.4049  1.4084  3.9875 
## 
## Coefficients:
##                        Estimate Std. Error t value Pr(>|t|)    
## (Intercept)           6.7424970  0.0922882  73.059  < 2e-16 ***
## emplrelEmployee      -0.4212041  0.0770035  -5.470 4.53e-08 ***
## emplrelSelf-employed -0.1823491  0.0818965  -2.227    0.026 *  
## agea                 -0.0035493  0.0006072  -5.846 5.08e-09 ***
## eduyrs                0.0736300  0.0027860  26.428  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.061 on 40073 degrees of freedom
## Multiple R-squared:  0.0232, Adjusted R-squared:  0.0231 
## F-statistic: 237.9 on 4 and 40073 DF,  p-value: < 2.2e-16

stargazer(bi_model, multi_model, type = "text", 
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("Bivariates Modell", "Multivariates Modell"), 
          single.row = TRUE, 
          omit.stat = c("f", "ser", "adj.rsq"), 
          digits = 2, digits.extra = 5,
          star.cutoffs = c(.05, .01, .001),
          notes = "Daten: ESS(2016), Standardfehler in Klammern", 
          title = "Der Effekt vom Anstellungsverhältnis auf die Lebenszufriedenheit",
          out = "bi_mult_model_2.doc")

## 
## Der Effekt vom Anstellungsverhältnis auf die Lebenszufriedenheit
## ==================================================================
##                                                                   
##                        Bivariates Modell     Multivariates Modell 
##                               (1)                    (2)          
## ------------------------------------------------------------------
## emplrelEmployee         -0.35*** (0.08)        -0.42*** (0.08)    
## emplrelSelf-employed      -0.11 (0.08)          -0.18* (0.08)     
## agea                                          -0.004*** (0.001)   
## eduyrs                                         0.07*** (0.003)    
## Constant                 7.47*** (0.08)         6.74*** (0.09)    
## ------------------------------------------------------------------
## Observations                 40,078                 40,078        
## R2                           0.002                   0.02         
## ==================================================================
## Note:                                *p<0.05; **p<0.01; ***p<0.001
##                       Daten: ESS(2016), Standardfehler in Klammern

Die weitere Formatierung und Aufbereitung zur Publikationswürdigkeit führen wir extern z.B. in MS-Word durch:

Koeffizienten und Teststatistik der UV verändern sich bei deren Bereinigung durch die Kontrollvariablen. So nimmt die Grösse beider Koeffizienten im Vergleich zum Basismodell zu, die bereinigten Mittelwertunterschiede sind also grösser als die unbereinigten. Insbesondere der bereinigte LZ-Unterschied zwischen Selbständigen und FamilienbetrieblerInnen erscheint substanziell: Er beträgt immerhin fast einen halben Skalenpunkt und somit mehr als ein Fünftel der Standardabweichung der Lebenszufriedenheit.

Zudem unterschreiten beide p-Werte im multivariaten Modell die Signifikanzgrenze. Entsprechend sind die Koeffizienten im Vergleich zwischen den beiden Modellen zwar nicht “robust”, aber die formulierte Hypothese (deren Prüfung, sofern sie in kausaler Terminologie formuliert wurde, immer auf Grundlage des Nettomodells erfolgt), lässt sich auf Grundlage des multiplen Regressionsmodells nun statistisch umfassender stützen.

Hintergrund der Koeffizientenveränderung der Selbständigen vom bi- zum multivariaten Modell dürfte deren überdurchschnittlich hohe Bildung sein, die zu einem bloss geringen unbereinigten Zufriedenheitsdefizit gegenüber den Familienbetrieblern führt bzw. den “wahren” negativen Selstständigeneffekt im bivariaten Modell kaschiert.

4. Visualisierung multipler Regressionsergebnisse mit kategorialer UV

Als Stütze für die Interpretation lassen sich die Ergebnisse der beiden Analysen auch visuell aufbereiten. Beim der bivariaten Analyse zum Zusammenhang von Geschlecht und Lebenszufriedenheit könnten wir die typische ggplot2 Funktionalität zur Visualisierung von Mittelwertvergleichen (siehe Statistik 1) nutzen. Bei der zweiten, multivariaten Analyse können wir die bereinigten Koeffizienten der Erwerbsform aus dem multiplen Regressionsmodell wiederum mit dem visreg()-Befehl darstellen.

4.1 Visualisierung - Trichotome UV

multiplot_2 <- visreg(multi_model, xvar = "emplrel", 
                      gg = TRUE, partial = FALSE, rug = FALSE) + 
  scale_y_continuous(breaks = seq(from = 6.2, to = 8.2, by = 0.2)) +
  coord_cartesian(ylim = c(6.2, 8.2)) + 
  labs(x = "Anstellungsverhältnis", y = "Lebenszufriedenheit",
       title = "Lebenszufriedenheit nach Anstellungsverhältnis",
       subtitle = "Vorhergesagte Werte unter Konstanthaltung des Alters und der Bildung im Mittelwert",
       caption = "Daten: ESS(2016), N = 40078 \n Lebenszufriedenheit: 0 = gar nicht zufrieden / 10 = sehr zufrieden") + 
  theme_bw()
multiplot_2

In beiden grafischen Aufbereitungen zeigen die blauen Geraden die vorhergesagten Zufriedenheitswerte unter Konstanthaltung der Drittvariablen (per visreg()-Default in den jeweiligen Mittelwerten bzw. Medianen) an. Die Unterschiede zwischen den blauen Geraden bilden folglich die bereinigten Mittelwertunterschiede bzw. den geschätzten Effekt der Erwerbsform ab. Die grauen Flächen wiederum indizieren das Konfidenzband.

Makel der Abbildung sind noch die Wertelabel der unabhängigen Variablen. visreg() ist in diesem Punkt sehr störrisch und bietet leider keine Option zur befehlsinternen Umbenennung der Wertelabel an. Um eine akzeptable Notation der x-Achse zu erhalten, müssen wir daher die Wertelabel der Ausgangsvariable vor der Regression und deren Darstellung modifizieren:

ess8_noNA$emplrel <- factor(ess8_noNA$emplrel, labels = c("Working for own family business" = "Familienbetrieb", "Employee" = "Angestellt", "Self-employed" = "Selbstständig"))
multi_model <- lm(stflife ~ emplrel + agea + eduyrs, data = ess8_noNA)

multiplot_2 <- visreg(multi_model, xvar = "emplrel", 
                      gg = TRUE, partial = FALSE, rug = FALSE) + 
  scale_y_continuous(breaks = seq(from = 6.2, to = 8.2, by = 0.2)) +
  coord_cartesian(ylim = c(6.2, 8.2)) + 
  labs(x = "Anstellungsverhältnis", y = "Lebenszufriedenheit",
       title = "Lebenszufriedenheit nach Anstellungsverhältnis",
       subtitle = "Vorhergesagte Werte unter Konstanthaltung des Alters und der Bildung im Mittelwert",
       caption = "Daten: ESS(2016), N = 40078 \n Lebenszufriedenheit: 0 = gar nicht zufrieden / 10 = sehr zufrieden, \n Achsenauschnitt entspricht etwas einer Standardabweichung") + 
  theme_bw()
multiplot_2

jpeg("multiplot_2.jpeg", width = 42, height = 30, units ="cm", res = 1200)
multiplot_2
dev.off()

Alternative Speicher/Export-Optionen für Grafiken: siehe Einheit Drittvariablen.

4.2 Visualisierung - Kategoriale Variabel als Kontrollvariable

Ist die kategoriale Variable nicht zentrales unabhängiges Merkmal sondern Kontrollvariable, verschiebt sich der Auswertungs- und Darstellungsfokus des multivariaten Modells: von der kategorialen zu der zentralen unabhängigen (metrischen) Variable.

Wäre die vorliegende Analyse nicht auf den Erwerbsstatus, sondern die Bildungsjahre als zentrales unabhängiges Merkmal gerichtet, müssten wir zunächst die Reihenfolge der Darstellung in der Tabelle ändern - einer Konvention folgend steht stets die hypothesenrelevante UV an der ersten Position. Entsprechend würde sich auch der Interpretationsfokus verändern: Unter Konstanthaltung des Erwerbsstatus und des Alters steigt mit jedem zusätzlichen Bildungsjahr die Lebenszufriedenheit um 0,07 Skalenpunkte - wäre nun das wichtigste Auswertungsstatement.

Weiterhin gelten alle Konventionen zur Darstellung und Kommunikation von multivariaten Regressionsergebnissen (siehe Einheit Drittvariablen):

Im Rahmen der Vorhersagen mit predict muss nun zusätzlich eine sinnvoll ausgewählte Kategorie der kategorialen Variable konstant gehalten werden.

predict(multi_model, data.frame(eduyrs=c(8,12,16), emplrel="Angestellt", agea=22))

##        1        2        3 
## 6.832247 7.126767 7.421288

bei Vorhersagen mit emmeans werden für alle nicht im Befehl spezifizierten Variablen per Default arithmetische Mittel in die Regressionsgleichung eingesetzt. Dies gilt auch für die kategorialen Variablen: ihr Vorhersagebeitrag ergibt sich aus den Mittelwerten der einzelnen Dummy-Variablen.

library(emmeans)
emmeans(multi_model,~ eduyrs, at = list(eduyrs=c(8,12,16)))

##  eduyrs emmean     SE    df lower.CL upper.CL
##       8   6.95 0.0309 40073     6.89     7.01
##      12   7.25 0.0277 40073     7.19     7.30
##      16   7.54 0.0287 40073     7.48     7.60
## 
## Results are averaged over the levels of: emplrel 
## Confidence level used: 0.95

Die Default-Konstanthaltungslogik von visreg unterscheidet sich davon in zweifacher Hinsicht. Bei (a) metrischen unabhängigen Variablen wird der Querschnitt aus dem Regressionskörper per Default nicht beim arithmetischen Mittel, sondern dem Median gesetzt (siehe Einheit Drittvariablen). Bei (b) kategorialen unabhängigen Variablen beruht die Vorhersagefunktion von visreg nicht auf den Mittelwerten der abgeleiteten Dummy-Variablen, sondern der Modalkategorie.

visreg(multi_model, xvar = "eduyrs", 
                      gg = TRUE, partial = FALSE, rug = FALSE) + 
  scale_y_continuous(breaks = seq(from = 5, to = 9, by = 0.5)) +
  scale_x_continuous(breaks = seq(from = 6, to = 20, by = 2)) +
  coord_cartesian(ylim = c(6.2, 8.2), xlim=c(5,20)) + 
  labs(       title = "Lebenszufriedenheit nach Bildung",
       subtitle = "Vorhergesagte Werte unter Konstanthaltung des Alters und des Angestelltenverhältnisses",
       caption = "Daten: ESS(2016), N = 40078 \n Lebenszufriedenheit: 0 = gar nicht zufrieden / 10 = sehr zufrieden, \n Achsenauschnitt entspricht etwa einer Standardabweichung") + 
  theme_bw()

Werden Vorhersagen und Effektplot in einen Interpretationszusammenhang gerückt, ist es daher sinnvoll, den Default von emmeans zu überschreiben und die Vorhersagen ebenfalls auf Mediane und Modalwerte der Drittmerkmale auszurichten.

median (ess8_noNA$agea)

## [1] 50

table (ess8_noNA$emplrel)

## 
## Familienbetrieb      Angestellt   Selbstständig 
##             732           34594            4752

emmeans(multi_model,~ eduyrs, at = list(eduyrs=c(8,12,16), agea=50, emplrel="Angestellt"))

##  eduyrs emmean      SE    df lower.CL upper.CL
##       8  6.733 0.01822 40073    6.697    6.769
##      12  7.027 0.01157 40073    7.005    7.050
##      16  7.322 0.01356 40073    7.295    7.348
## 
## Confidence level used: 0.95

Folgende Texte empfehlen wir euch zur Lektüre und Vertiefung:

BCP Kapitel 6.2-6.3: Beckerman, A.P., Childs, D.Z., Petchey, O.L. (2017): Getting Started with R. Oxford: University Press.

SRM Kapitel II.8.1: Marco R. Steenbergen (2016): Regression Analysis: A Primer for the Social Sciences

Hier gehts weiter zur Übung VI

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”