Drittvariablen

link zum pdf

# Statistik 2: R Tutorat
# Übungsskript zum Thema Multivariate Regression und Drittvariablen
# Datum: 18.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("ggplot2")
library(ggplot2)
#install.packages("visreg")
library(visreg)
#install.packages("table1")
library(table1)
#install.packages("stargazer")
library(stargazer)
#install.packages("labelled")
library(labelled)
#install.packages("rgl")
library(rgl)
#install.packages("forcats")
library(forcats)

Theoretische Vorüberlegungen

In dieser Einheit wirst du lernen, wie du in R multiple Regressionsanalysen durchführen und die Ergebnisse anschaulich aufbereiten kannst. Hierfür werden wir uns beispielhaft mit dem Zusammenhang zwischen dem Einkommen und der Migrationswertschätzung befassen.

Wir vermuten, dass das Einkommen auf verschiedenen Kanälen Toleranz, Anerkennung und Offenheit befördert und formulieren dementsprechend die Hypothese, dass Einkommen hat einen positiven Einfluss auf die Migrationswertschätzung hat (H1)

Wir könnten nun wie bisher eine einfache bivariate Regression durchführen und bei positivem Einkommenskoeffizienten unsere H1 als empirisch gestützt betrachten. Dieser Schluss wäre allerdings voreilig: Der bivariate Zusammenhang könnte Einflüsse möglicher Störmerkmale aufgreifen und daher eine blosse “Scheinkausalität” darstellen. Insbesondere könnte der positive Zusammenhang auf die höhere Bildung von Besserverdienenden zurückgehen: Personen mit hohen Einkommen haben deswegen positivere Migrationseinstellungen, weil sie besser gebildet sind und hohe Bildung zu positiveren Migrationseinstellungen führt. Der positive Zusammenhang zwischen Einkommen und Migrationswertschätzung würde somit einen Bildungseffekt transportieren. Folglich würden wir den Effekt des Einkommens auf Basis des Einkommenskoeffizienten der bivariaten Analyse überschätzen.

Die Regressionsanalyse erlaubt es uns, solche Störmerkmale - sofern sie als Variablen im Datensatz vorhanden sind - in die Analyse zu integrieren und ihre Effekte somit zu kontrollieren.

Rein formal ist ein Störmerkmal (wie die Bildung in diesem Beispiel) dadurch definiert, dass es einen Einfluss sowohl auf die UV als auch die AV hat. So wäre z.B. das Alter ein weiterer möglicher Störfaktor, da dieses Merkmal gut begründbar sowohl das Einkommen als auch die Einstellung zur Migration beeinflussen könnte. Wir beschränken unsere Analyse aus didaktischen Gründen vorerst auf diese beiden zusätzlichen Kontrollmerkmale (*Bildung und Alter), fügen dann später noch zusätzlich das Geschlecht hinzu - ein kategoriales Merkmal zwar, welches aber trotzdem im Rahmen der Regressionsanalyse kontrolliert werden kann.

1. Datenmanagement

1.1 Selektion und Inspektion

Wir beschränken den Datensatz auf die für uns relevanten Beobachtungen und Variablen. Für uns sind das die Variablen imueclt, hinctnta, eduyrs,agea und gndr für alle Beobachtungen aus der Schweiz.

ess8_CH <- filter(ess8, cntry == "CH")
ess8_CH_ss <- select(ess8_CH, idno, imueclt, hinctnta, eduyrs, agea, gndr)

Detailinfos zu den Variablen und Labels gibt’s per look_for():

look_for(ess8_CH, "hinctnta")
look_for(ess8_CH, "imueclt")
look_for(ess8_CH, "eduyrs")
look_for(ess8_CH, "agea")
look_for(ess8_CH, "gndr")
attributes(ess8_CH$gndr)

Im ESS wird das Einkommen in Dezilen gemessen. Auf der Website des ESS finden wir unter “Data” und “Documentation” mehr Informationen zu bestimmten Variablen, z.B. ein Dokument namens ESS Appendix A2 Income ed. 2.1, welches Informationen zu länderspezifischen Dezilsgrenzen enthält.

Wichtiger analysevorbereitender Schritt ist die Inspektion der Verteilungen der analyserelevanten Variablen:

summary(ess8_CH_ss)

##       idno         imueclt          hinctnta          eduyrs    
##  Min.   :   1   Min.   : 0.000   Min.   : 1.000   Min.   : 0.0  
##  1st Qu.:1017   1st Qu.: 5.000   1st Qu.: 3.000   1st Qu.: 9.0  
##  Median :1827   Median : 6.000   Median : 5.000   Median :10.0  
##  Mean   :1807   Mean   : 6.072   Mean   : 5.472   Mean   :11.3  
##  3rd Qu.:2645   3rd Qu.: 8.000   3rd Qu.: 7.000   3rd Qu.:13.0  
##  Max.   :3460   Max.   :10.000   Max.   :10.000   Max.   :26.0  
##                 NA's   :14       NA's   :274      NA's   :3     
##       agea            gndr      
##  Min.   :15.00   Min.   :1.000  
##  1st Qu.:32.00   1st Qu.:1.000  
##  Median :48.00   Median :1.000  
##  Mean   :47.83   Mean   :1.483  
##  3rd Qu.:62.00   3rd Qu.:2.000  
##  Max.   :94.00   Max.   :2.000  
##  NA's   :6

Die Geschlechtervariable steckt noch im falschen Variablenformat, deutlich an den numerischen Kategorien im summary(). Wir faktorisieren sie folglich mit as_factor().

ess8_CH_ss$gndr <- as_factor(ess8_CH_ss$gndr)
summary(ess8_CH_ss)

##       idno         imueclt          hinctnta          eduyrs    
##  Min.   :   1   Min.   : 0.000   Min.   : 1.000   Min.   : 0.0  
##  1st Qu.:1017   1st Qu.: 5.000   1st Qu.: 3.000   1st Qu.: 9.0  
##  Median :1827   Median : 6.000   Median : 5.000   Median :10.0  
##  Mean   :1807   Mean   : 6.072   Mean   : 5.472   Mean   :11.3  
##  3rd Qu.:2645   3rd Qu.: 8.000   3rd Qu.: 7.000   3rd Qu.:13.0  
##  Max.   :3460   Max.   :10.000   Max.   :10.000   Max.   :26.0  
##                 NA's   :14       NA's   :274      NA's   :3     
##       agea              gndr    
##  Min.   :15.00   Male     :788  
##  1st Qu.:32.00   Female   :737  
##  Median :48.00   No answer:  0  
##  Mean   :47.83                  
##  3rd Qu.:62.00                  
##  Max.   :94.00                  
##  NA's   :6

1.2 Stichprobenstatistik

Vor jeder multiplen Regressionsanalyse stellen wir die univariaten Verteilungen der beteiligten Variablen tabellarisch dar. Mit dem Befehl table1() aus dem gleichnamigen Paket können wir leicht die Rohversion einer solchen Stichprobenstatistik erstellen. Diese geht über die einfachen univariaten Darstellungen aus summary() hinaus und ist auf Publizierbarkeit ausgelegt. Achtung! Unsere Stichprobenstatistiken beziehen sich immer auf den vollständigen Datensatz (also inklusive der Personen mit unvollständigen Angaben) - Es ist gerade auch Sinn und Zweck der Stichprobenstatistik, über Ausfallmuster in den analyserelevanten Variablen zu informieren.

