Übung zur Regression mit kategorialen UV

Dies ist das Dokument mit Lösungen

Klick hier, um die Musterlösungen auszublenden.

Zur Vorbereitung

ESS8-Datensatz einlesen:

# Working directory setzen (z.B. "C:/daten")
setwd("mein_laufwerk/mein_datenverzeichnis")
# Daten einlesen
library(haven)
ess8 <- read_dta("ESS8e02_2.dta")

Ein neues Übungsskript erstellen:

# Statistik 2: R Tutorat
# Übungsskript zur Regression mit kategorialen Variablen
# Datum: XXX
# AutorIn: XXX

Aktiviere die Packages aus dem tidyverse, sowie zusätzlich table1, stargazer, summarytools, labelled und visreg:

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

I. Datenaufbereitung & Inspektion

1. Wir analysieren den Zusammenhang zwischen dem Zivilstand und der Lebenszufriedenheit. Gibt es einen kausalen Einfluss des Zivilstandes auf die Lebenszufriedenheit? Formuliere eine Hypothese dazu und begründe diese in 2 Sätzen.

Der Ehestand hat einen positiven Einfluss auf die Lebenszufriedenheit (H1)

Mögliche Begründungen:

Die enge Beziehung zu einem Partner oder einer Partnerin kann als Ausgleich zu einem stressigen oder isolierten Alltag dienen und so die Lebenszufriedenheit erhöhen.
Die Ehe bietet ein gewisses Mass an finanzieller Stabilität und damit Schutz vor wirtschaftlichen Risikofaktoren, was wiederum die Lebenszufriedenheit erhöhen kann.
Der Ehestand ist nach wie vor als soziale Norm breit etabliert und als Leitbild der Lebensführung positiv konnotiert, mithin zufriedenheitsstiftend.

2. Weshalb ist es sinnvoll, das Alter im Kontext des empirischen Tests der Hypothese als Störmerkmal zu betrachten?

Zum einen kann davon ausgegangen werden, dass das Alter einen Einfluss auf die Lebenszufriedenheit hat. Zum anderen ist plausibel, dass das Alter den Rahmen der individuellen Lebensgestaltung beeinflusst und damit auch die Wahrscheinlichkeit zur Realisierung bestimmter Zivilstände altersabhängig ist.

3. Datenmanagement: Erstelle einen neuen Teildatensatz des ESS, in dem nur die Variablen idno, maritalb, stflife, agea enthalten sind. Wie gross ist dieser neue Datensatz?

#####################
## Datenmanagement ##
#####################

# Variablenauswahl
ess8_samp <- select(ess8, idno, maritalb, stflife, agea)

# Stichprobenumfang: 44387 Personen (mit Messungen in 4 Variablen)
dim(ess8_samp)

## [1] 44387     4

4. Inspiziere die ausgewählten Variablen. Welcher Variablenklasse gehört maritalb an? Und: Sind Rekodierungen dieser Variable notwendig oder sinnvoll?

# Inspektion: Variablentypen
look_for(ess8_samp)

##  pos variable label                      col_type missing
##  1   idno     Respondent's identificati~ dbl      0      
##  2   maritalb Legal marital status, pos~ dbl+lbl  878    
##                                                          
##                                                          
##                                                          
##                                                          
##                                                          
##                                                          
##                                                          
##                                                          
##  3   stflife  How satisfied with life a~ dbl+lbl  187    
##                                                          
##                                                          
##                                                          
##                                                          
##                                                          
##                                                          
##                                                          
##                                                          
##                                                          
##                                                          
##                                                          
##                                                          
##                                                          
##  4   agea     Age of respondent, calcul~ dbl+lbl  155    
##  values                    
##                            
##  [1] Legally married       
##  [2] In a legally register~
##  [3] Legally separated     
##  [4] Legally divorced/Civi~
##  [5] Widowed/Civil partner~
##  [6] None of these (NEVER ~
##  [NA(b)] Refusal           
##  [NA(c)] Don't know        
##  [NA(d)] No answer         
##  [0] Extremely dissatisfied
##  [1] 1                     
##  [2] 2                     
##  [3] 3                     
##  [4] 4                     
##  [5] 5                     
##  [6] 6                     
##  [7] 7                     
##  [8] 8                     
##  [9] 9                     
##  [10] Extremely satisfied  
##  [NA(b)] Refusal           
##  [NA(c)] Don't know        
##  [NA(d)] No answer         
##  [NA(d)] Not available

