Übung zu Drittvariablen

Dies ist das Dokument mit Lösungen

Klick hier, um die Musterlösungen auszublenden.

Zur Vorbereitung

ESS8-Datensatz einlesen.

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

Ein neues Übungsskript erstellen.

# Statistik 2: R Tutorat
# Übungsskript zum Thema Multivariate Regression und Drittvariablen
# Datum: 25.05.2021
# AutorIn: XXX

Installiert die Packages dplyr, ggplot2, rgl, visreg, ggpubr, table1, stargazer und aktiviert sie. Ladet den ESS-Datensatz ein.

# install.packages("dplyr")
library(dplyr)
# install.packages("ggplot2")
library(ggplot2)
# install.packages("rgl")
library(rgl)
# install.packages("visreg")
library(visreg)
# install.packages("ggpubr")
library(ggpubr)
# install.packages("table1")
library(table1)
# install.packages("stargazer")
library(stargazer)

I. Datenaufbereitung & Inspektion

1. In dieser Übung analysierst Du den Einfluss von Bildung auf Religiösität. Formuliere eine Hypothese dazu und begründe diese kurz.

Hypothese (H1)
Bildung hat einen negativen Einfluss auf den Grad an Religiösität einer Person

Begründung (Beispiel)
Ein Ziel des Bildungssystems ist die Förderung kritischen Denkens. Je länger jemand im Bildungssystem verbleibt, desto stärker wird folglich die kritische Auseinandersetzung mit Regelsystemen und Autoritäten, wie der Religion.

2. Begründe das Geburtsjahr und die elterliche Bildung als Störmerkmale des Zusammenhangs von Bildung und Religiösität

Geburtsjahr
Ältere Geburtsjahrgänge sind stärker durch religiöse Institutionen und Konventionen geprägt als jüngere (z->y). Gleichzeitig ist das Bildungsniveau in älteren Generationen aufgrund der fortschreitenden Bildungsexpansion in modernen Gesellschaften tiefer als in jüngeren (z->x). Ein negativer Zusammenhang zwischen Bildung und Religiösität kann somit lediglich Ausdruck eines Generationenunterschieds in der Religiösität sein (der aufgrund der Bildungsungleicheiten zwischen Kohorten in bildungsspezifischen Religionsunterschieden transportiert wird).

Elterliche Bildung
Die elterliche Bildung kann als Störmerkmal betrachtet werden, da zum einen die eigene Bildung von der sozialen Herkunft beeinflusst wird (z->x). So lässt sich im Sinne Bourdieus argumentieren, dass besser gebildete Eltern ihren Kindern mehr kulturelles Kapital zuteil werden lassen, wodurch sich deren Bildungschancen erhöhen. Zum anderen befördert der elterliche Bildungskontext durch den intergenerationalen Wertetransfer direkt die kritische Auseinandersetzung mit Religionen und Autoritäten (z->y). Dass hochgebildete Personen weniger religiös sind, mag also gar nicht ursächlich an deren Bildungsniveau liegen, sondern am typischerweise gehobenen Bildungsumfeld in der Herkunftsfamilie.

3. Warum ist das individuelle Erwerbseinkommen für diesen Zusammenhang kein Störmerkmal?

Eine Variable ist als Störmerkmal qualifiziert, wenn sie sich theoretisch-logisch als Einflussfaktor sowohl der UV als auch der AV herleiten und begründen lässt. Dieses ist für den Lohn nicht der Fall, da dieser keinen gerichteten Einfluss auf das Bildungsniveau ausübt. Zwar liesse sich argumentieren, dass ein hohes Wohlstandniveau, als typische Folge hoher Bildung, den Bedarf nach spirituellem Trost eher drosselt (z->y) - dieser Prozess wäre aber (wegen x->z) nicht störend, sondern konsistent zur Kausalinterpretation des Zusammenhangs zwischen Bildung und Religiösität. Er beschreibt nämlich, eben im Einklang mit der kausal formulierten Hypothese, einen kausalen Mechanismus, über den die Bildung auf die Religiösität wirkt.

