invisible header

Das hier vorgestellte Analyseskript dient als vereinfachtes Muster einer empirischen Ausarbeitung der multiplen Regressionsanalyse auf BA-Niveau am Soziologischen Institut der UZH. Es verdeutlicht zudem unsere Standards für Analyse und Ergebnisdarstellung. Gewähr für die Funktionalität des Codes und Korrektheit der Ergebnisse wird nicht übernommen. Fehlermeldungen und Verbesserungsvorschläge bitte an den Autor.

Das Skript behandelt die regressionsbasierte Analyse einer binären abhängigen Variablen. Dieses Skript ist weniger ausführlich aufbereitet und wird seltener aktualisiert als das zur Analyse eines metrischen abhängigen Merkmals. Es dient daher insbesondere der Betonung von Unterschieden in der analytischen Aufbereitung einer metrischen vs. dichotomen AV. Beispielhaft widmen wir uns hier der Frage, inwieweit die Wohnungsgrösse die Neigung zu Übergewicht beeinflusst.

 

Vorab: Lineare vs. Logistische Regression

Im Falle einer binären AV muss die Wahl zwischen einer einfachen linearen Modellierung (Lineares Wahrscheinlichkeitsmodell) und einer komplexen Modellierung (Logistische Regression) getroffen werden. Wir finden grundsätzlich beide Varianten legitim. Galt die logistische Regression lange Zeit als Goldstandard der Analyse binärer abhängiger Variablen, weichen heute die Forschungsstandards in dieser Frage auf. Im Einzelfall mögen neben individuellen Präferenzen insbesondere Nuancen der Zusammenhangsannahme entscheiden: Ist es plausibel, von einem eng gefassten Transitionsbereich der unabhängigen Variablen auszugehen, in dem die Wahrscheinlichkeit für das durch die Kategorie 1 beschriebene Ereignis rasant zu- oder abnimmt? Dies wäre dann ein Argument für die logistische Modellierung. Wenn ihr es besonders gründlich machen wollt, führt beide Analyseformen durch. Setzt dann aber eine der beiden Modellierungsoptionen lediglich in den Anhang Eurer Arbeit und verweist auf diesen Anhang, bzw. die grundlegende Robustheit der Ergebnisse, in einer Fussnote zur Ergebnisauswertung.

Für den hier zu untersuchenden Zusammenhang scheint die funktionale Form der logistischen Regression einigermassen plausibel, es spricht a priori aber auch nichts gegen die Modellierung des durchschnittlichen linearen Effektes im Rahmen einer linearen Wahrscheinlichkeitsregression. Die lineare Modellierung einer dichotomen Variable hat grundsätzlich den Vorteil, dass sie klar und leicht zu interpretierende Koeffizienten erzeugt und zudem robustere Ergebnisse insbesondere bei der Analyse von Moderations- und Mediationseffekten generiert.

Untersucht wird der Effekt der Wohnungsgrösse auf Übergewichtigkeit auf Basis des SHP. Die Wohnungsgrösse wird, stark vereinfacht, über die Anzahl an Zimmern in der Wohnung gemessen und stammt aus dem Haushaltsdatensatz des SHP. Übergewicht wird als Quotient von Grösse und Gewicht, über den BMI, nach Standards der WHO operationalisiert. Die Analyse ist querschnittlich ausgerichtet, obgleich Fragestellung & Variationseigenschaften eine längsschnittliche Analyseperspektive eröffnen - und auch das SHP grundsätzlich auf längsschnittliche Analysen ausgerichtet ist.

Im querschnittlichen Ansatz erfolgt die empirische Effektabsicherung ausschliesslich über die Integration von Kontrollvariablen. Daher werden in der Analyse zusätzlich das Geschlecht, der Erwerbsstatus, Bildung, das Haushaltseinkommen, der Zivilstand sowie das Alter berücksichtigt. Das Geschlecht wird zudem als Moderator des Wohnungsgrösseneffektes modelliert, da das (hier unsichtbare) Erklärungsmodell unterschiedlich starke Wirkungsmechanismen für Männer und Frauen postuliert hat.

Der vorgelagerte Prozess des Datenmanagements ist nicht in diesem Skript dokumentiert, sondern auf der Seite zum Datenmanagement. Infos zum SHP gibt es auf der Homepage von FORS.

Los geht’s wie immer mit dem Kopf der Syntax:

# BAa
# SHP Analyse zur Übergewichtigkeit
# Datum: 23.03.2021
# AutorIn: XXX

Als nächstes lesen wir die zuvor produzierte analysebereite Datenmatrix ein und installieren/aktivieren die für die Analyse relevanten packages.

# Daten einlesen (z.B. von "c:\meine_BA")
library(haven)
datensatz <- read_dta("mein_laufwerk/mein_zielverzeichnis/shp_base2.rds")
# install.packages("package")
library(tidyverse)
library(ggplot2)
library(visreg)
library(summarytools)
library(stargazer)
library(fastDummies)
library(labelled)
library(ggridges)

 

1 Vorbereitungen und Inspektion

Bevor es losgeht verschaffen wir uns einen Überblick über den Datensatz und die Variablen:

dim (shp_base2)
## [1] 9338   28
tibble::view(shp_base2)
summary(shp_base2)
##       pid             w11101_2017         alter         educyears    
##  Min.   :5.101e+03   Min.   :   0.0   Min.   :18.00   Min.   : 9.00  
##  1st Qu.:7.915e+06   1st Qu.: 407.4   1st Qu.:38.00   1st Qu.:12.00  
##  Median :2.133e+07   Median : 592.0   Median :53.00   Median :12.00  
##  Mean   :3.103e+07   Mean   : 723.3   Mean   :52.26   Mean   :13.53  
##  3rd Qu.:6.268e+07   3rd Qu.: 816.6   3rd Qu.:67.00   3rd Qu.:15.00  
##  Max.   :1.054e+09   Max.   :5205.8   Max.   :99.00   Max.   :21.00  
##                                                       NA's   :15     
##    geschlecht      zivilstand        hhgr         numkids           erwerb     
##  Min.   :1.000   Min.   :1.00   Min.   : 1.0   Min.   :0.0000   Min.   :1.000  
##  1st Qu.:1.000   1st Qu.:1.00   1st Qu.: 2.0   1st Qu.:0.0000   1st Qu.:1.000  
##  Median :2.000   Median :1.00   Median : 2.0   Median :0.0000   Median :2.000  
##  Mean   :1.542   Mean   :1.55   Mean   : 2.6   Mean   :0.3183   Mean   :1.999  
##  3rd Qu.:2.000   3rd Qu.:2.00   3rd Qu.: 4.0   3rd Qu.:0.0000   3rd Qu.:3.000  
##  Max.   :2.000   Max.   :5.00   Max.   :10.0   Max.   :4.0000   Max.   :3.000  
##                  NA's   :1                                      NA's   :5      
##     hhnetto          idhous_2017           height          weight      
##  Min.   :   146.7   Min.   :      51   Min.   :1.180   Min.   : 36.00  
##  1st Qu.: 57085.7   1st Qu.:   79211   1st Qu.:1.640   1st Qu.: 62.00  
##  Median : 86623.7   Median :  213346   Median :1.700   Median : 71.00  
##  Mean   : 95507.9   Mean   :  313594   Mean   :1.709   Mean   : 72.58  
##  3rd Qu.:121032.6   3rd Qu.:  626831   3rd Qu.:1.780   3rd Qu.: 81.00  
##  Max.   :736158.6   Max.   :10539611   Max.   :2.060   Max.   :180.00  
##  NA's   :8                             NA's   :325     NA's   :419     
##     lifesat           health          kanton         zustand     
##  Min.   : 0.000   Min.   :1.000   Min.   : 1.00   Min.   :1.000  
##  1st Qu.: 8.000   1st Qu.:2.000   1st Qu.: 6.00   1st Qu.:2.000  
##  Median : 8.000   Median :2.000   Median :16.00   Median :2.000  
##  Mean   : 8.081   Mean   :1.994   Mean   :14.44   Mean   :2.162  
##  3rd Qu.: 9.000   3rd Qu.:2.000   3rd Qu.:23.00   3rd Qu.:2.000  
##  Max.   :10.000   Max.   :5.000   Max.   :26.00   Max.   :3.000  
##  NA's   :313      NA's   :13                      NA's   :28     
##    anz_zimmer         noise           sample   status_2017     
##  Min.   : 1.000   Min.   :1.000   Min.   :1   Min.   :0.00000  
##  1st Qu.: 3.000   1st Qu.:2.000   1st Qu.:1   1st Qu.:0.00000  
##  Median : 4.000   Median :2.000   Median :1   Median :0.00000  
##  Mean   : 4.479   Mean   :1.792   Mean   :1   Mean   :0.03309  
##  3rd Qu.: 5.000   3rd Qu.:2.000   3rd Qu.:1   3rd Qu.:0.00000  
##  Max.   :10.000   Max.   :2.000   Max.   :1   Max.   :1.00000  
##  NA's   :346      NA's   :19                                   
##      idpers              erwach       oecd_weight     oecd_inc_mon     
##  Min.   :5.101e+03   Min.   :1.000   Min.   :1.000   Min.   :   12.22  
##  1st Qu.:7.915e+06   1st Qu.:2.000   1st Qu.:1.500   1st Qu.: 3061.22  
##  Median :2.133e+07   Median :2.000   Median :1.500   Median : 4227.94  
##  Mean   :3.103e+07   Mean   :2.282   Mean   :1.736   Mean   : 4639.25  
##  3rd Qu.:6.268e+07   3rd Qu.:3.000   3rd Qu.:2.100   3rd Qu.: 5635.79  
##  Max.   :1.054e+09   Max.   :7.000   Max.   :4.700   Max.   :44435.40  
##                                                      NA's   :8         
##       bmi          overweight     zivilstand_class  
##  Min.   :13.06   Min.   :0.0000   Length:9338       
##  1st Qu.:21.80   1st Qu.:0.0000   Class :character  
##  Median :24.21   Median :0.0000   Mode  :character  
##  Mean   :24.77   Mean   :0.1061                     
##  3rd Qu.:26.99   3rd Qu.:0.0000                     
##  Max.   :66.79   Max.   :1.0000                     
##  NA's   :529     NA's   :529