freq(ess8_samp$maritalb)

## Frequencies  
## 
##                Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
## ----------- ------- --------- -------------- --------- --------------
##           1   21711     49.90          49.90     48.91          48.91
##           2     443      1.02          50.92      1.00          49.91
##           3     648      1.49          52.41      1.46          51.37
##           4    3912      8.99          61.40      8.81          60.18
##           5    3756      8.63          70.03      8.46          68.65
##           6   13039     29.97         100.00     29.38          98.02
##        <NA>     878                               1.98         100.00
##       Total   44387    100.00         100.00    100.00         100.00

Die kategoriale Variable zum Zivilstand ist als numerische Variable angelegt. Alle fehlenden Werte sind als NA repräsentiert. Die Variable enthält 6 substanzielle Kategorien, die mit Blick auf unser Erkenntnisinteresse rekodiert werden sollten. Unsere Vorschläge dazu, mit Begründungen (allerdings nicht in Stein gemeisselt und durchaus diskussionswürdig):

Kategorie 2 bezieht sich auf zertifizierte Gemeinschaften gleichgeschlechtlicher Paare. Auch weil unser Erklärungsmodell solche Partnerschaften mit einschliesst, stellen wir sie mit der Ehe gleich, binden sie also mit der Kategorie 1 zusammen.
Kategorien 3 und 4 beziehen sich auf unterschiedliche Zeitpunkte im Prozess der formalen Auflösung einer Ehe. Wir fassen sie daher zusammen.
Verwitwete Personen stellen einen Sonderfall mit Blick auf die Lebenszufriedenheit dar, wir schliessen sie daher aus der Analyse aus und rekodieren die Werte der Kategorie 5 folglich auf Missing.

5. Rekodierungen kategorialer Variablen gehen meistens leichter von der Hand, wenn sie vor der Faktorisierung durchgeführt werden. Fasse daher nun die Ausprägungen 1 und 2 sowie die Ausprägungen 3 und 4 der Variable “maritalb” in jeweils einer Kategorie zusammen und definieren die Ausprägung 5 (´Widowed´) als Missing.

Tipp: Hier - wie bei nahezu jeder Rekodierung - ist es hilfreich…

die Rekodierung in einer neuen Variable vorzunehmen, damit Ziel- und Bedingungsvariable nicht identisch sind.
direkt Kategoriennamen statt numerischer Codes zu vergeben, zumindest wenn es sich um kategoriale Variablen handelt.
aufgrund der Unwägbarkeiten von R nach jedem Rekodierungsschritt die Datenmatrix zu inspizieren.
nach der Recodierung die neue Variable umfassend zu prüfen.

# Rekodierung / Klassierung der UV "martialb" zu neuer Variable "marital_class"
ess8_samp$marital_class[ess8_samp$maritalb == 1 | ess8_samp$maritalb == 2] <- "Verheiratet/Verpartnert" 
ess8_samp$marital_class[ess8_samp$maritalb == 3 | ess8_samp$maritalb == 4] <- "Getrennt"
ess8_samp$marital_class[ess8_samp$maritalb == 6 ] <- "Single" 
ess8_samp$marital_class[ess8_samp$maritalb == 5 ] <- NA 

# Kontrolle (Ausgangs vs. Zielvariable)
freq (ess8_samp$maritalb)

## Frequencies  
## 
##                Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
## ----------- ------- --------- -------------- --------- --------------
##           1   21711     49.90          49.90     48.91          48.91
##           2     443      1.02          50.92      1.00          49.91
##           3     648      1.49          52.41      1.46          51.37
##           4    3912      8.99          61.40      8.81          60.18
##           5    3756      8.63          70.03      8.46          68.65
##           6   13039     29.97         100.00     29.38          98.02
##        <NA>     878                               1.98         100.00
##       Total   44387    100.00         100.00    100.00         100.00

freq (ess8_samp$marital_class)

## Frequencies  
## 
##                                  Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
## ----------------------------- ------- --------- -------------- --------- --------------
##                      Getrennt    4560     11.47          11.47     10.27          10.27
##                        Single   13039     32.80          44.27     29.38          39.65
##       Verheiratet/Verpartnert   22154     55.73         100.00     49.91          89.56
##                          <NA>    4634                              10.44         100.00
##                         Total   44387    100.00         100.00    100.00         100.00

