Linear Regression

1 Wczytanie i przygotowanie danych

1.1 Wczytanie zbioru danych z pliku płaskiego

autodane<-read.table("/usr/miswdm/autodane.csv", header = TRUE, sep = ",", dec = ".")

1.2 Zmienne pomocnicze

p<-2                # liczba predyktorów
no<-nrow(autodane)
indices<-seq(no)
df<-no-p-1          # liczba stopni swobody

1.3 Zmienne egzogeniczne

wiek<-autodane$Wiek
przebieg<-autodane$Przebieg

1.4 Zmienna endogeniczna

cena<-autodane$Cena

1.5 Nowa obserwacja

Do prognozowania ceny.

new.car<-data.frame(wiek=4, przebieg=145)

2 Obliczanie kowariancji

Kowariancja – miara wielkości zależności liniowej pomiędzy zmiennymi losowymi: \(\begin{aligned} Cov(X, Y)=\mathbb{E}[(X-\mathbb{E}[X])\cdot(Y-\mathbb{E}[Y])]=\mathbb{E}[XY]-\mathbb{E}[X]\cdot\mathbb{E}[Y] \end{aligned}\)

2.1 Wariancja wieku

\(\begin{aligned} Var(X)=\mathbb{E}[(X-\mathbb{E}[X])^2]=\mathbb{E}[X^2]-{\mathbb{E}[X]}^2 \end{aligned}\)

Zgodny, nieobciążony estymator z próbki:
\(\begin{aligned} Var(X)=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})^2 \end{aligned}\)

wiek.mean<-sum(wiek)/no
wiek.diff<-wiek-wiek.mean
wiek.diff.power<-wiek.diff^2
var.wiek<-sum(wiek.diff.power)/(no-1)

[1] 7.066667

2.2 Wariancja przebiegu

Zgodny, nieobciążony estymator z próbki.

przebieg.mean<-sum(przebieg)/no
przebieg.diff<-przebieg-przebieg.mean
przebieg.diff.power<-przebieg.diff^2
var.przebieg<-sum(przebieg.diff.power)/(no-1)

[1] 2540

2.3 Kowariancja zmiennych egzogenicznych

Zgodny, nieobciążony estymator z próbki:
\(\begin{aligned} Cov(X, Y)=\frac{1}{n-1}\sum_{i=1}^n(x_i-\bar{x})\cdot(y_i-\bar{y}) \end{aligned}\)

x.mult<-wiek.diff*przebieg.diff
cov.x<-sum(x.mult)/(no-1)

[1] 119.8889

Obliczenia pomocnicze.

cena.mean<-sum(cena)/no
cena.diff<-cena-cena.mean
wiek.cena.mult<-wiek.diff*cena.diff
przebieg.cena.mult<-przebieg.diff*cena.diff

Kowariancja wieku i ceny.

cov.wiek.cena<-sum(wiek.cena.mult)/(no-1)

[1] -24.66667

Kowariancja przebiegu i ceny.

cov.przebieg.cena<-sum(przebieg.cena.mult)/(no-1)

[1] -432.7778

3 Wyznaczanie wartości współczynników modelu

3.1 Macierz kowariancji zmiennych niezależnych

A<-matrix(c(var.wiek, cov.x, cov.x, var.przebieg), nrow=p)
dimnames(A)<-list(c('wiek', 'przebieg'), c('wiek', 'przebieg'))

Alternatywnie.

A<-cov(cbind(wiek, przebieg))

               wiek  przebieg
wiek       7.066667  119.8889
przebieg 119.888889 2540.0000

3.2 Wektor kowariancji wartości zmiennych niezależnych i zmiennych zależnych (obserwacji)

B<-c(cov.wiek.cena, cov.przebieg.cena)

Alternatywnie.

B<-cov(cbind(wiek, przebieg), cena)

[1]  -24.66667 -432.77778

3.3 Układ równań liniowych

7.0666667 * B₁ + 119.8888889 * B₂ = -24.6666667
119.8888889 * B₁ + 2540 * B₂ = -432.7777778

Rozwiązanie układu.

coef<-solve(A, B)
B1<-coef[1]; B2<-coef[2]

[1] "-3.01122024477396; -0.0282543025323229"

3.4 Wartość wyrazu wolnego B₀

Obliczenia pomocnicze.

aux<-B1*wiek+B2*przebieg
aux.numerator<-cena-aux

Wartość wyrazu wolnego B₀.

B0<-sum(aux.numerator)/no

Alternatywnie

B0<-mean(cena-as.matrix(cbind(wiek, przebieg)) %*% coef)

[1] 46.9232

4 Model regresji liniowej

4.1 Funkcja prognozy wartości zmiennej zależnej

est.fnc<-function(a, b) B0+B1*a+B2*b