Der Datensatz enthält 26 Variablen und 9338 Personen bzw. Beobachtungen. Die Verteilung der Altersvariable und die konstante Ausprägung der von uns angelegten sample-Variable lassen erkennen, dass wir die Stichprobe bereits im Zuge der Datenaufbereitung eingegrenzt haben: Unsere Grundgesamtheit bilden alle erwachsene Personen der Schweiz, die Stichprobe und Daten sind dementsprechend begrenzt auf alle Personen, die mindestens 18 Jahre alt sind und die ein gültiges Interview abgeliefert haben. Zur Sicherheit führen wir die entsprechende Fallauswahl noch mal aus (auch wenn sie hier den bestehenden Datensatz 1:1 reproduziert):

shp_sample<-filter(shp_base2, alter>=18 & sample==1)
dim (shp_base2)
## [1] 9338   28

In unsere Analyse gehen auch kategoriale unabhängige Variablen als sogenannte Faktoren ein. Damit R diese auch solche verarbeitet (und nicht etwa als metrische Variablen versteht), definieren wir sie vorab typkonsistent:

shp_sample$geschlecht<-as_factor(shp_sample$geschlecht)
shp_sample$erwerb<-as_factor(shp_sample$erwerb)
shp_sample$zivilstand<-as_factor(shp_sample$zivilstand)
shp_sample$overweight<-as_factor(shp_sample$overweight)

Achtung: Mögliche Probleme im Rahmen von Faktorisierungen, die gegebenenfalls Postfaktorisierungs-Bereinigungen im Datensatz erfordern, werden im Skript zur Analyse eines metrischen abhängigen Merkmals behandelt.

Alle relevanten kategorialen Variablen sind nun als Faktoren im Datensatz angelegt, wie sich leicht nachprüfen lässt:

sapply(shp_sample, class)

Univariate Illustrationen (Histogramme etc.) von UV und AV gehören zwar meist nicht in die Schriftfassung eurer Arbeit, sind aber nützlich für die interne Inspektion und Vergegenwärtigung der zentralen Variablen.

ggplot(shp_sample, aes(x = anz_zimmer)) +
  geom_bar(aes(y = (..count..)/sum(..count..)), fill="blue") +
  xlab("")+ ylab("Anzahl")+
  scale_x_continuous(breaks=c(0,1,2,3,4,5,6,7,8,9,10))+
  scale_y_continuous(labels=scales::percent,breaks=c(0,0.10,0.20,0.30, 0.40))+
  labs(title="Wohnungsgrösse in der Schweiz 2017",
          subtitle="Anzahl Zimmer in der Wohnung",
          caption="SHP v19, erwachsene Personen, n=10.936")

Die Verteilung der Wohnungsgrösse folgt näherungsweise einer Normalverteilung. Schön, wenn auch zur Einbindung in die Regressionsanalyse kein notwendiges Kriterium. Wichtig ist in erster Linie die erkennbare Existenz von Variation.

Die Visualisierung der Verteilung einer binären Variablen halten wir im Rahmen einer wissenschaftlichen Arbeit nicht für hilfreich.

 

2 Stichprobenstatistik

Wichtiges Element der empirischen Ausarbeitung und zentraler Bestandteil des Methodenteils ist die Übersichtstabelle zur Stichprobenstatistik. Durch Beschränkung des Datensatzes auf die analysebewährten Merkmale erleichtern wir die Erstellung der Tabelle (und auch der weiteren empirischen Bearbeitung):

shp_sample<-select(shp_sample, pid, erwerb, oecd_inc_mon, alter, educyears, geschlecht, anz_zimmer, zivilstand, overweight, bmi, w11101_2017)

Eine einfache, aber effektive Funktion zur Erstellung einer deskriptiven Übersichtstabelle bietet das Package table1. Praktischerweise verarbeitet dies Faktor-Variablen und weist für diese in der generierten Tabelle automatisch Häufigkeiten (statt Lage- und Streuungsmasse) aus.

Achtung: Achtet vor Anwendung des Tabellen-Kommandos darauf, dass alle kategorialen Variablen als Faktoren definiert (siehe oben) wurden. Ansonsten würden deren Ausprägungen vom Kommando als numerische Werte interpretiert und entsprechend unsinnige Tabelleneinträge generieren.

library(table1)
table1(~bmi+overweight+anz_zimmer+geschlecht+erwerb+zivilstand+oecd_inc_mon+alter+educyears, data=shp_sample)
Overall
(N=9338)
bmi
Mean (SD) 24.8 (4.30)
Median [Min, Max] 24.2 [13.1, 66.8]
Missing 529 (5.7%)
overweight
0 7874 (84.3%)
1 935 (10.0%)
Missing 529 (5.7%)
Number of rooms: not including the kitchen and the bathroom
Mean (SD) 4.48 (1.51)
Median [Min, Max] 4.00 [1.00, 10.0]
Missing 346 (3.7%)
Gender of Individual
.S (-3) 0 (0%)
.M (-2) 0 (0%)
.C (-1) 0 (0%)
Male 4274 (45.8%)
Female 5064 (54.2%)
Employment Level of Individual
.S (-3) 0 (0%)
.M (-2) 0 (0%)
.C (-1) 0 (0%)
Full Time 3452 (37.0%)
Part Time 2441 (26.1%)
Not Working 3440 (36.8%)
Missing 5 (0.1%)
Marital Status of Individual
.S (-3) 0 (0%)
.M (-2) 0 (0%)
.C (-1) 0 (0%)
Married or Living with a Partner 6299 (67.5%)
Single, not Living with a Partner 1788 (19.1%)
Widowed, not Living with a Partner 516 (5.5%)
Divorced, not Living with a Partner 617 (6.6%)
Separated (Legally Married), not Living with a Partner 117 (1.3%)
Not Living with a Partner (Individuals age 18 and older) 0 (0%)
Not Living with a Partner (Individuals under age 18) 0 (0%)
Missing 1 (0.0%)
oecd_inc_mon
Mean (SD) 4640 (2620)
Median [Min, Max] 4230 [12.2, 44400]
Missing 8 (0.1%)
Age of Individual
Mean (SD) 52.3 (18.3)
Median [Min, Max] 53.0 [18.0, 99.0]
Number of Years of Education
Mean (SD) 13.5 (3.10)
Median [Min, Max] 12.0 [9.00, 21.0]
Missing 15 (0.2%)


Es fehlen noch:

…umsetzbar am besten direkt in Word nach copy & paste.

Bei den in der Tabelle ausgewiesenen Missing-Kategorien (-3,-2,-1, other error etc.) handelt es sich um Relikte ursprünglicher Missing-Kodierungen. Diese Relikte sind nach wie vor in einigen Variable angelegt - obgleich wir die entsprechenden Ausprägungen längst in NAs umcodiert haben, wie auch die ausgewiesenen Belegungen in den Missing-Kategorien zeigen. Vor der Publikation der Tabelle sollten diese Kategorien-Relikte im Rahmen der manuellen Nachbearbeitung entfernt werden. Achtung: Im Skript zur Analyse eines metrischen abhängigen Merkmals werden im ersten Abschnitt Optionen für Faktorvariablen behandelt, mit denen das Problem solcher Phantomkategorien umgangen werden kann.

Unkonventionell am Output von table1 ist die zeilenübergreifende Organisation der statistischen Kennwerte, welche in der traditionellen Stichprobenstatistik die Spalten der Übersichtstabelle definieren. Insbesondere bei Analysen mit vielen Variablen bricht die durch table1 arrangierte Aufbereitung daher mit dem Leitbild darstellungsökonomischer Sparsamkeit und Übersichtlichkeit.