table (ess8_samp$marital_class, ess8_samp$maritalb, useNA = c("always"))

##                          
##                               1     2     3     4     5     6  <NA>
##   Getrennt                    0     0   648  3912     0     0     0
##   Single                      0     0     0     0     0 13039     0
##   Verheiratet/Verpartnert 21711   443     0     0     0     0     0
##   <NA>                        0     0     0     0  3756     0   878

class (ess8_samp$marital_class)

## [1] "character"

6. Bereite die 3 Variablen analysefähig auf. Achte darauf, dass…

kategoriale Variablen als Faktorvariable angelegt sind.
alle fehlende Werte korrekt als NA codiert sind.
Phantomkategorien eliminiert sind.

# Faktorisierung: "marital_class" als "character" angelegt, soll aber "factor" sein
ess8_samp$marital_class<-as_factor(ess8_samp$marital_class)

# Dateninspektion / Post-Faktorisierungs-Check
summary (ess8_samp)

##       idno              maritalb        stflife            agea       
##  Min.   :        1   Min.   :1.000   Min.   : 0.000   Min.   : 15.00  
##  1st Qu.:     1208   1st Qu.:1.000   1st Qu.: 6.000   1st Qu.: 34.00  
##  Median :     2589   Median :2.000   Median : 8.000   Median : 49.00  
##  Mean   : 31545782   Mean   :3.153   Mean   : 7.155   Mean   : 49.14  
##  3rd Qu.:    11058   3rd Qu.:6.000   3rd Qu.: 9.000   3rd Qu.: 64.00  
##  Max.   :551603139   Max.   :6.000   Max.   :10.000   Max.   :100.00  
##                      NA's   :878     NA's   :187      NA's   :155     
##                  marital_class  
##  Single                 :13039  
##  Verheiratet/Verpartnert:22154  
##  Getrennt               : 4560  
##  NA's                   : 4634  
##                                 
##                                 
##

freq (ess8_samp$marital_class)

## Frequencies  
## 
##                                  Freq   % Valid   % Valid Cum.   % Total   % Total Cum.
## ----------------------------- ------- --------- -------------- --------- --------------
##                        Single   13039     32.80          32.80     29.38          29.38
##       Verheiratet/Verpartnert   22154     55.73          88.53     49.91          79.29
##                      Getrennt    4560     11.47         100.00     10.27          89.56
##                          <NA>    4634                              10.44         100.00
##                         Total   44387    100.00         100.00    100.00         100.00

Fehlende Werte sind alle korrekt als NA kodiert, auch Phantomkategorien sind in dieser Anwendung keine vorhanden: Weder eine Postfaktorisierung noch die Entfernung von Phantomkategorien sind also vor dem Einstieg in die Analyse notwendig.

7. Nutze table1() um eine Tabelle mit der Stichprobenstatistik zu erstellen. Harmonisiere anschliessend die Stichprobe durch Ausschluss aller Beobachtungen mit fehlenden Werten in einer der Analysevariablen.

# Stichprobenstatistik
table1(~ stflife + marital_class + agea, data = ess8_samp)

	Overall (N=44387)
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%)
marital_class
Single	13039 (29.4%)
Verheiratet/Verpartnert	22154 (49.9%)
Getrennt	4560 (10.3%)
Missing	4634 (10.4%)
Age of respondent, calculated
Mean (SD)	49.1 (18.6)
Median [Min, Max]	49.0 [15.0, 100]
Missing	155 (0.3%)

ess8_samp_harm <- na.omit(ess8_samp)

Achtung: Vor einer Veröffentlichung (Paper, Präsentation) müsste die Tabelle nun zwingend noch weiter bearbeitet werden (z.B. in MS Word). Dabei sollte insb. auf knappe & prägnante Variablennamen geachtet werden. Auch die Reihenfolge der Kategorien der Zivilstand-Variable könnte überdacht werden.

# Stichprobenharmonisierung
ess8_samp_harm <- na.omit(ess8_samp)

8. Führe eine bivariate Testregression mit Familienstand als UV und Lebenszufriedenheit als AV durch. Was messen die Koeffizienten? Welches ist die Referenzkategorie der UV? Ist diese Default-Referenzkategorie sinnvoll? Führe ggf. eine Referenzmodifikation durch.

##################
## Datenanalyse ##
##################