4. Wir untersuchen die Hypothese auf Basis der Schweiz, operationalisieren Generationenzugehörigkeit über das Alter, Bildung über die Bildungsjahre, elterliche Bildung über den Bildungsabschluss des Vaters und den Grad an Religiösität über eine subjektive Selbsteinschätzung. Erstelle folglich einen neuen Teildatensatz namens “ess8_CH”, in dem:

Nur SchweizerInnen inkludiert sind
Ausschliesslich die Variablen *idno, eduyrs, rlgdgr, agea, eiscedf** enthalten sind

Wie viele Merkmalsträger sind in diesem neuen Datensatz?

ess8_CH <- filter(ess8, cntry == "CH")
ess8_CH <- select(ess8_CH, eduyrs, agea, rlgdgr, eiscedf)
dim (ess8_CH)

## [1] 1525    4

Es gibt noch 1525 Merkmalsträger.

5. Überprüfe, ob die Variablen rlgdgr und eiscedf richtig kodiert sind

attributes(ess8_CH$rlgdgr)

## $label
## [1] "How religious are you"
## 
## $format.stata
## [1] "%20.0g"
## 
## $labels
## Not at all religious                    1                    2 
##                    0                    1                    2 
##                    3                    4                    5 
##                    3                    4                    5 
##                    6                    7                    8 
##                    6                    7                    8 
##                    9       Very religious              Refusal 
##                    9                   10                   NA 
##           Don't know            No answer 
##                   NA                   NA 
## 
## $class
## [1] "haven_labelled" "vctrs_vctr"     "double"

attributes(ess8_CH$eiscedf)

## $label
## [1] "Father's highest level of education, ES - ISCED"
## 
## $format.stata
## [1] "%51.0g"
## 
## $labels
##             Not possible to harmonise into ES-ISCED 
##                                                   0 
##              ES-ISCED I , less than lower secondary 
##                                                   1 
##                        ES-ISCED II, lower secondary 
##                                                   2 
##           ES-ISCED IIIb, lower tier upper secondary 
##                                                   3 
##           ES-ISCED IIIa, upper tier upper secondary 
##                                                   4 
##        ES-ISCED IV, advanced vocational, sub-degree 
##                                                   5 
##     ES-ISCED V1, lower tertiary education, BA level 
##                                                   6 
## ES-ISCED V2, higher tertiary education, >= MA level 
##                                                   7 
##                                               Other 
##                                                  55 
##                                             Refusal 
##                                                  NA 
##                                          Don't know 
##                                                  NA 
##                                           No answer 
##                                                  NA 
## 
## $class
## [1] "haven_labelled" "vctrs_vctr"     "double"

attributes(ess8_CH$eduyrs)

## $label
## [1] "Years of full-time education completed"
## 
## $format.stata
## [1] "%12.0g"
## 
## $labels
##    Refusal Don't know  No answer 
##         NA         NA         NA 
## 
## $class
## [1] "haven_labelled" "vctrs_vctr"     "double"

Wir stellen fest, dass bei der Variable eiscedf die Ausprägung ‘Other’(=55) bestimmte fehlende Angaben repräsentiert.

5.1 Die väterliche Bildung liegt nicht über Bildungsjahre operationalisiert in den Daten vor: bei der Bildungsmessung nach “ISCED” handelt es sich um eine ordinale (also nicht-metrische) Variable, welche wir nicht einfach in ein lineares Regressionsmodell integrieren können, denn: die Abstände zwischen den einzelnen Messpunkte der ISCED-Skala sind nicht einheitlich, genau dies setzt allerdings die Berechnung und Interpretation eines Regressionskoeffizienten voraus. Wir müssen die Variable daher vorab metrisieren bzw. die einzelnen Ausprägungen in einer neuen Variable nach Bildungsjahren rekodieren.
Erweitere also den ess8_CH Datensatz um eine Variable für väterliche Bildung, in welcher die ISCED Werte nach Bildungsjahren umkodiert sind. Die Umwandlung kann in Anlehnung an den Umrechnungsschlüssel des bfs erfolgen: https://www.bfs.admin.ch/bfs/de/home/statistiken/bildung-wissenschaft/bildungssystem.html):

