Główną ideą regresji jest wyjaśnianie, przewidywanie, prognozowanie danych dla pewnej zmiennej (określanej jako zmienna zależna, objaśniana, wyjaśniana) na podstawie innych zmiennych (zmienne niezależne, predyktory, zmienne objaśniające, zmienne wyjaśniające).

Budowa liniowego modelu zależności między dwoma zmiennymi składa się z trzech etapów:

  1. Określenie czy między zmiennymi istnieje zależność.

    • wykonanie wykresu rozrzutu
    • współczynnik korelacji

  2. Budowa modelu.

  3. Analiza wynków modelu (w tym ocena dopasowania modelu).

Po zbudowaniu modelu oraz ocenie dopasowania modelu można go wykorzystać do wykonania predykcji.

dane <- readRDS("dane/dane2007.rds")

1 Określenie zależności między temperaturą gruntu na 5 cm i 10 cm.

Wykres rozrzutu

library(ggplot2)
ggplot(dane, aes(x=tg5cm, y=tg10cm)) + geom_point() +
  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'

Współczynnik korelacji liniowej

cor(dane$tg5cm, dane$tg10cm, use= "complete.obs", method = "pearson")
## [1] 0.9985322

2 Budowa modelu

W R do budowy modelu liniowego służy funkcja lm(). Obowiązkowo należy podać 2 argumenty:

  • formula - opisuje równanie modelu liniowego
  • data - wskazuje zbiór danych, zawierający dane wejściowe - zmienne objaśniające.
model_lm <- lm(formula = tg10cm ~ tg5cm, data = dane)

3 Wyniki modelu

Funkcja lm() dopasowuje model liniowy, wyznacza oceny współczynników \(\beta\) oraz wylicza wartości reszt. Wynikiem działania funkcji jest obiekt klasy lm, który będzie przechowywał informacje o dopasowanym modelu. Do ważniejszych informacji należą:

  • coefficients - wartości dopasowanych współczynników modelu
  • residuals - wektor reszt (różnica dla wartości obserwowanej i oszacowanej przez model)
  • fitted.values - wektor ocen modelu
  • call - formuła użyta do budowy modelu
  • model - ramka danych użyta do budowy modelu
str(model_lm)

3.1 Współczynniki modelu (coefficients)

model_lm$coefficients
## (Intercept)       tg5cm 
##   0.2291464   0.9604959
#to samo uzyskamy podając model_lm$coeff

W przypadku regresji linowej funkcja regresji przyjmuje postać funkcji liniowej:

\[y = ax + b\] gdzie,

  • \(y\) to zmienna, której wartość chcemy przewidywać (zmienna zależna, zmienna wyjaśniana, przewidywana)
  • \(x\) to zmienna objaśniająca (określane też jako zmienna niezależne, predyktory, zmienne wyjaśniające)
  • \(a\) to współczynnik kierunkowy (mówiący o nachyleniu prostej regresji)
  • \(b\) to wyraz wolny

W wyniku otrzymamy wartości dwóch współczynników regresji:

  • współczynnika mówiącego o nachyleniu prostej regresji - określa o ile jednostek przeciętnie wzrośnie (lub zmaleje gdy współczynnik ten ma wartość ujemną) wartość zmiennej zależnej (objaśnianej, przewidywanej), gdy wartość zmiennej niezależnej (objaśniającej, predyktorów) wzrośnie o jedną jednostkę.

    • W przykładzie - na każdy 1C wzrostu temperatury gruntu na 5 cm (tg5cm) przeciętnie przypada wzrost temperatury gruntu na 10 cm (tg10cm) o 0.96C.

  • wyraz wolny - Intercept

W efekcie otrzymamy równanie modelu regresji dla zależności między temperaturą gruntu na 5 i 10 cm.

\[tg10cm = 0.960496 *tg5cm + 0.229146\] Interpretacja: W przykładzie - na każdy 1C wzrostu temperatury gruntu na 5 cm (tg5cm) przeciętnie przypada wzrost temperatury gruntu na 10 cm o 0.96C.

3.2 Podsumowanie wyników modelu