# Testregression
summary (lm(stflife ~ marital_class, data = ess8_samp_harm))

## 
## Call:
## lm(formula = stflife ~ marital_class, data = ess8_samp_harm)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -7.388 -1.206  0.612  1.612  3.452 
## 
## Coefficients:
##                                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                           7.20576    0.01779 404.956  < 2e-16 ***
## marital_classVerheiratet/Verpartnert  0.18224    0.02242   8.127 4.53e-16 ***
## marital_classGetrennt                -0.65788    0.03494 -18.829  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.024 on 39476 degrees of freedom
## Multiple R-squared:  0.01621,    Adjusted R-squared:  0.01616 
## F-statistic: 325.2 on 2 and 39476 DF,  p-value: < 2.2e-16

Die Koeffizienten messen, wie sich die Lebenszufriedenheits-Mittelwerte der jeweils ausgewiesenen Familienstandskategorien von dem der Referenzkategorie unterscheiden. Diese Referenzkategorie ist nicht explizit ausgewiesen, kann aber im Ausschlussverfahren leicht hergeleitet werden: Es sind die Sinlges. Verheiratete/Verpartnerte Personen haben also im Mittel eine um 0.18 Skalenpunkte höhere, Getrennte eine um .66 Skalenpunkte niedrigere Lebenszufriedenheit als Singles.

Die hier automatisch von R gewählte Referenzkategorie ist nicht auf die Formulierung der Hypothese abgestimmt, in der explizit Unterschiede zur Gruppe der Verheirateten adressiert werden. Auf Auswertungsebene äussert sich diese Unstimmigkeit hier konkret in der Unmöglichkeit, auf Basis des dargestellten Ergebnisses den Mittelwertunterschied zwischen Verheirateten/Verpartnerten und Getrennten statistisch zu beurteilen. Daher ist in diesem Fall eine Referenzmodifikation sinnvoll.

# Referenzmodifikation entsprechend Hypothesenausrichtung
ess8_samp_harm$marital_class <- relevel(ess8_samp_harm$marital_class, ref = "Verheiratet/Verpartnert")

II. Regressionsanalyse

1. Wir lassen Management, Inspektion und Testanalyse hinter uns und starten in die Regressionsanalyse. Berechne zunächst eine einfache bivariate Regression. Speicher das Regressionsergebnis als Objekt ab.

# Bivariate Regressionsanalyse (unbereinigte Unterschiede)
bi_model <- lm(stflife ~ marital_class, data = ess8_samp_harm)
bi_model

## 
## Call:
## lm(formula = stflife ~ marital_class, data = ess8_samp_harm)
## 
## Coefficients:
##           (Intercept)    marital_classSingle  marital_classGetrennt  
##                7.3880                -0.1822                -0.8401

2. Integriere nun das Kontrollmerkmal Alter in die Analyse, Speichere auch dieses Ergebnis als Objekt.

# Multiple Regressionsanalyse (altersbereinigte Unterschiede)
tri_model <- lm(stflife ~ marital_class + agea, data = ess8_samp_harm)
tri_model

## 
## Call:
## lm(formula = stflife ~ marital_class + agea, data = ess8_samp_harm)
## 
## Coefficients:
##           (Intercept)    marital_classSingle  marital_classGetrennt  
##               7.90844               -0.38515               -0.83170  
##                  agea  
##              -0.00969

3. Mache dir die Funktionalität des stargazer()-Befehls zu Nutze und setze eine publikationswürdige Tabelle auf, welche beide Regressionsergebnisse darstellt.

# Regressionstabelle
stargazer(bi_model, tri_model, type = "text", 
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("Bivariates Modell", "Multivariates Modell"),
          omit.stat = c("f", "ser", "adj.rsq"), 
          digits = 2, digits.extra = 5,
          star.cutoffs = c(.05, .01, .001),
          notes = "Daten: ESS(2016), N = 39479, Standardfehler in Klammern", 
          title = "Determinanten der Lebenszufriedenheit",
          out = "bi_tri_model.doc")