#install.packages("table1")
library(table1)
table1(~imueclt + hinctnta + eduyrs + agea + gndr, data = ess8_CH_ss)

	Overall (N=1525)
Country's cultural life undermined or enriched by immigrants
Mean (SD)	6.07 (2.26)
Median [Min, Max]	6.00 [0, 10.0]
Missing	14 (0.9%)
Household's total net income, all sources
Mean (SD)	5.47 (2.57)
Median [Min, Max]	5.00 [1.00, 10.0]
Missing	274 (18.0%)
Years of full-time education completed
Mean (SD)	11.3 (3.50)
Median [Min, Max]	10.0 [0, 26.0]
Missing	3 (0.2%)
Age of respondent, calculated
Mean (SD)	47.8 (18.8)
Median [Min, Max]	48.0 [15.0, 94.0]
Missing	6 (0.4%)
Gender
Male	788 (51.7%)
Female	737 (48.3%)
No answer	0 (0%)

table1() stellt bei numerischen Variablen automatisch die wichtigsten Parameter der zusammenfassenden Statistik (Mittelwert, Standardabweichung etc.) dar, bei Faktorvariablen praktischerweise dagegen die Häufigkeitsverteilung.
NAs werden jeweils unter der Kategorie “Missing” aufgeführt. Bei der Geschlechtervariable taucht allerdings die Variablenkategorie “No Answer” auf. Im Rahmen von Faktorisierungen werden oftmals solche Phantomkategorien angelegt und dann in die verschiedene R-Anwendungen hineingespült - wie hier in den table1()-Output. Um dieses Darstellungsproblem frühzeitig zu umgehen entfernen wir die Phantomkategorien einer Variable einfach mit fct_drop() - einem Befehl stammt aus dem forcats Package, welches Teil des tidyverse ist:

ess8_CH_ss$gndr <- fct_drop(ess8_CH_ss$gndr)
table1(~imueclt + hinctnta + eduyrs + agea + gndr, data = ess8_CH_ss)

	Overall (N=1525)
Country's cultural life undermined or enriched by immigrants
Mean (SD)	6.07 (2.26)
Median [Min, Max]	6.00 [0, 10.0]
Missing	14 (0.9%)
Household's total net income, all sources
Mean (SD)	5.47 (2.57)
Median [Min, Max]	5.00 [1.00, 10.0]
Missing	274 (18.0%)
Years of full-time education completed
Mean (SD)	11.3 (3.50)
Median [Min, Max]	10.0 [0, 26.0]
Missing	3 (0.2%)
Age of respondent, calculated
Mean (SD)	47.8 (18.8)
Median [Min, Max]	48.0 [15.0, 94.0]
Missing	6 (0.4%)
Gender
Male	788 (51.7%)
Female	737 (48.3%)

Vor der Veröffentlichung müssen wir die Tabelle nun noch in ein Textverarbeitungsprogramm kopieren, um sie weiterzuverarbeiten und zu finalisieren (Überschrift, ggf. Sprache anpassen, Variablennamen modifizieren/vereinfachen etc.). Im Rahmen von Präsentationen sind Sprachdualismen Deutsch/Englisch meist ok, im Rahmen formeller textbasierter Publikationen dagegen eher nicht.

1.3 Stichprobenharmonisierung

Nur Einheiten mit vollständigen Messungen in den analyserelevanten Variablen gehen üblicherweise in die Regressionsanalyse ein. Diese Strategie, in R typischerweise umgesetzt mit na.omit(), folgt dem Leitbild des “paarweisen Ausschlusses” und führt zur Harmonisierung der Fallzahl über verschiedene Modellierungen einer Analyse.

ess8_noNA <- na.omit(ess8_CH_ss)

Durch Entfernung von Personen mit NAs in einer der analyserelevanten Variablen reduziert sich der Datensatz von 1525 auf 1242 Merkmalsträger. Achtung: Die Stichprobenharmonisierung unbedingt erst nach Erstellung der Stichprobenstatistik durchführen!

2. Regressionsanalyse

2.1 Bivariates Modell

Ausgangspunkt der Regressionsanalyse ist immer die einfache bivariate Regressionsanalyse:

bi_model <- lm(imueclt ~ hinctnta, data = ess8_noNA)
bi_model

## 
## Call:
## lm(formula = imueclt ~ hinctnta, data = ess8_noNA)
## 
## Coefficients:
## (Intercept)     hinctnta  
##     5.64297      0.09451

Der Koeffizient des Einkommens beträgt 0.095: mit jedem Dezilsprung im Haushaltseinkommen erhöht sich die Migrationswertschätzung im Schnitt um etwa einen zehntel Skalenpunkt. Unserer Argumentation oben folgend steht dieser Koeffizient gleichwohl im Verdacht, in Wirklichkeit keinen Einkommens-, sondern einen Bildungseffekt abzubilden. Daher integrieren wir im folgenden Analyseschritt eine Bildungsvariable zur Bereinigung des Einkommenskoeffizienten mit in die Analyse. Geometrische betrachtet erweitern wir damit unser Streudiagramm um eine zusätzliche Dimension (zu einem Streuraum) und ermitteln nun statt einer bestmöglichen Regressionsgerade die bestmögliche Regressionsebene.

2.3 Trivariates Modell

Wir integrieren die Bildungsvariable in die Analyse, indem wir sie mit + als zusätzliche unabhängige Variable in den Regressionsbefehl einbinden.

tri_model <- lm(imueclt ~ hinctnta + eduyrs, data = ess8_noNA)

Mit summary() bekommen wir wieder mehr Informationen zum Analyseergebnis:

summary(bi_model)

## 
## Call:
## lm(formula = imueclt ~ hinctnta, data = ess8_noNA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.3045 -1.3991  0.4119  1.6955  4.2625 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.64297    0.14781  38.177  < 2e-16 ***
## hinctnta     0.09451    0.02441   3.872 0.000114 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.213 on 1240 degrees of freedom
## Multiple R-squared:  0.01194,    Adjusted R-squared:  0.01115 
## F-statistic: 14.99 on 1 and 1240 DF,  p-value: 0.0001137

summary(tri_model)

## 
## Call:
## lm(formula = imueclt ~ hinctnta + eduyrs, data = ess8_noNA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7692 -1.5479  0.2602  1.4930  4.9479 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.93456    0.21889  17.975   <2e-16 ***
## hinctnta     0.02937    0.02429   1.209    0.227    
## eduyrs       0.18114    0.01767  10.254   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.125 on 1239 degrees of freedom
## Multiple R-squared:  0.08923,    Adjusted R-squared:  0.08776 
## F-statistic: 60.69 on 2 and 1239 DF,  p-value: < 2.2e-16

Allerdings ist der Vergleich beider Modelle so ziemlich mühsam. Ausserdem können wir die Ergebnisse in dieser Form nicht publizieren.

2.4 Regressionstabelle

