2 Lineare Regression
rm(list = ls())
#Einlesen der benötigten Pakete
library(plotly)
library(data.table)
library(knitr)
library(lmtest)
library(sandwich)
library(nlme)
library(tidyquant)
library(data.table)
#definiere Farben
the_grey <- "rgb(215, 215, 215)"
the_green <- "rgb(185, 205, 165)"
the_blue <- "rgb(135, 135, 201)"
the_red <- "rgb(200, 130, 130)"Um die Funktionsweise eines typischen Regressionsmodells besser zu verstehen, wollen wir im ersten Schritt Daten simulieren, die genau aus diesem Datengenerierenden Prozess stammen.
#Anzahl der Datenpunkte
n <- 100
#Standardabweichung der abhängigen Variablen
sigma <- 5
#Erzeuge fiktive deterministische erklärende, unabhängige Variablen zwischen 20 und 30
x <- seq(20, 30, length.out = 15)
x_new <- sample(x, size = n, replace = TRUE)
#Simuliere Störterme
eps <- rnorm(n, mean = 0, sd = sigma)
#wahre Werte der Regressionsgeraden
beta_0 <- 20
beta_1 <- 3
#abhängige Variable
y <- beta_0 + beta_1 * x_new + eps
plot_ly() %>%
add_markers(x = x_new, y = y, marker = list(color = the_green)) %>%
layout(
xaxis = list(title = 'Temperatur'),
yaxis = list(title = 'Eiskugeln je Stunde')
)2.1 Galton’s Datensatz - Haben große Väter große Kinder?
Sir Galton kann als der Ideengeber der linearen Regressionsanalyse angesehen werden. Seine Idee zur Regressionanalyse entstand aus seinen Beobachtungen zwischen der Größe der Kinder und der Größe von deren Vätern.
df <- data.table(read.csv(file = 'Data/Galton.csv', sep = ';'))
avg_heigt_per_cat <- df[ ,.(m_heigt = mean(Height)), by = .(Father)]
p1 <- plot_ly() %>%
add_markers(x = df$Father, y = df$Height, marker = list(color = the_green)) %>%
layout(xaxis = list(title = 'Größe des Vaters'),
yaxis = list(title = 'Größe des Kindes'))
p2 <- plot_ly() %>%
add_markers(x = avg_heigt_per_cat$Father, y = avg_heigt_per_cat$m_heigt, marker = list(color = the_green)) %>%
layout(xaxis = list(title = 'Größenkategorie des Vaters'), yaxis = list(title = 'Durchschnittsgröße der Kinder je Kategorie'))
subplot(p1, p2, nrows=1, shareX = T,shareY = T, titleX=T, titleY=T) %>% layout(showlegend=F)Schätzt man die Regressionsgerade scheint ein lineares Modell plausibel.
lin_reg <- lm(Height ~ Father, data = df)
beta0_hat <- lin_reg$coefficients[1]
beta1_hat <- lin_reg$coefficients[2]
plot_ly() %>%
add_markers(x = df$Father, y = df$Height, marker = list(color = the_green), name = 'Daten') %>%
add_lines(x = df$Father, y = fitted(lin_reg), name = 'Regressionslinie', line = list(color = the_blue)) %>%
layout(xaxis = list(title = 'Größe des Vaters'),
yaxis = list(title = 'Größe des Kindes'))2.2 Metrische Variablen mit Variablentransformation im linearen Regressionsmodell
In folgendem Beispiel besteht zwischen der Kovariablen und der abhängigen Variablen offensichtlich kein linearer, sondern ein quadratischer Zusammenhang. Transformiert man jedoch die Kovariable durch \(f(x) = x^2\) und nimmt die quadrierte Kovariable ins Modell auf, so kann dieser Zusammenhang durch das lineare Modell abgebildet werden.
#SP-Umfang
n <- 100
#Parameter der Regressionsgeraden
beta0 <- 20
beta1 <- 3
#Standardabweichung der abh. Variable
sigma <- 5
x <- rnorm(n, mean = 0, sd = 2)
x_squared <- x**2
y <- beta0 + beta1 * x_squared + rnorm(n, mean = 0, sd = sigma)
p1 <- plot_ly() %>%
add_markers(x = x, y = y, marker = list(color = the_green)) %>%
layout(
xaxis = list(title = TeX('$x$')),
yaxis = list(title = TeX('$y$'))
) %>%
config(last_plot(), mathjax = 'local')
p2 <- plot_ly() %>%
add_markers(x = x_squared, y = y, marker = list(color = the_green)) %>%
layout(
xaxis = list(title = TeX('$x^2$')),
yaxis = list(title = TeX('$y$'))
) %>%
config(last_plot(), mathjax = 'local')
subplot(p1, p2, nrows=1, shareX = T,shareY = T, titleX=T, titleY=T) %>% layout(showlegend=F) 2.3 Modell mit kategorialen Kovariablen in Dummykodierung
Angenommen wir haben drei Kategorien für das Wetter (1 = schön, 2 = neutral, 3 = schlecht). Wähle 2 als die Referenzkategorie und erzeuge zwei Dummyvariablen. Mit \(\beta_0 = 50, ~ \beta_{\text{schön}} = 30, ~ \beta_{\text{schlecht}} = -20\) würde dies implizieren, dass sich der Effekt bei schönem Wetter erhöht und bei schlechtem Wetter reduziert.
#SP Umfang
n <- 100
#Parameter der Regressionsgeraden
beta0 <- 50
beta1 <- 30
beta2 <- -20
#Standardabweichung der abh. Variable
sigma <- 5
x <- sample(c(1,2,3), size = n, replace = TRUE)
xGut <- (x == 1)
xSchlecht <- (x == 3)
y <- beta0 + beta1 * xGut + beta2 * xSchlecht + rnorm(n, mean = 0, sd = sigma)
plot_ly()%>%
add_markers(x = x, y = y, marker = list(color = the_green)) %>%
layout(xaxis = list(title = 'Wetterindex'), yaxis = list(title = 'Wohlfühlindex'))2.4 Modellschätzung des linearen Regressionsmodells
Verwenden wir das Beispiel des simulierten Modells. Hierbei geben wir uns Werte für \(\beta_0 = 20\), \(\beta_1 = 3\) und \(\sigma^2 = 25\) vor. In der Realität sind uns diese Werte unbekannt und wir müssen sie schätzen. Betrachten wir erst einmal die Koeffizienten der Regressionsgerade und was mit der Gerade passiert, wenn wir hierfür verschiedene Werte einsetzen.
#SP-Umfang
n <- 100
#Standardabw. der abh. Variable
sigma <- 5
#Deterministische Variable
x <- seq(20, 30, length.out = 15)
x_new <- sample(x, size = n, replace = TRUE)
#Parameter der Regressionsgeraden
beta0 <- 20
beta1 <- 3
y = beta0 + beta1 * x_new + rnorm(n, mean = 0, sd = sigma)
plot_ly()%>%
add_markers(x = x_new, y = y, marker = list(color = the_green)) %>%
add_lines(x = c(min(x_new), max(x_new)), y = c(20 + 3 * min(x_new), 20 + 3 * max(x_new)), line = list(color = the_blue)) %>%
add_lines(x = c(min(x_new), max(x_new)), y = c(-5 + 4 * min(x_new), -5 + 4 * max(x_new)), line = list(color = the_grey)) %>%
add_lines(x = c(min(x_new), max(x_new)), y = c(40 + 2.3 * min(x_new), 40 + 2.3 * max(x_new)), line = list(color = the_red)) %>%
layout(xaxis = list(title = 'Temperatur'), yaxis = list(title = 'Eiskugeln je Stunde'), showlegend = FALSE)2.5 Schätzen des linearen Regressionsmodells
In R ist es am einfachsten das lineare Regressionsmodell durch die im Basismodul vorhandene lm-Funktion zu schätzen. Wie zu den meisten anderen Regressionsmodellen in R, existieren sogenannte generische Funktionen, die spezifische Methoden für das Modell, wie z.B. eine Zusammenfassung oder die Prognose, verfügbar machen.
#Schätzen des Modells
lin_reg <- lm(Height ~ Father, data = df)
#Ausgabe der Zusammenfassung des Modells
summary(lin_reg)##
## Call:
## lm(formula = Height ~ Father, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.2683 -2.6689 -0.2092 2.6342 11.9329
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.11039 3.22706 12.120 <2e-16 ***
## Father 0.39938 0.04658 8.574 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.446 on 896 degrees of freedom
## Multiple R-squared: 0.07582, Adjusted R-squared: 0.07479
## F-statistic: 73.51 on 1 and 896 DF, p-value: < 2.2e-16## 1 2 3 4 5 6
## 70.46182 70.46182 70.46182 70.46182 69.26367 69.26367##
## Call:
## lm(formula = Height ~ Father, data = df)
##
## Coefficients:
## (Intercept) Father
## 39.1104 0.3994#eine generische Funktion ist z.B. predict mit der neue Prognosen für neue Daten erstellt werden können
new_data <- data.frame('Father' = sample(df$Father, 100, replace = TRUE))
pred <- predict(lin_reg, newdata = new_data)
#bzw. mit Prognoseintervall
p_pred <- predict(lin_reg, newdata = new_data, interval = 'prediction')
#oder Konfidenzintervall
c_pred <- predict(lin_reg, newdata = new_data, interval = 'confidence')
fig <- plot_ly() %>%
add_markers(x = new_data$Father, y = as.vector(p_pred[,1]),
marker = list(color = 'rgba(0,100,80,1)'),
showlegend = FALSE, name = 'Punktprognose') %>%
add_lines(x = new_data$Father, y = as.vector(p_pred[,2]),
line = list(color = 'rgba(0,100,80,1)'),
showlegend = FALSE, name = 'untere Grenze') %>%
add_lines(x = new_data$Father, y = as.vector(p_pred[,3]), fill = 'tonexty', fillcolor='rgba(0,100,80,0.2)',
line = list(color = 'rgba(0,100,80,1)'),
showlegend = FALSE, name = 'Low 2014') %>%
layout(title = "Prognosen inkl. Prognoseintervall",
xaxis = list(title = "Größe des Vaters"),
yaxis = list(title = "Größe des Kindes"))
figEin praktisch wichtiger Tipp ist, dass man die Werte innerhalb eines Objektes mittels der str-Funktion heraus finden und auf diese zugreifen kann.
## List of 12
## $ coefficients : Named num [1:2] 39.11 0.399
## ..- attr(*, "names")= chr [1:2] "(Intercept)" "Father"
## $ residuals : Named num [1:898] 2.74 -1.26 -1.46 -1.46 4.24 ...
## ..- attr(*, "names")= chr [1:898] "1" "2" "3" "4" ...
## $ effects : Named num [1:898] -2000.59 -29.55 -1.4 -1.4 4.25 ...
## ..- attr(*, "names")= chr [1:898] "(Intercept)" "Father" "" "" ...
## $ rank : int 2
## $ fitted.values: Named num [1:898] 70.5 70.5 70.5 70.5 69.3 ...
## ..- attr(*, "names")= chr [1:898] "1" "2" "3" "4" ...
## $ assign : int [1:2] 0 1
## $ qr :List of 5
## ..$ qr : num [1:898, 1:2] -29.9666 0.0334 0.0334 0.0334 0.0334 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:898] "1" "2" "3" "4" ...
## .. .. ..$ : chr [1:2] "(Intercept)" "Father"
## .. ..- attr(*, "assign")= int [1:2] 0 1
## ..$ qraux: num [1:2] 1.03 1.12
## ..$ pivot: int [1:2] 1 2
## ..$ tol : num 1e-07
## ..$ rank : int 2
## ..- attr(*, "class")= chr "qr"
## $ df.residual : int 896
## $ xlevels : Named list()
## $ call : language lm(formula = Height ~ Father, data = df)
## $ terms :Classes 'terms', 'formula' language Height ~ Father
## .. ..- attr(*, "variables")= language list(Height, Father)
## .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. ..$ : chr [1:2] "Height" "Father"
## .. .. .. ..$ : chr "Father"
## .. ..- attr(*, "term.labels")= chr "Father"
## .. ..- attr(*, "order")= int 1
## .. ..- attr(*, "intercept")= int 1
## .. ..- attr(*, "response")= int 1
## .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. ..- attr(*, "predvars")= language list(Height, Father)
## .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
## .. .. ..- attr(*, "names")= chr [1:2] "Height" "Father"
## $ model :'data.frame': 898 obs. of 2 variables:
## ..$ Height: num [1:898] 73.2 69.2 69 69 73.5 72.5 65.5 65.5 71 68 ...
## ..$ Father: num [1:898] 78.5 78.5 78.5 78.5 75.5 75.5 75.5 75.5 75 75 ...
## ..- attr(*, "terms")=Classes 'terms', 'formula' language Height ~ Father
## .. .. ..- attr(*, "variables")= language list(Height, Father)
## .. .. ..- attr(*, "factors")= int [1:2, 1] 0 1
## .. .. .. ..- attr(*, "dimnames")=List of 2
## .. .. .. .. ..$ : chr [1:2] "Height" "Father"
## .. .. .. .. ..$ : chr "Father"
## .. .. ..- attr(*, "term.labels")= chr "Father"
## .. .. ..- attr(*, "order")= int 1
## .. .. ..- attr(*, "intercept")= int 1
## .. .. ..- attr(*, "response")= int 1
## .. .. ..- attr(*, ".Environment")=<environment: R_GlobalEnv>
## .. .. ..- attr(*, "predvars")= language list(Height, Father)
## .. .. ..- attr(*, "dataClasses")= Named chr [1:2] "numeric" "numeric"
## .. .. .. ..- attr(*, "names")= chr [1:2] "Height" "Father"
## - attr(*, "class")= chr "lm"## (Intercept) Father
## 39.1103868 0.39938132.6 Interpretation des Outputs
Betrachten wir kurz die Zusammenfassung des Modells
##
## Call:
## lm(formula = Height ~ Father, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -10.2683 -2.6689 -0.2092 2.6342 11.9329
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.11039 3.22706 12.120 <2e-16 ***
## Father 0.39938 0.04658 8.574 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.446 on 896 degrees of freedom
## Multiple R-squared: 0.07582, Adjusted R-squared: 0.07479
## F-statistic: 73.51 on 1 and 896 DF, p-value: < 2.2e-16Mit den Resiuals sind die geschätzten Residuen, also gemäß unserer Notation \(\hat{\epsilon}\). Hier werden Minimum, Maximum, Median und das erste bzw. dritte Quartil angegeben, um einen Überblick zu bekommen, wie hoch die Bandbreite der Residuen ist.
Bei der Koeffizienten werden in den ersten beiden Spalten Punktschätzer und deren Standardabweichung angegeben. Der t-Wert ist: \(\frac{\text{Punktschätzer}}{\text{Standardabweichung}}\) und der p-value wir über den t-Wert und die Student t Verteilung bestimmt.
Der Residual standard error ist die geschätzte Standardabweichung der Residuen: \(\hat{\sigma}\)
2.7 Diagnose des Regressionsmodells
Nicht alle Daten sehen so schön aus, wie die aus dem Galton Datensatz. Oft kommt es zu Fehlern in den Daten in Form von Ausreissern, falsch eingegbenen Daten. Zudem sind die Annahmen des Regressionsmodells selten alle erfüllt. Die im folgenden Aspekte sind für Sie zunächst möglicherweise etwas viel, sollen Sie jedoch sensibilieren, wie theoretisch in der Praxis mit dem linearen Regressionsmodell umgegangen werden sollte.
2.7.1 Diagnose von ungewöhnlichen Datenpunkten
Ausreisser oder Datenfehler können die Schätzung des Regressionsmodells in Extremfällen enorm verändern. Daher gehört es eigenlich zur guten Praxis das Regressionsmodell auf Besonderheiten in den Daten zu untersuchen. Hierzu konstruieren wir uns drei Beispieldatensätze. Die Datensätze bestehen aus elf Beobachtungen, die Werte für \(\boldsymbol{x}\) sind: \(\boldsymbol{x} = (1,2,3,4,5,6,7,8,9,10,50)\). Die Gleichung des Regressionsmodell lautet:
\[ y_i = 1 + x_i + \epsilon_i \]
Der simulierte Datensatz sieht folgendermaßen aus:
set.seed(123)
x <- c(1:10, 50)
eps <- rnorm(11)
y <- 1 + x + eps
kable(data.frame(x = x, y = round(y, 2)))| x | y |
|---|---|
| 1 | 1.44 |
| 2 | 2.77 |
| 3 | 5.56 |
| 4 | 5.07 |
| 5 | 6.13 |
| 6 | 8.72 |
| 7 | 8.46 |
| 8 | 7.73 |
| 9 | 9.31 |
| 10 | 10.55 |
| 50 | 52.22 |
Der letzte Punkt \((50, 52.22)\) ist eindeutig ein extremer Punkt. Dies muss an sich kein Fehler sein. Folgende drei Datensätze werden erzeugt.
- Datensatz A folgt dem Modell
- Bei Datensatz B wurde der letzte Punkt falsch eingetragen mit: \((50, 5.22)\)
- Bei Datensatz C wurde der letzte Punkt falsch eingetragen mit: \((5, 52.22)\)
Die Datensätze mit der geschätzten Regressionslinie sehen folgendermaßen aus:
df_A <- data.frame(x = x, y = y)
df_B <- data.frame(x = x, y = c(y[1:10], 5.22))
df_C <- data.frame(x = c(x[1:10], 5.5), y = y)
t1 <- list(
text = "Datensatz A",
xref = "paper",
yref = "paper",
yanchor = "bottom",
xanchor = "center",
align = "center",
x = 0.5,
y = 1,
showarrow = FALSE
)
t2 <- list(
text = "Datensatz B",
xref = "paper",
yref = "paper",
yanchor = "bottom",
xanchor = "center",
align = "center",
x = 0.5,
y = 1,
showarrow = FALSE
)
t3 <- list(
text = "Datensatz C",
xref = "paper",
yref = "paper",
yanchor = "bottom",
xanchor = "center",
align = "center",
x = 0.5,
y = 1,
showarrow = FALSE
)
p1 <- plot_ly(data = df_A) %>%
add_markers(x = ~x, y = ~y, marker = list(color = the_green)) %>%
add_lines(x = x, y = lm(y ~ x, data = df_A)$coef[1] + lm(y ~ x, data = df_A)$coef[2] * x, line = list(color = the_grey)) %>%
layout(annotations = t1,
xaxis = list(title = 'x'),
yaxis = list(title = 'y', range = c(0,60)))
p2 <- plot_ly(data = df_B) %>%
add_markers(x = ~x, y = ~y, marker = list(color = the_green)) %>%
add_lines(x = x, y = lm(y ~ x, data = df_B)$coef[1] + lm(y ~ x, data = df_B)$coef[2] * x, line = list(color = the_grey)) %>%
layout(annotations = t2,
xaxis = list(title = 'x'),
yaxis = list(title = 'y', range = c(0,60)))
p3 <- plot_ly(data = df_C) %>%
add_markers(x = ~x, y = ~y, marker = list(color = the_green)) %>%
add_lines(x = x, y = lm(y ~ x, data = df_C)$coef[1] + lm(y ~ x, data = df_C)$coef[2] * x, line = list(color = the_grey)) %>%
layout(annotations = t3,
xaxis = list(title = 'x'),
yaxis = list(title = 'y', range = c(0,60)))
subplot(p1, p2, p3, nrows = 2,
titleX = TRUE, titleY = TRUE, margin = 0.1) %>% layout(showlegend = FALSE)Bei den Datensätzen B und C kann man gut erkennen, wie sehr der jeweils eine veränderte Datenpunkte die Schätzung des Modells verändert. Wichtig ist hierbei die Unterscheidung zwischen einer hohen aber an sich richtigen Realisierung wie in Datensatz A und den fehlerhaften Daten in den anderen beiden Datensätzen. Um derartige Punkte zu identifizieren kann man sich an drei Kennzahlen bedienen.
- Leverages
- Studentized resdiuals
- Cook’s Distance - Cook’s D
Leverages
Um zu verstehen was mit leverage bzw. leverages gemeint ist betrachten wir die Schätzung aus den Slides der Vorlesung:
\[ \hat{\boldsymbol{\beta}} = (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T \boldsymbol{y} \] Somit werden die die Werte für \(\boldsymbol{y}\) durch: \(\hat{\boldsymbol{y}} = \boldsymbol{X}\hat{\boldsymbol{\beta}}\) geschätzt. Durch Ausklammer sieht man,
\[ \hat{\boldsymbol{y}} = \underbrace{\left(\boldsymbol{X} (\boldsymbol{X}^T \boldsymbol{X})^{-1} \boldsymbol{X}^T \right)}_{\boldsymbol{H}} \boldsymbol{y} \]
dass die Prognosen \(\hat{\boldsymbol{y}}\) von den Realisierungen \(\boldsymbol{y}\) abhängen und hierbei mit den Werten in der Matrix \(\boldsymbol{H}\) gewichtet werden. Die Matrix \(\boldsymbol{H}\) wird auch als “Hat-Matrix” bezeichnet, wobei die Diagonalelemente \(H_{ii}\) als “leverage” bezeichnet werden. Je größer die leverage \(H_{ii}\), umso stärker hängt die Prognose für \(\hat{y}_i\) von ihrer eigenen Realisierung ab. Daraus folgt, dass der Schätzer für das Residum \(\hat{\epsilon}_i = Y_i - \hat{Y}_i\) zu kleinen Schätzungen für \(\epsilon_i\) führt und somit nicht gut ist. Mathematisch lässt sich zeigen, dass \(E(H_{ii}) = \frac{(p + 1)}{n}\) ist. Im allgemeinen werden Werte als zu groß angesehen für die \(H_{ii} \geq 2 \cdot \frac{(p + 1)}{n}\). \(p\) ist die die Anzahl der Prediktoren, \(n\) der Datenumfang. Wie bestimmen manuell die Werte für unser Beispiel.
#Für Datensatz A
X <- cbind(rep(1, 11), df_A$x)
y <- df_A$y
H_mat <- X %*% solve(t(X) %*% X) %*% t(X)
critical <- 2/11
diag(H_mat)## [1] 0.12969580 0.12114920 0.11366490 0.10724288 0.10188315 0.09758571
## [7] 0.09435056 0.09217769 0.09106712 0.09101883 0.96016417## [1] 0.7133269 0.6663206 0.6251569 0.5898358 0.5603573 0.5367214 0.5189281
## [8] 0.5069773 0.5008691 0.5006036 5.2809029#Für Datensatz B
X <- cbind(rep(1, 11), df_B$x)
y <- df_B$y
H_mat <- X %*% solve(t(X) %*% X) %*% t(X)
critical <- 2/11
diag(H_mat)## [1] 0.12969580 0.12114920 0.11366490 0.10724288 0.10188315 0.09758571
## [7] 0.09435056 0.09217769 0.09106712 0.09101883 0.96016417## [1] 0.7133269 0.6663206 0.6251569 0.5898358 0.5603573 0.5367214 0.5189281
## [8] 0.5069773 0.5008691 0.5006036 5.2809029#Für Datensatz C
X <- cbind(rep(1, 11), df_C$x)
y <- df_C$y
H_mat <- X %*% solve(t(X) %*% X) %*% t(X)
critical <- 2/11
diag(H_mat)## [1] 0.33636364 0.23939394 0.16666667 0.11818182 0.09393939 0.09393939
## [7] 0.11818182 0.16666667 0.23939394 0.33636364 0.09090909## [1] 1.8500000 1.3166667 0.9166667 0.6500000 0.5166667 0.5166667 0.6500000
## [8] 0.9166667 1.3166667 1.8500000 0.5000000Studentized Residuals
Im Allgemeinen sind die geschätzten Residuen \(\hat{\epsilon}_i = Y_i - \hat{Y}_i\) eine wichtige Informationsquelle für die Modelldiagnose des linearen Regressionsmodells. Unter den Annahmen entsprechenden Bedingungen gilt: \(\epsilon \sim N(0, \sigma_{\epsilon}^2)\), weshalb unter geschätzten Bedinungen Grenzen wie \(\mp 2 \hat{\sigma_{\epsilon}}\) oder \(\mp 2 \hat{\sigma_{\epsilon}}\) in Frage kommen, um nach extremen Werten zu suchen. Jedoch sind gerade beim Vorhandensein von Ausreissern und Datenfehlern diese Bedingungen nicht gegeben, insbesondere da die Schätzungen \(\hat{\epsilon}_i\) aufgrund der Ausreisser verzerrt sind. Um die geschätzten Residuen \(\hat{\epsilon}_i\) bzgl. dieser Verzerrung zu bereinigen, werden sie mit Hilfe ihres Standardfehlers \(\hat{\sigma_{\epsilon}} \sqrt{1 - H_{ii}}\) normalisiert. Da jedoch der Schätzer \(\hat{\sigma_{\epsilon}}\) ebenso durch die Ausreisser verzerrt wird, wird dieser noch einmal ohne den entsprechenden Punkt \(i\) geschätzt. Dieser Schätzer folgt der Notation: \(\hat{\sigma_{\epsilon, (-i)}}\). Das Resultat sind sogenannte “studentized” Residuals:
\[ \frac{\hat{\epsilon}_i}{\hat{\sigma_{\epsilon, (-i)}} \sqrt{1 - H_{ii}}} \]
Im Kern sind dies geschätzte Residuen, die so normalisiert werden, dass sie nicht von dem \(i\)-ten Ausreisser verzerrt werden. Das Schöne ist, dass wir diese in R nicht manuell bestimmen müssen, sondern dafür eine Funktion existiert.
## 1 2 3 4 5 6
## -0.592445809 -0.259080373 1.869200320 0.008633285 0.047066831 2.013312439
## 7 8 9 10 11
## 0.343548514 -1.662227126 -0.910792785 -0.660975327 0.980212019## 1 2 3 4 5 6
## -2.09118361 -1.35822503 -0.27922350 -0.45109984 -0.09149556 0.80410275
## 7 8 9 10 11
## 0.70398323 0.44076800 1.02016593 1.56781533 -8.62418247## 1 2 3 4 5 6 7
## -0.4146901 -0.3553723 -0.2044831 -0.2970066 -0.2832238 -0.1684848 -0.2524359
## 8 9 10 11
## -0.3787952 -0.3464609 -0.3449693 44.5591356Cook’s D
Keine der beiden vorherigen Methoden ist in der Lage, zu quantifizieren, wie stark das Modell durch den jeweiligen Datenpunkt beinflusst wird. Cook’s D hingegen misst den Einfluss. Hierfür wird das Modell \(i = 1, ..., n\) jeweils unter Wegnahme des \(i\)-ten Datenpunktes geschätzt. Der Schätzer \(\hat{Y}_j(-i)\) bezeichnet die Prognose für den \(j\)-Wert, wenn beim Schätzen des Modells der \(i\)-te Wert ignoriert wurde. Cook’s D für die \(i\)-te Beobachtung bestimmt sich durch:
\[ \frac{\sum_{j = 1}^n \left( \hat{Y}_j - \hat{Y}_j(-i) \right) }{(p + 1)s^2} \]
wobei \(s^2\) ein Schätzer für \(\sigma_{\epsilon}^2\) ist.
## 1 2 3 4 5 6
## 2.818555e-02 5.161409e-03 1.754223e-01 5.036233e-06 1.413194e-04 1.636457e-01
## 7 8 9 10 11
## 6.815890e-03 1.172965e-01 4.235876e-02 2.333335e-02 1.162992e+01## 1 2 3 4 5 6
## 0.237014517 0.116239994 0.005569842 0.013408935 0.000533627 0.036389110
## 7 8 9 10 11
## 0.027348141 0.010832957 0.051901512 0.105907861 97.930185136## 1 2 3 4 5 6
## 0.047996777 0.022011050 0.004679544 0.006577543 0.004631690 0.001649667
## 7 8 9 10 11
## 0.004765964 0.015857733 0.020937094 0.033431008 0.4481941462.7.2 Verletzung der Annahmen
Die Annahmen des linearen Regressionsmodell sind sehr restriktiv und in der Realität entsprechend oft verletzt. Insbesondere die Annahme homoskedastischer unabhängiger Fehlerterme muss oft für reale Vorgänge negiert werden. Die gute Nachricht ist, dass dies keine Auswirkung auf die Erwartungstreue der geschätzten Regressionsgeraden hat, also weitherhin \(E(\hat{\boldsymbol{\beta}}) = \boldsymbol{\beta}\) gilt. Die schlechte Nachricht ist, dass der OLS Schätzer nicht mehr effizient ist, was bedeutet, dass es effizientere, insbesondere weniger volatile Schätzer für \(\boldsymbol{\beta}\) gibt. Besonders negative Konsequenzen kann es haben, wenn man dennoch den OLS Schätzer des linearen Regressionsmodells verwendet, obwohl die Annahmen verletzt, da die Kovarianzmatrix des Schätzers der Regressionsgeraden verzerrt geschätzt wird. Dies hat unmittelbare Auswirkungen auf alle Größen, die auf der geschätzten Kovarianzmatrix basieren, insbondere dem t-Wert bzw. p-value der einzelnen Koeffizienten. Eine Folge hieraus kann eine Fehlspezifikation der Einschätzung bzgl. statistischer Signifikanz der Kovariablen sein.
Um Verletzungen der Annahmen des klassichen linearen Regressionsmodells bei der Schätzung zu berücksichtigen, genügt es bei dieser ein paar Anpassungen vorzunehmen. Die Annahme von homoskedastischen unkorrelierten Störtermen wird durch:
\[ \text{Cov}(\boldsymbol{\epsilon}) = \sigma^2 \boldsymbol{I} = \begin{pmatrix} \sigma^2 & 0 & ... & 0 \\ 0 & \sigma^2 & ... & 0 \\ \vdots & \ddots & & \vdots \\ 0 & ... & & \sigma^2 \end{pmatrix} \]
festgelegt. Bei Heteroskedastizität sind die Werte auf der Diagonalen nicht mehr identisch und bei korrelierten Störtermen sind die Elemente jenseits der Diagonlen von null verschieden (mindestens eins). Somit kann die Kovarianzmatrix durch:
\[ \text{Cov}(\boldsymbol{\epsilon}) = \sigma^2 \boldsymbol{W} = \begin{pmatrix} \sigma_{1}^2 & \sigma_{1,2} & ... & \sigma_{1,n} \\ \sigma_{2,1} & \sigma_{2}^2 & ... & \sigma_{2,n} \\ \vdots & \ddots & & \vdots \\ \sigma_{n, 1} & ... & & \sigma_{n}^2 \end{pmatrix} \]
beschrieben werden, wobei \(\boldsymbol{W}\) eine \(n\times n\) Gewichtsmatrix ist und \(\sigma_{i,j}\) die Kovarianz zwischen Störterm \(i\) und \(j\) ist. Wie bereits erwähnt gilt unter diesen Annahmen \(E(\hat{\boldsymbol{\beta}}) = \boldsymbol{\beta}\), jedoch ist verändert sich \(\text{Cov}(\boldsymbol{\hat{\beta}}) = \sigma^2 \left( \boldsymbol{X}^T \boldsymbol{X} \right)^{-1}\) zu:
\[ \text{Cov}(\boldsymbol{\hat{\beta}}) = \sigma^2 \left( \boldsymbol{X}^T \boldsymbol{X} \right)^{-1} \boldsymbol{X}^T \boldsymbol{W} \boldsymbol{X} \left( \boldsymbol{X}^T \boldsymbol{X} \right)^{-1} \] Um dies zu berücksichtigen wird die Schätzung einfach angepasst. Hierfür werden die ursprünglichen Fehlerterme sowie die unabhängige und abhängige Variable so transformiert, dass die Schätzung des klassischen Regressionsmodells wieder verwendet werden kann. Ohne zu sehr ins Detail zu gehen, erfolgt die Transformation auf Basis der Gewichtsmatrix \(\boldsymbol{W}\), so dass das transformierte Regressionsmodell fogendermaßen aussieht:
\[ \boldsymbol{W}^{-1/2} \boldsymbol{y} = \boldsymbol{W}^{-1/2} \boldsymbol{X} \boldsymbol{\beta} + \boldsymbol{W}^{-1/2} \boldsymbol{\epsilon} \]
Der neue, angepasste Schätzer der Parameter der Regressionsgerade lautet dann:
\[ \hat{\boldsymbol{\beta}} = \left( \boldsymbol{X}^T \boldsymbol{W}^{-1} \boldsymbol{X} \right)^{-1} \boldsymbol{X}^{T} \boldsymbol{W}^{-1} \boldsymbol{y} \] Der erwartungstreue Schätzer für \(\sigma^2\) lautet:
\[ \hat{\sigma}^2 = \frac{1}{n - p -1} \hat{\boldsymbol{\epsilon}}^T \boldsymbol{W}^{-1} \hat{\boldsymbol{\epsilon}} \]
Ein nicht ganz unwesentliches Detail ist jedoch, dass uns die Gewichtsmatrix in der Regel nicht bekannt ist. Wie wir hiermit umgehen, hängt von der Annahmenverletzung aus und wird für jeden einzelnen Fall besprochen.
2.7.2.1 Homoskedastizität
Zunächst befassen wir uns ausschließlich mit Homoskedastizität der Fehlerterme, dies bedeuet an Stelle von \(\text{Cov}(\boldsymbol{\epsilon)} = \sigma^2 \boldsymbol{I}\) die Annahme \(\text{Cov}(\boldsymbol{\epsilon)} = \sigma^2 \boldsymbol{W}\) mit \(\boldsymbol{W} = \text{diag}(w_1, ..., w_n)\) getroffen wird. Dies hat zur Folge, dass \(\text{Var}(\epsilon_i) = \sigma_i^2 = \sigma^2 w_i\) und entsprechend die abhängigen Variablen \(Y_i\) unterschiedlich volatil sind.
Ohne weiter auf die mathematische Handhabung dieses Problem einzugehen, betrachten wir, was im Rahmen einer Simulation mit der geschätzten Regressionsgeraden passiert, wenn die Störterme nicht homo-, sondern heteroskedastisch sind. Im Folgenden werden 100 Mal aus einem Regressionsmodell mit und ohne heteroskedastischen Fehlertermen simuliert und die Regressionsgerade geschätzt. In der folgenden Grafik repräsentieren die grauen Linien die jeweilgen geschätzten Regressionsgeraden, die grüne Linie ist die wahre Regressionsgerade.
x <- 1:10
sigma_1 <- rep(1, 10)
sigma_2 <- 1/3 * x
p1 <- plot_ly()
p2 <- plot_ly()
set.seed(42)
for(i in 1:100){
eps_1 <- rnorm(10, mean = 0, sd = sigma_1)
eps_2 <- rnorm(10, mean = 0, sd = sigma_2)
y_1 <- 1 + x + eps_1
y_2 <- 1 + x + eps_2
lm_1 <- lm(y_1 ~ x)
lm_2 <- lm(y_2 ~ x)
p1 <- p1 %>% add_lines(x = x, y = lm_1$coef[1] + lm_1$coef[2] * x, line = list(color = the_grey))
p2 <- p2 %>% add_lines(x = x, y = lm_2$coef[1] + lm_2$coef[2] * x, line = list(color = the_grey))
}
t1 <- list(
text = "Schätzungen der Regressionsgeraden - homoskedastisch",
xref = "paper",
yref = "paper",
yanchor = "bottom",
xanchor = "center",
align = "center",
x = 0.5,
y = 1,
showarrow = FALSE
)
t2 <- list(
text = "Schätzungen der Regressionsgeraden - hetereoskedastisch",
xref = "paper",
yref = "paper",
yanchor = "bottom",
xanchor = "center",
align = "center",
x = 0.5,
y = 1,
showarrow = FALSE
)
p1 <- p1 %>%
add_lines(x = x, y = 1 + 1 * x, line = list(color = the_green)) %>%
layout(annotations = t1, xaxis = list(title = 'x'), yaxis = list(title = 'y'))
p2 <- p2 %>%
add_lines(x = x, y = 1 + 1 * x, line = list(color = the_green)) %>%
layout(annotations = t2, xaxis = list(title = 'x'), yaxis = list(title = 'y'))
subplot(p1, p2, nrows = 2, titleX = TRUE, titleY = TRUE, margin = 0.1) %>% layout(showlegend = FALSE)Wir können erkennen, dass die Bandbreite der geschätzten Regressionsgerade mit heteroskedastischen Fehlertermen deutlich breiter ist, was wiederum die größere Volatilität der Schätzer repräsentiert. Nun stellen sich die Fragen, wie man Heteroskedastizität identifiert und wie man damit umgeht.
Identifikation von Heteroskedastizität
Hierzu existieren zwei Möglichkeiten:
- graphische Betrachtung der geschätzten Residuen
- Breusch-Pagan-Test
Betrachtet man die geschätzten studentisierten Residuen (nicht studentisierte Residuen sind heteroskedastisch, wenn sie unterschiedliche leverage aufweisen, da \(\hat{\sigma_{\epsilon, (-i)}} \sqrt{1 - H_{ii}}\)) bei Vorliegen von Heteroskedastizität, sollte man erkennen können, dass diese für verschiedene Bereiche unterschiedlich stark schwanken. Am besten plottet man hierzu die geschätzen Werte der abhängigen Variable \(\hat{Y}_i\) und/oder die Kovariablen gegen die studentisierten Residuen.
x <- 1:10
sigma_1 <- rep(1, 10)
sigma_2 <- 1/3 * x
set.seed(123)
eps_1 <- rnorm(10, mean = 0, sd = sigma_1)
eps_2 <- rnorm(10, mean = 0, sd = sigma_2)
y_1 <- 1 + x + eps_1
y_2 <- 1 + x + eps_2
lm_1 <- lm(y_1 ~ x)
lm_2 <- lm(y_2 ~ x)
res_1 <- rstudent(lm_1)
res_2 <- rstudent(lm_2)
t1 <- list(
text = "Homoskedastische Residuen",
xref = "paper",
yref = "paper",
yanchor = "bottom",
xanchor = "center",
align = "center",
x = 0.5,
y = 1,
showarrow = FALSE
)
t2 <- list(
text = "Hetereoskedastische Residuen",
xref = "paper",
yref = "paper",
yanchor = "bottom",
xanchor = "center",
align = "center",
x = 0.5,
y = 1,
showarrow = FALSE
)
p1 <- plot_ly() %>%
add_markers(x = fitted(lm_1), y = res_1, marker = list(color = the_green)) %>%
layout(xaxis = list(title = 'Geschätzte Werte für y'), yaxis = list(title = 'Studentisierte Residuen', range = c(-3, 3)), annotations = t1)
p2 <- plot_ly() %>%
add_markers(x = x, y = res_1, marker = list(color = the_green)) %>%
layout(xaxis = list(title = 'x'), yaxis = list(title = 'Studentisierte Residuen', range = c(-3, 3)), annotations = t1)
p3 <- plot_ly() %>%
add_markers(x = fitted(lm_2), y = res_2, marker = list(color = the_green)) %>%
layout(xaxis = list(title = 'Geschätzte Werte für y'), yaxis = list(title = 'Studentisierte Residuen', range = c(-3, 3)), annotations = t2)
p4 <- plot_ly() %>%
add_markers(x = x, y = res_2, marker = list(color = the_green)) %>%
layout(xaxis = list(title = 'x'), yaxis = list(title = 'Studentisierte Residuen', range = c(-3, 3)), annotations = t2)
subplot(p1, p2, p3, p4, nrows = 2, titleX = TRUE, titleY = TRUE, margin = 0.15) %>% layout(showlegend = FALSE)Zusätzlich kann man einen statistischen Test, den Breusch-Pagan-Test, zu Hilfe ziehen. Die Ablehnung der Nullhypothese dieses Tests spricht für das Vorliegen heteroskedastischer Fehlerterme.
##
## studentized Breusch-Pagan test
##
## data: y_1 ~ x
## BP = 0.24532, df = 1, p-value = 0.6204##
## studentized Breusch-Pagan test
##
## data: y_2 ~ x
## BP = 1.9472, df = 1, p-value = 0.1629Wie man erkennen kann schlägt der Test in unserem Beispiel nicht an. Dies sollte Ihnen vor Augen führen, dass man sich selbstverständlich auch mit der Verwendung des Tests nicht sicher sein kann.
Umgang mit Heteroskedastiztät
Nun stellt sich die Frage, was zu tun ist, wenn man den Verdacht hat, dass die zugrunde liegenden Daten heteroskedastisch sind. Zunächst noch einmal ein kleiner allgemeiner Exkurs. Eine Idee ist, die abhängige Variable, Kovariablen und die Störgrößen so zu transformieren, dass diese Werte wieder den Annahmen des klassischenn linearen Regressionmodell genügen. Multiplizieren wir die Residuen \(\epsilon_i\) mit \(1/\sqrt{w_i}\) haben die transformierten Residuen \(\epsilon_i^{*} = \epsilon_i / \sqrt{w_i}\) wieder konstante Varianz \(\text{Var}(\epsilon_i^{*}) = \sigma^2\). Die übrigen Variablen müssen entsprechend mit angepasst werden, also \(y_i^{*} = y_i / \sqrt{w_i}\) und \(x_{ij}^{*} = x_{ij} / \sqrt{w_i},~\text{für }j=1, ..., p\). Für die transformierten Variablen kann wie bereits erwähnt nun das klassischen Regressionmodell geschätzt werden, da die Annahme der Homoskedastiztät wieder erfüllt ist. Es lässt sich zeigen, dass der Schätzer \(\hat{\boldsymbol{\beta}}\) aus diesem Modell die Summe der gewichteten quadratischen Abweichungen minimiert:
\[ GKQ(\boldsymbol{\beta}) = \sum_{i = 1}^{n} \frac{1}{w_i} (y_i - \underbrace{\boldsymbol{x_i}^T \boldsymbol{\beta}}_{\mu_i})^2\]
Das bedeutet, dass bei der Schätzung Beobachtungen mit höherer Varianz (höheres \(w_i\)) ein kleineres Gewicht (\(1/w_i\)) bekommen. Man spricht hier von der gewichteten Methoden der kleinesten Quadrate oder im englischen auch “weighted least squared”, falls Sie einmal auf diesen Ausdruck stoßen, wissen Sie was gemeint ist.
Ein Problem ist jedoch, dass die Art der Heteroskedastizität und somit die Gewichte \(w_i\) normalerweise unbekannt sind. Daher ist eine Möglichkeit die Gewichte in einem beispielsweise zweistufigen Verfahren mit zu schätzen oder auf deren Schätzung zu verzichten und einfach die geschätze Kovarianz der Regressionskoeffizienten \(\text{Cov}(\hat{\boldsymbol{\beta}})\) anuzupassen, da dies ja das eigentliche Problem beim Vorliegen von Heteroskedastizität ist. Eine bekannte Art die geschätzte Kovarianz anzupassen ist der “White” Schätzer. Hierbei wird die Kovarianz an Stelle von \(\text{Cov}(\boldsymbol{\hat{\beta}}) = \sigma^2 \left( \boldsymbol{X}^T \boldsymbol{X} \right)^{-1}\) durch
\[ \widehat{\text{Cov}(\boldsymbol{\hat{\beta}})} = \left( \boldsymbol{X}^T \boldsymbol{X} \right)^{-1} \boldsymbol{X}^T \text{diag}(\epsilon_1^2, ..., \epsilon_n^2) \boldsymbol{X} \left( \boldsymbol{X}^T \boldsymbol{X} \right)^{-1} \]
geschätzt. In dem folgenden Beispiel kann man sehen, wie es zu einer unterschiedlichen Inferenz bzgl. des Einflusses der Kovariablen kommen kann, wenn die Heteroskedastizität bei Vorliegen nicht beim Schätzen berücksichtigt wird.
#Anzahl der Beobachtungen
n <- 10
#deterministische Kovariable
x <- 1:n
#heteroskedastische Standardabwweichung
sigma <- rep(c(1,2), n / 2)
#heteroskedastische Fehlerterme
set.seed(42)
eps <- rnorm(n, mean = 0, sd = sigma)
#abhängige Variable
y <- 1 + 0.23 * x + eps
#geschätztes Modell
lm_mod <- lm(y ~ x)
#Zusammenfassung des geschätzten Modells
sum_mod <- summary(lm_mod)
#wir rechnen einmal manuell
X <- cbind(rep(1, n), x)
print('Kovarianzschätzung des klassischen Regressionsmodells')## [1] "Kovarianzschätzung des klassischen Regressionsmodells"## x
## 0.50260112 -0.07180016
## x -0.07180016 0.01305457## [1] "Kovarianzschätzung mittels White-Schätzer"cov_heterosked <- solve(t(X) %*% X) %*% t(X) %*% diag(residuals(lm_mod)^2) %*% X %*% solve(t(X) %*% X)
cov_heterosked## x
## 0.51294108 -0.07639425
## x -0.07639425 0.01367127## (Intercept) x
## (Intercept) 0.50260112 -0.07180016
## x -0.07180016 0.01305457## (Intercept) x
## (Intercept) 0.51294108 -0.07639425
## x -0.07639425 0.01367127#Auswirkung auf die Standardfehler der Koeffzienten
sd_beta_0 <- sqrt(cov_homosked[1,1])
sd_beta_1 <- sqrt(cov_homosked[2,2])
t_beta_0 <- coef(lm_mod)[1] / sd_beta_0
t_beta_1 <- coef(lm_mod)[2] / sd_beta_1
p_val_beta_0 <- 2 * (1 - pt(t_beta_0, df = n - 2))
p_val_beta_1 <- 2 * (1 - pt(t_beta_1, df = n - 2))
sd_beta_0_w <- sqrt(cov_heterosked[1,1])
sd_beta_1_w <- sqrt(cov_heterosked[2,2])
t_beta_0_w <- coef(lm_mod)[1] / sd_beta_0_w
t_beta_1_w <- coef(lm_mod)[2] / sd_beta_1_w
p_val_beta_0_w <- 2 * (1 - pt(t_beta_0_w, df = n - 2))
p_val_beta_1_w <- 2 * (1 - pt(t_beta_1_w, df = n - 2))
print(paste('Der p-value von beta_1 bei klassicher Schätzung lautet: ', p_val_beta_1, sep = ' '))## [1] "Der p-value von beta_1 bei klassicher Schätzung lautet: 0.048361726640914"## [1] "Der p-value von beta_1 bei mit Anpassung gemäß dem White Schätzer lautet: 0.0525409728176802"2.7.2.2 Korrelation der Fehlerterme
Eine weitere Annahmen, die in der Realität und vor allem bei Zeitreihendaten oft verletzt ist, besteht in der Unabhängigkeit der Fehlerterme. Die Annahme der Unabhängigkeit \(\text{Cov}(\boldsymbol{\epsilon})=\sigma^2 \boldsymbol{I}\) impliziert, dass alle Elemente jenseits der Diagonalen null sind
\[ \text{Cov}(\boldsymbol{\epsilon}) = \begin{pmatrix} \sigma^2 & 0 & ... & 0 \\ 0 & \sigma^2 & ... & 0 \\ \vdots & \ddots & & \vdots \\ 0 & ... & & \sigma^2 \end{pmatrix} \]
Bestehen Abhängigkeiten unter den Fehlertermen, so ist mindestens ein Element jenseits der Diagonalen von Null verschieden. Für ein Model mit homoskedastischen, jedoch autokorrelierten Fehlertermen sieht die Kovarianzmatrix der Residuen somit folgendermapen aus:
\[ \text{Cov}(\boldsymbol{\epsilon}) = \sigma^2 \boldsymbol{W} = \begin{pmatrix} \sigma^2 & \sigma_{1,2} & ... & \sigma_{1,n} \\ \sigma_{2,1} & \sigma^2 & ... & \sigma_{2,n} \\ \vdots & \ddots & & \vdots \\ \sigma_{n, 1} & ... & & \sigma^2 \end{pmatrix} \]
Existieren derartige Abhängigkeiten in der Realität, kann man hierzu das Verhalten der geschätzten Residuen auf Korrelation analysieren. Beispielsweise graphisch. Bei einer Abhängigkeitsstruktur unter den Residuen sollte der Plot der Resdiuen einen systematischen Verlauf aufweisen. Darüber hinaus existieren statistische Tests, wie der Durbin-Watson-Test. Liegt Autokorrelation unter den Fehlertermen vor und wird dies beim Schätzen nicht berücksichtig, können wieder Fehleinschätzungen der Kovarianzmatrix der Regressionskoeffizienten resultieren.
Identfikation von Autokorrelation
Wir betrachten zunächst die Residuen, wenn die Fehlerterme tatsächlich korreliert sind.
#unabhängige Variable
x <- 1:100
#erzeuge korrelierte Störterme, die einem AR(1)-Prozess folgen
eps <- 0
u <- rnorm(100)
rho <- 0.7
for(i in 1:100){
eps <- c(eps, rho * eps[i] + u[i])
}
eps <- eps[-c(1)]
#abhängige Variable
y <- 1 + x + eps
plot_ly() %>%
add_lines(x = 1:100, y = residuals(lm(y ~x)), line = list(color = the_green)) %>%
layout(yaxis = list(title = 'Residuen'))Wie können relativ klar erkennen, dass die Realisierung des vorherigem Residums die Realisierung des aktuellen Residums beinflusst. Um dies noch deutlicher zu betrachten, können wir uns auch die empirische Autokorrelation der Residuen ansehen. Bei der empirischen Autokorrelation handelt es sich es sich um den geschätzten Korrelationskoeffizienten der Residuen mit ihren um \(j\) verzögerten Einheiten:
\[ \widehat{\text{ACF}}(j) = \frac{\widehat{\text{Cov}}(\epsilon_i, \epsilon_{i - j})}{\widehat{\text{Var}(\epsilon_i)}} = \frac{\frac{1}{n} \sum_{i = j+1}^n \hat{\epsilon_i} \hat{\epsilon_{i - j}}}{\widehat{\text{Var}(\epsilon_i)}} \]
Betrachten wir dies hier,
auto <- acf(residuals(lm(y ~ x)), plot = FALSE)
plot_ly() %>%
add_bars(x = as.vector(auto$lag), y = as.vector(auto$acf), width = 0.1, marker = list(color = the_grey), name = 'Autokorrelation') %>%
add_lines(x = c(0, 35), y = c(1.96/sqrt(auto$n.used), 1.96/sqrt(auto$n.used)), line = list(dash = 'dash', color = the_green), name = 'Kritischer Wert') %>%
add_lines(x = c(0, 35), y = c(-1.96/sqrt(auto$n.used), -1.96/sqrt(auto$n.used)), line = list(dash = 'dash', color = the_green), name = 'Kritischer Wert', showlegend = FALSE) %>%
layout(
title = 'Autokorrelation der Residuen',
xaxis = list(title = '')
)so sehen wir, dass die Autokorrelation für aufeinder folgende Realsierungen der Residuen mit steigendem Abstand abnimmt. Alternativ oder zusätzlich kann man den Durbin-Watson-Test verwenden, bei dem eine Ablehnung der Nullhypothese für die Existenz von autokorrelierter Störterme spricht.
##
## Durbin-Watson test
##
## data: y ~ x
## DW = 0.68621, p-value = 7.882e-12
## alternative hypothesis: true autocorrelation is greater than 0Umgang mit Autokorrelation
Um die Autokorrelation bei der Schätzung der Regressionsparameter zu berücksichtigen, muss diese mit modelliert werden. Typischerweise, indem die Fehlerterme von vergangenen Fehlertermen abhängen. Also ganz allgemein formuliert existiert ein funktionaler Zusammenhang zwischen aktuellem und vergangenen Fehlertermen,
\[ \epsilon_i = f(\epsilon_{i - t}, ..., \epsilon_i) + u_i \]
wobei \(u_i\) den ursprünglichen Annahmen des Fehlerterms folgt.
- \(E(u_i) = 0\)
- \(\text{Var}(u_i) = \sigma_{u_i}^2\) für alle \(i = 1, ..., n\)
- \(\text{Cov}(u_i, u_j) = 0\), wenn \(i \neq j\)
Wie wollen als Beispiel für den Fehlerterm ein einfaches Modell betrachten, einen AR(1)-Prozess:
\[ \epsilon_i = \rho \epsilon_{i - 1} + u_i \]
Aus diesen Annahmen lässt sich zeigen, dass \(E(\boldsymbol{\epsilon}) = 0\) und
\[ \text{Cov}(\boldsymbol{\epsilon}) = \sigma^2 \boldsymbol{W} = \frac{\sigma_u^2}{1 - \rho^2} \begin{pmatrix} 1 & \rho & ... & & \rho^{n - 1} \\ \rho & 1 & ... & & \rho^{n - 2} \\ \vdots & \ddots & & & \vdots \\ \rho^{n - 1} & \rho^{n - 2} & \rho^{n - 3} & ... & 1 \end{pmatrix} \]
Somit kann man bei Kenntnis von \(\boldsymbol{W}\) die transformierten Daten durch ein klassisches Regressionsmodell schätzen. Da die Gewichtsmatrix in der Regel unbekannt ist, muss diese mit geschätzt werden. Hierfür schätzt man das Modell iterativ über zwei Stufen.
- Schätze ein ungewichtetes klassisches Regressionsmodell.
- Schätze die Korrelation durch:
\[ \hat{\rho} = \frac{\sum_{i = 2}^n \hat{\epsilon_i} \hat{\epsilon_{i - 1}} }{ \sqrt{ \sum_{i = 2}^n \hat{\epsilon_i}^2 } \sqrt{ \sum_{i = 2}^n \hat{\epsilon_{i-1}^2} } } \] 3. Schätze die Gewichtsmatrix \(\boldsymbol{W}\) durch \(\hat{ \boldsymbol{W}}\), in dem \(\hat{\rho}\) für \(\rho\) verwendet wird und führe mit dieser eine gewichtete Schätzung des Regressionsmodells durch.
Die Schritte 2. und 3. sollten mehrfach hintereinander durchgeführt werden.
#1.Schritt
lm_start <- lm(y ~ x)
lm_current <- lm_start
iter <- 10
for(i in 1:iter){
#2.Schritt
rho_nom <- sum( tail(residuals(lm_current), n = 99) * head(residuals(lm_current), n = 99))
rho_denom <- sqrt(sum(tail(residuals(lm_current), n = 99)^2)) * sqrt(sum(head(residuals(lm_current), n = 99)^2))
#Schätze Korrelation zwischen den Residuen
rho_hat <- rho_nom / rho_denom
#3.Schritt
#geschätzte Gewichtsmatrix
rho_hoch <- c()
for(i in 1:99 ){
rho_hoch <- c(rho_hoch, rho_hat^(i))
}
offdiag <- c()
for(i in 99:1){
offdiag <- c(tail(rho_hoch, n = 99))
}
diag <- rep(1, length(y))
W_hat <- matrix(NA, ncol = length(diag), nrow = length(diag))
W_hat[lower.tri(W_hat)] <- offdiag
W_hat[upper.tri(W_hat)] <- t(W_hat)[upper.tri(t(W_hat))]
diag(W_hat) <- diag
#transformiere Werte, hierzu wird die negative Wurzel der geschätzten Gewichtsmatrix benötigt
#diese wird hier auf Basis der Eigenwertzerlegung bestimmt, ist jedoch für Sie hier nicht weiter wichtig
ev <- eigen(W_hat)
W_hat_neg_root <- ev$vectors %*% solve(diag(ev$values)) %*% t(ev$vectors)
#transformierte Werte
y_trans <- solve(W_hat_neg_root) %*% as.matrix(y)
X_trans <- solve(W_hat_neg_root) %*% cbind(rep(1, length(x)), x)
eps_trans <- solve(W_hat_neg_root) %*% as.matrix(residuals(lm_current))
#Regressionsmodell der transformierten Werte
lm_current <- lm(y_trans ~ X_trans[,2])
print(coef(lm_current))
print(rho_hat)
}## (Intercept) X_trans[, 2]
## 1.196198 1.012521
## [1] 0.6533764
## (Intercept) X_trans[, 2]
## 1.042533 1.011139
## [1] 0.5408649
## (Intercept) X_trans[, 2]
## 1.088345 1.011645
## [1] 0.5810557
## (Intercept) X_trans[, 2]
## 1.071797 1.011472
## [1] 0.5672809
## (Intercept) X_trans[, 2]
## 1.077461 1.011532
## [1] 0.5720884
## (Intercept) X_trans[, 2]
## 1.075483 1.011512
## [1] 0.5704206
## (Intercept) X_trans[, 2]
## 1.076169 1.011519
## [1] 0.5710004
## (Intercept) X_trans[, 2]
## 1.075931 1.011516
## [1] 0.570799
## (Intercept) X_trans[, 2]
## 1.076014 1.011517
## [1] 0.570869
## (Intercept) X_trans[, 2]
## 1.075985 1.011517
## [1] 0.5708447Man kann erkennen, dass sie die Werte \(\hat{\boldsymbol{\beta}}\) und \(\hat{\rho}\) ab der fünften Iteration einpendeln und sich nicht mehr stark verändern. Jedoch ist es für Sie praktikabler auf bestehende R-Funktionen, wie die gls-Funktion des nlme Pakets zurückzugreifen. Den Erebnissen kann an dieser Stelle auch mehr getraut werden, als die des eben druchgeführten Verfahren, welches numerisch doch etwas ineffizient ist. So einfach geht es dann:
## Generalized least squares fit by REML
## Model: y ~ x
## Data: NULL
## Log-restricted-likelihood: -148.1026
##
## Coefficients:
## (Intercept) x
## 0.5783329 1.0089387
##
## Correlation Structure: AR(1)
## Formula: ~1
## Parameter estimate(s):
## Phi
## 0.6890173
## Degrees of freedom: 100 total; 98 residual
## Residual standard error: 1.4227782.7.3 Verteilung der Störterme
Meist wird mit dem klassischen Regressionsmodell die Annahme normal verteilter Störterme getroffen. Dies ist zumindest für große Stichproben nicht zwingend erforderlich, da der OLS-Schätzer auch bei beliebiger Verteilung asymptotisch normal verteilt ist. Die Folge daraus ist, dass auf die gleiche Art und Weise statistische Tests durchgeführt werden können, wie unter dem Vorliegen normal verteilter Störterme. Dennoch ist hier ein wenig Augenmaß geboten. Bei extremen bzw. sehr schiefen Verteilungen ist der Schätzer der Kovarianzmatrix der geschätzen Parameter tendenziell auch für große Datenmengen verzerrt. Hier sollte man eher ein Modell mit Maximum-Likelihood-Schätzung unter alternativer Verteilungsannahme schätzen oder eine geeignete Transformation der Variablen vor Schätzung des Modells vornehmen. Das Vorliegen nicht normal verteilter Störterme ist oft auch durch Fehlspezifikation des Modells bedingt. So können beispielsweise sogenannte “heavy tails” Folge heteroskedastischer Störterme sein.
2.8 Spurious Regression
Zuletzt noch ein wichtiger Hinweis zur inhaltlichen Interpretation von Regressionsmodellen. Durch Regression wird ein kausaler Zusammenhang zwischen unabhängigen und unabhängiger Variable angenommen. Dennoch liefern Regressionsmodelle auch immer Implikationen eines kausalen Zusammenhangs, wenn unabhängige und abhängige Variable über einen weiteren gemeinsamen Faktor korreliert sind. Dies bedeutet einer der wichtigsten Grundsätze bei der inhaltlichen Auslegung ist:
Kausalität ist nicht gleich Korrelation
Man spricht bei derartige Phänomenen von “spurious correlation”. Ein Beispiel wäre z.B. die Volkswagen und die BMW Aktie. Da es sich um die gleichen Geschäftsfelder handelt, werden tendenziell das Geschäftfeld betreffende Entwicklungen zu gleichläufigen Veränderungen der Aktienkurse führen, wobei nicht zwingend Entscheidungen oder Entwicklungen in einem der beiden Unternehmen, wie beispielsweise die neue Besetzung des CEOs, direkten Einfluss auf die Entwicklung des anderen Unternehmens haben muss. Ein Weg mit diesem Effekt umzugehen, ist weitere Variablen im Modell mit aufzunehmen, die Entwicklungen des Geschäftsfeldes beider Unternehmen betreffen. Man sagt hierbei: es wird für die Variablen kontrolliert. Dennoch bleibt natürlich das Problem bestehen, dass man auch hier eine signifikante Variable ignoriert, weshalb “spurious correlation” weiter ein Thema bleibt.
Eine zum Teil sehr unterhaltsame Liste solcher Scheinkorrelation kann hier betrachtet werden.