## 
## Determinanten der Lebenszufriedenheit
## ==============================================================================
##                                                                               
##                            Bivariates Modell          Multivariates Modell    
##                                   (1)                          (2)            
## ------------------------------------------------------------------------------
## marital_classSingle             -0.18***                    -0.39***          
##                                  (0.02)                      (0.03)           
##                                                                               
## marital_classGetrennt           -0.84***                    -0.83***          
##                                  (0.03)                      (0.03)           
##                                                                               
## agea                                                        -0.01***          
##                                                              (0.001)          
##                                                                               
## Constant                        7.39***                      7.91***          
##                                  (0.01)                      (0.04)           
##                                                                               
## ------------------------------------------------------------------------------
## Observations                     39,479                      39,479           
## R2                                0.02                        0.02            
## ==============================================================================
## Note:                                            *p<0.05; **p<0.01; ***p<0.001
##                        Daten: ESS(2016), N = 39479, Standardfehler in Klammern

Achtung: Vor einer Veröffentlichung (Paper, Präsentation) müsste die Tabelle nun zwingend noch weiter bearbeitet werden (z.B. in MS Word). Dabei sollte insb. auf knappe & prägnante Variablennamen, einheitliche Sprache und Indikation der Referenzkategorie(n) geachtet werden.

4. Welche Schlüsse lassen die Ergebnisse der Regressionsanalyse auf die Hypothese (H1) zu?

Insbesondere unter Kontrolle des Alters weisen verheiratete bzw. verpartnerte Personen im Mittel eine deutlich höhere Lebenszufriedenheit auf als Singles und getrennte Personen: Der bereinigte Mittelwertunterschied umfasst mit 0,83 bzw. 0.39 Skalenpunkte immerhin 40% bzw. 20% der Standardabweichung der Lebenszufriedenheit (siehe auch Abbildung weiter unten). Beide Mittelwertunterschied sind zudem statistisch signifikant von 0 verschieden. Die oben forumlierte Hypothese wird durch die empirischen Ergebnisse also gestützt.

5. Interpretiere die Konstanten der beiden Analysen.

Im Rahmen des bivariaten Basismodells repräsentiert der Achsenabschnitt den Mittelwert der Gruppe der Verheirateten und Verpartnerten. Im multivariaten Modell misst die Konstante den vorhergesagten Wert für eine nulljährige Person, die verheiratet oder verpartnert ist.

6. Bei Konstanthaltung des Einkommens verdoppelt sich der Mittelwertunterschied zwischen Verheirateten und Singles. Dies zeigt der Vergleich der Koeffizienten von Modell 1 und Modell 2. Deute diese Veränderung.

Entsprechend eines typischen demographischen Lebensverlaufs sind Singles eher jüngere Personen. Jüngere Personen zeichnen sich, ausweislich der Ergebnisse des zweiten Modells, im Mittel durch eine relativ hohe Lebenszufriedenheit aus. Der dämpfende Einfluss des Single-Daseins auf die Lebenszufriedenheit wird folglich im ersten, bivariaten Modell durch den positiven Alterseffekt der (überwiegend jungen) Singles kaschiert bzw. unterschätzt. Erst, wenn der kaschierende Alterseffekt durch Alterskontrolle in Modell 2 abgespalten wird, zeigt sich der tatsächliche, ungeschönte negative Einfluss des Single-Daseins auf die Lebenszufriedenheit.

7. Stelle mit visreg() die bereinigten Mittelwertunterschiede graphisch dar

triplot <- visreg(tri_model, xvar = "marital_class", 
                    gg = T, partial = F, rug = F) +
  scale_y_continuous(breaks = seq(from = 6, to = 8, by = 0.5)) +
  coord_cartesian(ylim = c(6, 8)) + 
  labs(x = "Zivilstand", y = "",
       title = "Lebenszufriedenheit nach Zivilstand: Altersbereinigte Unterschiede",
       subtitle = "Lebenszufriedenheit (0 = unzufrieden / 10 = zufrieden)",
       caption = "Daten: ESS(2016), N = 39.479 \n Vorgergesagte Werte für Personen mit mittlerem Alter, \n Achsenauschnitt entspricht etwas einer Standardabweichung" ) +
  theme_bw()
triplot

Um Fehleinschätzungen zur Effektgrösse vorzubeugen wurde durch Wahl des Achsenausschnittes eine Standardabweichung der AV implizit als Vergleichsgrösse angelegt.

jpeg("triplot.jpeg", width = 14, height = 10, units ="cm", res = 1200)
triplot
dev.off()

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”

Übung zur Regression mit kategorialen UV

Dies ist das Dokument mit Lösungen

Zum European Social Survey

Zur Vorbereitung

I. Datenaufbereitung & Inspektion

II. Regressionsanalyse