Alternative Darstellungen bieten dann die Packages summarytools und stargazer. Allerdings verarbeiten diese keine Faktor-Variablen. Diese müssen folglich zunächst in Sets numerischer Dummy-Variablen rekodiert werden, wenn ihr summarytools und stargazer zur Produktion der Übersichtsstatistik nutzen wollt:

shp_sample_D<-dummy_cols(shp_sample, select_columns=c("overweight", "geschlecht", "erwerb", "zivilstand"), ignore_na = TRUE)

…so werden im (neu gebildeten) Datensatz neben den grundlegenden Faktorvariablen auch die davon abgeleiteten Dummies angelegt - und zwar formal als metrische Variablen (check: is.numeric(shp_sample_D$geschlecht_Male)). Durch die doppelte Variablenredefinition (metrisch→faktor→dummy) werden kategoriale Variablen nun auch von den alternativen Kommandos zur Stichprobenstatistik korrekt aufgegriffen:

stats<-descr(shp_sample_D, stats = c("N.Valid", "Pct.Valid", "mean", "sd", "min", "max"), transpose=T)
print(stats, method="browser", file="stats.htm")
statsdf<-data.frame(stats)
stats

Zugegeben: Direkt veröffentlichungswürdig ist der Output so noch nicht, ABER: der Weg dahin ist nun relativ kurz - der “print”-Befehl legt eine html-Tabelle im aktuellen Working Directory an, die leicht z.B. nach Word kopiert & dort weiterbearbeitet werden kann. Durch

…entsteht binnen 5 Minuten eine veröffentlichungswürdige Tabelle zur Stichprobenstatistik (siehe die beispielhafte Ausarbeitung der Übersichtsstatistik auf der Nachbarseite).

Alternativ hätten wir die oben angelegte descr-Output-Matrix stats_df zur Weiterverarbeitung in Word oder Excel verwenden oder auch das Package stargazer nutzen können, welches auch ein gut übertragbares html-Dokument kreiert.

stargazer(as.data.frame(shp_sample_D), type="text", out="table1.html",
          title="Übersichtstabelle: Univariate Statistik", digits=2)

Allerdings müssen hier im Rahmen der Nachbearbeitung noch die Anteile der NAs per Hand berechnet und ergänzt werden.

Achtung Stichprobenfresser!

Die Übersicht zur Stichprobenstatistik, die typischerweise vor der Reduktion des Datensatzes auf die Modellstichprobe erstellt wird, ist zwar ein wichtiges Instrument der Datenkommunikation, dient aber auch der inneren Vergegenwärtigung möglicher Probleme der Variablenkomposition. Merke: Die besetzungsschwächste Drittvariable prägt die Grösse der Modellstichprobe; daher muss die Kontrollleistung einer Variable immer gegenüber den mit ihrer Integration verbundenen Stichprobeneinbussen abgewogen werden. Die vorliegende Stichprobenstatistik weisst die Übergewichtigkeit als ausfallsensibelste Variable aus: da dies eine der beiden zentralen Variablen der Analyse ist, sind die entsprechenden Ausfälle hier alternativlos.

In vielen anderen denkbaren Szenarios könntet ihr aber auf Grundlage eines Nebendarstellers einen Grossteil der Stichprobe verlieren, was in vielen Fällen inakzeptabel wäre. Ihr würdet dann stattdessen entweder (a) auf die Kontrolle des entsprechenden Stichprobenfressers verzichten, (b) eine andere Variable, die ein ähnliches Konzept abbildet, verwenden, oder (c) eine modellierbare Zusatzkategorie “Fehlend” zu den Ausprägungen des Stichprobenfressers integrieren, in welche die nicht-modellierbaren Missings verschoben werden.

 

3 Bivariate Inspektion und Visualisierung

Die deskriptive Aufarbeitung des bivariaten Zusammenhangs von UV und AV dient in erster Linie der Datenexploration. Gleichwohl wertet eine gelungene Visualisierung des bivariaten Zusammenhangs jede Forschungsarbeit auf. Das grundlegende bivariate Visualisierungsinstrument ist das einfache Streudiagramm:

shp_sample_scatter<-filter(shp_sample, !is.na(overweight))
scatterplot <- ggplot(shp_sample_scatter, 
       aes(x = anz_zimmer, y = overweight)) + 
  geom_point() +
  scale_x_continuous(breaks = seq(0,10,1)) +
  scale_y_discrete(breaks = seq(0,1,1)) +
    labs(x = "Anzahl Zimmer", 
       y = "Übergewichtigkeit",
       title = "Wohungsgrösse und Übergewicht in der Schweiz",
       caption = "Quelle: SHP v.19, n = 8809") + 
  theme_classic()
scatterplot

Das für grosse Datensätze typische Problem Overplotting ist in dieser Anwendung noch durch die Dichotomie der abhängigen Variablen verstärkt, das Streudiagramm daher völlig unbrauchbar. Die folgenden Varianten versuchen die Anteilswerte auf verschiedene Arten graphisch herausarbeiten. Empfehlung: In der konkreten empirischen Situation immer nur eine der Varianten - nämlich die, die den Charakter des Zusammenhangs am besten abzubilden vermag - in den Bericht (die Präsentation, das Poster) integrieren

Die geom_jitter-Funktion fügt den Werten eine Streuungskomponente zu, zudem lässt das alpha-Argument die Punkte transparent erscheinen, so dass Überlagerungen besser sichtbar werden. Die jeweilige Gruppengrösse wird so deutlich dargestellt und einigermassen gut kontrastiert.

scatterplot <- ggplot(shp_sample_scatter, 
       aes(x = anz_zimmer, y = overweight)) +  
  geom_jitter(size=2, alpha=0.1) +
  scale_x_continuous(breaks = seq(0,10,1)) +
  scale_y_discrete(breaks = seq(0,1,1)) +
    labs(x = "Anzahl Zimmer", 
       y = "Übergewichtigkeit",
       title = "Wohungsgrösse und Übergewicht in der Schweiz",
       caption = "Quelle: SHP v.19, n = 8809") + 
  theme_classic()
scatterplot

Die Heatmap vermag in diesem Beispiel ebenfalls den Charakter des Zusammenhangs zu visualisieren. Dieses Visualisierungsinstrument ist insbesondere dann hilfreich, wenn wir es mit grossen Datensätzen zu tun haben. Andererseits ist ein konservativer Umgang mit dieser Darstellungsoption geboten, da zum einen das zentrale Repräsentationsniveau (1 Messung = 1 Darstellungspunkt) aufgegeben wird und zum anderen die implizite Klassierung Zusammenhänge verschleiern oder überbetonen kann. Wichtig: Kalibriere über das bins-Argument die Geographie bzw. Klassenbildung der Map so, dass keine strukturbedingten Lücken in der Abbildung entstehen. Nur dann können abgebildete Lücken eindeutig, nämlich als geringe bzw. fehlende Belegung, interpretiert werden.

heatmap <- ggplot(shp_sample_scatter, 
       aes(x = anz_zimmer, y = overweight)) +
  geom_bin2d(bins = c(9, 4)) +
  scale_x_continuous(breaks = seq(0,10,1)) +
  scale_y_discrete(breaks = seq(0,1,1)) +
    scale_fill_gradient(low = "azure", high="darkblue") +
    labs(x = "Anzahl Zimmer", 
       y = "Übergewichtigkeit",
       title = "Wohungsgrösse und Übergewicht in der Schweiz",
       caption = "Quelle: SHP v.19, n = 8809") + 
  theme_classic()
heatmap

Der Bubbleplot, realisiert über die geom_count - Funktionalität, wird oft zur Visualisierung bivariater Zusammenhänge verwendet. Die “Bubble”-Grösse repräsentiert dabei die Anzahl an Personen in der jeweiligen Merkmalskombination.

scatterplot <- ggplot(shp_sample_scatter, 
       aes(x = anz_zimmer, y = overweight)) +  
  geom_count() +
  scale_size_area(max_size = 15)+
  scale_x_continuous(breaks = seq(0,10,1)) +
  scale_y_discrete(breaks = seq(0,1,1)) +
    labs(x = "Anzahl Zimmer", 
       y = "Übergewichtigkeit",
       title = "Wohungsgrösse und Übergewicht in der Schweiz",
       caption = "Quelle: SHP v.19, n = 8809") + 
  theme_classic()
scatterplot

Am sinnvollsten erscheint in diesem konkreten Beispiel allerdings die Darstellung des Zusammenhangs über ein Balkendiagram, welches konditionale Anteilswerte abbildet:

ggplot(shp_sample_scatter, aes(x = anz_zimmer, fill = overweight, na.rm= TRUE )) + 
  geom_bar(position = "fill", na.rm=TRUE ) + 
  labs(title = "Übergewicht nach Anzahl Zimmern",
       x = "Anzahl Zimmer", y = "Prozent", 
       caption="Quelle: SHP") + 
  scale_y_continuous(labels = scales::percent_format()) +
  scale_x_continuous(breaks=seq(1,10,1)) +
  theme_bw()

