Regresja prosta

Regresja jest to podstawowy model opisujący współzależność zjawisk. Zajmiemy sie w tym rozdziale regresją prostą czyli liniowym opisem zależności dwóch zmiennych. Przeanalizujemy znów przykład na podstawie danych mtcars. Funkcja, która dopasowuje model regresji liniowej do danych w R to lm. Zachęcam do dokładnego przeczytrania opisu działania tej funkcji ?lm.

Zbadajmy liniową zależność między spalaniem, objetością silnika.

data(mtcars)

#dopasowujemy model liniowy zależności spalania (zmienna objaśniana) od objętości silnika
fit.mpg.disp<-lm(mpg~disp, data=mtcars)

#podsumowanie wyników estymacji
summary(fit.mpg.disp)
## 
## Call:
## lm(formula = mpg ~ disp, data = mtcars)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.8922 -2.2022 -0.9631  1.6272  7.2305 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 29.599855   1.229720  24.070  < 2e-16 ***
## disp        -0.041215   0.004712  -8.747 9.38e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.251 on 30 degrees of freedom
## Multiple R-squared:  0.7183, Adjusted R-squared:  0.709 
## F-statistic: 76.51 on 1 and 30 DF,  p-value: 9.38e-10

W podsumowaniu mamy większość interesujących wartości. Najważniejsze z praktycznego punktu widzenia są parametry modelu. W naszym modelu mamy dwa parametry. Parametr \(\beta_0\) (wyraz wolny, intercept) interpretujemy jako wartość oczekiwaną zmiennej objaśnianej w przypadku, gry zmienna objaśniająca jest równa 0. W naszym przypadku ta interpretacja jest raczej mało sensowna ponieważ mówi o tym, że gdyby samochód miał silnik o pojemności 0 cu.in. (cale sześcienne) to przejechałby 29.6 mil na galonie. Parametr \(\beta_1\) (disp) mówi nam ile średnio wzrośnie wartość zmiennej objaśnianej przy wzroście zmiennej objaśniającej o jedną jednostkę. W naszym przypadku, gdy samochód ma pojemność większo o 1 cu.in. to przejeżdża na galonie o średnio 0.04 mili.

Nie mniej ważne od wartości parametrów są istotności parametrów. W naszym przypadku, oba parametry są istotne (ostatnia kolumna tabeli z parametrami, Pr(>|t|))

Problemy w regresji prostej.

Jest wiele rzeczy na które należy zwrócić uwagę przy budowie modelu. Regresja prosta jest dobrym modelem, aby zrozumieć i skwantyfikować zależność między dwoma zmiennymi. Należy jednak pamiętać, że zalezność między zmiennymi musi być liniowa (w przybliżeniu), aby interpretacja wyników była prawdziwa. Innym problemem są zależności między zmiennymi nie uwzględnionymi w modelu. Bardziej technicznym problemem jest homoskedastyczność, która mówi o tym, iż wariancja jest równa dla każdego poziomu zmiennej objaśniającej. Dla naszego przykładu znaczyłoby to iż odchylenie od średniej spalania dla samochodu z silnikiem o pojemności 160 cu.in. jest takie samo jak dla samochodu z silnikiem 360 cu.in. To założenie w wielu przypadkach nie jest spełnione.

Pierwszym badaniem, czy model regresji liniowej jest odpowiedni jest analiza wykresów. Po pierwsze stwórzmy wykres zależności mpgod disp oraz narysujmy prostą wyestymowanej regresji.

#dla wygody przypiszemy odpowiednie zmienne do x i y
x<-mtcars$disp
y<-mtcars$mpg
#model regresji dla x i y
fit.lm<-lm(y~x)

#współczynniki
beta0<-coef(fit.lm)[1]
beta1<-coef(fit.lm)[2]

#prosta regresji
curve(beta0+beta1*x, from = min(x), to = max(x), col='red', ylim=c(min(y), max(y)),
      xlab='pojemność silnika (cu.in.)', ylab='spalanie (mpg)')
points(x,y)

Aby łatwiej zdiagnozować ewentualne problemy z modelem regresji należy stworzyć wykres reszt z modelu. Reszty z modelu to róznica między zobserwowaną wartościa a wartościa średnią wynikająca z modelu. W przypadku modelu liniowego to odległość w lini pionowej punktu od prostej regresji (punkty pod prostą mają ujemne reszty). Aby stworzyć wktor reszt użyjemy funkcji resid.

#tworzymy wektor reszt z modelu
res.lm<-resid(fit.lm)

#wykres
plot(x, res.lm,  xlab='pojemność silnika (cu.in.)', ylab='reszta', pch=19, col='red')
abline(h=0)

Test modelu

W przypadku modelu regresji prostej, istotność parametrów mówi nam o tym, czy model jest odpowiedni. Dla większej ilości zmiennych objasniających moglibyśmy się zastanawiać czy dodatkowe zmienne wprowadzają dodatkowe informacje do modelu.

W najprostszym przypadku możemy przetestować czy zmienna objaśniana jest rzeczywiście zalezna od zmiennej objaśniającej. Drugim testem jest test czy może wyraz wolny (\(\beta_0\)) jest istotnie różny id zera. Do testowania istotności modelu w stosunku do innego modelu (muszą to być modele zagnieżdżone) służy funkcja anova.

Zobaczmy jak działa test anova w dwóch najprostszych przypadkach. Najpierw regresja liniowa względem modelu bez zmiennej objaśniającej

#dopasowujemy model bez zmiennej x
fit.0<-lm(y~1)

#przeprowadzamy trst anova dla fit.lm i fit.0
anova(fit.lm, fit.0)
## Analysis of Variance Table
## 
## Model 1: y ~ x
## Model 2: y ~ 1
##   Res.Df     RSS Df Sum of Sq      F   Pr(>F)    
## 1     30  317.16                                 
## 2     31 1126.05 -1   -808.89 76.513 9.38e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Zauważmy, że p-value tego modelu jest dokładnie takie jak p-value dla o istotności parametru disp w modelu fit.lm.

Drugim przypadkiem jest regresja liniowa względem regresji liniowej o wyrazie wolnym równy 0.

#dopasowujemy model z beta0=0
fit.1<-lm(y~x-1)

#przeprowadzamy trst anova dla fit.lm i fit.1
anova(fit.lm, fit.1)
## Analysis of Variance Table
## 
## Model 1: y ~ x
## Model 2: y ~ x - 1
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     30  317.2                                  
## 2     31 6442.4 -1   -6125.2 579.38 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1