Wir bereiten unsere Analyseergebnisse statt dessen im Rahmen einer integrierten Regressionstabelle auf, um sie anschaulich darzustellen. Wir erstellen diese Tabelle mithilfe des Befehls stargazer() aus dem gleichnamigen Package. Innerhalb des Befehls spezifizieren wir die Analyseergebnisse, die in unserer Tabelle dargestellt werden sollen. Damit die Tabelle in unserer Konsole erscheint, wählen wir zusätzlich type = “text”.

#install.packages("stargazer")
library(stargazer)
stargazer(bi_model, tri_model, 
          type = "text")

## 
## =====================================================================
##                                    Dependent variable:               
##                     -------------------------------------------------
##                                          imueclt                     
##                               (1)                      (2)           
## ---------------------------------------------------------------------
## hinctnta                    0.095***                  0.029          
##                             (0.024)                  (0.024)         
##                                                                      
## eduyrs                                               0.181***        
##                                                      (0.018)         
##                                                                      
## Constant                    5.643***                 3.935***        
##                             (0.148)                  (0.219)         
##                                                                      
## ---------------------------------------------------------------------
## Observations                 1,242                    1,242          
## R2                           0.012                    0.089          
## Adjusted R2                  0.011                    0.088          
## Residual Std. Error    2.213 (df = 1240)        2.125 (df = 1239)    
## F Statistic         14.990*** (df = 1; 1240) 60.694*** (df = 2; 1239)
## =====================================================================
## Note:                                     *p<0.1; **p<0.05; ***p<0.01

…sieht schon ganz gut aus, allerdings müssen wir noch einige Dinge anpassen, bevor wir die Tabelle veröffentlichen können. Dies tun wir, indem wir im stargazer()-Befehl die verschiedenen Unteroptionen spezifizieren - so wir es zuvor auch immer für die graphischen Darstellungen getan haben. Eine Übersicht dieser Optionen erhältst du, wie auch bei anderen Befehlen, über den Aufruf der R-interne Hilfe per ?stargazer. Achtung: Bedingung zur Nutzung dieser R-internen Hilfsoption ist, dass das befehltragende Package aktiviert wurde.

Spaltenüberschriften entfernen

Entfernen wir nun zuerst die Überschrift “Dependent Variable” aus der Tabelle. Dies machen wir über das Argument dep.var.caption = ““.

stargazer(bi_model, tri_model, 
          type = "text", 
          dep.var.caption = "")

## 
## =====================================================================
##                                          imueclt                     
##                               (1)                      (2)           
## ---------------------------------------------------------------------
## hinctnta                    0.095***                  0.029          
##                             (0.024)                  (0.024)         
##                                                                      
## eduyrs                                               0.181***        
##                                                      (0.018)         
##                                                                      
## Constant                    5.643***                 3.935***        
##                             (0.148)                  (0.219)         
##                                                                      
## ---------------------------------------------------------------------
## Observations                 1,242                    1,242          
## R2                           0.012                    0.089          
## Adjusted R2                  0.011                    0.088          
## Residual Std. Error    2.213 (df = 1240)        2.125 (df = 1239)    
## F Statistic         14.990*** (df = 1; 1240) 60.694*** (df = 2; 1239)
## =====================================================================
## Note:                                     *p<0.1; **p<0.05; ***p<0.01

In unserer Tabelle können wir den Namen der AV über das Argument dep.var.labels anpassen. Da wir jedoch nur eine AV haben, entfernen wir den Namen ganz aus der Tabelle.

stargazer(bi_model, tri_model, 
          type = "text", 
          dep.var.caption = "",
          dep.var.labels = "")

## 
## =====================================================================
##                                                                      
##                               (1)                      (2)           
## ---------------------------------------------------------------------
## hinctnta                    0.095***                  0.029          
##                             (0.024)                  (0.024)         
##                                                                      
## eduyrs                                               0.181***        
##                                                      (0.018)         
##                                                                      
## Constant                    5.643***                 3.935***        
##                             (0.148)                  (0.219)         
##                                                                      
## ---------------------------------------------------------------------
## Observations                 1,242                    1,242          
## R2                           0.012                    0.089          
## Adjusted R2                  0.011                    0.088          
## Residual Std. Error    2.213 (df = 1240)        2.125 (df = 1239)    
## F Statistic         14.990*** (df = 1; 1240) 60.694*** (df = 2; 1239)
## =====================================================================
## Note:                                     *p<0.1; **p<0.05; ***p<0.01

Spaltennamen anpassen

Über das Argument column.labels können wir die Spalten benennen. Wir haben zwei Spalten, also müssen wir zwei Namen angeben. Hierfür geben wir einen Vektor mit zwei Elementen an. Das heisst, wir listen die beiden Namen auf und “verpacken” sie als Vektor mit der Funktion c(). Hinweis: Das erste Element des Vektors ist der Name des ersten Modells, das zweite Element ist der Name des zweiten Modells, etc.

stargazer(bi_model, tri_model, 
          type = "text", 
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("Bivariates Modell", "Nettomodell 1"))

## 
## =====================================================================
##                                                                      
##                        Bivariates Modell          Nettomodell 1      
##                               (1)                      (2)           
## ---------------------------------------------------------------------
## hinctnta                    0.095***                  0.029          
##                             (0.024)                  (0.024)         
##                                                                      
## eduyrs                                               0.181***        
##                                                      (0.018)         
##                                                                      
## Constant                    5.643***                 3.935***        
##                             (0.148)                  (0.219)         
##                                                                      
## ---------------------------------------------------------------------
## Observations                 1,242                    1,242          
## R2                           0.012                    0.089          
## Adjusted R2                  0.011                    0.088          
## Residual Std. Error    2.213 (df = 1240)        2.125 (df = 1239)    
## F Statistic         14.990*** (df = 1; 1240) 60.694*** (df = 2; 1239)
## =====================================================================
## Note:                                     *p<0.1; **p<0.05; ***p<0.01

Position des Standardfehlers verändern

Per Default werden die Standardfehler immer unter dem Koeffizienten abgebildet. Dies können wir über das Argument single.row anpassen. In der Dokumentation steht, dass single.row “a logical value” ist. Das bedeutet, wir können nur angeben, ob single.row = TRUE oder ob single.row = FALSE ist. Hinweis: Statt TRUE oder FALSE reicht es auch aus, wenn wir T bzw. F angeben

stargazer(bi_model, tri_model, 
          type = "text", 
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("Bivariates Modell", "Nettomodell 1"), 
          single.row = T)

## 
## =====================================================================
##                                                                      
##                        Bivariates Modell          Nettomodell 1      
##                               (1)                      (2)           
## ---------------------------------------------------------------------
## hinctnta                0.095*** (0.024)          0.029 (0.024)      
## eduyrs                                           0.181*** (0.018)    
## Constant                5.643*** (0.148)         3.935*** (0.219)    
## ---------------------------------------------------------------------
## Observations                 1,242                    1,242          
## R2                           0.012                    0.089          
## Adjusted R2                  0.011                    0.088          
## Residual Std. Error    2.213 (df = 1240)        2.125 (df = 1239)    
## F Statistic         14.990*** (df = 1; 1240) 60.694*** (df = 2; 1239)
## =====================================================================
## Note:                                     *p<0.1; **p<0.05; ***p<0.01

Unnötige Statistiken entfernen