Fazit: Die Abbildungen vermögen alle nicht 100%ig zu überzeugen, was insbesondere den vielen Datenpunkten bzw. der grossen Stichprobe geschuldet ist. Der Jitterplot und der Bubbleplot funktionieren einigermassen, wären allerdings auch kaum veröffentlichungswürdig. Sinnvoll ist in diesem Beispiel eine einfache Darstellung des Zusammenhangs über bedingte Anteilswerte, in der wir Anteilswerte im Übergewicht zwischen Bewohnern verschiedener Wohnungsgrössen miteinander vergleichen. Diese Darstellung könnte entweder als einfache Kreuztabelle oder, wie oben angedeutet, ein Balkendiagramm erfolgen.

 

4 Linearitätsdiagnose

Die Linearität bzw. die Form des Zusammenhangs sollte vor jeder Regressionsanalyse mit metrischer unabhängiger Variable theoretisch begründet und empirisch abgestützt werden. Obgleich nicht jede Linearitätsabweichung direkt zu einer Transformation führen muss, können so die Analyseergebnisse später besser eingeordnet werden. Die Diagramme oben lassen im Rahmen der visuellen Zusammenhangsinspektion Zweifel an der Linearität des Zusammenhangs aufkommen. Wir nutzen nun das Instrument der Multigruppenanalyse, um unseren Verdacht empirisch zu fundieren.

Die Multigruppenanalyse vergleicht die Steigungskoeffizienten zwischen verschiedenen Abschnitten der UV:

shp_sample_g1<-filter(shp_sample, anz_zimmer<=4)
shp_sample_g2<-filter(shp_sample, anz_zimmer>=4)
lm(overweight~anz_zimmer, shp_sample_g1)
## Warning in model.response(mf, "numeric"): Benutzung von type="numeric" wird für
## eine faktorielle Zielgröße ignoriert
## Warning in Ops.factor(y, z$residuals): '-' ist nicht sinnvoll für Faktoren
## 
## Call:
## lm(formula = overweight ~ anz_zimmer, data = shp_sample_g1)
## 
## Coefficients:
## (Intercept)   anz_zimmer  
##     1.08237      0.01241
lm(overweight~anz_zimmer, shp_sample_g2)
## Warning in model.response(mf, "numeric"): Benutzung von type="numeric" wird für
## eine faktorielle Zielgröße ignoriert

## Warning in model.response(mf, "numeric"): '-' ist nicht sinnvoll für Faktoren
## 
## Call:
## lm(formula = overweight ~ anz_zimmer, data = shp_sample_g2)
## 
## Coefficients:
## (Intercept)   anz_zimmer  
##     1.17476     -0.01382
Unterhalb des Medians ist der Zusammenhang positiv, darüber ist er negativ. Dieser Befund zeugt von einem übergreifenden Zusammenhang, der zumindest graduell nicht linear ist. Diese Einschätzung entspricht auch dem im Balkendiagramm oben abgebildeten Verlauf. Wir bleiben zwar bei der linearen bzw. stetigen Modellierung (zumal wir hinter dem positiven Zusammenhang im unteren Bereich der UV einen Alterseffekt vermuten und die angedeutete Wellenförmige Form des Zusammenhangs auch mit der logistischen Regression nicht aufgegriffen würde), sollten aber diesen Befund im Rahmen der Ergebnisdarstellung berichten und als Einschränkung bzw. Limitation diskutieren.

 

5 Modellierung: Lineare Wahrscheinlichkeitsregression

Fallzahlharmonisierung

Bevor wir die multivariate regressionsbasierte Analysereihe etablieren, schränken wir die Stichprobe auf alle Merkmalsträger mit gültigen Messungen in allen aufgegriffenen Variablen ein, um einer Konvention folgend die Fallzahl über die Modellreihe konstant zu halten. Achtung: Diesen Schritt bitte erst dann durchführen, wenn alle nicht modellierten Variablen ausgeschlossen sind und Ihr Euch zudem vergewissert habt, dass keine unnötigen “Stichprobenfresser” (siehe oben) im Modell sind.

shp_model<-drop_na(shp_sample)

Basis- und Nettomodell

Los geht’s mit der regressionsbasierten Analyse des einfachen bivariaten Zusammenhangs, quasi als Brücke von der deskriptiven, bivariaten Visualisierung in die Modellreihe. Da lm offenbar die Formulierung eines linearen Wahrscheinlichkeitsmodells verweigert, stellen wir vorab die abhängige Variable als numerische Variable dar.

shp_model$overweight<-as.numeric(shp_model$overweight)
shp_model$overweight[shp_model$overweight==1]<-0
shp_model$overweight[shp_model$overweight==2]<-1
basis<-lm(overweight~anz_zimmer, shp_model)
summary (basis)
## 
## Call:
## lm(formula = overweight ~ anz_zimmer, data = shp_model)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.13764 -0.11168 -0.10303 -0.09438  0.94023 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.146293   0.010566  13.846  < 2e-16 ***
## anz_zimmer  -0.008652   0.002231  -3.879 0.000106 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3094 on 8468 degrees of freedom
## Multiple R-squared:  0.001774,   Adjusted R-squared:  0.001656 
## F-statistic: 15.05 on 1 and 8468 DF,  p-value: 0.0001057

Ad-hoc Interpretation: Mit jedem zusätzlichen Zimmer in der Wohnung nimmt die Wahrscheinlichkeit für Übergewicht im Mittel um etwa 0.9 Prozentpunkte ab.