Cena = 46.9232017 - 3.0112202 * Wiek - 0.0282543 * Przebieg

4.2 Prognozowana cena nowej obserwacji

result<-est.fnc(new.car$wiek, new.car$przebieg)[[1]]

[1] 30.78145

4.3 Wariancja resztowa

Obliczenia pomocnicze.

cena.est<-est.fnc(wiek, przebieg)
cena.est.diff<-cena-cena.est
cena.est.diff.power<-cena.est.diff^2

Suma kwadratów resztowych (suma kwadratów błędów).

RSS<-sum(cena.est.diff.power)

[1] 93.4586

Nieobciążony estymator wariancji błędu obserwacji \(\sigma^2\):
\(\begin{aligned} S^2=\frac{\sum_{i=1}^n(y_i-\hat{y_i})^2}{df}=\frac{\sum_{i=1}^n(y_i-\hat{y_i})^2}{n-k}=\frac{RSS}{n-p-1}={RSE}^2 \end{aligned}\)
\(df\) - liczba stopni swobody
\(n\) - liczba obserwacji
\(k\) - liczba estymowanych parametrów; \(k=p+1\)

S2<-RSS/df

[1] 13.35123

5 Dopasowanie modelu

5.1 Współczynnik R²

Całkowita suma kwadratów.

cena.diff.power<-cena.diff^2
TSS<-sum(cena.diff.power)

[1] 872

\(\begin{aligned} R^2=\frac{Regresyjna\ suma\ kwadratów}{Całkowita\ suma\ kwadratów} \end{aligned}\)

R2<-(TSS-RSS)/TSS

[1] 0.8928227

5.2 Współczynnik R² dopasowany

R2.adj<-1-(1-R2)*(no-1)/df

[1] 0.8622006

5.3 Statystyka F

F<-df*(TSS-RSS)/(p*RSS)
F.pValue<-pf(F, p, df, lower.tail=FALSE)

[1] 0.0004030507

5.4 Test t-Studenta

5.4.1 Macierz zmiennych niezależnych

Pierwszą kolumną jest wartość 1 dla wyznaczenia wartości błędu parametru wyrazu wolnego.

X<-as.matrix(cbind(rep(1, no), wiek, przebieg))

5.4.2 Estymator macierzy kowariancji estymatora metody najmniejszych kwadratów

Macierz kowariancji parametrów modelu.

Sigma<-S2*solve(t(X)%*%X)  #to samo co Sigma<-S2*solve(t(X)%*%X, diag(2))
colnames(Sigma)[1]<-rownames(Sigma)[1]<-'Intercept'

5.4.3 Błędy współczynników

Błąd współczynnika B₀.

B0.err<-sqrt(diag(Sigma)['Intercept'])  # to samo co B0.err<-sqrt(Sigma[1,1])

Intercept 
 2.962759 

Błąd współczynnika B₁.

B1.err<-sqrt(diag(Sigma)['wiek'])       # to samo co B1.err<-sqrt(Sigma[2,2])

    wiek 
1.026498 

Błąd współczynnika B₂.

B2.err<-sqrt(diag(Sigma)['przebieg'])   # to samo co B2.err<-sqrt(Sigma[3,3])

  przebieg 
0.05414379 

5.4.4 Wartości statystyki t

Dla współczynników B₀, B₁, B₂.

B0.t<-B0/B0.err
B1.t<-B1/B1.err
B2.t<-B2/B2.err

[1] "Intercept: 15.8376706620683; wiek: -2.93348798741743; przebieg: -0.521838285707935"

H₀: nie ma zależności pomiędzy predyktorem a zmienną objaśnianą; B_i = 0.
Przyjęty poziom istotności: 5%.

alpha<-1-.05/2                          # test obustronny
tValue<-qt(alpha, df, lower.tail=TRUE)  # kwatyl rozkładu t-Studenta

Moduł wartości brzegowej obszaru krytycznego wynosi zatem:

[1] 2.364624

p-wartości

Współczynnik	W obszarze krytycznym 5%?	p-wartość
B0 Intercept	TRUE	Pr(>|15.84|)=9.7e-07
B1 wiek	TRUE	Pr(>|-2.933|)=0.0219
B2 cena	FALSE	Pr(>|-0.5218|)=0.618

5.5 Przedziały ufności współczynników modelu

Poziom ufności ciągle 95%.

tValues<-qt(c(1-alpha, alpha), df)
B0.conf<-B0+B0.err*tValues
B1.conf<-B1+B1.err*tValues
B2.conf<-B2+B2.err*tValues

Intercept: [39.917389890249; 53.929013589753]
wiek: [-5.43850286812919; -0.583937621418716]
przebieg: [-0.156284021335215; 0.099775416270569]

5.6 Przedział ufności wyestymowanych wartości modelu

Poziom istotności 5%

