Co to znaczy, że model jest dobrze dopasowany?

1 Założenia modelu

Założenia modelu dotyczą typu zależności (czy jest to zależność liniowa) oraz tzw. reszt z modelu. Reszty z modelu to różnica między wartością obserwowaną \(y\) oraz wartością oszacowaną z modelu \(\hat{y}\).

\[reszty = y - \hat{y}\]

  1. Liniowość - związek między zmiennymi x oraz y ma charakter liniowy
  2. Normalność - reszty (błędy) mają rozkład normalny, ze średnią 0
  3. Homoscedastyczność - wariancja reszt jest taka sama dla wszystkich obserwacji
  4. Niezależność - błędy są niezależne

Do sprawdzenia założeń modelu wykorzystuje się wykresy diagnostyczne.

2 Budowa modelu

dane <- readRDS("dane/dane2007.rds")
model_lm <- lm(formula = tg10cm ~ tg5cm, data = dane)   
smr_model_lm <-summary(model_lm)
smr_model_lm
## 
## Call:
## lm(formula = tg10cm ~ tg5cm, data = dane)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.46323 -0.23488 -0.01063  0.26094  0.94909 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 0.229146   0.036519   6.275 9.98e-10 ***
## tg5cm       0.960496   0.002734 351.257  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4108 on 363 degrees of freedom
## Multiple R-squared:  0.9971, Adjusted R-squared:  0.9971 
## F-statistic: 1.234e+05 on 1 and 363 DF,  p-value: < 2.2e-16

3 Metody oceny dopasowania modelu

Istnieje kilka metod oceny dopasowania modelu:

  • współczynnik determinacji R2 - określa jaki procent zmienności zmiennej zależnej (objaśnianej) jest uwarunkowany zmiennością zmiennej niezależnej (objaśniającej); im wyższa wartość tym lepiej dopasowany model.
  • analiza reszt - reszty z modelu to różnica wartości obserwowanej i oszacowanej przez model. W dobrze dopasowanym modelu reszty będą miały rozkład normalny, ze średnią równą 0. Do analizy reszt, poza statystykami reszt, wykorzystuje się tzw. wykresy diagnostyczne.

3.1 Potencjalne błędy modelu:

  • model niedoszacowuje wartości

    • wartość obserwowana \(y\) jest wyższa od wartości wyliczonej z modelu (\(\hat{y}\), fitted)
    • przykład: wartość obserwowana temperatury: 10C, wartość wyliczona z modelu 7C. Reszta z modelu: 10C-7C=3C. Model przewidział o 3 stopnie za niską temperaturę.
    • średnia wartość reszt dodatnia

  • model przeszacowuje wartości

    • wartość obserwowana \(y\) jest niższa od wartości wyliczonej z modelu (\(\hat{y}\), fitted)
    • przykład: wartość obserwowana temperatury: 10C, wartość wyliczona z modelu 12C. Reszta z modelu: 10C-12C=-2C. Model przewidział o 2 stopnie za wysoką temperaturę.
    • średnia wartość reszt ujemna

“Dobry model” to taki, który nie będzie regularnie przeszacowywał lub niedoszacowywał przewidywanych wartości. W “dobrym modelu” reszty dodatnie i ujemne będą się równoważyć - średnia reszt będzie zbliżona do 0 (model czasami przewidzi niższą wartość, czasami wyższą, ale w większości przypadków różnice między tym co obserwowane a tym co wyliczone z modelu będą oscylowały koło 0).

3.2 Współczynnik determinacji (Multiple R-Squared) oraz zmododyfikowany współczynnik determinacji (Adjusted R-Squared)

  • Przy ocenie dopasowania modelu należy zwrócić uwagę na wartość R2 przedstawiającą procent wariancji wyjaśnionej przez model, tzn. jaki procent zmienności zmiennej zależnej (objaśnianej) jest uwarunkowany zmiennością zmiennej niezależnej (objaśniającej).

  • Im wyższa wartość współczynnika R2 (maksymalnie to 1) tym lepsze dopasowanie modelu do danych. Niestety również im więcej zmiennych w modelu tym wyższa jest wartość współczynnika R2.

  • Aby uwzględnić liczbę zmiennych w modelu stosuje się modyfikację współczynnika R2, tzw. Zmodyfikowany R2 (ang. Adjusted R2)

c( R2 = smr_model_lm$r.squared, R2_adj = smr_model_lm$adj.r.squared)
##        R2    R2_adj 
## 0.9970665 0.9970585

Interpretacja

99,71% zmienności zmiennej tg10cm jest wyjaśniona przez zmienność zmiennej tg5cm.