Über das Argument omit.stat können wir bestimmte Statistiken entfernen. Jede Statistik ist durch ein bestimmtes Kürzel gekennzeichnet. Eine Übersicht der Kürzel findet man, wenn man über die Hilfe (z.B. ?stargazer) unter omit.stat auf “list of statistic codes” klickt. Wir entfernen das adjustierte R-Quadrat, den residualen Std.-Fehler und die F-Statistik, da diese für uns irrelevant sind.

stargazer(bi_model, tri_model, 
          type = "text", 
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("Bivariates Modell", "Nettomodell 1"), 
          single.row = T, 
          omit.stat = c("f", "ser", "adj.rsq"))

## 
## ===============================================
##                                                
##              Bivariates Modell  Nettomodell 1  
##                     (1)              (2)       
## -----------------------------------------------
## hinctnta     0.095*** (0.024)   0.029 (0.024)  
## eduyrs                         0.181*** (0.018)
## Constant     5.643*** (0.148)  3.935*** (0.219)
## -----------------------------------------------
## Observations       1,242            1,242      
## R2                 0.012            0.089      
## ===============================================
## Note:               *p<0.1; **p<0.05; ***p<0.01

Nachkommastellen verändern

Wir wollen grundsätzlich zwei (bis drei) Nachkommastellen angeben, aber sicherstellen, dass unsere Zahlen immer mindestens eine substanzielle Stelle haben (also eine Stelle ungleich 0). Die Anzahl an Nachkommastellen können wir mit dem Argument digits verändern. Mit dem Argument digits extra können wir bestimmen, wie viele Nachkommastellen maximal angegeben werden, bis die Zahl eine substantielle Stelle aufweist. Dieses Argument wird derzeit von stargazer noch nicht stringent umgesetzt (mit digits.extra=5 bekommen in einigen Stargazer-Versionen ALLE Zahlen 5 Nachommastellen, sobald eine Zahl mit 4 führenden Nullen im Output ist).

stargazer(bi_model, tri_model, 
          type = "text", 
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("Bivariates Modell", "Nettomodell 1"), 
          single.row = T, 
          omit.stat = c("f", "ser", "adj.rsq"), 
          digits = 2,
          digits.extra = 3)

## 
## =============================================
##                                              
##              Bivariates Modell Nettomodell 1 
##                     (1)             (2)      
## ---------------------------------------------
## hinctnta      0.09*** (0.02)    0.03 (0.02)  
## eduyrs                         0.18*** (0.02)
## Constant      5.64*** (0.15)   3.93*** (0.22)
## ---------------------------------------------
## Observations       1,242           1,242     
## R2                 0.01             0.09     
## =============================================
## Note:             *p<0.1; **p<0.05; ***p<0.01

Signifikanzlevel anpassen

Standardmässig ist das niedrigste dargestellte Signifikanzlevel bei Stargazer auf 0.1 gesetzt. Wir wollen uns aber an die herrschende Konventionen halten und setzen das minimale Signifikanzlevel auf 0.05. Dazu verwende wir den Befehl star.cutoffs und geben ein, in welchen Intervallen wir die 1 bis 3 Sternchen setzen wollen.

stargazer(bi_model, tri_model, 
          type = "text", 
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("Bivariates Modell", "Nettomodell 1"), 
          single.row = T, 
          omit.stat = c("f", "ser", "adj.rsq"), 
          digits = 2,
          digits.extra = 5,
          star.cutoffs = c(0.05, 0.01, 0.001))

## 
## =============================================
##                                              
##              Bivariates Modell Nettomodell 1 
##                     (1)             (2)      
## ---------------------------------------------
## hinctnta      0.09*** (0.02)    0.03 (0.02)  
## eduyrs                         0.18*** (0.02)
## Constant      5.64*** (0.15)   3.93*** (0.22)
## ---------------------------------------------
## Observations       1,242           1,242     
## R2                 0.01             0.09     
## =============================================
## Note:           *p<0.05; **p<0.01; ***p<0.001

Fussnote hinzufügen

Über das Argument notes können wir die Fussnote unserer Tabelle bearbeiten. Wir benutzen die Fussnote insbesondere, um die Datenquelle anzugeben.

stargazer(bi_model, tri_model, 
          type = "text", 
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("Bivariates Modell", "Nettomodell 1"), 
          single.row = T, 
          omit.stat = c("f", "ser", "adj.rsq"), 
          digits = 2,
          digits.extra = 5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = "Daten: ESS(2016), Teilstichprobe CH. Standardfehler in Klammern")

## 
## =============================================================================
##                                                                              
##                      Bivariates Modell                 Nettomodell 1         
##                             (1)                             (2)              
## -----------------------------------------------------------------------------
## hinctnta              0.09*** (0.02)                    0.03 (0.02)          
## eduyrs                                                 0.18*** (0.02)        
## Constant              5.64*** (0.15)                   3.93*** (0.22)        
## -----------------------------------------------------------------------------
## Observations               1,242                           1,242             
## R2                         0.01                             0.09             
## =============================================================================
## Note:                                           *p<0.05; **p<0.01; ***p<0.001
##               Daten: ESS(2016), Teilstichprobe CH. Standardfehler in Klammern

Titel hinzufügen

Mit dem Argument title können wir nun unserer Tabelle einen passenden Titel geben.

stargazer(bi_model, tri_model, 
          type = "text", 
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("Bivariates Modell", "Nettomodell 1"), 
          single.row = T, 
          omit.stat = c("f", "ser", "adj.rsq"), 
          digits = 2,
          digits.extra = 5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = "Daten: ESS(2016), Teilstichprobe CH. Standardfehler in Klammern", 
          title = "Determinanten der Migrationswertschätzung")

## 
## Determinanten der Migrationswertschätzung
## =============================================================================
##                                                                              
##                      Bivariates Modell                 Nettomodell 1         
##                             (1)                             (2)              
## -----------------------------------------------------------------------------
## hinctnta              0.09*** (0.02)                    0.03 (0.02)          
## eduyrs                                                 0.18*** (0.02)        
## Constant              5.64*** (0.15)                   3.93*** (0.22)        
## -----------------------------------------------------------------------------
## Observations               1,242                           1,242             
## R2                         0.01                             0.09             
## =============================================================================
## Note:                                           *p<0.05; **p<0.01; ***p<0.001
##               Daten: ESS(2016), Teilstichprobe CH. Standardfehler in Klammern

Tabellen abspeichern, kopieren und weiterverarbeiten

Inzwischen haben wir unsere Tabelle so angepasst, dass der Inhalt publikationswürdig ist. Nun muss die Tabelle noch aus R exportiert werden (in die Arbeit, in das Poster, in die Präsentation…). Zunächst spezifizieren wir im Befehl type = “html”. Dies stellt sicher, dass Ränder und Linien der Tabelle richtig formatiert werden. Im zweiten Schritt müssen wir über out=““ Dateiformat und Namen des Output-Dokuments festlegen. Z.B. in Word können wir unsere Tabelle dann weiterbearbeiten und beispielsweise die Namen der Zeilen anpassen oder den p-Wert Schlüssel als Fussnote ergänzen bzw. korrigieren.
Hinweis: Unsere Tabelle wurde nun direkt in unsere working directory gespeichert. Wenn du nicht mehr weisst, welchen Dateipfad du dafür festgelegt hattest, kannst du dies mit dem Befehl getwd() herausfinden. Alternativ kannst du über das Files ‘Tab’ (im Hilfe-Fenster unten rechts) deine working directory finden, indem du dort auf “more” klickst und anschliessend auf “Go To Working Directory”.