Das zentrales Analysemodell, welches die Effektinterpretation anvisiert und somit zur Hypothesenprüfung eingesetzt wird, ist üblicherweise das Nettomodell (manchmal auch bereinigtes Modell, Modell 2 oder, etwas euphemistisch, Kausalmodell.genannt), welches den Effekt der als Störmerkmale identifizierten Variablen durch deren Integration abkoppelt. Die Spezifikationsstrategie des Nettomodells sollte im Methodenabschnitt erläutert und begründet werden (siehe das Beispiel mit metrischer Variable.

options(scipen=2, digits=4)
netto<-lm(overweight~anz_zimmer+geschlecht+erwerb+zivilstand+oecd_inc_mon+alter, shp_model)
summary (netto)
## 
## Call:
## lm(formula = overweight ~ anz_zimmer + geschlecht + erwerb + 
##     zivilstand + oecd_inc_mon + alter, data = shp_model)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.1935 -0.1225 -0.1044 -0.0803  1.0726 
## 
## Coefficients:
##                                                                     Estimate
## (Intercept)                                                       0.14712498
## anz_zimmer                                                       -0.00784054
## geschlechtFemale                                                 -0.00979329
## erwerbPart Time                                                  -0.00963161
## erwerbNot Working                                                -0.00108456
## zivilstandSingle, not Living with a Partner                      -0.03133971
## zivilstandWidowed, not Living with a Partner                      0.02114666
## zivilstandDivorced, not Living with a Partner                    -0.01153485
## zivilstandSeparated (Legally Married), not Living with a Partner  0.01087821
## oecd_inc_mon                                                     -0.00000503
## alter                                                             0.00061713
##                                                                   Std. Error
## (Intercept)                                                       0.01872146
## anz_zimmer                                                        0.00230533
## geschlechtFemale                                                  0.00730355
## erwerbPart Time                                                   0.00917841
## erwerbNot Working                                                 0.00973185
## zivilstandSingle, not Living with a Partner                       0.01037691
## zivilstandWidowed, not Living with a Partner                      0.01562216
## zivilstandDivorced, not Living with a Partner                     0.01401596
## zivilstandSeparated (Legally Married), not Living with a Partner  0.03230363
## oecd_inc_mon                                                      0.00000132
## alter                                                             0.00026494
##                                                                  t value
## (Intercept)                                                         7.86
## anz_zimmer                                                         -3.40
## geschlechtFemale                                                   -1.34
## erwerbPart Time                                                    -1.05
## erwerbNot Working                                                  -0.11
## zivilstandSingle, not Living with a Partner                        -3.02
## zivilstandWidowed, not Living with a Partner                        1.35
## zivilstandDivorced, not Living with a Partner                      -0.82
## zivilstandSeparated (Legally Married), not Living with a Partner    0.34
## oecd_inc_mon                                                       -3.82
## alter                                                               2.33
##                                                                  Pr(>|t|)    
## (Intercept)                                                       4.4e-15 ***
## anz_zimmer                                                        0.00067 ***
## geschlechtFemale                                                  0.17999    
## erwerbPart Time                                                   0.29403    
## erwerbNot Working                                                 0.91127    
## zivilstandSingle, not Living with a Partner                       0.00253 ** 
## zivilstandWidowed, not Living with a Partner                      0.17589    
## zivilstandDivorced, not Living with a Partner                     0.41054    
## zivilstandSeparated (Legally Married), not Living with a Partner  0.73631    
## oecd_inc_mon                                                      0.00013 ***
## alter                                                             0.01987 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.309 on 8459 degrees of freedom
## Multiple R-squared:  0.00869,    Adjusted R-squared:  0.00752 
## F-statistic: 7.41 on 10 and 8459 DF,  p-value: 7.96e-12

Hinweis: Der erste Befehl im Chunk oben sorgt dafür, dass alle Zahlenwerte im Output im Dezimalsystem dargestellt werden und so besser lesbar sind.

Ad-hoc Interpretation: Der bereinigte Koeffizient hat sich gegenüber der bivariaten Regressionsschätzung geringfügig verändert; ein kleiner Teil des negativen bivariaten Zusammenhangs ist offenbar den Effekten der nun als Störmerkmale integrierten Drittvariablen geschuldet. Gleichwohl: auch der bereinigte Koeffizient weisst einen inhaltlich wie statistisch bedeutsamen Einfluss der Wohnungsgrösse auf die Neigung zu Übergewicht aus. Die Nullhypothese kann also abgelehnt werden.

Moderationsmodell

Falls Erklärungsmodell und Hypothese eine Wechselwirkung benennen, spezifizieren wir zusätzlich ein Moderationsmodell. Hier, in unserer Beispielanalyse, integrieren wir einen Interaktionsterm mit der Erwerbsvariablen in den Regressionsbefehl - nehmen also an, dass im theoretischen Abschnitt nicht nur für einen Einfluss der Wohnungsgrösse argumentiert wurde, sondern auch dafür, dass die Wohnungsgrösse für nicht-erwerbstätige Personen einen höhere Bedeutung hat als für Erwerbstätige.

moderation<-lm(overweight~geschlecht+anz_zimmer*erwerb+zivilstand+oecd_inc_mon+alter, shp_model)
summary (moderation)
## 
## Call:
## lm(formula = overweight ~ geschlecht + anz_zimmer * erwerb + 
##     zivilstand + oecd_inc_mon + alter, data = shp_model)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.2071 -0.1216 -0.1039 -0.0789  1.0655 
## 
## Coefficients:
##                                                                     Estimate
## (Intercept)                                                       0.12694520
## geschlechtFemale                                                 -0.00933366
## anz_zimmer                                                       -0.00306550
## erwerbPart Time                                                   0.01181177
## erwerbNot Working                                                 0.04189468
## zivilstandSingle, not Living with a Partner                      -0.03134403
## zivilstandWidowed, not Living with a Partner                      0.01985654
## zivilstandDivorced, not Living with a Partner                    -0.01251327
## zivilstandSeparated (Legally Married), not Living with a Partner  0.01051180
## oecd_inc_mon                                                     -0.00000484
## alter                                                             0.00057160
## anz_zimmer:erwerbPart Time                                       -0.00481855
## anz_zimmer:erwerbNot Working                                     -0.00943758
##                                                                   Std. Error
## (Intercept)                                                       0.02273236
## geschlechtFemale                                                  0.00731278
## anz_zimmer                                                        0.00374499
## erwerbPart Time                                                   0.02794576
## erwerbNot Working                                                 0.02569825
## zivilstandSingle, not Living with a Partner                       0.01038001
## zivilstandWidowed, not Living with a Partner                      0.01564139
## zivilstandDivorced, not Living with a Partner                     0.01402560
## zivilstandSeparated (Legally Married), not Living with a Partner  0.03230186
## oecd_inc_mon                                                      0.00000132
## alter                                                             0.00026633
## anz_zimmer:erwerbPart Time                                        0.00580284
## anz_zimmer:erwerbNot Working                                      0.00522232
##                                                                  t value
## (Intercept)                                                         5.58
## geschlechtFemale                                                   -1.28
## anz_zimmer                                                         -0.82
## erwerbPart Time                                                     0.42
## erwerbNot Working                                                   1.63
## zivilstandSingle, not Living with a Partner                        -3.02
## zivilstandWidowed, not Living with a Partner                        1.27
## zivilstandDivorced, not Living with a Partner                      -0.89
## zivilstandSeparated (Legally Married), not Living with a Partner    0.33
## oecd_inc_mon                                                       -3.67
## alter                                                               2.15
## anz_zimmer:erwerbPart Time                                         -0.83
## anz_zimmer:erwerbNot Working                                       -1.81
##                                                                  Pr(>|t|)    
## (Intercept)                                                       2.4e-08 ***
## geschlechtFemale                                                  0.20187    
## anz_zimmer                                                        0.41306    
## erwerbPart Time                                                   0.67255    
## erwerbNot Working                                                 0.10309    
## zivilstandSingle, not Living with a Partner                       0.00254 ** 
## zivilstandWidowed, not Living with a Partner                      0.20430    
## zivilstandDivorced, not Living with a Partner                     0.37233    
## zivilstandSeparated (Legally Married), not Living with a Partner  0.74487    
## oecd_inc_mon                                                      0.00025 ***
## alter                                                             0.03189 *  
## anz_zimmer:erwerbPart Time                                        0.40635    
## anz_zimmer:erwerbNot Working                                      0.07077 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.308 on 8457 degrees of freedom
## Multiple R-squared:  0.00907,    Adjusted R-squared:  0.00766 
## F-statistic: 6.45 on 12 and 8457 DF,  p-value: 1.47e-11

Der Koeffizient der Produktterme misst hier den Einfluss des Moderators auf den Effekt der unabhängigen Variable. Ad-hoc Interpretation: Für Teilzeitbeschäftigte und Nichterwerbstätige ist der Einfluss der Wohnungsgrösse auf die Neigung zu Übergewicht stärker als für Erwerbstätige: Ein zusätzlicher Raum verringert für eine teilzeitbeschäftigte Person die Übergewichtswahrscheinlichkeit um etwa 0.5 Prozentpunkte stärker, für eine nichterwerbstätige Person sogar um fast einen ganzen Prozentpunkt stärker als für eine erwerbstätige Person. Die Moderationshypothese wird durch diese Ergebnisse folglich gestützt, zumal letzterer Unterschied auch (schwach) statistisch signifikant ist.

Achtung: Wie fast immer bei Spezifikation eines Interaktionsterms mit metrischer unabhängiger Variable kann der Koeffizient des Haupteffektes der Moderatorvariable - hier: Erwerbstätigkeit - nicht sinnvoll interpretiert werden, da dieser sich, wie sich leicht algebraisch zeigen liesse, nun ausschliesslich auf den Nullpunkt der moderierten Variable - also Personen in Wohnungen mit Null Zimmern - bezieht.

 

5.1 Tabellarische Darstellung

Die Regressionstabelle ist obligatorisches Element eines jeden Forschungsbeitrags mit regressionsbasierter Analyse. Insbesondere das package stargazer erleichtert die Formatierung; ganz ohne manuelles Nachhelfen geht es allerdings auch hier nicht.

stargazer(basis, netto, moderation, type="text", 
          title="Lineares Wahrscheinlichkeitsmodell: Determinanten des Übergewichtigkeit",
          notes=c("in Klammern: Standardfehler", 
                  "SHP v19, n=8.470. Eigene Berechnungen"),
          dep.var.labels="",  dep.var.caption = "", 
          column.labels =c("Basis", "Netto", "Moderation"),
          single.row = F,
          digits=3,
          keep.stat=c("n", "rsq"),
          digits.extra=5,
          out="regtab2.htm")
## 
## Lineares Wahrscheinlichkeitsmodell: Determinanten des Übergewichtigkeit
## =======================================================================
##                                                                        
##                                     Basis        Netto      Moderation 
##                                      (1)          (2)          (3)     
## -----------------------------------------------------------------------
## Anzahl Zimmer                     -0.009***    -0.008***      -0.003   
##                                    (0.002)      (0.002)      (0.004)   
##                                                                        
## Weiblich                                         -0.010       -0.009   
##                                                 (0.007)      (0.007)   
##                                                                        
## Teilzeit                                         -0.010       0.012    
##                                                 (0.009)      (0.028)   
##                                                                        
## Nichterwerbstätig                                -0.001       0.042    
##                                                 (0.010)      (0.026)   
##                                                                        
## Single                                         -0.031***    -0.031***  
##                                                 (0.010)      (0.010)   
##                                                                        
## Verwittwet                                       0.021        0.020    
##                                                 (0.016)      (0.016)   
##                                                                        
## Geschieden                                       -0.012       -0.013   
##                                                 (0.014)      (0.014)   
##                                                                        
## Getrennt                                         0.011        0.011    
##                                                 (0.032)      (0.032)   
##                                                                        
## Einkommen                                     -0.00001***  -0.000005***
##                                                (0.000001)   (0.000001) 
##                                                                        
## Alter                                           0.001**      0.001**   
##                                                 (0.0003)     (0.0003)  
##                                                                        
## Teilzeit*Anazhl Zimmer                                        -0.005   
##                                                              (0.006)   
##                                                                        
## Nichterwerbstätig*Anazhl Zimmer                              -0.009*   
##                                                              (0.005)   
##                                                                        
## Konstante                         0.146***      0.147***     0.127***  
##                                    (0.011)      (0.019)      (0.023)   
##                                                                        
## -----------------------------------------------------------------------
## Observations                        8,470        8,470        8,470    
## R2                                  0.002        0.009        0.009    
## =======================================================================
## Note:                                       *p<0.1; **p<0.05; ***p<0.01
##                                             in Klammern: Standardfehler
##                                   SHP v19, n=8.470. Eigene Berechnungen

In der Mustersyntax zur metrischen abhängigen Variable ist die Weiterverarbeitung bzw. notwendige Formatierung der Regressionstabelle beispielhaft dargestellt.

 

5.2 Ausreisserdiagnostik

Der Test deiner Analyseergebnisse auf Robustheit gegenüber Ausreissern ist ein Element der klassischen Regressionsdiagnostik, welches wir für hinreichend relevant halten, um eine Integration dieses Tests in Eure Abschlussarbeit zu empfehlen. Die Ausarbeitung sollte allerdings nicht im Hauptteil, sondern im Anhang der Arbeit dokumentiert werden. Im Hauptteil reicht in der Regel ein kurzer Verweis auf die Robustheit der Ergebnisse gegenüber Ausreisserausschluss; bei fehlender Robustheit sollte dieser Befund als Limitation in der Diskussion kritisch diskutiert und eingeordnet werden.

Die Ausreisserbestimmung beruht auf Berechnung des standardisierten Einfluss einzelner Messungen auf den Regressionskoeffizienten, genannt DFbeta. Zunächst legen wir den Grenzwert fest, ab dem wir eine Messung als “einflussreich” bezeichnen wollen. Wir orientieren uns an dem von Belsley et al. (1980) vorgeschlagene Schwellenwert (2/Wurzel(n), hier |0,021|). Anschliessend ermitteln wir, ob es nach dieser Definition einflussreiche Fälle auf den Einkommenskoeffizienten gibt.

summary(dfbetas(netto))
##   (Intercept)         anz_zimmer       geschlechtFemale   erwerbPart Time   
##  Min.   :-0.13475   Min.   :-0.13196   Min.   :-0.06969   Min.   :-0.06510  
##  1st Qu.:-0.00314   1st Qu.:-0.00229   1st Qu.:-0.00330   1st Qu.:-0.00345  
##  Median :-0.00046   Median : 0.00001   Median :-0.00105   Median : 0.00095  
##  Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.00000   Mean   : 0.00000  
##  3rd Qu.: 0.00196   3rd Qu.: 0.00231   3rd Qu.: 0.00351   3rd Qu.: 0.00306  
##  Max.   : 0.10375   Max.   : 0.12147   Max.   : 0.05805   Max.   : 0.09749  
##  erwerbNot Working  zivilstandSingle, not Living with a Partner
##  Min.   :-0.08104   Min.   :-0.06713                           
##  1st Qu.:-0.00311   1st Qu.:-0.00176                           
##  Median : 0.00029   Median : 0.00061                           
##  Mean   : 0.00000   Mean   : 0.00000                           
##  3rd Qu.: 0.00289   3rd Qu.: 0.00229                           
##  Max.   : 0.11351   Max.   : 0.10592                           
##  zivilstandWidowed, not Living with a Partner
##  Min.   :-0.03981                            
##  1st Qu.:-0.00050                            
##  Median : 0.00045                            
##  Mean   : 0.00000                            
##  3rd Qu.: 0.00176                            
##  Max.   : 0.14922                            
##  zivilstandDivorced, not Living with a Partner
##  Min.   :-0.02471                             
##  1st Qu.:-0.00018                             
##  Median : 0.00089                             
##  Mean   : 0.00000                             
##  3rd Qu.: 0.00164                             
##  Max.   : 0.12914                             
##  zivilstandSeparated (Legally Married), not Living with a Partner
##  Min.   :-0.05827                                                
##  1st Qu.: 0.00002                                                
##  Median : 0.00041                                                
##  Mean   : 0.00000                                                
##  3rd Qu.: 0.00068                                                
##  Max.   : 0.29723                                                
##   oecd_inc_mon         alter         
##  Min.   :-0.0552   Min.   :-0.09435  
##  1st Qu.:-0.0017   1st Qu.:-0.00275  
##  Median : 0.0003   Median :-0.00010  
##  Mean   : 0.0000   Mean   : 0.00000  
##  3rd Qu.: 0.0024   3rd Qu.: 0.00237  
##  Max.   : 0.5585   Max.   : 0.10411

Die Übersichtsstatistik zu den Einflusswerten zeigt, dass es Personen gibt, deren Einfluss auf den Regressionskoeffizienten oberhalb des Grenzwertes liegt. Wir schliessen diese Personen nun für den Robustheitstest aus der Analyse aus, um zu prüfen, ob das Regressionsergebnis durch diese einflussreiche Fälle unbotmässig stark geprägt wird.

dfbetas_matrix <- data.frame(betas=dfbetas(netto))
shp_model_betas <- cbind(shp_model, dfbetas_matrix)
shp_model_oA <- filter(shp_model_betas, betas.anz_zimmer<0.021 & betas.anz_zimmer>-0.021 )
netto2<-lm(overweight~anz_zimmer+geschlecht+erwerb+zivilstand+oecd_inc_mon+alter, shp_model_oA)
stargazer (netto, netto2, type="text", 
           column.labels =c("Originalnetto", "ohne Ausreisser"))
## 
## =====================================================================
##                                    Dependent variable:               
##                     -------------------------------------------------
##                                        overweight                    
##                          Originalnetto           ohne Ausreisser     
##                               (1)                      (2)           
## ---------------------------------------------------------------------
## Anzahl Zimmer              -0.008***                -0.008***        
##                             (0.002)                  (0.002)         
##                                                                      
## Weiblich                     -0.010                   0.001          
##                             (0.007)                  (0.006)         
##                                                                      
## Teilzeit                     -0.010                   -0.012         
##                             (0.009)                  (0.008)         
##                                                                      
## Nichterwerbstätig            -0.001                   -0.005         
##                             (0.010)                  (0.008)         
##                                                                      
## Single                     -0.031***                -0.027***        
##                             (0.010)                  (0.009)         
##                                                                      
## Verwittwet                   0.021                    0.011          
##                             (0.016)                  (0.013)         
##                                                                      
## Geschieden                   -0.012                  -0.0005         
##                             (0.014)                  (0.012)         
##                                                                      
## Getrennt                     0.011                    0.033          
##                             (0.032)                  (0.027)         
##                                                                      
## Einkommen                 -0.00001***              -0.00000***       
##                            (0.00000)                (0.00000)        
##                                                                      
## Alter                       0.001**                  0.001**         
##                             (0.0003)                 (0.0002)        
##                                                                      
## Konstante                   0.147***                 0.103***        
##                             (0.019)                  (0.016)         
##                                                                      
## ---------------------------------------------------------------------
## Observations                 8,470                    8,112          
## R2                           0.009                    0.009          
## Adjusted R2                  0.008                    0.008          
## Residual Std. Error    0.309 (df = 8459)        0.251 (df = 8101)    
## F Statistic         7.414*** (df = 10; 8459) 7.504*** (df = 10; 8101)
## =====================================================================
## Note:                                     *p<0.1; **p<0.05; ***p<0.01

Der zentrale, hypothesenbelastende Koeffizient verändert sich nicht bei Ausschluss einflussreicher Fälle. Ergo: Das Originalergebnis ist robust gegenüber einem Ausreisserausschluss nach Belsley et al (1980).

 

5.3 Gewichtete Ergebnisse

Das Schweizer Haushaltspanel ist eine grundsätzlich repräsentativ ausgerichtete Studie. Allerdings weichen die rohen Merkmalsverteilungen in den Daten von der Repräsentativität ab; bestimmte Gruppen sind hier über- bzw. unterrepräsentiert. Der Grad an Über- bzw. Unterrepräsentation bestimmter Merkmalsträger wird in den Daten über die bereitgestellten Gewichtungsvariablen aufgegriffen. Werden diese Gewichtungsvariablen als Gewichtungsanweisung in die Analysekommandos integriert, wird der Einfluss einzelner Merkmalsträger auf die Analyseergebnisse entsprechend seines Repräsentationsfaktors modifiziert. Gewichtete Ergebnisse nicht-repräsentativer Daten bilden also Ergebnisse repräsentativer Daten nach.

Zusammenhänge zwischen Variablen sind, anders als Anteilswerte, in der Regel robust gegenüber Abweichungen von der Repräsentativität. Da die in unseren Forschungsberichten untersuchen Fragen in der Regel auf die Analyse von Zusammenhängen ausgerichtet ist, halten wir die Verwendung von Gewichten in Euren Analyse daher nicht für essentiell. Wir stellen Euch folglich frei, in Euren Analysen Gewichte zu verwenden oder im Appendix zusätzlich gewichtete Ergebnisse als Robustheitstest zu präsentieren (und in einer Fussnote kurz zu referenzieren). Auch in unserem Beispiel liegen die Ergebnisse einer zusätzlichen gewichteten Analyse sehr eng an den ungewichteten Ergebnissen:

shp_model_weight<-filter(shp_model, w11101_2017>0)
netto3<-lm(overweight~anz_zimmer+geschlecht+erwerb+zivilstand+oecd_inc_mon+alter, shp_model_weight, weights=w11101_2017)
stargazer (netto, netto3, type="text", column.labels =c("Originalnetto", "Gewichtete Ergebnisse"))
## 
## ======================================================================
##                                    Dependent variable:                
##                     --------------------------------------------------
##                                         overweight                    
##                          Originalnetto         Gewichtete Ergebnisse  
##                               (1)                       (2)           
## ----------------------------------------------------------------------
## Anzahl Zimmer              -0.008***                 -0.009***        
##                             (0.002)                   (0.002)         
##                                                                       
## Weiblich                     -0.010                  -0.019**         
##                             (0.007)                   (0.008)         
##                                                                       
## Teilzeit                     -0.010                   -0.008          
##                             (0.009)                   (0.010)         
##                                                                       
## Nichterwerbstätig            -0.001                   0.021**         
##                             (0.010)                   (0.010)         
##                                                                       
## Single                     -0.031***                 -0.035***        
##                             (0.010)                   (0.010)         
##                                                                       
## Verwittwet                   0.021                     0.015          
##                             (0.016)                   (0.018)         
##                                                                       
## Geschieden                   -0.012                   -0.002          
##                             (0.014)                   (0.015)         
##                                                                       
## Getrennt                     0.011                     0.013          
##                             (0.032)                   (0.039)         
##                                                                       
## Einkommen                 -0.00001***               -0.00001***       
##                            (0.00000)                 (0.00000)        
##                                                                       
## Alter                       0.001**                   0.0004          
##                             (0.0003)                 (0.0003)         
##                                                                       
## Konstante                   0.147***                 0.182***         
##                             (0.019)                   (0.019)         
##                                                                       
## ----------------------------------------------------------------------
## Observations                 8,470                     8,299          
## R2                           0.009                     0.013          
## Adjusted R2                  0.008                     0.012          
## Residual Std. Error    0.309 (df = 8459)         8.893 (df = 8288)    
## F Statistic         7.414*** (df = 10; 8459) 10.790*** (df = 10; 8288)
## ======================================================================
## Note:                                      *p<0.1; **p<0.05; ***p<0.01

 

5.4 Visualisierung: Haupteffekt

Der Conditional Effect Plot bildet den Verlauf der bereinigten Regressionsfunktion (aus dem Nettomodell) ab - empfehlen wir für jeden Forschungsbeitrag mit regressionsbasierter Analyse. Auch hier gibt es unterschiedliche Workflows; Package unserer Wahl ist derzeit visreg, welches den in diesem Beispiel wichtigen Vorteil hat, dass es per Default die logarithmische Kurve über der linearen Skala (statt der linearen Gerade über der logarithmierten Skala) darstellt.

visreg hält bei der Darstellung der Vorhersagefunktion automatisch die Kovariaten in den Mittelwerten konstant (daher auch conditional effect plot), wodurch der bereinigte Effekt der “durchschnittlichen” Person abgebildet wird. Zudem kann bei Wahl von gg=TRUE die gewohnte ggplot-Funktionalität für die Formatierung der Abbildung verwendet werden:

vis1<-visreg(netto, "anz_zimmer", partial=F, rug=F, gg=T) 
vis1

Die Charakteristika des Zusammenhangs und somit die praktischen Implikationen des ermittelten Koeffizienten werden so schön deutlich. Anreichern liesse sich die Darstellung durch Ableitung bzw. Verdeutlichung charakteristischer vorhergesagter Werte:

predict(netto, data.frame(anz_zimmer=2, oecd_inc_mon=5000, zivilstand="Married or Living with a Partner", alter=35, geschlecht="Female", erwerb="Full Time")) 
##      1 
## 0.1181
predict(netto, data.frame(anz_zimmer=4, oecd_inc_mon=5000, zivilstand="Married or Living with a Partner", alter=35, geschlecht="Female", erwerb="Full Time"))  
##      1 
## 0.1024
predict(netto, data.frame(anz_zimmer=6, oecd_inc_mon=5000, zivilstand="Married or Living with a Partner", alter=35, geschlecht="Female", erwerb="Full Time"))  
##       1 
## 0.08673

Für eine 35-jährige, vollzeiterwerbstätige und verheiratete Frau mit Haushaltseinkommen von 5000 CHF sinkt die Wahrscheinlich für Übergewichtigkeit von etwa 11,8% auf 8,6%, wenn sich die Wohnungsgrösse von 2 auf 4 Räume erhöht. Solche eine inhaltliche Interpretationsstrategie, die konditionale Vorhersagen auf Basis der Regressionsschätzung vergleichend einordnet, funktioniert häufig sehr gut im Rahmen von Forschungsberichten oder Präsentationen zur Darstellung der substanziellen Grösse des Koeffizienten bzw. des Effect-Size.

Als nächstes trimmen wir die durch visreg erstellte Grafik zur publikationsreife:

vis1+labs(title="Regressionsfunktion mit 95%-CI: Wohnungsgrösse und Übergewichtigkeit",
       caption="Vorhergesagte Werte unter Konstanthaltung von Standarddemographie, Einkommen und Erwerbstätigkeit (Modell 'Netto') \nSHP v19, n=8470")+
       ylab("p(Übergewichtigkeit)")+
  xlab("Anzahl Zimmer")+
  scale_x_continuous(breaks=seq(1,10, 1))

Gutes Beispiel: # https://www.science.org/doi/full/10.1126/science.aay0214?adobe_mc=MCMID%3D58848560922295680121103752198191409082%7CMCORGID%3D242B6472541199F70A4C98A6%2540AdobeOrg%7CTS%3D1656657585

 

5.5 Visualisierung: Interaktion

visreg eignet sich auch zur Visualisierung von Regressionsanalysen mit Interaktionsterm:

vis2<-visreg(moderation, "anz_zimmer", by="erwerb", overlay=T, partial=F, rug=F, gg=T) 
vis2

Auch hier gilt es, die Darstellung durch Spezifikation des Befehls zu verbessern:

vis2+
  scale_color_manual(values=c("#00C1C9", "#D63EFF", "#FF4E37")) +
  scale_fill_manual(values=c("#00C1C910", "#D63EFF10", "#FF4E3710"))+
  labs(title="Regressionsfunktion: Wohnungsgrösse und Übergewichtigkeit nach Erwerbsstatus",
       caption="Vorhergesagte Werte unter Konstanthaltung der Kontrollvariablen (Modell ´moderation´)\nSHP v19, n=8470")+
       ylab("p(Übergewichtigkeit)")+
  xlab("Anzahl Zimmer")+
  scale_x_continuous(breaks=seq(1,10, 1))

… noch nicht perfekt (Legende nicht gut formatiert, Titel abgeschnitten etc), aber eine gute Basis um weiterzuarbeiten. Möglicherweise hilfreiche Anregungen zum Troubleshooting von visreg findet ihr in dem Skript zur Analyse eines metrischen abhängigen Merkmals.

 

5.8 Interpretation

Die tabellarische und grafische Aufbereitung stützt die Zugänglichkeit der Ergebnisse, ersetzt aber nicht die prosaische Interpretation. Es reicht dabei nicht aus, die statistische Signifikanz des Ergebnisses zu berichten und an die Hypothese rückzubinden: die inhaltliche Interpretation muss vielmehr den Koeffizienten bzw. seine Grösse einordnen und damit greifbar machen. Einige Anhaltspunkte findet ihr oben in den jeweiligen Erläuterungen und ad-hoc Interpretationen, aber grundsätzlich stellt jede Analyse ihre eigenen Anforderungen, so dass es für die Interpretation kein pauschales Rezept gibt. Vielmehr sind Geschick und statistische Intuition gefragt. Ein paar Tipps dazu:

  • Kernaufgabe ist die Auswertung des Koeffizienten der unabhängigen Variable aus dem Nettomodell, da dieser die zentrale empirische Aussage der Studie enthält und ausserdem der Belastung der Hypothese dient
  • Die Standardinterpretation des Regressionskoeffizienten wird oftmals durch Reskalierung der Bezugseinheit (Interpretation z.B. in 10er oder 1000er Einheiten der UV) deutlich zugänglicher - z.B. bei Einkommens- oder Altersvariablen
  • Grade bei abstrakten Skaleneinheiten ist es häufig hilfreich, in der Interpretation auch mit Standardisierungen zu arbeiten (siehe Beispiele aus unserer Statistik 2 - Veranstaltung)
  • Effekteigenschaften lassen sich, wie oben demonstriert, oft auch gut durch geeignete Vorhersagepaare verdeutlichen.
  • Versucht auch, die Veränderung des zentralen Koeffizienten über die Modelle zu erklären
  • Bitte nur am Rande: Auswertung der Koeffizienten der Kontrollvariablen.

     

6 Alternative Modellierung: Logistische Regression

Schliesslich widmen wir uns der alternativen, logistischen Modellierungsoption.

logit_basis<-glm(data=shp_model, family=binomial(link = "logit"), overweight~anz_zimmer)
logit_netto<-glm(data=shp_model, family=binomial(link = "logit"), overweight~anz_zimmer+geschlecht+erwerb+zivilstand+oecd_inc_mon+alter)
logit_moderation<-glm(data=shp_model, family=binomial(link = "logit"), overweight~geschlecht+anz_zimmer*erwerb+erwerb+zivilstand+oecd_inc_mon+alter)

stargazer(logit_basis, logit_netto, logit_moderation, type="text", 
          title="Logistisches Regressionsmodell: Determinanten des Übergewichtigkeit",
          notes=c("in Klammern: Standardfehler", 
                  "SHP v19, n=8.470. Eigene Berechnungen"),
          dep.var.labels="",  dep.var.caption = "", 
          column.labels =c("Basis", "Netto", "Moderation"),
          single.row = F,
          digits=3,
          keep.stat=c("n", "rsq"),
          digits.extra=5,
          out="regtab2.htm")
## 
## Logistisches Regressionsmodell: Determinanten des Übergewichtigkeit (Odds Ratios)
## ==============================================================
##                                                               
##                            Basis        Netto      Moderation 
##                             (1)          (2)          (3)     
## --------------------------------------------------------------
## Anzahl Zimmer            -0.093***    -0.087***      -0.037   
##                           (0.024)      (0.026)      (0.041)   
##                                                               
## Weiblich                                -0.098       -0.093   
##                                        (0.077)      (0.077)   
##                                                               
## Teilzeit                                -0.119       0.159    
##                                        (0.101)      (0.311)   
##                                                               
## Nichterwerbstätig                       -0.051       0.342    
##                                        (0.106)      (0.273)   
##                                                               
## Single                                -0.394***    -0.393***  
##                                        (0.119)      (0.119)   
##                                                               
## Verwittwet                              0.143        0.132    
##                                        (0.146)      (0.146)   
##                                                               
## Geschieden                              -0.138       -0.147   
##                                        (0.145)      (0.145)   
##                                                               
## Getrennt                                0.084        0.080    
##                                        (0.314)      (0.314)   
##                                                               
## Einkommen                             -0.0001***   -0.0001*** 
##                                       (0.00002)    (0.00002)  
##                                                               
## Alter                                  0.007**      0.007**   
##                                        (0.003)      (0.003)   
##                                                               
## Teilzeit*Anazhl Zimmer                               -0.063   
##                                                     (0.066)   
##                                                               
## No Job*Anazhl Zimmer                                 -0.089   
##                                                     (0.057)   
##                                                               
## Konstante                -1.706***    -1.658***    -1.871***  
##                           (0.110)      (0.202)      (0.247)   
##                                                               
## --------------------------------------------------------------
## Observations               8,470        8,470        8,470    
## ==============================================================
## Note:                              *p<0.1; **p<0.05; ***p<0.01
##                                    in Klammern: Standardfehler
##                          SHP v19, n=8.470. Eigene Berechnungen

Ausgewiesen wurden hier die Originalkoeffizienten der logistischen Regression, die sich hinsichtlich der Wahrscheinlichkeitsveränderung lediglich in ihrer Tendenz interpretieren lassen:

Alternativ können die Koeffizienten exponiert und dann als multiplikative Chancenveränderung (“Odds Rations”) interpretiert werden - obgleich der tatsächliche Gewinn an Aussagekraft Gegenstand kritischer Diskussionen ist: Einerseits verheisst das Konzept der Chancenveränderung semantische Zugänglichkeit, andererseits ist der substanzielle Gehalt einer konkreten Chancenveränderung nicht besonders gut greifbar. Es steht zu befürchten, dass der Rückgriff auf den in der Umgangssprache gebräuchlichen, aber in der Statistik mit komplexer Definition belegten Begriff der “Chance” eine Scheinzugänglichkeit eigentlich nicht intuitiv prozessierbarer statistischer Grössen herstellt.

Stagazer bietet die Option, alle Koeffizienten direkt exponiert zur Basis e auszuweisen. Allerdings zerschiesst diese Option die Teststatistik (p-Werte, Sternchen…), welche dann nicht mehr interpretierbar, bzw. sogar ungültig ist. Daher empfehlen wir die Darlegung originaler Koeffizienten und deren manuelle Exponierung im Rahmen der Auswertung, wenn dezidiertes Interesse an der Interpretation als multiplikative Chancenveränderung besteht:

Interpretationen, die auf den Koeffizienten der logistischen Regression basieren, sind grundsätzlich sperrig und schwer greifbar, wie diese Interpretationsversuche verdeutlichen. Daher ist eine gute Visualisierung der Ergebnisse als Wahrscheinlichkeitsfunktion aus unserer Sicht unerlässlich:

vis3<-visreg(logit_netto, "anz_zimmer", scale="response", partial=F, rug=F, gg=T) 
vis3

vis3+labs(title="Logistische Regressionsvorhersage mit 95%-CI",
       caption="Vorhergesagte Werte unter Konstanthaltung von Standarddemographie, Einkommen und Erwerbstätigkeit (Modell 'Netto') \nSHP v19, n=8470")+
       ylab("p(Übergewichtigkeit)")+
         ylab("p(Übergewichtigkeit)")+
  xlab("Anzahl Zimmer")+
  scale_x_continuous(breaks=seq(1,10, 1))

Der näherungsweise lineare Verlauf der logistischen Vorhersagekurve und die Ähnlichkeit der Vorhersagen der beiden Modellklassen (siehe unten) liefern in diesem Fall eine gute a posteriori Begründung dafür, die Ergebnisse einer linearen Schätzung zu berichten und zu interpretieren. In einer Fussnote könnte dann auf die Ähnlichkeit der Vorhersagen und grundsätzliche Robustheit der Ergebnis im Vergleich mit einer (im Anhang ggf. dokumentierten) logistischen Regressionsanalyse verwiesen werden.

Für die Interpretation logistischer Regressionsergebnisse gelten grundsätzlich die gleichen Standards wie für lineare Regressionen (siehe oben). Allerdings sind hier viele der üblichen Ansätze (Reskalierung, Standardisierung) nicht umsetzbar bzw. nicht sinnvoll. So kommt insbesondere der grafisch gestützten vorhersagebasierten Interpretation in der Regel eine wichtige Rolle zu. Vorhergesagte Werte lassen sich wie gehabt mit dem predict-Befehl - allerdings mit dem Zusatz response, damit analog zur visuellen Darstellung Wahrscheinlichkeiten ausgegeben werden - erzeugen:

predict(logit_netto, data.frame(anz_zimmer=2, oecd_inc_mon=5000, zivilstand="Married or Living with a Partner", alter=35, geschlecht="Female", erwerb="Full Time"), type=c("response")) 
##      1 
## 0.1185
predict(logit_netto, data.frame(anz_zimmer=4, oecd_inc_mon=5000, zivilstand="Married or Living with a Partner", alter=35, geschlecht="Female", erwerb="Full Time"), type=c("response"))   
##      1 
## 0.1014
predict(logit_netto, data.frame(anz_zimmer=6, oecd_inc_mon=5000, zivilstand="Married or Living with a Partner", alter=35, geschlecht="Female", erwerb="Full Time"), type=c("response"))   
##       1 
## 0.08658

Diese Vorhersagewerte können nun sehr gut die visuelle Darstellung der Regressionsvorhersage flankieren (bzw. umgekehrt: die visuelle Darstellung kann die Vorhersagewerte flankieren): Für eine 35-jährige, vollzeiterwerbstätige und verheiratete Frau mit einem Haushaltseinkommen von 5000 CHF sinkt die Wahrscheinlich für Übergewichtigkeit von etwa 11,9% auf 8,7%, wenn sich die Wohnungsgrösse von 2 auf 4 Räume erhöht.

Der Vollständigkeit halber nun noch die Visualisierung des logistischen Moderationsmodells:

vis4<-visreg(logit_moderation, "anz_zimmer", scale="response", by="erwerb", overlay=T, partial=F, rug=F, gg=T) 
vis4

Auch hier gilt, es, die Darstellung durch weitere Spezifikation des Befehls zu verbessern:

vis4+
  scale_color_manual(values=c("#00C1C9", "#D63EFF", "#FF4E37")) +
  scale_fill_manual(values=c("#00C1C910", "#D63EFF10", "#FF4E3710"))+
  labs(title="Logistische Regressionsvorhersage mit 95%-CI",
       caption="Vorhergesagte Werte unter Konstanthaltung von Kontrollvariablen (Modell 'Moderation') \nSHP v19, n=8470")+
       ylab("p(Übergewichtigkeit)")+
  xlab("Anzahl Zimmer")+
  scale_x_continuous(breaks=seq(1,10, 1))

Anregungen zum Troubleshooting von visreg findet ihr im Skript zur Analyse eines metrischen abhängigen Merkmals.

 

logo.utf8