table(ess8_CH$eiscedf) # Ausprägungen vor Umkodierung

## 
##   1   2   3   4   5   6   7  55 
## 158 278 580  62 168  33 141   5

ess8_CH$eiscedf <- as.numeric(ess8_CH$eiscedf)  # Umwandlung in numeric, da sonst recode() nicht funktioniert
ess8_CH$eduyrsf <- recode(ess8_CH$eiscedf,
                          '1'=8, 
                          '2'=11, 
                          '3'=15, 
                          '4'=15, 
                          '5'=15, 
                          '6'=18, 
                          '7'=20)

table(ess8_CH$eduyrsf) # Ausprägungen nach Umkodierung

## 
##   8  11  15  18  20  55 
## 158 278 810  33 141   5

Für die metrisierte Variable müssen wir die Ausprägung 55 in ein NA umwandeln. Schliesslich vergeben wir noch ein Label für die Variable, damit sie auch in den Rohoutputs gut lesbar ist:

ess8_CH$eduyrsf[ess8_CH$eduyrsf==55] <- NA
library(labelled)
var_label(ess8_CH$eduyrsf) <- "Father years of full-time education"

6. Erstelle nun mit table1() die Stichprobenübersicht zu den univariaten Statistiken der vier Variablen

table1(~eduyrs+agea+rlgdgr+eduyrsf, data = ess8_CH)

	Overall (N=1525)
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%)
How religious are you
Mean (SD)	4.79 (3.03)
Median [Min, Max]	5.00 [0, 10.0]
Missing	14 (0.9%)
Father years of full-time education
Mean (SD)	14.0 (3.21)
Median [Min, Max]	15.0 [8.00, 20.0]
Missing	105 (6.9%)

7. Stichprobenharmonisierung: Erstelle nun einen neuen Datensatz “ess8_noNA”, aus dem alle Personen mit NA in einer der analyserelevanten Variablen ausgeschlossen sind

ess8_noNA <- na.omit(ess8_CH)
dim (ess8_noNA)

## [1] 1405    5

II. Bivariate Analyse

1. Analysiere den Zusammenhang zwischen Bildung und Religiösität auf Basis einer bivariaten Regressionsanalyse. Weise dem Modell den Namen “bi_model” zu. Interpretiere den Regressionskoeffizienten.

bi_model <- lm(rlgdgr~eduyrs, data = ess8_noNA )
summary(bi_model)

## 
## Call:
## lm(formula = rlgdgr ~ eduyrs, data = ess8_noNA)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -5.3006 -2.4145  0.0932  2.2901  6.2748 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  5.89141    0.27403  21.499  < 2e-16 ***
## eduyrs      -0.09846    0.02304  -4.274 2.05e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.009 on 1403 degrees of freedom
## Multiple R-squared:  0.01285,    Adjusted R-squared:  0.01215 
## F-statistic: 18.27 on 1 and 1403 DF,  p-value: 2.048e-05

summary(ess8_noNA $rlgdgr)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   2.000   5.000   4.772   7.000  10.000

sd(ess8_noNA $rlgdgr)

## [1] 3.027611

Mit jedem zusätzlichen Bildungsjahr nimmt die Religiosität im Schnitt um einen zehntel Skalenpunkt ab, dies entspricht etwa einem Dreissigstel des typischen Unterschieds in der Religiösität zwischen zwei Personen. Um zumindest eine halbe Standardabweichung im Grad der Religiösität zu überspringen, braucht es im Schnitt folglich etwa 15 zusätzliche Bildungsjahre - diese entsprechen etwas mehr als dem Unterschied zwischen einem Hauptschulsabschluss und einer Promotion. Bei der Verbindung von Bildung und Religiösität handelt sich folglich um keinen besonders grossen, aber auch nicht um einen marginalen Zusammenhang.

2. Visualisiere den bivariaten Zusammenhang. Achte auf angemessene Beschriftung und Formatierung deines Plots.