3.3 Określenie dowolnego przedziału ufności dla współczynników modelu

confint(model_lm, level = 0.99)
##                 0.5 %    99.5 %
## (Intercept) 0.1345839 0.3237090
## tg5cm       0.9534151 0.9675766

W powyższym przykładzie został obliczony 99% przedział ufności. Objaśnienie wyniku:

  • 0,5% - dolna granica przedziału ufności
  • 99,5% - górna granica przedziału ufności

Interpretacja

Obliczona wartość wyrazu wolnego (model_lm$coeff) równa jest 0.229146, a z 99% prawdopodobieństwem będzie się mieścić w przedziale 0.1345839 a 0.3237090

3.4 Identyfikacja wartości odstających - Bonferroni outlier test

Test identyfikuje najbardziej odstające obserwacje. Podaje numer przypadku (w danych 324 i 127)

library(car)
## Loading required package: carData
outlierTest(model_lm)
##      rstudent unadjusted p-value Bonferroni p
## 324 -6.324241         7.5041e-10   2.7390e-07
## 127 -5.153857         4.2053e-07   1.5349e-04

Możemy się tym obserwacjom bliżej przyjrzeć:

dane[c(324, 127), ]
##           data mce  tps tmin tmax tg5cm tg10cm tg20cm tg50cm tg100cm TwP TwMP
## 324 11/20/2007  11 -0.3 -1.3  0.8   4.2    1.8   2.07   3.50     6.5   3  3.0
## 127   5/7/2007   5 13.5  2.4 21.5  17.4   14.9  14.93  13.23    12.1  11 10.3
library(ggplot2)
ggplot(dane, aes(x=tg5cm, y=tg10cm)) + 
  geom_point() +
  geom_point(data = dane[c(324, 127),], aes(x=tg5cm, y=tg10cm), color = "red", size = 4) +
  coord_fixed(ratio = 1) + 
  xlab("tg5cm: Temperatura gruntu na 5cm [C]") + 
  ylab("tg10cm: Temperatura gruntu na 10cm [C]") + 
  geom_smooth() + 
  geom_abline(color = 'grey') + 
  theme_bw()
## `geom_smooth()` using method = 'loess' and formula 'y ~ x'

3.5 Podstawowe wykresy diagnostyczne

Po dopasowaniu modelu powinniśmy wpierw sprawdzić czy spełnione są przyjęte założenia. Jeśli zalożenia modelu są spełnione, to zakłócenia losowe (błędy) powinny mieć rozkład normalny o równych wariancjach. Założenia modelu możemy badać weryfikując właściwości reszt. Do badania reszt służą tzw. wykresy diagnostyczne.

  • Residuals vs Fitted (Reszty vs. wartości dopasowane przez model)

    • używany do sprawdzenia czy jest liniowa zależność między zmiennymi (Założenie 1. Liniowość).
    • Na osi poziomej (x) przedstawione są wartości dopasowane przez model \(\hat{y_i}\) (fitted values), a na osi pionowej (y) wartości residuów-reszt (residuals) \(\hat{\varepsilon_i} = y_i - \hat{y_i}\).
    • Linia czerwona o przebiegu poziomym dla y = 0 wskazuje na zależność liniową.

  • Normal Q-Q (wykres kwantylowy dla rozkładu normalnego)

    • używany do sprawdzenia czy reszty mają rozkład normalny (Założenie 2. Normalność)
    • punkty na wykresie powinne układać się po przekątnej przerywanej linii
    • odstępstwa od tej linii sugerują brak normalności oraz upoważniają do zastosowania transformacji nieliniowej badanych zmiennych
    • dla dużych prób nie należy się przejmować niewielkimi odstępstwem od normalności residuów

  • Scale-Location

    • używany do sprawdzenia, czy wariancja reszt jest taka sama dla wszystkich obserwacji (Założenie 3. Homoscedastyczność)
    • na osi poziomej (x) przedstawione są wartości dopasowane przez model \(\hat{y_i}\) (fitted values), a na osi pionowej (y) pierwiastek ze standaryzowanych residuów-reszt (residuals) \(\hat{\varepsilon_i} = y_i - \hat{y_i}\). Standaryzowane reszty to wartości reszt podzielone przez ich błąd standardowy (estimated standard error).
    • czerwona pozioma linia wraz z równomiernie rozłożonymi punktami na całej długości x jest dobrym wskaźnikiem, że założenie jest spełnione.

  • Residuals vs Leverage

    • używany do identyfikacji obserwacji nietypowych - tj. odstających wartości, które mogą mieć wpływ na przebieg linii regresji. Nie każda wartość odstająca jest nietypowa (wpływowa); niektóre wartości odstające nie zmienią przebiegu linii regresji.
    • Punkty znajdujące się w prawym lub lewym górnym rogu wykresu mogą mieć wpływ na przebieg linii regresji. Trzeba się takim punktom bliżej przyjrzeć oraz ewentualnie wyeliminować je z obliczeń.
    • Na osi pionowej przedstawione są standaryzowane reszty a na osi poziomej tzw. siły dźwigni \(h_i\) (miary wpływu tej obserwacji na ocenę współczynników modelu regresji, ang. levarege).
    • Standaryzowane reszty to wartości reszt podzielone przez ich błąd standardowy (estimated standard error). Standaryzowane reszty mogą być interpretowane jako liczba błędów standardowych od linii regresji. Punkt, który ma wartośc standaryzowanych reszt > 3 potencjalnie może być wartością odstającą (trzeba się takiej obserwacji bliżej przyjrzeć)
    • Siła dźwigni określa jaki wpływ na ocenę współczynników modelu regresji ma obserwacja \(X_i\), czyli jak bardzo różnią się oceny współczynników dla modelu z tą obserwacją i dla modelu bez tej obserwacji.
    • Obserwacje nietypowe identyfikuje się używając tzw. odległości Cooka (Cook distance). Przyjmuje się, że obserwacja będzie miała duży wpływ na model jeśli przekroczy wartość 4/(n - p - 1), gdzie n to liczba obserwacji w zbiorze danych, a p to liczba predyktorów (zmiennych objaśniających).

