Regresja wielokrotna
Budowa modelu liniowego prognozującego ceny mieszkań we Wrocławiu.
Anna Dmowska (dmowska@amu.edu.pl)
1 Dane
W ćwiczeniu zostaną wykorzystane dane dołączone do podręcznika “Przewodnik po pakiecie R” Przemysława Biecka, wydanie IV, 2017. Dane dotyczą 200 mieszkań w 3 dzielnicach Wrocławia (cena, liczba pokoi, powierzchnia, typ zabudowy, dzielnica).
- Instalacja pakietu zawierającego dane używane w książce “Przewodnik po pakiecie R”:
#install.packages("Przewodnik")- Wczytanie biblioteki zawierajacej dane:
library(Przewodnik)- Wczytanie danych mieszkania:
dane <- mieszkania - Podstawowe informacje o danych:
summary(dane)## cena pokoi powierzchnia dzielnica
## Min. : 83280 Min. :1.00 Min. :17.00 Biskupin :65
## 1st Qu.:143304 1st Qu.:2.00 1st Qu.:31.15 Krzyki :79
## Median :174935 Median :3.00 Median :43.70 Srodmiescie:56
## Mean :175934 Mean :2.55 Mean :46.20
## 3rd Qu.:208741 3rd Qu.:3.00 3rd Qu.:61.40
## Max. :295762 Max. :4.00 Max. :87.70
## typ.budynku
## kamienica :61
## niski blok:63
## wiezowiec :76
##
##
## 2 Model liniowy z jednym predyktorem
Najprostszy model regresji liniowej będzie miał postać:
zmObjaśniana ~ zmObjaśniajaca
Przykładem takiego modelu będzie zależność ceny mieszkań od powierzchni lub ceny mieszkań od liczby pokoi.
Przykład 1: Zależność ceny mieszkania od powierzchni mieszkania
model_1 <- lm(formula = cena ~ powierzchnia, data = dane)
summary_m1 <- summary(model_1)summary_m1##
## Call:
## lm(formula = cena ~ powierzchnia, data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -39522 -9337 -1207 9443 35652
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 82238.27 2503.36 32.85 <2e-16 ***
## powierzchnia 2028.18 49.72 40.79 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14080 on 198 degrees of freedom
## Multiple R-squared: 0.8937, Adjusted R-squared: 0.8931
## F-statistic: 1664 on 1 and 198 DF, p-value: < 2.2e-16Przykład 2: Zależność ceny mieszkania od liczby pokoi
model_2 <- lm(formula = cena ~ pokoi, data = dane)
summary_m2 <- summary(model_2)summary_m2##
## Call:
## lm(formula = cena ~ pokoi, data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -56655 -11952 -167 13649 68582
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 85812 3589 23.91 <2e-16 ***
## pokoi 35342 1296 27.27 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 19800 on 198 degrees of freedom
## Multiple R-squared: 0.7897, Adjusted R-squared: 0.7887
## F-statistic: 743.7 on 1 and 198 DF, p-value: < 2.2e-163 Uwzględnienie kilku zmiennych objaśniających w modelu
W modelu regresji można uwzględnić kilka zmiennych objaśniających. Wtedy model liniowy będzie mieć postać:
zmObjasniana ~ zmObjasniajaca1 + zmObjasniajaca2 + … + zmObjasniajacaN
Przykład 3: Zależność ceny mieszkania od powierzchni mieszkania oraz liczby pokoi
W poniższym przykładzie stworzymy model zależności ceny mieszkań od liczby pokoi oraz powierzchni mieszkania:
model_3 <- lm(formula = cena ~ powierzchnia + pokoi, data = dane)
summary_m3 <- summary(model_3)summary_m3##
## Call:
## lm(formula = cena ~ powierzchnia + pokoi, data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -39705 -9386 -863 9454 35097
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 82407.1 2569.9 32.066 <2e-16 ***
## powierzchnia 2070.9 149.2 13.883 <2e-16 ***
## pokoi -840.1 2765.1 -0.304 0.762
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 14110 on 197 degrees of freedom
## Multiple R-squared: 0.8937, Adjusted R-squared: 0.8926
## F-statistic: 828.3 on 2 and 197 DF, p-value: < 2.2e-164 Uwzględnienie zmiennych jakościowych w modelu
- W modelu regresji liniowej zmienna objaśniana musi być zmienną ciągłą.
- Zmienne objaśniające (predyktory) mogą być zmiennymi jakościowymi.
- Gdy chcemy w modelu umieścić zmienne jakościowe, musimy zamienić poziomy zmiennej na wartości liczbowe.
- Domyślnie przyjmuje się poziom 1 za poziom referencyjny.
Odczytanie poziomów zmiennej dzielnica:
#str(dane)
levels(dane$dzielnica)## [1] "Biskupin" "Krzyki" "Srodmiescie"Dzielnica Biskupin będzie traktowana jako poziom referencyjny.
Przykład 4: Budowa modelu zależności ceny mieszkania od powierzchni, liczby pokoi oraz dzielnicy w której mieszkanie się znajduje
model_4 <- lm(cena ~ powierzchnia + pokoi + dzielnica, dane)
summary_m4 <- summary(model_4)summary_m4##
## Call:
## lm(formula = cena ~ powierzchnia + pokoi + dzielnica, data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30501 -8480 -144 7346 35730
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 94222.02 2320.36 40.607 < 2e-16 ***
## powierzchnia 2022.99 116.31 17.393 < 2e-16 ***
## pokoi 34.36 2157.52 0.016 0.987
## dzielnicaKrzyki -20934.86 1842.79 -11.360 < 2e-16 ***
## dzielnicaSrodmiescie -12722.60 2008.03 -6.336 1.6e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 11000 on 195 degrees of freedom
## Multiple R-squared: 0.9362, Adjusted R-squared: 0.9348
## F-statistic: 714.8 on 4 and 195 DF, p-value: < 2.2e-16Interpretacja wyników modelu:
W kolumnie Estimate wartość oceny współczynnika dla dzielnicy Krzyki wynosi -20934. Oznacza to, że o tyle średnio tańsze są mieszkania w dzielnicy Krzyki od podobnego mieszkania w dzielnicy Biskupin.
zmienna objaśniająca pokoi jest nie istotna statystycznie, ponieważ zmienna ta jest skorelowana ze zmienną powierzchnia
5 Rozbudowywanie modelu
Celem analizy regresji jest zbudowanie modelu w jak najlepszym stopniu wykonującego predykcje.
Dąży się do budowy modeli, które będą dawały jak najlepszą predykcję a równocześnie nie będą zbyt skomplikowane (uzależnione od dużej ilości zmiennych).
Przykład 5: Zależność ceny mieszkania od powierzchni mieszkania, liczby pokoi, dzielnicy oraz typu budynku
Poniższy model uwzględnia wszystkie zmienne zawarte w zbiorze danych
model_5 <- lm(cena ~ powierzchnia + pokoi + dzielnica + typ.budynku, dane)
summary_m5 <- summary(model_5)summary_m5##
## Call:
## lm(formula = cena ~ powierzchnia + pokoi + dzielnica + typ.budynku,
## data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -31406.2 -6056.1 378.6 5709.4 29052.7
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 92267.0 2295.4 40.197 < 2e-16 ***
## powierzchnia 1977.8 104.0 19.023 < 2e-16 ***
## pokoi 401.1 1917.5 0.209 0.835
## dzielnicaKrzyki -21473.4 1642.5 -13.073 < 2e-16 ***
## dzielnicaSrodmiescie -14302.6 1794.2 -7.972 1.34e-13 ***
## typ.budynkuniski blok 11453.1 1766.2 6.485 7.27e-10 ***
## typ.budynkuwiezowiec 407.2 1706.4 0.239 0.812
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9728 on 193 degrees of freedom
## Multiple R-squared: 0.9505, Adjusted R-squared: 0.949
## F-statistic: 618.3 on 6 and 193 DF, p-value: < 2.2e-165.1 Zestawienie wyników dotychczas stworzonych modeli:
## model r2 std_error
## 1 cena ~ powierzchnia 0.8931359 14082.444
## 2 cena ~ pokoi 0.7886867 19802.752
## 3 cena ~ powierzchnia + pokoi 0.8926437 14114.835
## 4 cena ~ powierzchnia + pokoi + dzielnica 0.9348432 10996.183
## 5 cena ~ powierzchnia + pokoi + dzielnica + typ.budynku 0.9490082 9727.745
## min_residuals max_residuals mean_residuals median_residuals
## 1 -39521.97 35652.26 -3.582157e-13 -1207.3956
## 2 -56654.62 68582.04 -2.605294e-13 -166.6661
## 3 -39705.05 35097.49 6.132872e-13 -863.4924
## 4 -30501.42 35729.61 3.020695e-13 -144.1032
## 5 -31406.16 29052.69 -1.820766e-14 378.55955.2 Do jakiego momentu opłaca się komplikowanie modelu?
Pytanie: Do jakiego momentu “opłaca się” komplikować model - tj. dodawać kolejne zmienne wyjaśniajace?
Kolejne komplikowanie postaci modelu skutkuje coraz mniejszym przyrostem jakości dopasowania, dlatego ważne jest ustalenie kryterium które umożliwiłoby w obiektywny i powtarzalny sposób zachowanie proporcji między dążeniem do prostoty modelu a jakością modelu.
Taką miarą jest kryterium informacyjne Akaike (AIC). Wartość AIC spada tak długo jak siła predykcyjna modelu wzrasta. Preferowany jest model o najmniejszej wartości AIC.
res_AIC = c(M1 = AIC(model_1),
M2 = AIC(model_2),
M3 = AIC(model_3),
M4 = AIC(model_4),
M5 = AIC(model_5))
res_AIC## M1 M2 M3 M4 M5
## 4392.639 4528.996 4394.545 4296.633 4249.545Który model ma najmniejszą wartość AIC?
which.min(res_AIC)## M5
## 5zestawienie$AIC <- res_AICzestawienie## model r2 std_error
## 1 cena ~ powierzchnia 0.8931359 14082.444
## 2 cena ~ pokoi 0.7886867 19802.752
## 3 cena ~ powierzchnia + pokoi 0.8926437 14114.835
## 4 cena ~ powierzchnia + pokoi + dzielnica 0.9348432 10996.183
## 5 cena ~ powierzchnia + pokoi + dzielnica + typ.budynku 0.9490082 9727.745
## min_residuals max_residuals mean_residuals median_residuals AIC
## 1 -39521.97 35652.26 -3.582157e-13 -1207.3956 4392.639
## 2 -56654.62 68582.04 -2.605294e-13 -166.6661 4528.996
## 3 -39705.05 35097.49 6.132872e-13 -863.4924 4394.545
## 4 -30501.42 35729.61 3.020695e-13 -144.1032 4296.633
## 5 -31406.16 29052.69 -1.820766e-14 378.5595 4249.5456 Regresja z krokową eliminacją zmiennych: automatyczny dobór zmiennych do modelu
Dodawanie lub usuwanie zmiennych z modelu można wykonać ręcznie (przykład powyżej) lub też automatycznie.
Dodając ręcznie kolejne zmienne należy z modelu usuwać te, które nie są istotne statystycznie (zmienna pokoi w powyższym przykładzie).
Regresja krokowa jest odmianą analizy regresji (nie tylko regresji liniowej), w której do modelu wprowadzane są jedynie istotne statystycznie zmienne, predyktory, które rzeczywiście „poprawiają” zbudowany model.
Metoda regresji krokowej stanowi bardzo użyteczne narzędzie w eksploracji dużej ilości zebranych danych.
Procedura regresji krokowej opiera się na krokowym (jedna zmienna – jeden krok) wprowadzaniu bądź usuwaniu zmiennych z modelu. Metody krokowe można podzielić zatem na dwa rodzaje:
- metody postępujące (forward selection), w których na początku w modelu nie ma żadnego predyktora (rozpoczynamy od modelu tylko z wyrazem wolnym) i co krok do modelu wprowadzane są kolejne, najbardziej istotne statystycznie predyktory. Procedura powtarzana jest tak długo, aż nie będzie już żadnych istotnych zmiennych do dodania.
- metody wsteczne (backward selection), w których na początku wprowadzane są wszystkie analizowane predyktory (budujemy model zależny od wszystykich zmiennych) i z każdym krokiem kolejne nieistotne predyktory są usuwane z modelu.
Niektóre algorytmy (np. funkcja step()) pozwalają także na wykonanie regresji krokowej w obu kierunkach. Procedura taka składa się z następujących kroków:
- rozpoczynamy od modelu z jednym składnikiem (wyrazem wolnym).
- dodajemy do modelu zmienną z najmniejszą p-wartością o ile jest ona istotna (istotność oznacza poprawę kryterium AIC),
- usuwamy z modelu zmienną z największą p-wartością o ile jest ona nieistotna,
- krok 2 i 3 powtarzamy tak długo aż model przestanie się zmieniać.
Do automatycznego wyboru istotnych zmiennych objaśniających wykorzystuje się funkcję step(). Funkcja step() uwzględnia 3 “kierunki” działania:
- regresję postępującą (direction = Forward selection)
- regresja wsteczna (direction = Backward selection)
- regresja w obu kierunkach (direction = both)
6.1 Regresja postępująca (direction = Forward selection)
Pzy budowaniu modelu za pomocą regresji krokowej wykorzystywane jest kryterium informacyjne Akaika. Preferowany jest model o najmniejszej wartości AIC.
start_m <- lm(cena ~ 1, data = dane) #model tylko z wyrazem wolnym
end_m <- lm(cena ~ ., data = dane) #model ze wszystkimi zmiennymi w zbiorze danych
model_koniec_1 <- step(start_m, scope=list(lower=start_m, upper=end_m), direction="forward")
#model_koniec_1 <- step(lm(cena ~ ., data = dane), direction = "forward")
model_koniec_1
Wyniki regresji z krokową eliminacją zmiennych - regresja postępująca
6.2 Regresja wsteczna (direction = Backward selection)
W regresji krokowej wstecznej buduje się możliwie najbardziej rozbudowany model i sprawdza jak usunięcie pojedyńczej zmiennej wpłynie na obniżenie AIC.
model_koniec_2 <- step(lm(cena ~ ., data = dane), direction = "backward")## Start: AIC=3679.97
## cena ~ pokoi + powierzchnia + dzielnica + typ.budynku
##
## Df Sum of Sq RSS AIC
## - pokoi 1 4.1403e+06 1.8268e+10 3678.0
## <none> 1.8263e+10 3680.0
## - typ.budynku 2 5.3152e+09 2.3579e+10 3727.1
## - dzielnica 2 1.6401e+10 3.4664e+10 3804.1
## - powierzchnia 1 3.4242e+10 5.2506e+10 3889.2
##
## Step: AIC=3678.01
## cena ~ powierzchnia + dzielnica + typ.budynku
##
## Df Sum of Sq RSS AIC
## <none> 1.8268e+10 3678.0
## - typ.budynku 2 5.3111e+09 2.3579e+10 3725.1
## - dzielnica 2 1.6398e+10 3.4665e+10 3802.1
## - powierzchnia 1 3.1032e+11 3.2858e+11 4253.9model_koniec_2##
## Call:
## lm(formula = cena ~ powierzchnia + dzielnica + typ.budynku, data = dane)
##
## Coefficients:
## (Intercept) powierzchnia dzielnicaKrzyki
## 92326.3 1998.3 -21463.6
## dzielnicaSrodmiescie typ.budynkuniski blok typ.budynkuwiezowiec
## -14310.0 11462.1 439.86.3 Regresja w obu kierunkach (direction = both)
model_koniec_3 <- step(lm(cena ~ ., data = dane), direction = "both")## Start: AIC=3679.97
## cena ~ pokoi + powierzchnia + dzielnica + typ.budynku
##
## Df Sum of Sq RSS AIC
## - pokoi 1 4.1403e+06 1.8268e+10 3678.0
## <none> 1.8263e+10 3680.0
## - typ.budynku 2 5.3152e+09 2.3579e+10 3727.1
## - dzielnica 2 1.6401e+10 3.4664e+10 3804.1
## - powierzchnia 1 3.4242e+10 5.2506e+10 3889.2
##
## Step: AIC=3678.01
## cena ~ powierzchnia + dzielnica + typ.budynku
##
## Df Sum of Sq RSS AIC
## <none> 1.8268e+10 3678.0
## + pokoi 1 4.1403e+06 1.8263e+10 3680.0
## - typ.budynku 2 5.3111e+09 2.3579e+10 3725.1
## - dzielnica 2 1.6398e+10 3.4665e+10 3802.1
## - powierzchnia 1 3.1032e+11 3.2858e+11 4253.9model_koniec_3##
## Call:
## lm(formula = cena ~ powierzchnia + dzielnica + typ.budynku, data = dane)
##
## Coefficients:
## (Intercept) powierzchnia dzielnicaKrzyki
## 92326.3 1998.3 -21463.6
## dzielnicaSrodmiescie typ.budynkuniski blok typ.budynkuwiezowiec
## -14310.0 11462.1 439.87 Model końcowy
model <- lm(cena ~ powierzchnia + dzielnica, dane)
summary_model <- summary(model)
summary_model##
## Call:
## lm(formula = cena ~ powierzchnia + dzielnica, data = dane)
##
## Residuals:
## Min 1Q Median 3Q Max
## -30510 -8484 -151 7367 35716
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 94228.68 2276.43 41.393 < 2e-16 ***
## powierzchnia 2024.73 38.79 52.194 < 2e-16 ***
## dzielnicaKrzyki -20933.74 1836.73 -11.397 < 2e-16 ***
## dzielnicaSrodmiescie -12723.11 2002.64 -6.353 1.44e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10970 on 196 degrees of freedom
## Multiple R-squared: 0.9362, Adjusted R-squared: 0.9352
## F-statistic: 957.9 on 3 and 196 DF, p-value: < 2.2e-16AIC(model)## [1] 4294.633Budując model regresji zależy nam na uzyskaniu modelu, który jak najlepiej będzie przewidywał wartości zmiennej objaśnianej (wysokie R2). W przypadku, gdy model zbudowany za pomocą regresji z krokową eliminacją zmiennych, bedzie posiadał zmienne o niskim poziomie istotności (lub nieistotne statystycznie) można je z modelu ręcznie usunąć, upraszczając w ten sposób model (im model posiada mniej zmiennych objaśniajacych tym łatwiej go praktycznie wykorzystać)
8 Predykcja na podstawie modelu
Mając dopasowany model liniowy możemy użyć otrzymanych ocen współczynników modelu w celu wyznaczenia predykcji - wyznaczenie oczekiwanej wartości zmiennej y dla zadanej wartości X. W przypadku mieszkań będzie to średnia cena mieszkania o pewnych parametrach
ndane <- data.frame(powierzchnia = c(40, 60), dzielnica = c("Biskupin", "Krzyki"))
pred <- predict(model, newdata = ndane) #podac jakic tam model z dzielnica i pow
pred## 1 2
## 175218.1 194779.0Na podstawie stworzonego modelu mieszkanie o powierzchni 40m2 w dzielnicy Biskupin będzie kosztować 175218.1 a cena mieszkania o powierzchni 60m2 w dzielnicy Krzyki to 194779.0.
9 Zadanie do samodzielnego rozwiązania
Celem zadania jest stworzenie optymalnego modelu przewidującego temperaturę wody w Parsęcie. Zbiór danych dane_tg.csv zawiera dane dotyczące temperatury wody w Parsęcie oraz temperatur gruntu na 5, 10, 20, 50 cm.
Używając opcji regresji z krokową eliminacją zmiennych (opcja backward, forward, both) proszę zbudować model dla TwP. W jakim stopniu utworzony model wyjaśnia temperaturę wody w Parsęcie?
Czy można ten model jeszcze uprościć, przy zachowaniu podobnego poziomu wyjaśnienia zmiennej TwP?