ggplot(ess8_noNA, aes(eduyrs, rlgdgr)) +
  geom_bin2d(binwidth = c(3, 1)) +
  scale_fill_gradient(low = "white", high = "deepskyblue") +
  geom_smooth(method = "lm") +
  scale_x_continuous(breaks = seq(0, 24, 3)) +
  scale_y_continuous(breaks = seq(0, 10, 1)) +
  coord_cartesian(xlim = c(3, 21), ylim = c(0, 10)) +
  labs(
    x = "Bildungsjahre",
    y = "Religiosität",
    title = "Bildung und Religiosität in der Schweiz",
    caption = "Daten: ESS (2016), Teilstichprobe CH (N=1405)\n0-'gar nicht religiös'; 10-'sehr religiös'",
    fill = "" 
  )

III. Multivariate Analyse 1

1. Integriere nun das Alter, als Indikator für das Störmerkmal “Generation”, in die Regressionsanalyse als Kontrollvariable. Speichere das Regressionsergebnis, z.B. unter dem Namen “multimodel1”.

multimodel1 <- lm(rlgdgr~eduyrs+agea, data = ess8_noNA )

2. Erstelle mit stargazer() eine Regressionstabelle mit den Koeffizienten beider Analysen

stargazer(bi_model,multimodel1, type="text", 
          title="Regressionsanalyse: Determinanten der Religiösität",
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("M1", "M2"),
          notes = c("Daten: ESS8 (European Social Survey)","Standardfehler in Klammern"),
          notes.align = "c",
          single.row = T,
          omit.stat = c("f", "ser", "adj.rsq"), 
          digits = 2,
          digits.extra = 5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          model.numbers= FALSE)

## 
## Regressionsanalyse: Determinanten der Religiösität
## ==================================================
##                                                   
##                      M1                 M2        
## --------------------------------------------------
## eduyrs        -0.10*** (0.02)    -0.08*** (0.02)  
## agea                             0.03*** (0.004)  
## Constant       5.89*** (0.27)     4.00*** (0.35)  
## --------------------------------------------------
## Observations       1,405              1,405       
## R2                  0.01               0.06       
## ==================================================
## Note:            *p<0.05; **p<0.01; ***p<0.001    
##              Daten: ESS8 (European Social Survey) 
##                   Standardfehler in Klammern

3. Kopiere und Finalisiere die Tabelle in einem separaten Textfile

stargazer(bi_model,multimodel1, type="html", 
          title="Regressionsanalyse: Determinanten der Religiösität",
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("M1", "M2"),
          notes = c("Daten: ESS8 (European Social Survey)","Standardfehler in Klammern"),
          notes.align = "c",
          single.row = T,
          omit.stat = c("f", "ser", "adj.rsq"), 
          digits = 2,
          digits.extra = 5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          model.numbers= FALSE,
          out = "reg_table.doc")

4. Werte deine (vorläufige) Regressionstabelle aus:

Interpretiere den Koeffizienten der Bildungsjahre im trivariaten Modell
Wie verändert sich der Koeffizient der Bildungsvariable? Wie erklärst du diese Veränderung?

a) Unter Kontrolle des Alters sinkt die Religiösität mit jedem zusätzlichen Bildungsjahr im Schnitt um etwa 0,08 Skalenpunkte. Zwar verringert sich der Bildungskoeffizient bei Integration der Kontrollvariable, aber er bleibt statistisch signifikant und ähnlich gross.

b) Das Alter hat wie erwartet einen positiven Einfluss auf die Religiösität (10 Jahre entsprechen immerhin einem drittel Skalenpunkt Religiösität): Ältere Menschen sind mehr, jüngere weniger religiös. Da offenbar junge Menschen im Mittel eine höhere Bildung haben, wird nun ihre altersbedingt geringe religiöse Neigung im Bildungskoeffizienten des ersten, bivariaten Modells mit transportiert. Dieser bivariate Koeffizient überschätzt somit die religiösitätshemmende Wirkung von Bildung. In der zweiten, trivariaten Analyse ergibt sich durch Abkopplung des störenden Alterseffekts folglich ein weniger akzentuierter Bildungsgradient der Religiösität.