stargazer(bi_model, tri_model, 
          type = "html", 
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("Bivariates Modell", "Nettomodell 1"), 
          single.row = T, 
          omit.stat = c("f", "ser", "adj.rsq"), 
          digits = 2,
          digits.extra = 5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = "Daten: ESS(2016), Teilstichprobe CH. Standardfehler in Klammern", 
          title = "Determinanten der Migrationswertschätzung",
          out = "reg_table.doc")

2.5 Interpretation

Jetzt können wir anhand unserer Regressionstabelle die Einkommenskoeffizienten vergleichen und stellen fest:

Unter Kontrolle der Bildung steigt die Migrationswertschätzung pro Dezilstufe im Schnitt um 0.03 Skalenpunkte.
Unter Kontrolle der Bildung wird kein signifikanter Einkommenseffekt auf die Migrationseinstellung ausgewiesen.
Der unbereinigte Einkommensschätzer im bivariaten Modell transportiert also zu einem grossen Anteil einen Bildungseffekt und bildet den kausalen Einfluss des Einkommens verzerrt ab.

2.6 Multiple Regressionsanalyse mit mehr als 2 UV

Wir möchten nun in unserer Analyse nicht nur die Bildung, sondern auch das Alter kontrollieren. Folglich integrieren wir zusätzlich die Altersvariable in die Regressionsanalyse. Mathematisch und mit Blick auf die Interpretation ist diese Erweiterung trivial. Auf grafischer bzw. geometrischer Ebene kommen wir mit der vierten Variable in der Regression allerdings an die Grenzen des Vorstellbaren: Wir erweitern unseren 3D-Raum um eine weitere Dimension. Wir haben also einen 4-dimensionalen Raum, in dem sich unsere Merkmalsträger als Punkte verteilen. Jeder Punkt auf diesen Raum hat eine Einkommens-, Bildungs-, Alters- und Migrationseinstellungskoordinate.

Analog zur bi- und trivariaten Regression bestimmen wir im quatrivariaten Fall den geometrischen Körper, welcher die quadrierten Abstände zu den einzelnen Punkten minimiert. Dadurch erhalten wir 3 Koeffizienten, welche den jeweilig eigenständigen Effekt der unabhängigen Variablen auf die abhängige Variable messen.
Nach dieser Logik können wir beliebig viele weitere UV in die Regressionsanalyse integrieren. Jede zusätzliche UV wird dabei als eigenständige, zusätzliche Dimension angelegt, wodurch ihr Effekte auf die AV von den Koeffizienten der anderen UVs entkoppelt wird.

Klingt kompliziert, die Berechnung eines multiplen Regressionsmodells mit mehr als 2 UV bleibt allerdings trivial. Wir erweitern dazu einfach die Variablenliste in lm():

multi_model <- lm(imueclt ~ hinctnta + eduyrs + agea, data = ess8_noNA)
summary(multi_model)

## 
## Call:
## lm(formula = imueclt ~ hinctnta + eduyrs + agea, data = ess8_noNA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -6.7041 -1.5242  0.2662  1.4653  4.9025 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 3.695148   0.302824  12.202   <2e-16 ***
## hinctnta    0.033515   0.024557   1.365    0.173    
## eduyrs      0.183239   0.017759  10.318   <2e-16 ***
## agea        0.003946   0.003449   1.144    0.253    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.125 on 1238 degrees of freedom
## Multiple R-squared:  0.09019,    Adjusted R-squared:  0.08799 
## F-statistic: 40.91 on 3 and 1238 DF,  p-value: < 2.2e-16

Mit stargazer() können wir unsere Ergebnisse wiederum tabellarisch darstellen:

stargazer(bi_model, tri_model, multi_model,
          type = "text", 
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("Bivariates Modell", "Nettomodell 1", "Nettomodell 2"), 
          single.row = T, 
          omit.stat = c("f", "ser", "adj.rsq"), 
          digits = 2,
          digits.extra = 5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          notes = "Daten: ESS(2016), Teilstichprobe CH. Standardfehler in Klammern", 
          title = "Determinanten der Migrationswertschätzung")

## 
## Determinanten der Migrationswertschätzung
## ==============================================================================
##                                                                               
##                 Bivariates Modell       Nettomodell 1        Nettomodell 2    
##                        (1)                   (2)                  (3)         
## ------------------------------------------------------------------------------
## hinctnta         0.09*** (0.02)          0.03 (0.02)          0.03 (0.02)     
## eduyrs                                  0.18*** (0.02)       0.18*** (0.02)   
## agea                                                         0.004 (0.003)    
## Constant         5.64*** (0.15)         3.93*** (0.22)       3.70*** (0.30)   
## ------------------------------------------------------------------------------
## Observations          1,242                 1,242                1,242        
## R2                    0.01                   0.09                 0.09        
## ==============================================================================
## Note:                                            *p<0.05; **p<0.01; ***p<0.001
##                Daten: ESS(2016), Teilstichprobe CH. Standardfehler in Klammern

Achtung: In einer publizierten Regressionstabelle wird meist nur noch das bivariate Modell und das „vollständige” Nettomodell dargestellt. Und: Natürlich gelten auch für diese Tabelle die gleichen Standards zur Weiterverarbeitung und Formatierung wie oben unter 2.4 bzw. in den Folien zum Tutorat dargestellt.

Sofern es sich um eine kausale Hypothese handelt, die der Regressionsanalyse zugrunde liegt, erfolgt die empirische Bewertung der Hypothese immer auf Basis des vollständig spezifizierten Nettomodells - weil dieses ja durch die Abkopplung der Effekte von Störmerkmalen versucht, den kausalen Kern des Zusammenhangs zwischen UV und AV zu isolieren. In der vorliegenden Analyse beurteilen wir die Hypothese folglich als empirisch nicht gestützt, da unter Kontrolle relevanter Störmerkmale kein statistisch signifikanter Einfluss des Einkommens auf die Migrationswertschätzung geschätzt wird.

2.7 Vorhersagen

Vorhersagen lassen sich auch in der multivariaten Regression durch einfaches Einsetzen bestimmter Werte in die Regressionsgleichung ableiten. So ergibt sich für eine (a) 22-jährige Person im (b) 5. Einkommensdezil und mit (c) 14 Bildungsjahren - auf Basis der Regressionskoeffizienten in der Regressionstabelle oben - ein vorhergesagter Wert der Migrationswertschätzung von

3,7 + 0,03°5 + 0,18°14 + 0,004°22 = 6,5.

Der predict()-Befehl kann uns Schreib- bzw. Tipparbeit abnehmen:

predict(multi_model, data.frame(hinctnta=5, eduyrs=14, agea=22))

##        1 
## 6.514878

Wie im bivariaten Fall können vorhergesagte Werte auch (bzw. insbesondere) bei multivariaten Analysen dabei helfen, den Verlauf der Regressionsfunktion zu verdeutlichen und die substantiellen Implikationen eines Regressionskoeffizienten anschaulich zu kommunizieren. Dazu braucht es, wie auch im bivariaten Fall, vorhergesagte Werte zu unterschiedlichen Ausprägungen der unabhängigen Variable, bei deren Berechnung - und das ist neu im multivariaten Fall - die Kontrollvariablen stets dieselbe Ausprägungen einnehmen bzw. konstant gehalten werden sollen:

