Regresja liniowa

Przykładowy zbiór danych: Boston z biblioteki MASS (zrodlo:http://lib.stat.cmu.edu/datasets/boston)

library(MASS)
head(Boston)
##      crim zn indus chas   nox    rm  age    dis rad tax ptratio  black
## 1 0.00632 18  2.31    0 0.538 6.575 65.2 4.0900   1 296    15.3 396.90
## 2 0.02731  0  7.07    0 0.469 6.421 78.9 4.9671   2 242    17.8 396.90
## 3 0.02729  0  7.07    0 0.469 7.185 61.1 4.9671   2 242    17.8 392.83
## 4 0.03237  0  2.18    0 0.458 6.998 45.8 6.0622   3 222    18.7 394.63
## 5 0.06905  0  2.18    0 0.458 7.147 54.2 6.0622   3 222    18.7 396.90
## 6 0.02985  0  2.18    0 0.458 6.430 58.7 6.0622   3 222    18.7 394.12
##   lstat medv
## 1  4.98 24.0
## 2  9.14 21.6
## 3  4.03 34.7
## 4  2.94 33.4
## 5  5.33 36.2
## 6  5.21 28.7

Założenia regresji liniowej

Przypuśćmy, ze wpływy każdej rozpatrywanej zmiennej objaśniającej na zmienną objaśnianą jest liniowy i nie zależy od wartości innych zmiennych. Wartość \(y_i\) traktujemy jako wartość zmiennej losowej \(Y_i\).

Liniowy model regresji wielokrotnej:

\(Y_i=\alpha_0+\alpha_1 * x_{1,i} + \alpha_2 * x_{2,i} + \dots + \alpha_{m,i} * x_{m,i} + \epsilon_i\),

gdzie:
\(\alpha_0, \dots, \alpha_m\) nieznane parametry,
\(\epsilon_i\) dla \(i=1,\dots, n\) błędy losowe o takim samym rozkładzie mającym wartość średnią 0 i nieznaną, stałą wariancję \(\sigma^2\)

\(n\) liczba obserwacji
\(m\) liczba zmiennych niezależnych
\(p=m+1\) liczba parametrów modelu

Budowa modelu

Podzial zbioru na ciag uczący i testowy

train.idx <-sample (1:nrow(Boston), 0.7*nrow(Boston), replace = FALSE)
train<-Boston[train.idx, ]
test<-Boston[-train.idx, ]

nrow(train)
## [1] 354
nrow(test)
## [1] 152
#
#model$fitted.values
#plot(model$fitted.values, df$medv)
#plot(model)
#model.matrix.default()

Model regresji liniowej

model<-lm(medv~age, data=train)
summary(model)
## 
## Call:
## lm(formula = medv ~ age, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.083  -5.116  -1.774   2.160  31.341 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 30.93089    1.19615  25.859  < 2e-16 ***
## age         -0.12272    0.01601  -7.665 1.76e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 8.35 on 352 degrees of freedom
## Multiple R-squared:  0.143,  Adjusted R-squared:  0.1406 
## F-statistic: 58.75 on 1 and 352 DF,  p-value: 1.757e-13

Model:
* Wartości współczynników (estymatory)
* Błąd standardowy współczynników
* Wartość statystyki t
* p-value dla współczynników

Ocena jakości modelu
* Istotność parametrów
* Wartość współczynnika \(R^2\)

Rozkład całkowitej zmienności zmiennej objaśnianej

  • Całkowita suma kwadratów (total sum of squares)

\(SST=\sum_{i=1}^{n}(y_i-\overline{y})^2\)

y_mean <-mean(train$medv)
SST <- sum((train$medv - y_mean)^2)
SST
## [1] 28635.53

\(SST = SSE + SSR\)
Równość powyższa zachodzi w sytuacji, gdy \(y_i\) oznaczają wartości przewidywane obliczne na podstawie prostej MNK i nie musi być spełniona dla innej prostej przybliżającej chmurę punktów.

  • Suma kwadratów błędów (error sum of squares)
    \(SSE=\sum_{i=1}^{n}(y_i-\hat{y_i})^2\)

Indeks zmienności wartości \(y_i\) nie wyjaśnionej przez zależność liniową.

SSE <- sum ((train$medv - model$fitted.values )^2)
SSE
## [1] 24539.58
SSE <-sum(model$residuals^2)
SSE
## [1] 24539.58
  • Regresyjna suma kwadratów (regression sum of squares)
    \(SSR=\sum_{i=1}^{n}(\hat{y_i}-\overline{y})^2\)
SSR <- sum (( model$fitted.values - y_mean )^2)
SSR
## [1] 4095.951
SSE + SSR 
## [1] 28635.53
  • Współczynnik determinacji \(R^2\)
    \(R^2=\frac{SSR}{SST}=1-\frac{SSE}{SST}\)
R2<- SSR/SST
R2
## [1] 0.1430374
  • Skorygowany współczynnik determinacji - uwzględnia liczbę zmiennych w modelu
n<-nrow(train)
p<-2
R2_adj<- 1 - (1-R2)*(n-1)/(n-p)
R2_adj
## [1] 0.1406028
  • MSR

\(MSR=\frac{SSR}{m}\)

MSR<- SSR/(1)
MSR
## [1] 4095.951
  • Błąd średniokwadratowy, MSE (Mean Square Error), SE (Standard Error)

Nieobciążony estymator wariancji \(\sigma^2\)

\(MSE=\frac{SSE}{n-p}\)

MSE<- SSE/(n-p)
MSE
## [1] 69.71472
  • Statystyka F

Badanie liniowej zależności pomiędzy zmienną celu a zbiorem zmiennych objaśniających traktyowanych jako całość. Hipoteza zerowa: \(H_{0}: \alpha_1=\alpha_2=\dots=\alpha_m=0\)
Hipoteza alternatywna: \(H_{1}:\) istnieje przynajmniej jeden współczynnik \(\alpha_i \ne 0\)

Statystyka F:

\(F=\frac{MSR}{MSE}\)

ma rozkład F-Snedeckora z parametrami p-1, n-m.

F<-MSR/MSE
F
## [1] 58.75303

Rozkład gęstości i dystrybuanta:

F<-MSR/MSE
F
## [1] 58.75303
par(mfrow=c(1,2))

plot(seq(0,20,0.1), df(seq(0,20,0.1),p-1,n-p), xlab="X", ylab="f(x)" )

plot(seq(0,20,0.1), pf(seq(0,20,0.1),p-1,n-p),xlab="X", ylab="F(x)" )

Obszar krytyczny: \({F: F \ge f_{1-\alpha, p-1, n-p}}\)

na poziomie istotności 0.05:

qf(0.95, p-1, n-p)
## [1] 3.868012

Ponieważ F nie wpada do przedziału krytycznego, zatem hipotezę zerową należy odrzucić.

  • p-value modelu

Jaki jest graniczny poziom istotności, przy którym hipoteza zerowa zostałaba odrzucona ?

pf(F, p-1, n-p,     lower.tail=FALSE)
## [1] 1.757375e-13

Weryfikacja hipotezy o istotności współczynników

Współczynniki i błędy standardowe współczynników:

alfa<-model$coefficients
alfa
## (Intercept)         age 
##  30.9308895  -0.1227233
alfa_se<-sqrt(diag(vcov(model)))
alfa_se
## (Intercept)         age 
##  1.19615087  0.01601076

Statystyka t

\(t_i=\frac{\alpha_i}{SE_{\alpha_{i}}}\)

t<-alfa/alfa_se
t
## (Intercept)         age 
##   25.858686   -7.665052

Hipoteza zerowa

Hipoteza zerowa: \(H_{0}: \alpha_i=0\)

Hipoteza alternatywna: \(H_{1}: \alpha_i \neq 0\)

Statystyka testowa t przy ma przy prawdziwości hipotezy H0 rozkład \(T_{n-2}\)

plot(seq(-5,5,0.1), dt(seq(-5,5,0.1), df=n-p), xlab="x", ylab="f(x)", main="Gęstość rozkładu T-Studenta")

Obszar krytyczny:

\(\{{t: |{t}| \ge t_{{1-\frac{\alpha}{2}},n-p}}\}\)

Obszar krytyczny na poziomie istotności 0.05 (granica):

t_gr<-qt(1-0.05/2, df=n-p)
t_gr
## [1] 1.966726

Ponieważ wartość statystyki t wpada do obszru krytycznego, hipotezę zerową należy odrzucić. Oznacza, to że nie ma podstaw do wnioskowania, że parametr jest nieistotny.

p-value

p_value<-2*pt(-abs(t), df=n-p, lower.tail=TRUE)
p_value
##  (Intercept)          age 
## 2.223914e-83 1.757375e-13

Weryfikacja założeń regresji liniowej

Wartości resztowe (residuals)

  • Różnica między wartością zmiennej objaśniającej Y dla i-tej obserwacji i odpowiednią wartoscią przewidywaną:

\(e_i=y_i-\hat{y_i}\)

#wartości prognozowane oraz reszty 

y_pred<- model$fitted.values
head(y_pred)
##      461      110      246      150      334       68 
## 19.88579 19.73853 22.31571 19.28445 26.25513 28.30461
#Reszty 
head(model$residuals)
##        461        110        246        150        334         68 
## -3.4857935 -0.3385256 -3.8157146 -3.8844494 -4.0551322 -6.3046111
y<-train$medv
e<- y - y_pred

head(e)
##        461        110        246        150        334         68 
## -3.4857935 -0.3385256 -3.8157146 -3.8844494 -4.0551322 -6.3046111

Obserwacje odstające (outliers, leverage)

Wpływ (leverage, hat-value):

\(h_i=\frac{1}{n} + \frac{ ( x_i-\overline{x})^2}{\sum_{j=1}^{n} (x_j - \overline{x})^2 }\)

Hat-Matrix H - projekcja Y na wartości przewidywane

\(\hat{y_i}=h_{1,j}* y_1 + h_{2,j}*y_2+ \dots + h_{n,j} *y_n = \sum_{i=1}^n h_{i,j}*y_i\)

age_mean<-mean(train$age)
n<-nrow(train)

h <- 1/n + (train$age - age_mean)^2/sum((train$age - age_mean)^2)
head(h)
## [1] 0.004388676 0.004575964 0.002827347 0.005220099 0.006422025 0.011288801
plot(h)
abline(2/n,0)

Obserwacje odstające:

\(h_i > \frac{2*p}{n}\)

h<-hatvalues(model)
head(h)
##         461         110         246         150         334          68 
## 0.004388676 0.004575964 0.002827347 0.005220099 0.006422025 0.011288801
plot(h, rstandard(model), xlab ="Leverage", ylab="Standardized residuals", main="Residuals vs Leverage")

Obserwacje wpływowe (influential observations)

  • Obserwacja wpływowa to taka, której usunięcie ze zbioru danych powoduje dużą zmianę wektora estymatorow regresji.
  • Obserwacje odstające mogą być, ale nie muszą wpływowymi.
  • Odleglość Cooke’a - wpływ, jaki na prognozę znanych wartości zmiennej objaśnianej ma usunięcie ze zbioru danych i-tej obserwacji:

\(D_i=\frac{\sum_{j=1}^{n}(\hat{Y_j}-\hat{Y_{j(i)}}) }{m*S^2}\),
gdzie
\(\hat{Y_{j(i)}}\) jest wartością przewidywaną dla j-tej obserwacji, obliczoną na podstawie danych z usuniętą obserwacją i-tą.

Diagram Cooke’a

d<-cooks.distance(model)
head(d)
##          461          110          246          150          334 
## 3.858360e-04 3.795725e-06 2.969174e-04 5.708594e-04 7.672266e-04 
##           68 
## 3.292085e-03
plot(d)

Homoskedastyczność czynnika losowego

Rozrzut punktów wokół wartości prognozowanych, sprawdzenie stałości wariancji składnika losowego (wykres scale-location)

plot(model$fitted.values, model$residuals, xlab="Fitted values", ylab=" SQRT standardized residuals")

plot( model$fitted.values, sqrt(rstandard(model)), xlab="Fitted values", ylab="Residuals")
## Warning in sqrt(rstandard(model)): wyprodukowano wartości NaN

Wykresy kwantylowe - weryfikacja rozkładu normalnego składnika losowego

  • Rezyduum ustandaryzowane (Standardized Residuals)

\(e_{i}^{'}=\frac{e_i}{MSE*\sqrt{1-h_i}}\)

e_std <- rstandard(model)
head(e_std)
##         461         110         246         150         334          68 
## -0.41840272 -0.04063733 -0.45764473 -0.46644853 -0.48723873 -0.75938382
e_std <- e/(MSE*sqrt(1-h))
head(e_std)
##          461          110          246          150          334 
## -0.050110907 -0.004867018 -0.054810811 -0.055865217 -0.058355200 
##           68 
## -0.090949245
eq<-quantile(e_std, seq(0.01,0.99, 0.01))
tq<-qnorm(seq(0.01,0.99, 0.01))


plot(tq, eq, xlab="Theoretical quantiles", ylab="Standardized residuals quantiles")

Diagnostyka modelu

par(mfrow=c(2,2))
plot(model)

Regresja wieloraka

Model

model<-lm(medv~age+lstat+rad+crim, data=train)
summary(model)
## 
## Call:
## lm(formula = medv ~ age + lstat + rad + crim, data = train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -10.291  -3.947  -1.444   1.897  24.906 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 33.27788    0.90548  36.752   <2e-16 ***
## age          0.03196    0.01470   2.174   0.0304 *  
## lstat       -0.96199    0.06042 -15.923   <2e-16 ***
## rad         -0.07221    0.05013  -1.440   0.1507    
## crim        -0.05193    0.05227  -0.993   0.3212    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.062 on 349 degrees of freedom
## Multiple R-squared:  0.5522, Adjusted R-squared:  0.547 
## F-statistic: 107.6 on 4 and 349 DF,  p-value: < 2.2e-16

Ocena jakości modelu na zbiorze testowym

y_pred_test<- predict(model, test)

plot(y_pred_test, test$medv, xlab="Actual", ylab="Predicted")
abline(0,1)

*MAE - Mean Absolute Error

\(MAE = \frac{1}{n}{\sum_{j=1}^{n}|y_j-\hat{y_j}|}\)

  • RMSE Root Mean Square Error
    Mean Squared Prediction Error
MAE <- (sum(abs(y_pred_test -  test$medv)) / nrow(test))
MAE 
## [1] 4.591432

\(RMSE = \sqrt {\frac{1}{n}\sum_{j=1}^{n}(y_j-\hat{y_j})^2}\)

RMSE <- sqrt (sum((y_pred_test -  test$medv)^2) / nrow(test))
RMSE 
## [1] 6.404864