Genutzt wird der Datensatz Duncan aus dem Paket car. Siehe ?Duncan für Informationen zum Datensatz.
library(car)
head(Duncan)
type income education prestige
accountant prof 62 86 82
pilot prof 72 76 83
architect prof 75 92 90
author prof 55 90 76
chemist prof 64 86 90
minister prof 21 84 87
\[\widehat{prestige} = \beta_0 + \beta_1 income\]
m1 <- lm(prestige ~ income, data=Duncan)
coef(summary(m1))
Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.456574 5.1900585 0.473323 6.383769e-01
income 1.080390 0.1073688 10.062420 7.144356e-13
\[\widehat{prestige} = \beta_0 + \beta_1 income + \beta_2 education\]
m2 <- lm(prestige ~ income + education, data=Duncan)
coef(summary(m2))
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.0646629 4.27194117 -1.419650 1.630896e-01
income 0.5987328 0.11966735 5.003310 1.053184e-05
education 0.5458339 0.09825264 5.555412 1.727192e-06
Der Effekt von income in der multiplen Regression ist der Netto-Effekt wenn education kontrolliert wird. Anders ausgedrückt: Wenn education bereits seine Wirkung entfaltet hat, welchen zusätzlichen, von education linear unabhängigen Einfluss hat dann noch income? Um dies zu veranschaulichen können zwei Hilfsregressionen gerechnet werden:
- \(prestige = \beta_0 + \beta_1 education + e_{prestige}\)
- \(income = \beta_0 + \beta_1 education + e_{income}\)
m3 <- lm(prestige ~ education, data=Duncan)
e.prestige <- residuals(m3) # residuals of prestige
m4 <- lm(income ~ education, data=Duncan)
e.income <- residuals(m4) # residuals of income
Die Residuen der ersten Regression \(e_{prestige}\) sind jene Werte, die durch education nicht linear erklärt werden können. In diesen Werte ist der Einfluss von education entfernt oder kontrolliert.
In der zweiten Regression wird der Einfluss von education aus dem Prädiktor income entfernt. Die Residuen \(e_{income}\) enthalten nun Werte, die nicht durch education linear erklärt werden können.
Um nun die Frage zu beantworten, welchen Nettoeinfluss income auf prestige hat wird eine dritte Regression gerechnet. Die Residuen der ersten Hilfsregression werden auf die Residuen der zweiten Hilfsregression regridiert. Eine Regressionskonstante ist diesmal nicht nötig, da der Mittelwert der Residuen jeweils Null ergibnt und die Gerade somit bei \(P(0,0)\) durch den Ursprung läuft.
\[e_{prestige} = \beta_1 e_{income} + e\]
Das heisst, dass nun das um education bereinigte income das um education bereinigte prestige voraussagen soll. Der Koeffizient dieser Regression wird partieller Regressionkoeffzient genannt, da er nur jene Teile (Parts) der Beziehung, an denen education keinen Anteil mehr hat. Sein Wert ist identisch mit dem Koeffizienten in der oben gerechneten multiplen Regression.
m5 <- lm(e.prestige ~ e.income -1) # Kein Intercept hier
coef(summary(m5))
Estimate Std. Error t value Pr(>|t|)
e.income 0.5987328 0.116916 5.121051 6.45802e-06
Weiterhin besteht ein enger Zusammenhang zur Partialkorrelation. Während der partielle Regressionskoeffizient die Regression von \(e_{prestige}\) auf \(e_{income}\) ist, so ist der partielle Korrelationskoeffizient die Korrelation zwischen diesen.
\[
r_{(prestige\, income) \cdot education} = r_{e_{prestige} \; e_{income} }
\]
cor(e.prestige, e.income) # Korrelation zwischen den Residuen der Hilfsregressionen
[1] 0.6111006
Zur Kontrolle berechen wir die Psartialkorrelation mit Hilfe von pcor.test aus dem Paket ppcor.
library(ppcor) # dasselbe mit der Funktion für Partialkorrelation
attach(Duncan)
pcor.test(prestige, income, education)
estimate p.value statistic n gp Method
1 0.6111006 5.635427e-07 5.00331 45 1 pearson