predict(multi_model, data.frame(hinctnta=c(1,5,10), eduyrs=14, agea=22))

##        1        2        3 
## 6.380818 6.514878 6.682453

Der Vektor c(…) beauftragt R hier mit der Berechnung von Vorhersagen für mehrere, unterschiedliche Ausprägungen der unabhängigen Variable. Da gleichzeitig für die Kontrollvariablen weiterhin jeweils nur ein Wert spezifiziert wird, sind die Kontrollvariablen für alle drei Vorhersagen konstant bzw. in einem bestimmten Wert fixiert. Alle drei Vorhersagen beziehen sich also auf eine 22-jährige Person mit 14 Bildungsjahren: befindet sich diese Person im ersten Einkommensdezil, wird für sie eine Migrationswertschätzung von 6,4 vorhergesagt. Befindet Sie sich dagegen im 5. (bzw. 10.) Einkommensdezil, werden Werte von 6,5 (bzw. 6,7) vorhergesagt.

Es gilt für solche regressionsbasierten Vorhersagesequenzen der Standard, die Kontrollvariablen nicht in einem beliebigen Wert, sondern in ihrem jeweiligen Mittelwert konstant zu halten. Eine entsprechende Umsetzung ist durch vorheriges Ermitteln der Mittel- bzw. Durchschnittswerte und anschliessender Befehlsintegration nach predict möglich:

mean (ess8_noNA$eduyrs)

## [1] 11.40258

mean (ess8_noNA$agea)

## [1] 48.85427

predict(multi_model, data.frame(hinctnta=c(1,5,10), eduyrs=11.4, agea=48.9))

##        1        2        3 
## 6.010542 6.144602 6.312177

Einer Person mit durchschnittlichen Ausprägungen in den Kontrollvariablen wird, wenn sie im ersten Einkommensdezil liegt, eine Migrationswertschätzung von 6,0 vorhergesagt. Befindet sich ihr Einkommen dagegen im 10. Dezil, steigt ihr Vorhersagewert auf 6.3.

Etwas mühsam wird dieses Vorgehen, wenn weitere Kontrollvariablen in die Analyse integriert werden. In diesem Fall ist die Verwendung alternativer Vorhersagekommandos sinnvoll (z.B. emmeans), die per Default alle Kovariaten im Mittelwert konstant halten:

# install.packages("emmeans")
library(emmeans)
emmeans(multi_model,~ hinctnta, at = list(hinctnta=c(1,5,10)))

##  hinctnta emmean     SE   df lower.CL upper.CL
##         1   6.01 0.1255 1238     5.76     6.26
##         5   6.14 0.0615 1238     6.02     6.27
##        10   6.31 0.1263 1238     6.06     6.56
## 
## Confidence level used: 0.95

Werden regressionsbasierte Vorhersagen nicht zur Illustration von Effektgrössen, sondern zur fallbezogenen Vorhersage für eine konkrete Person verwendet, ist es sinnvoll, zusätzlich zum Vorhersagewert auch das Vorhersageintervall (im multivariaten Fall eigentlich: die Länge der Vorhersagehülle) zu ermitteln und zu kommunizieren:

predict(multi_model, data.frame(hinctnta=5, eduyrs=14, agea=22), interval="predict",level = 0.5)

##        fit      lwr      upr
## 1 6.514878 5.078768 7.950988

# Kann "emmeans" leider nicht.

Für eine 22-jährige Person im 5. Einkommensdezil mit 14 Bildungsjahren liegt die Migrationswertschätzung mit 50% Sicherheit zwischen 5.1 und 7.9.

3 Visualisierung

3.1 Conditional Effect Plot - Grundgedanke

Im Falle einer trivariaten Analyse erweitert sich das 2D Streudiagramm zu einem 3D Streuraum. Zudem wird die bestmögliche Regressionsgerade durch die bestmögliche Regressionsebene abgelöst. Die Neigungskoeffizienten dieser Regressionsebene repräsentieren dabei unsere multiplen und somit bereinigten Regressionskoeffizienten. Wie diese darstellen?

Eine 3D-Abbildung veranschaulicht zwar didaktisch wertvoll das Konzept der Regressionsebene, ist in der wissenschaftlichen Kommunikation aber eher untypisch. Stattdessen visualisieren wir eine 2D-Scheibe aus dem Streuraum bzw. aus der Regressionsebene, welche die Neigung der Ebene in der x/y-Dimension - und somit den Koeffizienten der zentralen UV aus der trivariaten Regression - darstellt.

Du kannst dir vorstellen, dass wir, um eine solche Scheibe zu erhalten, uns auf der Achse der Drittvariable “Bildungsjahre” hin und her bewegen und bei einem bestimmten Punkt einen Querschnitt durch den Streuraum vornehmen. Die Gerade, die wir dann erhalten, bildet den Neigungsgradienten der bestmöglichen Ebene in der x/y Dimension und somit den bereinigten Effekt der zentralen UV ab. Auf welcher Höhe diese Ausschnittsgerade verläuft, ist bedingt dadurch, bei wie vielen Bildungsjahren wir den Schnitt ansetzen. Deswegen nennen wir einen solchen Plot auch “Conditional Effect Plot”. Einer Konvention folgend setzen wir den Schnitt beim Mittelwert der Drittvariable(n) an (also in unserem Fall beim Mittelwert der Bildungsvariable).

3.2 Conditional Effect Plot - Umsetzung