5. Erstelle nun einen “Conditional Effect Plot” für dem Einfluss der Bildung im trivariaten Modell

visreg(multimodel1, xvar="eduyrs", gg = T, partial=F, rug=F)+
  labs(title = "Konditionaler Effekt der Bildung auf die Religiosität", 
       x="Bildungsjahre", 
       y="Religiosität", 
       subtitle = "0=gar nicht religiös, 10=sehr religiös", 
       caption = "Daten: ESS (2016), Teilstichprobe CH (n=1405) \n Vorhergesagte Werte unter Konstanthaltung des Alters im Mittelwert")+
  theme_bw()

Hinweis: Achte darauf, in der Beschriftung des Plots anzugeben, welche Variable konstant gehalten wird.

IV. Multivariate Analyse 2

1. Berechne eine weiteres multivariates Modell mit dem Alter und der väterlichen Bildung als Kontrollvariablen. Benenne es z.B. “multimodel2”.

multimodel2 <- lm(rlgdgr~eduyrs+agea+eduyrsf, data = ess8_noNA )

2. Integriere nun das multimodel2 in deine Regressionstabelle

stargazer(bi_model,multimodel1,  multimodel2, type="text", 
          title="Regressionsanalyse: Determinanten der Religiösität",
          dep.var.caption = "", 
          dep.var.labels = "", 
          column.labels = c("M1", "M2", "M3"),
          notes = c("Daten: ESS8 (European Social Survey)","Standardfehler in Klammern"),
          notes.align = "c",
          single.row = T,
          omit.stat = c("f", "ser", "adj.rsq"), 
          digits = 2,
          digits.extra = 5,
          star.cutoffs = c(0.05, 0.01, 0.001),
          model.numbers= FALSE)

## 
## Regressionsanalyse: Determinanten der Religiösität
## ============================================================
##                                                             
##                    M1              M2              M3       
## ------------------------------------------------------------
## eduyrs       -0.10*** (0.02) -0.08*** (0.02)  -0.06* (0.02) 
## agea                         0.03*** (0.004) 0.03*** (0.004)
## eduyrsf                                      -0.08** (0.03) 
## Constant     5.89*** (0.27)  4.00*** (0.35)  4.99*** (0.48) 
## ------------------------------------------------------------
## Observations      1,405           1,405           1,405     
## R2                0.01            0.06            0.06      
## ============================================================
## Note:                 *p<0.05; **p<0.01; ***p<0.001         
##                   Daten: ESS8 (European Social Survey)      
##                        Standardfehler in Klammern

3. Interpretiere den Bildungskoeffizienten der dritten Analyse im Kontext der gesamten Modellreihe

Unter Kontrolle des Alter und der väterlichen Bildung hat die eigene Bildung (weiterhin) einen signifikanten Einfluss auf die Religiosität. Die Grösse des Bildungskoeffizienten (und auch sein Signifikanzniveau) verringert sich allerdings im Vergleich zur Analyse ohne Herkunftskontrolle nochmals. Der unbereinigten Bildungskoeffizient im bivariaten Modell transportiert also zu einem substanziellen Anteil die Einflüsse der beiden integrierten Kontrollmerkmale und bildet somit den Bildungseffekt verzerrt ab.

Die Hypothese, dass Bildung einen negativen Einfluss auf die Religiösität hat, wird gleichwohl auch unter Kontrolle der Störmerkmale gestützt. In rein statistisch-formaler Hinsicht handelt es sich also um einen Zusammenhang, der robust ist gegenüber der Kontrolle durch Drittmerkmale. Inhaltlich betrachtet gilt dies allerdings nicht: Auf Basis der bivariaten Analyse überschätzen wir recht deutlich die religiösitätshemmende Wirkung von Bildung. Wir landen schliesslich bei einem multivariaten Schätzer für den Bildungseffekt, der für 5 zusätzliche Bildungsjahre lediglich ein Absinken der Religiösität um etwa 0,3 Skalenpunkte (bzw. um eine zehntel Standardabweichung) voraussagt. Unter Kontrolle von Störeffekten stellt sich der Einfluss der Bildung auf die Religiösität also eher als klein dar - zumindest wenn Effekte elterlicher Bildung ausgeklammert werden.

