Część I

Przygotowanie danych

Ustalam wymagany seed, ładuję biblioteki

rm(list = ls())
set.seed(82338)
library(car)
if(!require("olsrr")) install.packages("olsrr")
if(!require("L1pack")) install.packages("L1pack")

Generuję kolejno 200 elementow \(X \sim N(3+8, 8+2)\) oraz \(e_1 \sim N(0, 3)\)

x <- rnorm(200, mean = (3+8), sd = (8+2))
e_1 <- rnorm(200, 0, 3)

Na ich podstawie tworzę y o wzorze \(y=0.7*X + e_1\)

y <- 3 + 0.7*x + e_1
plot(y~x,col="grey", pch=20)
abline(lm(y~x), col="blue")

X i Y “surowe”

Podstawiam nowe wartosci dla 75, 100 i 125 obserwacji

x1 <- x
x1[c(75, 100, 125)] <- c(-30, 0, 30)
plot(y~x1,col="grey", pch=20)
abline(lm(y~x1), col="blue")

X i Y “zmanipulowane”

Na wykresie są widoczne szczególnie dwie obserwacje podejrzane o bycie odstającymi. Przeprowadzam formalne wnioskowanie dwoma metodami z wykładu.

Część II - Analiza

QQ plot

Najpierw szacuję regresję liniową KMNK

lm1 <- lm(y~x1)
summary(lm1)
## 
## Call:
## lm(formula = y ~ x1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.4277  -2.0392   0.0317   1.8545  28.6552 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.04595    0.40749   9.929   <2e-16 ***
## x1           0.62343    0.02952  21.120   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.192 on 198 degrees of freedom
## Multiple R-squared:  0.6926, Adjusted R-squared:  0.691 
## F-statistic: 446.1 on 1 and 198 DF,  p-value: < 2.2e-16

Widać, że oszcowania parametrów znacząco odbiegają od ich realnych wartosci \(\hat{\beta_1} = 0.6234\), \(\hat{\beta_0} = 4.04595\) A powinny miec kolejno 0.7 oraz 3

Wykres kwantyl-kwantyl

qqPlot(lm1, pch=20, type = "o")

## [1]  75 100
QQ plot

Wyraźnie widać że reszty z oszacowania dla obserwacji 100 i 75 znacząco wykraczają poza przedział dla wartości teoretycznych

Studentyzowane reszty na tle dźwigni

ols_plot_resid_lev(lm1)

Na wykresie Studentyzowanych reszt na tle dźwigni widać, że wszystkie “zmanipulowane” obserwacje w istocie wykraczają poza przedział \([-2;2]\), Co ciekawe równiez obserwacja 65 wyracza poza niego, chociaż na nią nie wpływałem. Spośród obserwacji odstąjacych, najbardziej wpływowa jest obserwacja 75.

Część III - Inna metoda estymacji uwzględniająca możliwość występowania obserwacji odstąjacych

LAD - Least absolute diviations

Polega ona na dopasowaniu parametrów do mediany rozkładu, a nie średniej jak w przypadku KMNK

lad1 <-lad(y~x1)

Porównanie wyników LAD i KMNK

summary(lad1)
## Call:
## lad(formula = y ~ x1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -19.080  -2.266   0.000   1.527  31.550 
## 
## Coefficients:
##              Estimate Std.Error Z value p-value
## (Intercept)  3.4247   0.5228    6.5508  0.0000
## x1           0.6992   0.0379   18.4640  0.0000
## 
## Degrees of freedom: 200 total; 198 residual
## Scale estimate: 3.802652
## Log-likelihood: -536.4545 on 3 degrees of freedom
summary(lm1)
## 
## Call:
## lm(formula = y ~ x1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -17.4277  -2.0392   0.0317   1.8545  28.6552 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.04595    0.40749   9.929   <2e-16 ***
## x1           0.62343    0.02952  21.120   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.192 on 198 degrees of freedom
## Multiple R-squared:  0.6926, Adjusted R-squared:  0.691 
## F-statistic: 446.1 on 1 and 198 DF,  p-value: < 2.2e-16

Oszacowanie metodą LAD jest dużo bliższe rzeczywistości. Parametr dla \(x_1\) jest prawie doskonałym dopasowaniem, a oszacowanie wyrazu wolnego jest co prawda błędne punktowo, ale jest w odleglości już jednego odchylenia standardowego od prawdziwej wartości.