cena.est.conf<-data.frame(cena.est)
for (i in indices) {
  C<-c(1, wiek[i], przebieg[i])
  delta<-sqrt(drop(t(C)%*%Sigma%*%C))
  interval<-cena.est[i]+tValues*delta
  cena.est.conf$lwr[i]<-interval[1]
  cena.est.conf$upr[i]<-interval[2]
}

6 Prognoza punktowa i błąd ex-ante

6.1 Wariancja prognozy

\(\begin{aligned} S^2_\tau=S^2+S^2\cdot\left(\frac{1}{n}+\frac{(x_\tau-\bar{x})^2}{\sum_{i=1}^n(x_i-\bar{x})^2}\right)= \end{aligned}\\ \begin{aligned} =S^2\cdot\left(1+\frac{1}{n}+\frac{(x_\tau-\bar{x})^2}{\sum_{i=1}^n(x_i-\bar{x})^2}\right)= \end{aligned}\\ \begin{aligned} =S^2\cdot\left(1+\frac{\sum_{i=1}^nx_i^2+nx_\tau^2-2x_\tau\sum_{i=1}^nx_i}{n\sum_{i=1}^nx_i^2-(\sum_{i=1}^nx_i)^2}\right) \end{aligned}\)
\(\hat{\Sigma}\) - estymator macierzy kowariancji estymatora metody najmniejszych kwadratów
\(C\) - wektor wartości nowej obserwacji (do kombinacji liniowej z modelem regresji)
Standardowy błąd obserwacji:
\(\begin{aligned} \delta_\tau=S\sqrt{\frac{1}{n}+\frac{(x_\tau-\bar{x})^2}{\sum_{i=1}^n(x_i-\bar{x})^2}}=\sqrt{C^T\hat{\Sigma}C}=\sqrt{\delta_\tau^2} \end{aligned}\)

C<-c(1, new.car$wiek, new.car$przebieg)
est.var<-drop(t(C)%*%Sigma%*%C)          # to samo co est.var<-drop(C%*%Sigma%*%C)

[1] 24.9955

6.2 Błąd ex-ante

Błąd prognozy punktowej:
\(\begin{aligned} S_\tau=\sqrt{S^2+\delta_\tau^2}=\sqrt{S^2+C^T\hat{\Sigma}C} \end{aligned}\)

exAnte<-sqrt(S2+est.var)

[1] 6.192473

6.3 Względny błąd ex-ante

Względny błąd prognozy punktowej:
\(\begin{aligned} \eta_\tau=\frac{S_\tau}{|\hat{y_\tau}|}\cdot100\% \end{aligned}\)

exAnte.rel<-exAnte/abs(est.fnc(C[2], C[3]))

     wiek 
0.2011755 

6.4 Przedział ufności prognozy

result.conf<-result+exAnte*tValues

[1] 16.13857 45.42432

7 Wykorzystanie pakietu stats

Funkcja lm z pakietu stats (będąca przeciążeniem funkcji glm z tego samego pakietu) pozwala na stworzenie modelu regresji liniowej z zastosowaniem formuł.

est<-lm(cena~wiek+przebieg)

Współczynniki modelu.

est$coefficients  #to samo co coefficients(est)

(Intercept)        wiek    przebieg 
 46.9232017  -3.0112202  -0.0282543 

Podsumowanie modelu.

summary(est)

Call:
lm(formula = cena ~ wiek + przebieg)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.0194 -0.9033 -0.1010  1.4341  6.1219 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 46.92320    2.96276  15.838  9.7e-07 ***
wiek        -3.01122    1.02650  -2.933   0.0219 *  
przebieg    -0.02825    0.05414  -0.522   0.6179    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.654 on 7 degrees of freedom
Multiple R-squared:  0.8928,    Adjusted R-squared:  0.8622 
F-statistic: 29.16 on 2 and 7 DF,  p-value: 0.0004031

Przedziały ufności współczynników modelu (ang. confidence intervals).

confint(est)

                2.5 %      97.5 %
(Intercept) 39.917390 53.92901359
wiek        -5.438503 -0.58393762
przebieg    -0.156284  0.09977542

Przedziały ufności wyestymowanych wartości modelu (ang. narrow intervals).

predict(est, interval="confidence", level=.95)

        fit       lwr      upr
1  11.72522  5.731752 17.71870
2  27.01935 22.767096 31.27161
3  19.44292 15.775111 23.11074
4  33.46561 30.008315 36.92290
5  18.59529 13.701139 23.48945
6  23.86686 18.207126 29.52659
7  29.04167 24.593028 33.49031
8  26.87808 23.092696 30.66346
9  36.47683 32.178454 40.77520
10 43.48817 37.490825 49.48551

Predykcja z zadanym poziomem tolerancji (ang. wide intervals).

