Hinweise
Eine R-Nutzerin, die bereits vorher einen R-Kurs belegt hat, bewertete dieses Tutorial insgesamt mit einer Schwierigkeit von 5 (0=sehr leicht, 10=sehr schwer). Sie brauchte für dieses Tutorial ungefähr 1h50min. Klicke auf “Nächstes Kapitel” und es geht los.
Einführung
Regression ist das Allheilmittel der Psychologie, denn sie ist nicht nur von selbst sehr nützlich, auch kann man mit ihr die am häufigsten verwendeten Verfahren wie ANOVA und \(t\)-Test berechnen. Wenn also ein Methode der Statistik wirklich wichtig für Psychologen ist, so ist das die Regression. Deshalb schauen wir sie uns heute etwas genauer an. Wir fangen mit einem einfachen Beispiel an, bei dem wir die Bewertung von Filmen vorhersagen möchten. Dann versuchen wir die Äquivalenz zum \(t\)-Test und zur ANOVA aufzuzeigen. Anschließend kommt der Höhepunkt: mit einer linearen Regression lassen sich auch nicht-lineare Zusammenhänge modellieren. Linear ist nicht gleich linear! Zum Schluss schauen wir uns zwei einfache Wege an zu standardisierten Regressionskoeffizienten zu kommen, die R
standardmäßig nicht ausgibt.
Filmbewertungen vorhersagen
Vielleicht erinnerst Du Dich noch an die Filme-Datenbank aus Tag 2. Stell Dir vor Du kennst eine großartige Schauspielerin, die nur in großartigen Filmen mitspielen will. Die Schauspielerin beauftragt Dich ein Modell zu entwickeln mit dem sie entscheiden kann, welche Rollenangebote sie annimmt und welche sie ablehnt.
Laden wir uns erst mal die Daten:
imdb <- readRDS("data/imdb.rds")
imdb
Da die Regressions-Analyse Intervallskalenniveau voraussetzt, schaust Du Dir die entsprechenden Variablen an und nimmst explorativ zunächst das Budget des Films, die Filmlänge und das Alter ins Modell. Das ist nicht völlig atheoretisch; mehr Geld sollte zu besseren Filmen führen, zu kurze oder zu lange Filme könnten schlechter bewerten werden und ältere Filme polarisieren vermutlich.
Bevor wir das Modell bauen, müssen wir das Alter des Filmes aus dem Erscheinungsjahr berechnen:
current_year <- as.numeric(format(Sys.Date(), "%Y"))
imdb$age <- current_year - imdb$title_year
Das sieht etwas kryptisch aus, aber Sys.Date
gibt einfach nur das aktuelle Datum zurück und über format
holen wir uns dann das aktuelle Jahr. Das ist ein Detail, was nicht ganz so wichtig ist. Wir wollen einfach nur das aktuelle Alter des Films berechnen.
Und nun das Modell:
model1 <- lm(imdb_score ~ budget + duration + age, data = imdb)
summary(model1)
##
## Call:
## lm(formula = imdb_score ~ budget + duration + age, data = imdb)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.7897 -0.5797 0.1058 0.6984 3.1809
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.385e+00 7.584e-02 57.822 <2e-16 ***
## budget 7.748e-11 7.427e-11 1.043 0.297
## duration 1.588e-02 6.776e-04 23.439 <2e-16 ***
## age 1.478e-02 1.238e-03 11.941 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.028 on 4534 degrees of freedom
## (505 observations deleted due to missingness)
## Multiple R-squared: 0.1493, Adjusted R-squared: 0.1487
## F-statistic: 265.2 on 3 and 4534 DF, p-value: < 2.2e-16
Mit lm
machen wir ein lineares Modell, was äquivalent zu einer Regression ist. Die Syntax bedient sich hier der Formelschreibweise von R
. Zunächst kommt die AV, dann die Tilde ~
und dann die UVs mit einem +
voneinander getrennt. Genau wie beim Boxplot mit mehreren UVs (siehe Tag 3). Schließlich sagen wir lm
noch in welchem Datensatz sich diese Variablen befinden. Das ist sehr praktisch, ansonsten müssten wir immer wieder data$
vor die Variablen schreiben (data$imdb_score ~ data$action + ...
), was viel mehr Tipparbeit wäre.
Der Output erfolgt dann am einfachsten über summary
, wobei dies zugegebenerweise etwas unübersichtlich ist. Wir können uns hier eines schönen Pakets von Daniel Lüdecke bedienen:
install.packages("sjPlot")
library(sjPlot)
Und nun:
tab_model(model1)
imdb_score | |||
---|---|---|---|
Predictors | Estimates | CI | p |
(Intercept) | 4.38 | 4.24 – 4.53 | <0.001 |
budget | 0.00 | -0.00 – 0.00 | 0.297 |
duration | 0.02 | 0.01 – 0.02 | <0.001 |
age | 0.01 | 0.01 – 0.02 | <0.001 |
Observations | 4538 | ||
R2 / R2 adjusted | 0.149 / 0.149 |
Das ist mehr als nur gut und lässt Deinen Abschluss-Arbeits-Betreuer, der nur SPSS kann, ganz schön staunen!
Welche Bewertung bekommt ein Film, der 0 Jahre alt ist und 120 Minuten geht? (Ignoriere hierfür den Effekt des Budgets. Runde auf 2 Dezimalstellen)
Überraschenderweise hat das Budget keinen signifikanten Einfluss auf die Bewertung, aber durchaus die Länge des Films. Wir erinnern uns, dass bei IMDB auch Serien enthalten sind, die eventuell generell schlechter bewertet werden als Filme. Aber die Variable Film/Serie ist ja nicht intervallskaliert, also können wir sie nicht in die Regression aufnehmen. Falsch gedacht! Es spricht nichts dagegen dichotome Variablen in eine Regression aufzunehmen:
imdb$film <- imdb$duration > 75
model2 <- lm(imdb_score ~ budget + duration + age + film, data = imdb)
tab_model(model2)
imdb_score | |||
---|---|---|---|
Predictors | Estimates | CI | p |
(Intercept) | 5.29 | 5.00 – 5.58 | <0.001 |
budget | 0.00 | -0.00 – 0.00 | 0.310 |
duration | 0.02 | 0.02 – 0.02 | <0.001 |
age | 0.01 | 0.01 – 0.02 | <0.001 |
filmTRUE | -1.02 | -1.31 – -0.74 | <0.001 |
Observations | 4538 | ||
R2 / R2 adjusted | 0.158 / 0.158 |
Der Effekt der Serie lässt sich hier einfach als Mittelwertsunterschied lesen. Wenn film == 0
, dann ist es eine Serie und wir sind bei einer Schätzung von 5.33 (wenn budget, duration und age 0 wären). Ist es ein Film, dann müssen wir 1.02 abziehen. Heißt also einfach nur, dass Serien generell besser bewertet werden – überraschend, denn spontan hätte ich das Gegenteil erwartet. Aber vielleicht ist der Grund für die bessere Bewertung von Serien ein Selektionsprozess: Personen, denen die Serie gefällt (Fans also) machen sich die Mühe eine Bewertung abzugeben, wohingegen jemand der die Serie nicht mag, sie einfach nicht schaut und keine Bewertung abgibt. Bei Filmen könnte dieser Prozess anders ablaufen.
Trotz des riesigen Effekts von Film/Serie bleibt der Einfluss der Dauer erhalten (immer noch 0.02). Da unsere theoretischen Überlegungen hierzu eher primitiver Natur sind, wollen wir erst mal nicht tiefer graben. Stattdessen überlegst Dir nun, dass das Genre sicher auch eine Rolle spielt. Action-Filme werden Deiner Ansicht nach besser bewertet als andere Filme.
Wir erstellen uns also eine Variable für Action-Film/kein Action-Film. Das stellt uns schon vor eine Herausforderung, denn üblicherweise wird ein Film in mehrere Kategorien sortiert (z. B. Action und Komödie). Wir nutzen hier reguläre Ausdrücke über grepl
um die Variable zu erstellen. Reguläre Ausdrücke sind ein Thema über welches man einen separaten Kurs halten könnte. Hier sollst Du nur lernen, dass es so etwas gibt und wie man es im einfachsten Fall verwendet.
imdb$action <- grepl("Action", imdb$genres)
grepl
durchsucht alle Strings und gibt 1 zurück, falls der Target-String vorkommt, ansonsten 0. Wir bekommen also eine Dummy-Variable. Und nun können wir das Modell bauen:
model2 <- lm(imdb_score ~ budget + duration + age + film + action, data = imdb)
tab_model(model2)
imdb_score | |||
---|---|---|---|
Predictors | Estimates | CI | p |
(Intercept) | 5.28 | 4.99 – 5.57 | <0.001 |
budget | 0.00 | -0.00 – 0.00 | 0.121 |
duration | 0.02 | 0.02 – 0.02 | <0.001 |
age | 0.01 | 0.01 – 0.02 | <0.001 |
filmTRUE | -0.98 | -1.27 – -0.70 | <0.001 |
actionTRUE | -0.31 | -0.38 – -0.24 | <0.001 |
Observations | 4538 | ||
R2 / R2 adjusted | 0.173 / 0.172 |
Okay, Action-Filme werden tatsächlich eher schlechter bewertet, vielleicht weil sie nicht gerade deep sind.
Du bist schon recht zufrieden mit dem Modell und empfiehlst der Schauspielerin in Serien zu spielen, die wenig Action enthalten und bald ausgestrahlt werden, denn alte Produktionen werden tendenziell schlechter bewertet. Aber Deine Auftraggeberin ist noch nicht vollends zufrieden. Sie schlägt vor ins Modell Regisseure und Schauspieler einzubeziehen. Mache zunächst Dummy-Variablen für die Regisseure Christopher Nolan, Stanley Kubrick und Guy Ritchie (Variable director_name
) sowie den Schauspieler Christian Bale in der Hauptrolle (Variable actor_1_name
), rechne anschließend die Regression (die bisherigen Variablen sollen drin bleiben).
imdb$nolan <- grepl("Christopher Nolan", imdb$director_name)
imdb$kubrick <- grepl("Stanley Kubrick", imdb$director_name)
imdb$ritchie <- grepl("Guy Ritchie", imdb$director_name)
imdb$bale <- grepl("Christian Bale", imdb$actor_1_name)
model1 <- lm(imdb_score ~ budget + duration + age + film + action + nolan + kubrick + ritchie + bale, data = imdb)
tab_model(model1)
Wenn Die Schauspielerin die Chance hat neben Christian Bale zu spielen oder bei einem Film von Guy Ritchie oder Christopher Nolan, kann sie mit einer deutlich überdurchschnittlichen Bewertung des Films rechnen (bei Stanley Kubrick natürlich auch, aber dieser ist 1999 gestorben). Allerdings widerspricht das etwas der Idee, dass sie in Serien spielen soll, denn die genannten Regisseure machen eigentlich nur Filme. Und da kommen wir schon schnell zu einem richtig großen Problem.
Generell sollte Dir bei diesem Beispiel auffallen, dass wir ganz schön im Schlamm wühlen, weil wir wenig Theorie mitgebracht haben und hier natürlich kein Experiment vorliegt. Die Variablen sind also vermutlich alle untereinander korreliert und es bleibt unklar was Ursache und Wirkung ist. Beispielsweise spielt Christian Bale in Filmen von Christopher Nolan mit (Batman-Triologie), die alle ausgezeichnete Bewertungen haben. Unser Modell mag einen praktischen Nutzen haben, aber richtige Wissenschaft funktioniert anders.
Die Regression lässt sich aber durchaus auch in experimentellen Settings verwenden. Wenn wir die Möglichkeit hätten die Dauer des Films und den Regisseur zufällig zu Filmen zuzuweisen, könnten wir die Beziehungen viel klarer analysieren. Das ist jedoch praktisch kaum möglich und zeigt die notwendige Kompromissbereitschaft in der Psychologie zwischen dem was ideal wäre (ein Experiment) und dem was möglich ist (korrelative Studien des Gegebenen). Schauen wir uns jetzt mal den Idealfall an: die Anwendung der Regression bei Experimenten. Normalerweise würden wir diese mit dem \(t\)-Test oder einer ANOVA analysieren, aber es geht auch mit einer Regression. Der Vorteil liegt auf der Hand, statt verschiedene Methoden zu lernen, konzentrierst Du dich auf eine und kannst diese dann auf viele Fragestellungen anwenden.
Äquivalenz zu t-Test
In Methodenlehre I haben wir einen \(t\)-Test für folgende Tabelle gerechnet:
EG | KG |
---|---|
-0.09 | -0.26 |
0.01 | -0.26 |
0.03 | -0.21 |
-0.31 | -0.08 |
0.34 |
Die Experimentalgruppe bekam ein Medikament, die Kontrollgruppe ein Placebo. Die AV ist die Flugleistung in einem Simulator. Die Probanden waren erfahrene Piloten.
An Tag 4 haben wir die Berechnung in R
durchgeführt:
t.test(c(-.09, .01, .03, -.31, .34), c(-.26, -.26, -.21, -.08), var.equal = TRUE)
##
## Two Sample t-test
##
## data: c(-0.09, 0.01, 0.03, -0.31, 0.34) and c(-0.26, -0.26, -0.21, -0.08)
## t = 1.5902, df = 7, p-value = 0.1558
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.09667611 0.49367611
## sample estimates:
## mean of x mean of y
## -0.0040 -0.2025
Wir hätten stattdessen auch eine Regression rechnen können:
av <- c(-.09, .01, .03, -.31, .34, -.26, -.26, -.21, -.08)
uv <- c(rep(1, 5), rep(0, 4))
tab_model(lm(av ~ uv))
Dependent variable | |||
---|---|---|---|
Predictors | Estimates | CI | p |
(Intercept) | -0.20 | -0.42 – 0.02 | 0.066 |
uv | 0.20 | -0.10 – 0.49 | 0.156 |
Observations | 9 | ||
R2 / R2 adjusted | 0.265 / 0.160 |
Du siehst, der \(p\)-Wert von 0.156 ist identisch zum \(p\)-Wert beim normalen \(t\)-Test. Außerdem kannst Du die Durchschnittswerte für die Gruppen ablesen. Wenn die UV 0 ist (KG), dann ist nur der Intercept relevant und es kommt -0.20 heraus. Wenn die UV 1 ist (EG), dann wird zum Intercept 0.2 addiert und wir sind bei 0. Der Unterschied beträgt also 0.2 in der Flugleistung (für mehr Kontext, schau Dir die Übungsaufgabe an).
Vergleiche die Bewertung von Filmen auf imdb nach der Variable color
(schwarz-weiß oder Farbe) mit einem \(t-\)Test und einer Regression. Achte auf die Schreibweise der Ausprägungen und entferne zunächst die Filme ohne Angabe zur Variable color
(""
als Ausprägung).
imdb <- imdb[imdb$color != "",]
t.test(imdb$imdb_score[imdb$color == " Black and White"],
imdb$imdb_score[imdb$color == "Color"],
var.equal = TRUE)
summary(lm(imdb_score ~ color, data = imdb))
Äquivalenz zu ANOVA
Äquivalent ist die Regression auch zur ANOVA, allerdings nur für den einfaktoriellen Fall. Schauen wir uns ein Beispiel aus der Vorlesung an (Folie 5):
vlf5 <- data.frame(av = matrix(c(18, 18, 20, 17, 9, 16, 25, 24, 16)),
gruppe = rep(c("placebo", "einfache Dosis", "doppelte Dosis"), each = 3))
vlf5
Die ANOVA ist:
vlf5$wid <- 1:9
ezANOVA(vlf5, dv = av, between = gruppe, wid = wid, detailed = TRUE)
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd SSn SSd F p p<.05 ges
## 1 gruppe 2 6 89.55556 89.33333 3.007463 0.1245347 0.5006211
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 2 6 9.555556 69.33333 0.4134615 0.6788577
Die Regression:
summary(lm(av ~ gruppe, data = vlf5))
##
## Call:
## lm(formula = av ~ gruppe, data = vlf5)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.6667 -0.6667 1.3333 2.3333 3.3333
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.667 2.228 9.726 6.79e-05 ***
## gruppeeinfache Dosis -7.667 3.151 -2.433 0.0509 .
## gruppeplacebo -3.000 3.151 -0.952 0.3778
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.859 on 6 degrees of freedom
## Multiple R-squared: 0.5006, Adjusted R-squared: 0.3342
## F-statistic: 3.007 on 2 and 6 DF, p-value: 0.1245
Der gleiche \(F\)-Wert, der gleiche \(p\)-Wert und die gleiche Effektgröße.
Nun Du. In Methodenlehre I haben wir untersucht, ob die Art des Unterrichts einen Einfluss auf die Lernleistung hat:
av <- matrix(c(2, 1, 3, 3, 1, 3, 4, 3, 5, 0, 6, 8, 7, 4, 10, 5, 5, 5, 3, 2))
unterricht <- data.frame(av = av,
uv = rep(c("Gruppenunterricht", "Selbststudium",
"Einzelunterricht","ELearning"), each = 5),
wid = 1:20)
unterricht
Rechne eine ANOVA und Regression und zeige, dass beide Verfahren identisch sind:
ezANOVA(unterricht, dv = av, between = uv, wid = wid)
summary(lm(av ~ uv, data = unterricht))
Für die mehrfaktorielle ANOVA wird die Sache schon schwieriger. An Tag 4 haben wir eine ANOVA (aus der Vorlesung auf Folie 12) nachgerechnet:
vlf12 <- data.frame(av = c(18, 18, 20, 13, 15, 9, 17, 9, 16, 15, 17, 22, 25,
24, 16, 17, 12, 18),
wid = 1:18,
geschlecht = factor(rep(c(c("m", "m", "m"), c("f", "f", "f")), 3)),
medikament = factor(rep(c("p", "1", "2"), each = 6)))
ezANOVA(data = vlf12, dv = av, between = c(geschlecht, medikament), wid = wid, detailed = T)
## Warning: Converting "wid" to factor for ANOVA.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd SSn SSd F p p<.05
## 1 geschlecht 1 12 34.72222 154.6667 2.693966 0.12665828
## 2 medikament 2 12 34.77778 154.6667 1.349138 0.29613457
## 3 geschlecht:medikament 2 12 103.44444 154.6667 4.012931 0.04629565 *
## ges
## 1 0.1833382
## 2 0.1835777
## 3 0.4007749
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 5 12 9.777778 104 0.225641 0.9442807
In diesem Fall kommt nicht das gleiche heraus wie bei einer Regression:
summary(lm(av ~ geschlecht * medikament, data = vlf12))
##
## Call:
## lm(formula = av ~ geschlecht * medikament, data = vlf12)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.667 -2.500 1.000 2.333 4.000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.000 2.073 8.684 1.61e-06 ***
## geschlechtm -4.000 2.931 -1.365 0.1974
## medikament2 -2.333 2.931 -0.796 0.4415
## medikamentp -5.667 2.931 -1.933 0.0772 .
## geschlechtm:medikament2 10.000 4.146 2.412 0.0328 *
## geschlechtm:medikamentp 10.333 4.146 2.493 0.0283 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.59 on 12 degrees of freedom
## Multiple R-squared: 0.5279, Adjusted R-squared: 0.3312
## F-statistic: 2.684 on 5 and 12 DF, p-value: 0.07484
Man beachte hier die Spezifikation der Interaktion über geschlecht * medikament
, R
macht dann ganz automatisch die Dummykodierung für uns.
In der ANOVA gibt es kein Äquivalent zum \(F\)-Wert für das Gesamtmodell der Regression. Aber das können wir schnell beheben. Der \(F\)-Wert der Regression bezieht sich auf die Gruppeneinteilung der Regression, sie unterscheidet also nicht zwischen den einzelnen Faktoren (es gibt nur eine Gruppenvariable, nicht zwei!). Das Ganze können wir auch mit einer ANOVA nachbauen, indem wir so tun als gäbe es nur eine Gruppenvariable:
vlf12$group <- as.factor(rep(1:6, each = 3))
vlf12
Nun würden wir erwarten, dass der gleiche \(F\)-Wert herauskommt wie bei der Regression (2.684), wenn wir group
als between-Faktor verwenden (statt medikament
und geschlecht
)
ezANOVA(data = vlf12, dv = av, between = c(group), wid = wid, detailed = T)
## Warning: Converting "wid" to factor for ANOVA.
## Coefficient covariances computed by hccm()
## $ANOVA
## Effect DFn DFd SSn SSd F p p<.05 ges
## 1 group 5 12 172.9444 154.6667 2.683621 0.07484044 0.5278955
##
## $`Levene's Test for Homogeneity of Variance`
## DFn DFd SSn SSd F p p<.05
## 1 5 12 9.777778 104 0.225641 0.9442807
Passt. Obwohl also die mehrfaktorielle ANOVA nicht ganz identisch zur Regression ist, lassen sie sich durchaus ineinander überführen, sofern man die Details kennt. Die Ideen hinter den Verfahren sind aber grundverschieden: Die ANOVA vergleicht Varianzen, die Regression minimiert den Fehler bei einer Schätzung. Wird der Fehler kleiner, werden Varianzverhältnisse größer.
Linear ist nicht gleich linear
Die Linearität der linearen Regression bezieht sich auf die Parameter, nicht auf die Zusammenhänge. Man kann also sehr wohl nicht-lineare Zusammenhänge mit einer linearen Regression modellieren. Man sollte es sogar!
In R
gibt es den Datensatz mtcars
, der Informationen zu 32 Automobilen von 1974 enthält. Wir vermuten, dass die Beziehung zwischen Verbrauch (Miles/(US) gallon) und Gewicht (1000 lbs) nicht-linear ist. Machen wir zunächst eine Abbildung:
plot(mtcars$wt, mtcars$mpg)
lines(lowess(x = mtcars$wt, y = mtcars$mpg))
Eine leichte nicht-Linearität ist zu erkennen. Wir können laut Potenzleiter entweder \(x\) oder \(y\) mit einem Exponenten der kleiner 1 ist versehen. Üblicherweise transformiert man eher \(x\), da die Interpretation für ein transformiertes \(y\) schwierig ist. Aber wir müssen nicht unbedingt transformieren, sondern können einfach die gewünschte Beziehung mit lm
spezifizieren. Vergleichen wir mal einen linearen Zusammenhang mit einem Exponenten von 0.5.
tab_model(lm(mpg ~ wt, data = mtcars), # model 1, linearer Zusammenhang
lm(mpg ~ I(wt^0.5), data = mtcars)) # model 2, nicht-linearer Zusammenhang
mpg | mpg | |||||
---|---|---|---|---|---|---|
Predictors | Estimates | CI | p | Estimates | CI | p |
(Intercept) | 37.29 | 33.45 – 41.12 | <0.001 | 54.75 | 48.02 – 61.49 | <0.001 |
wt | -5.34 | -6.49 – -4.20 | <0.001 | |||
wt^0 5 | -19.55 | -23.30 – -15.79 | <0.001 | |||
Observations | 32 | 32 | ||||
R2 / R2 adjusted | 0.753 / 0.745 | 0.790 / 0.783 |
Zu Beachten ist hier die Syntax I(wt^0.5)
, was die direkte Umwandlung der Variable ermöglicht, ohne, dass wir sie vorher transformieren müssen. Wenn wir die Modelle aufschreiben würden, wären sie:
\(\mathrm{mpg} = b_0 + b_1 \cdot \mathrm{wt}+ e\)
\(\mathrm{mpg} = b_0 + b_1 \cdot \mathrm{wt}^{0.5} + e\)
Und die Tabelle gibt uns nun die Koeffizienten:
\(\mathrm{mpg} = 37.29 - 5.34 \mathrm{wt} + e\)
\(\mathrm{mpg} = 54.75 -19.55 \cdot \mathrm{wt}^{0.5} + e\)
Und welches Modell ist besser? Man sollte hier zum einen theoretische Argumente anführen, zum anderen statistische. Zu den statistischen gehört \(R^2\), welches für das nicht-lineare Modell etwas höher ist. Vermutlich könnten wir auch einen besseren Exponenten finden, im Idealfall aus der Theorie abgeleitet!
Füge nun noch zu dem existierenden Modell die Variable hp
(Horsepower) hinzu. Wähle für den Exponenten 1 (linearer Zusammenhang), 0.5 und 0.1 und vergleiche die Modelle miteinander.
tab_model(lm(mpg ~ I(wt^0.5), data = mtcars),
lm(mpg ~ I(wt^0.5) + hp, data = mtcars),
lm(mpg ~ I(wt^0.5) + I(hp^0.5), data = mtcars),
lm(mpg ~ I(wt^0.5) + I(hp^0.1), data = mtcars))
Erstelle nun noch eine Abbildung für den Zusammenhang zwischen Horsepower und Miles per Gallon. Füge eine lowess-Kurve hinzu.
plot(mtcars$hp, mtcars$mpg)
lines(lowess(x = mtcars$hp, y = mtcars$mpg))
Bevor wir das Thema abschließen, worauf bezieht sich denn nun die Linearität bei linearen Regressionen? Schau Dir hierzu die nächste Aufgabe an:
Na, hattest Du alles richtig? Wie findet man heraus ob ein Modell linear in den Parametern ist? Man bildet einfach die Ableitung nach den Parametern für \(y\):
\(\frac{\partial y}{\partial b_i}\)
und schaut ob diese von anderen Parametern abhängen. Wenn Du das machst, wirst Du sehen, dass nur eine Gleichung nicht-linear ist, die vierte: \(\frac{\partial y}{\partial b_1} = e^{b_2x}\) und \(\frac{\partial y}{\partial b_2} = b_1x e^{b_2x}\).
Beide Ableitungen sind abhängig von einem anderen Parameter (\(b_2\), \(b_1\)), also ist das Modell nicht-linear. Mehr dazu hier: https://youtu.be/40w50jAOBsY?feature=shared
Standardisierung
Standardisierte Regressionskoeffizienten bekommen wir nicht automatisch in R
. Du hast hier zwei Möglichkeiten. Die einfache ist:
tab_model(lm(mpg ~ wt + hp, data = mtcars), show.std = T)
mpg | |||||
---|---|---|---|---|---|
Predictors | Estimates | std. Beta | CI | standardized CI | p |
(Intercept) | 37.23 | 0.00 | 33.96 – 40.50 | -0.16 – 0.16 | <0.001 |
wt | -3.88 | -0.63 | -5.17 – -2.58 | -0.84 – -0.42 | <0.001 |
hp | -0.03 | -0.36 | -0.05 – -0.01 | -0.57 – -0.15 | 0.001 |
Observations | 32 | ||||
R2 / R2 adjusted | 0.827 / 0.815 |
Wir benutzen hier also einfach tab_model
und sagen, dass wir standardisierte Koeffizienten wollen (show.std = TRUE
). Den Rest übernimmt das Paket.
Alternativ kannst Du auch die Variablen vorher selbst \(z\)-standardisieren:
mtcars_standardisiert <- as.data.frame(scale(mtcars))
mtcars_standardisiert
Und nun nochmal das Modell rechnen:
tab_model(lm(mpg ~ wt + hp, data = mtcars_standardisiert))
mpg | |||
---|---|---|---|
Predictors | Estimates | CI | p |
(Intercept) | 0.00 | -0.16 – 0.16 | 1.000 |
wt | -0.63 | -0.84 – -0.42 | <0.001 |
hp | -0.36 | -0.57 – -0.15 | 0.001 |
Observations | 32 | ||
R2 / R2 adjusted | 0.827 / 0.815 |
Du siehst, es kommt das gleiche heraus.
Probier Du es nun. Mache ein Modell, bei dem Du die Zeit für 1/4 Meile vorhersagst (qsec
) aus den PS (hp
) und dem Gewicht (wt
). Berechne die standardisierten Koeffizienten.
mtcars_standardisiert <- as.data.frame(scale(mtcars))
tab_model(lm(qsec ~ hp + wt, data = mtcars_standardisiert))
# ODER
tab_model(lm(qsec ~ hp + wt, data = mtcars), show.std = TRUE)
Übungsaufgabe Tag 6
Heute sollst Du mal weniger R
wiederholen und Dich stattdessen intensiver mit der Regression beschäftigen. Deine Aufgabe besteht darin in R
eine eigene Funktion für die Berechnung der einfachen Regression zu schreiben. Vielleicht musst Du hierfür ein paar Dinge wiederholen oder auch neu lernen. Du wirst auf jeden Fall die Formelsammlung benötigen um die Koeffizienten zu berechnen.
Prüfe Deine Funktion anhand des folgenden Beispiels:
lm(qsec ~ hp, data = mtcars)
##
## Call:
## lm(formula = qsec ~ hp, data = mtcars)
##
## Coefficients:
## (Intercept) hp
## 20.55635 -0.01846
meine_regression <- function(x, y){
slope <- cor(x, y) * sd(y) / sd(x)
intercept <- mean(y) - slope * mean(x)
c(slope, intercept)
}
meine_regression(mtcars$hp, mtcars$qsec)
Wenn Du immer noch nicht genug hast, dann versuche eine eigene Funktion für die multiple Regression zu schreiben.