4. Wie religiös ist laut unserem Regressionsfit eine Person mit 15 Bildungsjahren im Alter von 24 Jahren, deren Vater 13 Bildungsjahre aufweist? In welchem Intervall liegt der Religiösitätsscore dieser Person mit 50% Sicherheit? Wo liegt der für Sie vorhergesagte Religiösitätsscore? Wo der Ihrer Mutter? Finden Sie sich in diesen Vorhersagen wieder?

predict(multimodel2, data.frame(eduyrs=15, agea = 24, eduyrsf=13))

##        1 
## 3.884648

predict(multimodel2, data.frame(eduyrs=15, agea = 24, eduyrsf=13), interval = "predict", level=0.9)

##        fit        lwr      upr
## 1 3.884648 -0.9508817 8.720178

predict(multimodel2, data.frame(eduyrs=15, agea = 24, eduyrsf=13), interval = "predict", level=0.5)

##        fit      lwr      upr
## 1 3.884648 1.902584 5.866712

mean(ess8_noNA$rlgdgr)

## [1] 4.77153

Auf Basis des Regressionsfits wird der oben charakterisierten Person ein unterdurchschnittlicher Religiösitätswert von 3,9 vorhergesagt. Gleichwohl gibt es viel Streuung der realisierten Werte um die Regressionsgerade bzw. -vorhersage. Legen wir einen hohen Sicherheitsstandard von 90% an, kann die Vorhersage “nur” auf das sehr grobe Intervall zwischen 0 und 8,8 begrenzt werden. Erst eine Absenkung des Sicherheitsstandards führt zu einem brauchbaren Vorhersageintervall: Mit 50% Sicherheit liegt der tatsächliche Religiösitätsscore der o.g. Person zwischen 3.4 und 7.4.

6. Visualisiere den Koeffizienten des finalen Modells, also den letztendlichen hypothesenbelastenden Effektschätzer deiner Regressionsanalyse

b_con_plot <-visreg(multimodel2, xvar="eduyrs", gg = T, partial = F, rug = F)+
  labs(title = "Konditionaler Effekt der Bildung auf die Religiosität", 
       x="Bildungsjahre", 
       y="Religiosität", 
       subtitle = "0=gar nicht religiös, 10=sehr religiös", 
       caption = "Daten: ESS (2016), Teilstichprobe CH (n=1405) \n Vorhergesagte Werte unter Konstanthaltung des Alters und väterlicher Bildung, jeweils im Mittelwert")+
  theme_bw()
b_con_plot

7. Zur Verdeutlichung der Kontrollwirkung der integrierten Drittvariablen integriere zusätzlich die Regressionsgerade des unbereinigten Zusammenhangs

b_con_plot+
  geom_smooth(method = "lm", se=FALSE, data = ess8_noNA , aes(eduyrs, rlgdgr), size = 0.5, color="red", linetype = "dashed")

## `geom_smooth()` using formula = 'y ~ x'

8. Speichere nun deinen konditionalen Effektplot über die Syntax als jpeg mit Länge 10cm, Breite 15cm und einer Auflösung von 1200ppi in deine working directrory und kopiere ihn beliebig in Folien, Poster, Dokumente etc.

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

BONUS. Visualisiere nun die Ergebnisse der trivariaten Regression (Religiosität ~ Bildung + Alter) in einem 3D-Plot. Bilde die Punktwolke, wie auch Regressionsebene ab. Beschrifte die Achsen auf passende Art und Weise. Vergegenwärtige dir nochmal die Bedeutung der Regressionsebene und ihrer beiden Neigungen, sowie überhaupt die multivariate Kontrolllogik der multiplen Regression.

multi_model <- lm(rlgdgr~eduyrs+agea, data = ess8_noNA )
plot3d(multi_model, 
       xlab = "Bildung", 
       ylab = "Alter", 
       zlab = "Religiositeat", 
       plane.col = "red")

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”