Czego możemy dowiedzieć się z wykresów diagnostycznych dla modelu zależności między tg10cm a tg5 cm?

3.5.1 Domyślne wykresy diagnostyczne

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

3.5.2 Pakiet ggfortify

library(ggfortify)
autoplot(model_lm) + theme_classic()

3.5.3 Pakiet ggResidpanel

Pakiet ggResidpanel dostarcza 5 funkcji, które pozwalają na wykonanie wykresów diagnostycznych zgodnych ze stylizacją pakietu ggplot() lub też interaktywnych wykresów diagnostycznych (wykorzystując bibliotekę plotly). Pakiet pozwala na wykonanie 9 różnych wykresów (argument plots definiuje jakie wykresy będą wykonane). Opis możliwości pakietu znajduje się w dokumentacji, na stronie .

Pakiet składa się z 5 funkcji:

  • resid_panel: tworzy panel z wykresami diagnostycznymi na podstawie obiektu modelu.
  • resid_interact: tworzy interaktywny panel z wykresami diagnostycznymi.
  • resid_xpanel: Creates a panel of diagnostic plots of the predictor variables
  • resid_compare: tworzy panel z wykresami diagnostycznymi dla kilku modeli.
  • resid_auxpanel: Creates a panel of diagnostic plots based on the predicted values and residuals from models not included in the package
library(ggResidpanel)

Wykresy diagnostyczne - przykłady

resid_panel(model_lm)

resid_panel(model_lm, plots = "R")

resid_panel(model_lm, plots = "all")

resid_panel(model_lm, plots = c("resid", "hist", "yvp"))

Interaktywne wykresy diagnostyczne

resid_interact(model_lm)
## Warning in geom2trace.default(dots[[1L]][[1L]], dots[[2L]][[1L]], dots[[3L]][[1L]]): geom_GeomFunction() has yet to be implemented in plotly.
##   If you'd like to see this geom implemented,
##   Please open an issue with your example code at
##   https://github.com/ropensci/plotly/issues
resid_interact(model_lm, plots = "R")

4 Zadania do samodzielnego rozwiązania

Używając danych z pliku dane2007.rds proszę wykonać zadania 1-2. Wyniki proszę przedstawić w postaci krótkiego raportu (wraz z rycinami do 2 stron) zawierającego odpowiednie wykresy oraz interpretację wyników.

Zadanie 1.

Proszę zbudować model zależności temperatury wody w Parsęcie (TwP) od temperatury powietrza (tps). Proszę wykonać ocenę dopasowania modelu obejmującą:

  • określenie współczynnika determinacji.

  • identyfikację wartości odstających

  • analizę reszt:

    • statystyki reszt
    • wykresy diagnostyczne

Zadanie 2.

Proszę zbudować model zależności temperatury wody w Parsęcie (TwP) od temperatury gruntu na 50 cm. Proszę wykonać ocenę dopasowania modelu obejmującą:

  • określenie współczynnika determinacji.

  • identyfikację wartości odstających

  • analizę reszt:

    • statystyki reszt
    • wykresy diagnostyczne