Wir können mit dem Befehl visreg() aus dem gleichnamigen Package Conditional Effect Plots herstellen. Die Syntax sieht wie folgt aus: Zunächst geben wir an, von welchem Modell wir einen solchen Plot erzeugen wollen (tri_model). Danach geben wir über xvar= an, welche Variable auf der x-Achse liegt, sprich welchen Koeffizienten wir abbilden möchten. An dieser Stelle ist darauf zu achten, dass wir den Variablennamen als “string” angeben, d.h. in Anführungszeichen setzen. Mit partial = F und rug = F entfernen wir die Punktwolke aus dem Plot. Mit gg = T stellen wir sicher, dass unser Plot ein ggplot-Objekt ist. Dies bewirkt, dass wir den Plot mit unseren bekannten Befehlen wie labs() oder scale_y_continuous() zusätzlich anpassen können. Der visreg()-Befehl erstellt uns dabei per Default einen Conditional Effect Plot, dessen Schnitt durch die Ebene wie üblich beim Mittelwert unserer Bildungsvariablen ansetzt. Achtung: Anders als oben im emmeans-Kommando wird im durch visreg erstellen Plot per Default nicht das arithmetische Mittel, sondern jeweils der Median der Drittmerkmale (bzw. der Modalwert bei kategorialen Variablen, siehe nächste Einheit zugrunde gelegt.

tri_plot<- visreg(tri_model, 
                  xvar = "hinctnta", 
                  gg = T, 
                  partial = F, 
                  rug = F) +
            scale_x_continuous(breaks = seq(from = 1, to = 10, by = 1))+
            scale_y_continuous(breaks = seq(from = 0, to = 10, by = 1))+
            coord_cartesian(ylim = c(0, 10))+
            labs(x = "Einkommen (Dezile)",
                 y = "Migrationswertschätzung",
            title = "Conditional Effect Plot: Einkommen und Migrationswertschätzung",          
            subtitle = "Vorhergesagte Werte unter Konstanthaltung der Bildung im Mittelwert", 
            caption = "Daten: ESS (2016), Teilstichprobe CH (N= 1242), 'Migration ist gut(10) oder schlecht(0) für das kulturelle Leben?'") +
            theme_bw()

tri_plot

Bei Conditional Effect Plots ist sinnvoll, auf der Y-Achse jeweils den Bereich anzuzeigen, der etwa eine Standardabweichung der AV umfasst. Dabei nimmt man den Mittelwert der abhängigen Variablen und fügt jeweils eine halbe Standardabweichung nach oben und nach unten zum Mittelwert hinzu. Dadurch wird die Geradenneigung in Referenz zu einem «typischen» Unterschied gesetzt und kann entsprechend sinnvoll eingeordnet werden. In unserem Fall bedeutet dies, dass wir mit dem Befehl coord_cartesian() die Y-Achse auf die Werte zwischen 5 und 7 einschränken.

mean(ess8_noNA$imueclt, na.rm = T)

## [1] 6.161031

sd(ess8_noNA$imueclt, na.rm = T)

## [1] 2.225374

tri_plot<- visreg(tri_model, 
                  xvar = "hinctnta", 
                  gg = T, 
                  partial = F, 
                  rug = F) +
            scale_x_continuous(breaks = seq(from = 1, to = 10, by = 1))+
            scale_y_continuous(breaks = seq(from = 0, to = 10, by = 0.5))+
            coord_cartesian(ylim = c(5, 7))+
            labs(x = "Einkommen (Dezile)",
                 y = "Migrationswertschätzung",
            title = "Conditional Effect Plot: Einkommen und Migrationswertschätzung",          
            subtitle = "Vorhergesagte Werte unter Konstanthaltung der Bildung im Mittelwert", 
            caption = "Daten: ESS (2016), Teilstichprobe CH (N= 1242), 'Migration ist gut(10) oder schlecht(0) für das kulturelle Leben?'") +
            theme_bw()

tri_plot

Letztlich ist der Bezugspunkt in der Drittvariable für die Visualisierung nicht zentral, da die x/y-Neigung der Regressionsebene für alle Werte von z gleichbleibend ist. Manchmal sind für die Kommentierung eines Effectplots konkrete ganzzahlige Bezugswerte der Drittvariable allerdings etwas dankbarer als «krumme» Mittelwerte. So zeigt unser Plot nun den Einkommensgradienten für eine Person mit mittlerer Bildung - falls Euch aber jemand fragt, wen genau wir uns unter einer Person mit 11.31 Bildungsjahren vorstellen sollen, kommen wir möglicherweise etwas ins Schwitzen. Daher kann es manchmal sinnvoll sein, die Vorhersagefunktion des Effectplots auf eine Person mit anderweitig typischen Merkmalen zu beziehen - hier z.B. eine Person mit BA-Abschluss und entsprechend typischerweise 15 absolvierten Bildungsjahren. Umsetzen können wir dies über die Suboption cond = list() - Achtung: Oft muss der Ankerpunkt auf der y-Achse entsprechend der Höhenverschiebung des Ebenenausschnitts gleich mit angepasst werden.

tri_plot2<- visreg(tri_model, 
                  cond = list(eduyrs=15),
                  xvar = "hinctnta", 
                  gg = T, 
                  partial = F, 
                  rug = F) +
            scale_x_continuous(breaks = seq(from = 1, to = 10, by = 1))+
            scale_y_continuous(breaks = seq(from = 0, to = 10, by = 0.5))+
            coord_cartesian(ylim = c(5.5, 7.5))+
            labs(x = "Einkommen (Dezile)",
                 y = "Migrationswertschätzung",
            title = "Conditional Effect Plot: Einkommen und Migrationswertschätzung",          
            subtitle = "Vorhergesagte Werte für eine Person mit 12 Bildungsjahren", 
            caption = "Daten: ESS (2016), Teilstichprobe CH (N= 1242), 'Migration ist gut(10) oder schlecht(0) für das kulturelle Leben?'") +
            theme_bw()

tri_plot2

Im übrigen können wir die Schriftgrösse von Titeln und Achsenbeschriftungen anpassen, indem wir unserem ggplot-Objekt um einen theme() Befehl ergänzen. Im theme()-Befehl können wir einzelne Elemente unseres Plots angeben (wie axis.text) und dann über element_text(size=Zahl) die Textgrösse ändern (siehe Abnschnitt 4 unten).

3.3 Conditional Effect Plot mit mehr als 2 UV

Bei multivariaten Regressionen mit mehr als zwei unabhängigen Variablen können wir ebenfalls ganz einfach einen Conditional Effect Plot erstellen. Während wir im trivariaten Fall nur für die Bildungsvariable festlegen, wo wir den Schnitt durch den Regressionsraum vornehmen, müssen wir dies im multivariaten Fall für alle Kontrollvariablen festlegen. R bzw. visreg bezieht sich beim Conditional Effect Plot automatisch auf die jeweiligen Mittelwerte der Drittvariablen, was meistens auch sinnvoll und ok ist. Die praktische Umsetzung ist identisch mit dem trivariaten Fall. Der einzige Unterschied ist, dass wir im visreg() den Ergebnisbezug auf das jeweilige multivariate Regressionsergebnis setzen müssen.

multi_plot<- visreg(multi_model, 
                  xvar = "hinctnta", 
                  gg = T, 
                  partial = F, 
                  rug = F) +
            scale_x_continuous(breaks = seq(from = 1, to = 10, by = 1))+
            scale_y_continuous(breaks = seq(from = 0, to = 10, by = 0.5))+
            coord_cartesian(ylim = c(5, 7))+
            labs(x = "Einkommen (Dezile)",
                 y = "Migrationswertschätzung",
            title = "Conditional Effect Plot: Einkommen und Migrationswertschätzung",          
            subtitle = "Vorhergesagte Werte unter Konstanthaltung von Bildung und Alter im Mittelwert", 
            caption = "Daten: ESS (2016), Teilstichprobe CH (N= 1242), 'Migration ist gut(10) oder schlecht(0) für das kulturelle Leben?'") +
            theme_bw()

multi_plot

Konditionalle Effektplots bilden vorhergesagte Werte für Realisationen der UV unter Konstanthaltung der Kontrollvariablen im jeweiligen Mittwelwert ab, allerdings stellen wir schnell fest, dass sich die Vorhersagen von emmeans und die visreg nicht exakt entsprechen:

# install.packages("emmeans")
library(emmeans)
emmeans(multi_model,~ hinctnta, at = list(hinctnta=c(1,5,10)))

##  hinctnta emmean     SE   df lower.CL upper.CL
##         1   6.01 0.1255 1238     5.76     6.26
##         5   6.14 0.0615 1238     6.02     6.27
##        10   6.31 0.1263 1238     6.06     6.56
## 
## Confidence level used: 0.95

Grund dafür ist, dass visreg die Kovariaten im Median konstant hält, emmeans dagegen auf arithmetische Mittel abstellt. Werden Vorhersagen und Effektplot in einen Interpretationszusammenhang gerückt, ist es daher sinnvoll, durch Überschreiben einer der beiden Defaults die beiden Referenzlogiken aneinander anzugleichen. Leichter ist es dabei (gerade mit Blick auf die Einbindung kategorialer Variablen, siehe nächste Einheit), bei emmeans (oder predict) nun ebenfalls Mediane als Grundlage der bedingten Vorhersage zu wählen:

median(ess8_noNA$eduyrs)

## [1] 10

median(ess8_noNA$agea)

## [1] 49

emmeans(multi_model,~ hinctnta, at = list(hinctnta=c(1,5,10), eduyrs=10, agea=49))

##  hinctnta emmean     SE   df lower.CL upper.CL
##         1   5.75 0.1226 1238     5.51     5.99
##         5   5.89 0.0652 1238     5.76     6.02
##        10   6.06 0.1339 1238     5.79     6.32
## 
## Confidence level used: 0.95

…nun passen die abgeleiteten Vorhersagen von emmeans mit den durch visreg dargestellten zusammen.

4 Exportieren von Grafiken

Für gut kalibrierbare und replizierbare Grafiken, seien es wie hier Effektplots oder andere ggplot-basierte Abbildungen, empfehlen wir folgende Schritte:

Nutzt die Optionen zur Definition der Schriftgrössen innerhalb das theme()-Arguments. So könnt ihr ein sinnvolles und über alle Auflösungen konstantes Verhältnis von Plotgrösse und Schriftgrösse herstellen. Empfehlung zum size: Titel etwa 18 (Wenn es seine Länge erlaubt; in diesem Beispiel nicht der Fall), Achsentitel etwas kleiner, Wertelabel nochmal etwas kleiner - Achtung: Es gibt keine Standardlösung, die immer funktioniert, geht also nicht ohne ausprobieren.
Speichern der Abbildung z.B. mit ggsave. Dieser Befehl nimmt automatisch Bezug auf die letzte angelegte Grafik. Die Abbildungsdatei befindet sich dann automatisch in eurer Working Directory. Dies ist per Default der Ordner, in dem der Syntax-file gespeichert ist, über den ihr R steuert.
Spezifikation der Abbildungsattribute und Export mit ggsave:
- Höhe/Breite: Das relative Verhältnis von Höhe zu Breite hängt stark von den Eigenschaften Eurer Abbildung sowie der beteiligten Variablen ab. Die absolute Grösse der beiden Seiten sollte zwischen 10 und 20 Zentimetern liegen. Sie muss ggf. auf die Schriftgrösse abgestimmt werden. Auch hier gilt: ohne Probieren geht es nicht!
- Auflösung: Gegeben der oben vorgeschlagenen Grösse seit ihr für Präsentationen bei 1000 dpi auf der sicheren Seite. Für Poster braucht es ggf. noch etwas mehr.
- Format: .png hat sich bewährt, manche empfehlen auch .pdf

multi_plot + 
    theme(
    plot.title = element_text(size = 16),  # Increase title size
    axis.title.x = element_text(size = 14),  # Increase X axis label size
    axis.title.y = element_text(size = 14),  # Increase Y axis label size
    axis.text.x = element_text(size = 13),   # Increase X axis text size
    axis.text.y = element_text(size = 13)    # Increase Y axis text size
  )

ggsave(filename = "effectplot.png", width = 20, height = 13, dpi=1000, units = "cm")

(schneidet in der Online-Version in einigen Browsern in diesem Beispiel zwar die Überschrift ab, erzeugte aber eine 1a-Grafik im Speicher)

3D-Plots

Disclaimer: Im folgenden Abschnitt wird das Herstellen von 3D-Plots erklärt. An dieser Stelle wird aber nochmals betont, dass 3D-Plots in der wissenschaftliche Kommunikation nicht gebräuchlich sind und nicht in eure Arbeiten eingebunden werden sollten. Allerdings eignen sich 3D-Plots gut für didaktische Zwecke, sowie auch die Vermittlung ausserhalb der Wissenschaft.

Wir können 3D-Scatterplots mit dem Befehl plot3d() aus dem Package rgl herstellen. Jedes Mal wenn wir über rgl-Funktionen einen Plot herstellen, wird in R ein zusätzliches Fenster geöffnet.

plot3d(x=ess8_noNA$hinctnta,
           z=ess8_noNA$imueclt,
           y=ess8_noNA$eduyrs)

Achsen beschriften.

plot3d(x=ess8_noNA$hinctnta,
       z=ess8_noNA$imueclt,
       y=ess8_noNA$eduyrs, 
       xlab = "Einkommen", 
       ylab="Bildungsjahre", 
       zlab = "Migrationseinstellung")

Grösse der Punkte verändern.

plot3d(x=ess8_noNA$hinctnta,
       z=ess8_noNA$imueclt,
       y=ess8_noNA$eduyrs, 
       xlab = "Einkommen", 
       ylab="Bildungsjahre", 
       zlab = "Migrationseinstellung", 
       size =5)

Durchsichtigkeit der Punkte verändern.

plot3d(x=ess8_noNA$hinctnta,
       z=ess8_noNA$imueclt,
       y=ess8_noNA$eduyrs, 
       xlab = "Einkommen", 
       ylab="Bildungsjahre", 
       zlab = "Migrationseinstellung", 
       size =5, 
       alpha = 0.25)

Plot automatisch rotieren.

plot3d(x=ess8_noNA$hinctnta,
       z=ess8_noNA$imueclt,
       y=ess8_noNA$eduyrs, 
       xlab = "Einkommen", 
       ylab="Bildungsjahre", 
       zlab = "Migrationseinstellung", 
       size =5, 
       alpha = 0.25)

play3d(spin3d( axis = c(0, 0, 1), rpm = 2), duration = 30)

Plot speichern (statisch).

plot3d(x=ess8_noNA$hinctnta,
       z=ess8_noNA$imueclt,
       y=ess8_noNA$eduyrs, 
       xlab = "Einkommen", 
       ylab="Bildungsjahre", 
       zlab = "Migrationseinstellung", 
       size =5, 
       alpha = 0.25)

rgl.snapshot(filename="3d_scatter.png", fmt = "png")

Plot speichern (gif).

plot3d(x=ess8_noNA$hinctnta,
       z=ess8_noNA$imueclt,
       y=ess8_noNA$eduyrs, 
       xlab = "Einkommen", 
       ylab="Bildungsjahre", 
       zlab = "Migrationseinstellung", 
       size =5, 
       alpha = 0.25)

movie3d(movie ="scatter_rotation", 
  spin3d( axis = c(0, 0, 1), rpm = 2),
  duration = 30,
  dir = getwd(), 
  type = "gif", 
  clean = T)

Regressionsebene einfügen.

plot3d(lm(imueclt~hinctnta+eduyrs, data=ess8_noNA),
       xlab = "Einkommen", 
       ylab="Bildungsjahre", 
       zlab = "Migrationseinstellung", 
       plane.col = "brown1")

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

IK Kapitel 6.2.1 & 6.2.2: Chester Ismay und Albert Y. Kim (2021): Statistical Inference via Data Science

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

Hier gehts weiter zur Übung V

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”