summary(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

Podsumowanie wyników modelu regresji składa się z kilku elementów:

  • formuły użytej do budowy modelu (Call)
  • statystyk reszt (Residuals) - reszty wyznacza się jako różnicę wartości obserwowanej oraz wartości oszacowanej z modelu.
  • współczynników modelu (Coefficients)
  • wyników mających na celu ocenę jakości uzyskanego modelu (Residual standard error, Multiple R-square, Adjusted R-squared)

Formuła użyta do budowania modelu

model_lm$call
Call:
lm(formula = tg10cm ~ tg5cm, data = dane)

Reszty z modelu

Do oceny poprawności dopasowania modelu wykorzystuje się reszty (tzw.residua). Reszty wyznacza się jako różnicę wartości obserwowanej Y a wartości oszacowanej z modelu Ŷ.

Statystyki reszt – Min – minimum, 1Q – pierwszy kwartyl, Median – mediana, 3Q – trzeci kwartyl, Max – maksimum.

## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.46323 -0.23488 -0.01063  0.26094  0.94909

Interpretacja wyników:

  • 50% odchyleń mieści się w zakresie -0.23 do 0.26
  • Mediana ma wartość -0.01063. Wartości tg10cm mieści się w zakresie 0 do 27.33 a zatem -0.01 to niewielkie odchylenie.
  • Duże odchylenia mediany od 0 mogą świadczyć o zależności nieliniowej.
  • Średnia wartość błędu powinna być równa 0 - oznacza to, że błędy dodatnie i ujemne się równoważą
mean(model_lm$residuals)
## [1] 2.839496e-17

Współczynniki

smr = summary(model_lm)
smr$coefficients
## 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

Tabela współczynniki zawiera następujące informacje dla każdego współczynnika:

  • Estimate – oceny wartości współczynników równania regresji
  • Std.Error – błąd standardowy współczynnika (wartość 1 odchylenia standardowego; dla rozkładu normalnego w przedziale średnia±1 odchylenie standardowe mieści się 66% wartości danych)
  • t-value – wartość statystyki t
  • p-value – poziom istotności dla współczynnika

Interpretacja:

  • Błąd standardowy współczynnika wskazuje, że z 66% prawdopodobieństwem możemy stwierdzić, że wartość współczynnika wyrazu wolnego wynosi 0.229146 ± 0.036519 (tj. 1 błąd standardowy ) a współczynnika kierunkowego 0.960496±0.002734.
  • Oba współczynniki są bardzo wysoce istotne statystycznie. Innymi słowy prawdopodobieństwo, że wartość współczynnika jest przypadkowa wynosi 9.98e-10 dla współczynnika kierunkowego oraz < 2.e-16 dla współczynnika wyrazu wolnego, a zatem jest znikomo mała.

Residual standard error (błąd standardowy reszt; standardowy błąd oceny)

  • Wartość 1 odchylenia standardowego; dla rozkładu normalnego w przedziale średnia±1 odchylenie standardowe mieści się 66% wartości danych.
  • Pokazuje o ile przeciętnie wartości empiryczne odchylają się od wartości teoretycznych (wynikających ze zbudowanego modelu)
Residual standard error: 0.4108 on 363 degrees of freedom

Z 66% prawdopodobieństwem możemy stwierdzić, że wartość rzeczywista nie będzie się różnić od prognozowanej o więcej niż ±0.4108 C.

Multiple R-Squared, 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)

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

Multiple R-squared:  0.9971,    Adjusted R-squared:  0.9971

p-value

– wartość poziomu istotności całego modelu Model jest wysoko istotny

## F-statistic: 1.234e+05 on 1 and 363 DF,  p-value: < 2.2e-16

3.3 Informacje zawarte w obiekcie klasy lm

Stworzenie ramki danych zawierającej wyniki modelu

#ramka danych zawierająca dane użyte do budowy modelu
model_df <- model_lm$model
#reszty - różnica między warrtością obserwowaną Y a wartością oszacowaną z modelu Ŷ. 
model_df$reszty <- model_lm$residuals
#wartości oszacowane przez model
model_df$fitted <-  model_lm$fitted.values
head(model_df)
##   tg10cm tg5cm     reszty   fitted
## 1    5.0   5.1 -0.1276753 5.127675
## 2    4.5   4.1  0.3328206 4.167179
## 3    4.5   4.2  0.2367710 4.263229
## 4    4.8   4.6  0.1525727 4.647427
## 5    5.5   5.2  0.2762752 5.223725
## 6    5.6   5.2  0.3762752 5.223725

4 Predykcja

Mając dopasowany model liniowy możemy zastosować go do predykcji wartości \(y\) dla zadanych wartości \(X\). W poniższym przykładzie zastosujemy model liniowy do predykcji wartości tg10cm dla roku 2009.

Obiekt model_lm zawiera liniowy model zależności pomiędzy temperaturą gruntu na 5cm i temperaturą gruntu na 10 cm, zbudowany na podstawie danych z 2007 roku. Równanie regresji ma postać:

\[tg10cm = 0.960496 *tg5cm - 0.229146\] Znając zatem wartości temperatury gruntu na 5 cm w roku 2009, możemy oszacować wartości temperatury gruntu na 10 cm w 2009 roku.

Do wykonania predykcji w R służy funkcja predict(). W funkcji tej musimy podać 2 argumenty: 1) obiekt modelu; 2) zbiór danych zawierający zmienne zależne (objaśniające). Nazwy kolumn muszą być zgodne z nazwami zmiennych zależnych użytych do budowy modelu.

dane2009 <- read.csv("dane/tg5cm_2009.csv")
head(dane2009)
##   ID tg5cm
## 1  1  -2.0
## 2  2  -3.9
## 3  3  -2.2
## 4  4  -1.1
## 5  5  -3.2
## 6  6  -3.8

Kolumna tg5cm zawiera wartości temperatury gruntu na 5 cm. Wartości te zostaną użyte do oszacowania wartości temperatury gruntu na 10cm.

p <- predict(model_lm, dane2009)
head(p)
##         1         2         3         4         5         6 
## -1.691845 -3.516787 -1.883944 -0.827399 -2.844440 -3.420738

Predykcja może także zostać wykonana dla wybranych wartości.

#new to data.frame zawierajaca zmienne niezależne (objasniajace)
new <- data.frame(tg5cm = c(-3, 4, 10))
p <- predict(model_lm, new)
p
##         1         2         3 
## -2.652341  4.071130  9.834105