pred<-predict(est, newdata=new.car, se.fit=TRUE, level=.95, interval="prediction")  # błąd.obserwacji=TRUE

$fit
       fit      lwr      upr
1 30.78145 16.13857 45.42432

$se.fit
[1] 4.99955

$df
[1] 7

$residual.scale
[1] 3.653933

Błąd ex-ante (wariancja modelu + wariancja obserwacji).

signif(sqrt(pred$residual.scale^2+pred$se.fit**2), 4)

[1] 6.192

8 Zmienne skorelowane

summary(lm(cena~wiek*przebieg))

Alternatywnie.

summary(lm(cena~wiek+przebieg+wiek:przebieg))

Call:
lm(formula = cena ~ wiek * przebieg)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.0207 -0.9115 -0.0935  1.4333  6.1194 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)   46.9001050  4.9745134   9.428  8.1e-05 ***
wiek          -3.0095093  1.1440731  -2.631    0.039 *  
przebieg      -0.0276481  0.1158114  -0.239    0.819    
wiek:przebieg -0.0000646  0.0106521  -0.006    0.995    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.947 on 6 degrees of freedom
Multiple R-squared:  0.8928,    Adjusted R-squared:  0.8392 
F-statistic: 16.66 on 3 and 6 DF,  p-value: 0.002582

9 Wizualizacja

library(scatterplot3d)
s3d<-scatterplot3d(x=wiek, y=przebieg, z=cena, angle=55, highlight.3d=T, scale.y=0.8)
s3d$plane3d(est)
residual3d<-function(x, y, z, residual){
  s3d$points3d(x=x, y=y, z=z-residual, type='p', lwd=1, pch=5, col=3)
  s3d$points3d(x=c(x, x), y=c(y, y), z=c(z-residual, z), type='l', lwd=1, col='dark green')
}
for (i in indices) {
  residual3d(wiek[i], przebieg[i], cena[i], est$residuals[i])
}
# Predykcja z zadanym poziomem tolerancji (ang. wide intervals).
s3d$points3d(x=new.car$wiek, y=new.car$przebieg, z=pred$fit[,'fit'], lwd=1, pch=21, bg=8, col=2)
s3d$points3d(x=rep(new.car$wiek, 2), y=rep(new.car$przebieg, 2), z=pred$fit[,2:3], lwd=1, type='l', col=2)

9.1 Przedziały ufności (osobny przykład)

x<-1:7
y<-c(8, 13, 14, 17, 18, 20, 22)
est<-lm(y~x)

Call:
lm(formula = y ~ x)

Residuals:
      1       2       3       4       5       6       7 
-1.5714  1.2857  0.1429  1.0000 -0.1429 -0.2857 -0.4286 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   7.4286     0.8806   8.436 0.000384 ***
x             2.1429     0.1969  10.882 0.000114 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.042 on 5 degrees of freedom
Multiple R-squared:  0.9595,    Adjusted R-squared:  0.9514 
F-statistic: 118.4 on 1 and 5 DF,  p-value: 0.0001138

Przedziały ufności.

confint(est)  # to samo co confint(est, level=.95))

               2.5 %   97.5 %
(Intercept) 5.164838 9.692304
x           1.636671 2.649043

t-statystyka, stopni swobody: 5.

res<-signif(residuals(est), 3)
pre<-predict(est, level=.99, interval="confidence")
plot(x, y)
abline(est)
segments(x, y, x, pre[,"fit"], col="red")
library(MASS)
library(calibrate)
textxy(x, y, res, cex=0.8)
# Przedziały ufności (ang. narrow intervals)
# 99%
lines(x, pre[,2], col=8)
lines(x, pre[,3], col=8)
# 95%
pre<-predict(est, level=0.95, interval="confidence")
lines(est$model[["x"]], pre[,2], col=5)
lines(est$model[["x"]], pre[,3], col=5)
# Predykcja z zadanym poziomem tolerancji (ang. wide intervals)
# 99%
result<-predict(est, newdata=data.frame(x=c(4.5)), level=.99, interval="prediction")
segments(4.5, result[,"lwr"], 4.5, result[,"upr"], col=8, lwd=5)
textxy(4.5, result[,"lwr"], round(result[,"lwr"], 2), cex=0.85, col=8)
textxy(4.5, result[,"upr"], round(result[,"upr"], 2), cex=0.85, col=8)
# 95%
result<-predict(est, newdata=data.frame(x=c(4.5)), level=.95, interval="prediction")
segments(4.5, result[,"lwr"], 4.5, result[,"upr"], col=5, lwd=2)
textxy(4.5, result[,"lwr"], round(result[,"lwr"], 2), cex=0.85, col=5)
textxy(4.5, result[,"upr"], round(result[,"upr"], 2), cex=0.85, col=5)
points(4.5, result[,"fit"], pch=23, col="darkgreen", bg="green")