Cvičenie 11.1 Jednoduchá regresia

Zadanie: Urobte lineárny regresný model pre premennú dist z datasetu cars, pričom prediktorom je premenná speed. Overte predpoklady pre tento model a otestujte jeho štatistickú významnosť. Vizualizujte vzťah reálnych a predikovaných hodnôt pre target. Je medzi týmito premennými skutočne lineárny vzťah alebo sa dá nájsť iný, ktorý bude dávať lepšie výsledky?

Na internete som našiel, že údaje udávajú rýchlosť (mph) áut a vzdialenosti (ft), ktoré prešli na zastavenie. Údaje boli zaznamenané v 20. rokoch 20. storočia. Pre ľahšie vnímanie prevedieme míle na kilometre a a stopy na metre.

data1 <- datasets::cars
data1$speed <- data1$speed * 1.61
data1$dist <- data1$dist * 0.305
plot(data1$speed,data1$dist, main ="Zavislosť medzi rýchlosťou a brzdnou dráhou", xlab = "Rýchlosť auta [km/h]", ylab = "Brzdná dráha [m]")

Vysvetľovaná (závislá) premenná Y je brzdná dráha a vysvetľujúca (nezávislá) je rýchlosť auta. Parametre regresnej priamky odhadujeme metódou najmenších štvorcov. A dostaneme:

xm <- mean(data1$speed)
ym <- mean(data1$dist)

a <- sum( (data1$speed-xm)*(data1$dist-ym) ) / sum( (data1$speed-xm)^2 )
b <- ym - a*xm

plot(data1$speed,data1$dist, main ="Zavislosť medzi rýchlosťou a brzdnou dráhou", xlab = "Rýchlosť auta [km/h]", ylab = "Brzdná dráha [m]")
abline(a=b, b=a)

Model je teda v tvare Y= 0.745X - 5.3616

Model vzťah medzi premennými vždy popisuje s určitou chybou ε_i. Tá by mala byť vzhľadom na ďalšie testy normálne rozdelená so strednou hodnotou 0 a konštantným rozptylom σ^2. Chyby by mali by byť navzájom nezávislé. Budeme hovoriť o rezíduách, sú to rozdiely medzi nameranými hodnotami a hodnotami z modelu a sú to bodové odhady náhodných chýb ε_i modelu.

Model

linMod <- lm(data1$dist ~ data1$speed)
print(linMod)
## 
## Call:
## lm(formula = data1$dist ~ data1$speed)
## 
## Coefficients:
## (Intercept)  data1$speed  
##      -5.362        0.745
summary(linMod)
## 
## Call:
## lm(formula = data1$dist ~ data1$speed)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.8661 -2.9052 -0.6929  2.8105 13.1764 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.36162    2.06132  -2.601   0.0123 *  
## data1$speed  0.74496    0.07872   9.464 1.49e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.691 on 48 degrees of freedom
## Multiple R-squared:  0.6511, Adjusted R-squared:  0.6438 
## F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12

Prvú časť výstupu tvoria charakteristiky týkajúce sa rezíduí. Z týchto hodnôt sa zdá, že by mohli byť normálne rozdelené.

Ďalšia časť je odhad koeficientov. Hodnota koefiecientu pre X: 0.745 vyjadruje, že brzdná dráha auta sa zvýši o 0.745m, ak sa jeho rýchlosť zvýši o jeden kilometer za hodinu.

V stĺpci Pr(>|t|) sú príslušné p-hodnoty k testom nulovosti regresného koeficienta. H0: ai=0 H1: ai!=0. A teda vhodný je model, kde sú koeficienty štatisticky významne odlišné od 0, teda H0 zamietame. (Predpoklad vhodne použitého testu: Chyby musia byť z rozdelenia N(0,σ^2).) Oba koeficienty sú štatisticky významné. Majú svoje miesto v modeli.

Ďalej RSE, výberový reziduálny rozptyl s^2- keďže ide o rozptyl bodov okolo regresnej priamky (teda vlastne o rozptyl chýb), vyberáme model, kde je táto charakteristika minimálna. Tomuto hodnota 4.691 vyhovuje. Skutočná brzdná dráha sa odchyľuje od odhadnutých hodnôt ležiacich na regresnej priamke približne o ± 4.691 metra.

Ďalšími dôležitými výstupmi sú (adjusted) R-squared. Ide o (korigovaný) koeficient determinácie a popisujú, koľko % variability dát model popisuje. Tu je vhodný model, ktorý popisuje najviac variability, ideálne blízky 1. Na základe R^2 môžeme povedať tiež, že model vhodne popisuje vzťah medzi X a Y. Popisuje takmer 65% celkovej variability dát.

Posledná časť výstupu je test vhodnosti celého modelu. Je to vlastne analýza rozptylu (ANOVA)

H0: Všetky koeficienty sú nulové, resp. model nemá zmysel. H1: Model má zmysel.

(Opäť musí platiť predpoklad o normalite chýb.) Na základe p-hodnoty, ktorá je menšia ako hladina významnosti 0.05 H0 zamietame a model je významný.

Informačné kritériá AIC a BIC, majú zmysel hlavne, keď porovnávame viacero modelov. Volíme model s nižšou hodnotou.

linAIC<-AIC(linMod)
linBIC<-BIC(linMod)
linAIC
## [1] 300.4125
linBIC
## [1] 306.1486

Analýza rezíduí

linMod$residuals 
##          1          2          3          4          5          6          7 
##  1.1740853  3.6140853 -1.8140688  3.6759312  0.6465466 -2.3828381 -1.1422228 
##          8          9         10         11         12         13         14 
##  1.2977772  3.7377772 -2.6466074  0.7083926 -4.7609921 -2.9309921 -1.7109921 
##         15         16         17         18         19         20         21 
## -0.4909921 -2.3003768  0.1396232  0.1396232  3.7996232 -3.4997615 -0.4497615 
##         22         23         24         25         26         27         28 
##  6.8702385 12.9702385 -6.5291461 -4.6991461  3.8408539 -4.0685308 -1.6285308 
##         29         30         31         32         33         34         35 
## -5.2679155 -2.8279155  0.2220845 -3.4173001  0.8526999  6.9526999  9.3926999 
##         36         37         38         39         40         41         42 
## -6.4466848 -3.3966848  3.3133152 -8.8660695 -3.9860695 -2.7660695 -1.5460695 
##         43         44         45         46         47         48         49 
##  0.8939305 -0.8948388 -5.7542235 -2.0736082  4.6363918  4.9413918 13.1763918 
##         50 
##  1.3020072

Sú normálne rozdelené? Shapiro-Wilk test

shapiro.test(linMod$residuals) 
## 
##  Shapiro-Wilk normality test
## 
## data:  linMod$residuals
## W = 0.94509, p-value = 0.02152

Rezíduá nie je normalne rozdelené na hladine významnosti 0.05, zdá sa, že sme vytvorili nesprávny model.

Je stredná hodnota nula? t-test

t.test(linMod$residuals, mu=0)
## 
##  One Sample t-test
## 
## data:  linMod$residuals
## t = 2.2768e-17, df = 49, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -1.31943  1.31943
## sample estimates:
##    mean of x 
## 1.494898e-17

Áno, stredná hodnota je nula.

Sú náhodné? Runs test a turning point test

library(randtests)
runs.test(linMod$residuals)
## 
##  Runs Test
## 
## data:  linMod$residuals
## statistic = -0.28577, runs = 25, n1 = 25, n2 = 25, n = 50, p-value =
## 0.7751
## alternative hypothesis: nonrandomness
turning.point.test(linMod$residuals)
## 
##  Turning Point Test
## 
## data:  linMod$residuals
## statistic = -1.4961, n = 49, p-value = 0.1346
## alternative hypothesis: non randomness

Rezíduá sú náhodné.

Platí homoskedasticita? Golfeld-Quandtov test

𝐻0: homoskedasticita

𝐻1: heteroskedasticita

library(lmtest)
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
gqtest(linMod)
## 
##  Goldfeld-Quandt test
## 
## data:  linMod
## GQ = 1.5512, df1 = 23, df2 = 23, p-value = 0.1498
## alternative hypothesis: variance increases from segment 1 to 2

Rezíduá maju konštantne rozptyly.

Sú nekorelované? Resp. neexistuje autokorelácia? Durbin-Watsonov test

𝐻0:𝜌= 0 (nekorelované)

𝐻1:𝜌!= 0

library(car)
## Loading required package: carData
durbinWatsonTest(linMod)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.1604322      1.676225   0.206
##  Alternative hypothesis: rho != 0

Rezíduá sú nekorelované

Použitie modelu na predikcie

Odhad brzdnej dráhy pre auto pri rýchlosti 50 km/h (X=50)

0.745*50 - 5.3616
## [1] 31.8884

Očakávame, že náhodne vybrané auto bude mať pri rýchlosti 50 km/h brzdnú dráhu 31,89 m.

Môžeme odhadnúť aj interval spoľahlivosti pre regresné koeficienty.

confint(linMod, level = 0.95)
##                 2.5 %     97.5 %
## (Intercept) -9.506194 -1.2170538
## data1$speed  0.586692  0.9032268

Porovnanie modelov

Na to aby sme vybrali najlepší model popisujúci závislosť medzi znakmi, bežne odhadneme viacero modelov a porovnáme ich. Na základe vyššie popísaných charakteristík vyberáme najvhodnejší model (najčastejšie berieme do uváhy kombináciu viacerých aj z fyziky vieme, že vzťah medzi rýchlosťou auta a jeho brzdnou dráhou musí byť kvadratický).

Keby sme skúsili vzťah medzi rýchlosťou a brzdnou dráhou popísať kvadratickou funkciou:

linMod2 <- lm(data1$dist ~ data1$speed + I(data1$speed^2))  
print(linMod2)
## 
## Call:
## lm(formula = data1$dist ~ data1$speed + I(data1$speed^2))
## 
## Coefficients:
##      (Intercept)       data1$speed  I(data1$speed^2)  
##          0.75339           0.17301           0.01176
summary(linMod2)
## 
## Call:
## lm(formula = data1$dist ~ data1$speed + I(data1$speed^2))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.7595 -2.8012 -0.9723  1.4116 13.7713 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)
## (Intercept)      0.753392   4.519235   0.167    0.868
## data1$speed      0.173014   0.385365   0.449    0.656
## I(data1$speed^2) 0.011762   0.007762   1.515    0.136
## 
## Residual standard error: 4.629 on 47 degrees of freedom
## Multiple R-squared:  0.6673, Adjusted R-squared:  0.6532 
## F-statistic: 47.14 on 2 and 47 DF,  p-value: 5.852e-12
lin2AIC<-AIC(linMod2)
lin2BIC<-BIC(linMod2)

linAIC
## [1] 300.4125
lin2AIC
## [1] 300.0277
linBIC
## [1] 306.1486
lin2BIC
## [1] 307.6758

Keď porovnáme hodnoty týchto dvoch modelov, zdá sa, že kvadratická funkcia je trochu vhodnejšia. Ak rozdiely v informačných kritériách nie sú príliš veľké pre dva modely, vo všeobecnosti sa snažíme zvoliť jednoduchší model. Teda zvolime linearný model.

Vieme to aj testovať pomocou ANOVY.

anova(linMod, linMod2)
## Analysis of Variance Table
## 
## Model 1: data1$dist ~ data1$speed
## Model 2: data1$dist ~ data1$speed + I(data1$speed^2)
##   Res.Df    RSS Df Sum of Sq     F Pr(>F)
## 1     48 1056.2                          
## 2     47 1007.0  1    49.192 2.296 0.1364

Cvičenie 11.2 Násobná regresia – generované dáta

Zadanie: Nájdite najvhodnejší lineárny regresný model pre nižšie uvedené náhodne generované dáta. Overte predpoklady. Zdôvodnite prečo je vybraný model najvhodnejší. Testujte štatistickú významnosť modelu a parametrov.

n <- 100  
size <- rnorm(n, mean = 2000, sd = 500) # House size  
bedrooms <- sample(2:5, n, replace = TRUE) # Number of bedrooms  
bathrooms <- sample(1:4, n, replace = TRUE) # Number of bathrooms  
age <- rpois(n, lambda = 10) # Age of the house  
price <- 50000 + size * 150 + bedrooms * 10000 - age * 2000 + rnorm(n, mean = 0, sd = 30000) # Create target variable 
house_data <- data.frame( Size = size, Bedrooms = bedrooms, Bathrooms = bathrooms, Age = age, Price = price ) # Create data frame
house_data
##          Size Bedrooms Bathrooms Age    Price
## 1   2070.3349        2         2  11 341574.1
## 2   2007.0240        2         1   6 361494.2
## 3   2623.4083        3         3  13 475299.5
## 4   1511.6169        5         1  10 328028.9
## 5   2602.5957        4         3  13 465996.7
## 6   1736.0763        3         2   8 325702.6
## 7   1936.3165        3         4   6 376000.3
## 8   1802.0652        4         2  11 337253.8
## 9   2890.3861        5         2  10 468926.8
## 10  2542.0847        2         4  11 415052.0
## 11  2776.0463        2         1   5 482030.6
## 12  1998.5960        2         4  12 312637.1
## 13  1837.8262        3         2  12 351482.0
## 14  1068.8008        3         2  11 208379.9
## 15  2101.8926        2         3  10 391533.2
## 16  1989.8118        2         3  11 305024.3
## 17  2696.6803        5         3  15 413730.4
## 18  2647.9479        3         4   9 526095.3
## 19  2406.3160        3         3  10 433077.3
## 20  2500.3337        4         1  11 508457.7
## 21  2426.9029        2         4  11 386200.5
## 22  1843.6768        4         3   8 395811.5
## 23  2520.3100        2         1  10 438425.7
## 24  1983.0292        2         4   5 361476.0
## 25  2279.1819        5         4  11 423192.6
## 26  1528.8887        2         1   9 250762.0
## 27  2730.6837        2         2  13 471588.1
## 28  1773.3746        2         3  15 306060.3
## 29  1426.5910        3         1  11 203369.9
## 30  1562.6342        2         3   7 262940.9
## 31  1953.2279        3         3   9 392427.2
## 32  1634.6591        3         3  11 339449.6
## 33  2349.0240        4         1  12 385644.2
## 34  1298.5248        5         1   9 314466.3
## 35  1822.5467        2         1   9 340379.9
## 36  2224.0628        5         4   5 412982.5
## 37  2092.5156        2         1   5 340654.7
## 38  1758.9268        3         3  13 338240.0
## 39  1842.9384        4         3   6 375063.8
## 40  1661.4436        5         4   9 310686.0
## 41  2135.1579        5         2   9 397685.0
## 42  1944.3653        5         4   7 407197.4
## 43  2117.3157        2         1  10 394030.9
## 44  3088.2693        3         1   7 523600.3
## 45  1678.5543        3         2   9 286034.0
## 46  1183.7315        4         1  11 252114.4
## 47  1486.7707        5         4   7 342397.3
## 48  1787.8517        5         2   3 366047.1
## 49   833.8945        2         3   9 206174.4
## 50  2629.3783        2         2  12 441220.3
## 51  2718.1907        4         4  19 456885.0
## 52  1530.5758        4         4  17 303791.8
## 53  1398.1705        4         3  10 257526.8
## 54  1791.4690        5         2   9 355223.7
## 55  1953.1133        2         4   6 393562.8
## 56  2736.6161        4         3   3 499468.2
## 57  2801.5687        3         1  12 465708.7
## 58  1874.9392        4         3  15 350858.6
## 59  2393.6240        2         1  10 478321.0
## 60  1377.7502        2         2   8 228203.3
## 61  2338.7413        5         2  10 427278.1
## 62  2059.1355        2         3   6 381969.4
## 63  1516.4716        4         3  17 282234.7
## 64  1739.0277        4         1  12 351980.8
## 65  1663.0484        3         1  12 327579.7
## 66  1729.4104        5         1   6 305926.3
## 67  1533.3999        3         1   9 270166.0
## 68  1692.2507        3         1  10 322174.2
## 69  2047.5894        5         3  11 358116.3
## 70  1338.4133        5         2  10 314569.0
## 71  1076.4833        3         2   8 233620.6
## 72  1470.2414        4         3  10 290016.4
## 73   782.1324        5         3  11 207221.7
## 74  1724.0549        4         1  16 343065.0
## 75  1309.8860        5         2   8 295693.9
## 76  2267.8515        2         1  11 400895.5
## 77  1716.2269        2         1   6 293007.5
## 78  1707.5533        3         3   5 358373.4
## 79  2634.8761        2         4  24 423466.1
## 80  1306.8145        5         3   7 303785.4
## 81  1851.6418        2         3   8 320197.8
## 82  1808.1816        2         2  12 263063.1
## 83  1956.8325        5         1  11 310994.5
## 84  1965.5953        4         1   8 377472.4
## 85  1956.9714        3         2   9 395171.0
## 86  2276.5158        4         1  10 382865.6
## 87  1997.0250        2         1   7 363122.9
## 88  1919.2264        5         2   9 380758.5
## 89  1660.4392        3         4   7 330657.9
## 90  2069.4058        4         3  11 393801.5
## 91  1908.1227        3         2  12 414829.9
## 92   905.6552        4         3   6 226036.6
## 93  2094.5244        5         1   7 368960.3
## 94  2334.3592        4         3  15 462832.7
## 95  2607.9179        2         1  12 415716.1
## 96  2079.7282        4         2   8 372204.6
## 97  1593.8333        3         1  10 318228.7
## 98  2408.1987        2         4   6 404974.3
## 99  2204.0545        4         1  11 432967.8
## 100 1405.7917        2         4  10 246796.0
plot(Price ~ ., data = house_data)

Model:

linMod1 <- lm(house_data$Price ~ house_data$Size + house_data$Bedrooms + house_data$Bathrooms + house_data$Age)
print(linMod1)
## 
## Call:
## lm(formula = house_data$Price ~ house_data$Size + house_data$Bedrooms + 
##     house_data$Bathrooms + house_data$Age)
## 
## Coefficients:
##          (Intercept)       house_data$Size   house_data$Bedrooms  
##              51316.6                 147.6               10285.7  
## house_data$Bathrooms        house_data$Age  
##               2569.3               -1965.9
summary(linMod1)
## 
## Call:
## lm(formula = house_data$Price ~ house_data$Size + house_data$Bedrooms + 
##     house_data$Bathrooms + house_data$Age)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -70296 -20203   3060  17305  70252 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)          51316.58   18085.87   2.837 0.005561 ** 
## house_data$Size        147.59       6.19  23.844  < 2e-16 ***
## house_data$Bedrooms  10285.74    2559.84   4.018 0.000117 ***
## house_data$Bathrooms  2569.33    2629.83   0.977 0.331052    
## house_data$Age       -1965.90     892.09  -2.204 0.029961 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28920 on 95 degrees of freedom
## Multiple R-squared:  0.8576, Adjusted R-squared:  0.8516 
## F-statistic: 143.1 on 4 and 95 DF,  p-value: < 2.2e-16

Prvú časť výstupu tvoria charakteristiky týkajúce sa rezíduí. Z týchto hodnôt sa zdá, že by mohli byť normálne rozdelené.

Ďalšia časť je odhad koeficientov. Hodnoty koefiecientov pre X vyjadruju, ak ktorúkoľvek premennú zvýšime o 1, tak sa cena domu zvýši o hodnotu zodpovedajúcej premennej.

V stĺpci Pr(>|t|) sú príslušné p-hodnoty k testom nulovosti regresného koeficienta. H0: ai=0 H1: ai!=0. A teda vhodný je model, kde sú koeficienty štatisticky významne odlišné od 0, teda H0 zamietame. (Predpoklad vhodne použitého testu: Chyby musia byť z rozdelenia N(0,σ^2).) Nie všetky koeficienty sú štatisticky významné. Vidíme, že počet kúpeľní nemá vplyv na cenu domu.

Ďalej RSE, výberový reziduálny rozptyl s^2- keďže ide o rozptyl bodov okolo regresnej priamky (teda vlastne o rozptyl chýb), vyberáme model, kde je táto charakteristika minimálna. Tomuto hodnota 31400 vyhovuje. Skutočná cena sa odchyľuje od odhadnutých hodnôt ležiacich na regresnej priamke približne o ± 31400 .

Ďalšími dôležitými výstupmi sú (adjusted) R-squared. Ide o (korigovaný) koeficient determinácie a popisujú, koľko % variability dát model popisuje. Tu je vhodný model, ktorý popisuje najviac variability, ideálne blízky 1. Na základe R^2 môžeme povedať tiež, že model vhodne popisuje vzťah medzi X a Y. Popisuje takmer 83,8% celkovej variability dát.

Posledná časť výstupu je test vhodnosti celého modelu. Je to vlastne analýza rozptylu (ANOVA)

H0: Všetky koeficienty sú nulové, resp. model nemá zmysel. H1: Model má zmysel.

(Opäť musí platiť predpoklad o normalite chýb.) Na základe p-hodnoty, ktorá je menšia ako hladina významnosti 0.05 H0 zamietame a model je významný.

Informačné kritériá AIC a BIC, majú zmysel hlavne, keď porovnávame viacero modelov. Volíme model s nižšou hodnotou.

linAIC1<-AIC(linMod1)
linBIC1<-BIC(linMod1)
linAIC1
## [1] 2345.097
linBIC1
## [1] 2360.728

Analýza rezíduí

linMod1$residuals 
##           1           2           3           4           5           6 
## -19384.3864   2619.4303  23790.7002  19276.7231   7273.8240  -2106.7630 
##           7           8           9          10          11          12 
##   9567.3712  -4682.7392 -45884.5282 -20669.8198   7691.4136 -40906.3306 
##          13          14          15          16          17          18 
##  16519.2398 -15049.8764  21381.9823 -46619.2905 -65232.1422  60531.7813 
##          19          20          21          22          23          24 
##   7711.0126  66034.3876 -32521.8408  39266.4723  11659.6664  -3531.2929 
##          25          26          27          28          29          30 
##  -4585.0683 -31648.0034  17101.8434  -5776.1368 -70296.1160 -33520.0082 
##          31          32          33          34          35          36 
##  31965.4263  29936.5456 -32481.6677  35198.1024  14629.5493 -18455.6968 
##          37          38          39          40          41          42 
## -32803.4523  14318.3984  14695.9599 -29852.7412  -7629.6520  20970.9828 
##          43          44          45          46          47          48 
##  26742.0283  -3173.2029 -31319.8622   4005.8533  23706.3892    195.2923 
##          49          50          51          52          53          54 
##  21198.5521   -280.3888  -9672.2019   8580.5741 -29334.9612    633.4694 
##          55          56          57          58          59          60 
##  34936.6106   1306.4957  -8921.6452   3460.9408  70252.3074 -36435.6962 
##          61          62          63          64          65          66 
##  -6117.0861  10265.0021  -8325.5840  23883.0427  20981.3240 -42833.2666 
##          67          68          69          70          71          72 
## -23195.4794   7334.0844 -32911.7881  28810.2227   3159.2905  -7482.1734 
##          73          74          75          76          77          78 
##   2959.9357  25040.6997  10213.5978  13355.2916 -22949.1068  26306.6408 
##          79          80          81          82          83          84 
##   -393.8374  14223.2714 -16951.2436 -57238.8168 -61500.3058   8072.4011 
##          85          86          87          88          89          90 
##  36726.0619 -28490.7782   7689.7643   7312.7520   6907.0891   9839.3325 
##          91          92          93          94          95          96 
##  69492.1227   4000.5549 -31719.7481  47630.1569 -20048.0146 -16609.4594 
##          97          98          99         100 
##  17913.8841 -20816.9897  34271.7825 -23188.4426

Sú normálne rozdelené? Shapiro-Wilk test

shapiro.test(linMod1$residuals) 
## 
##  Shapiro-Wilk normality test
## 
## data:  linMod1$residuals
## W = 0.98579, p-value = 0.3613

Rezíduá je normalne rozdelené na hladine významnosti 0.05.

Je stredná hodnota nula? t-test

t.test(linMod1$residuals, mu=0)
## 
##  One Sample t-test
## 
## data:  linMod1$residuals
## t = -3.1417e-16, df = 99, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -5620.687  5620.687
## sample estimates:
##     mean of x 
## -8.899548e-13

Áno, stredná hodnota je nula.

Sú náhodné? Runs test a turning point test

library(randtests)
runs.test(linMod1$residuals)
## 
##  Runs Test
## 
## data:  linMod1$residuals
## statistic = 1.6081, runs = 59, n1 = 50, n2 = 50, n = 100, p-value =
## 0.1078
## alternative hypothesis: nonrandomness
turning.point.test(linMod1$residuals)
## 
##  Turning Point Test
## 
## data:  linMod1$residuals
## statistic = -0.55848, n = 100, p-value = 0.5765
## alternative hypothesis: non randomness

Rezíduá sú náhodné.

Platí homoskedasticita? Golfeld-Quandtov test

𝐻0: homoskedasticita

𝐻1: heteroskedasticita

library(lmtest)
gqtest(linMod1)
## 
##  Goldfeld-Quandt test
## 
## data:  linMod1
## GQ = 0.90465, df1 = 45, df2 = 45, p-value = 0.6309
## alternative hypothesis: variance increases from segment 1 to 2

Rezíduá maju konštantne rozptyly.

Sú nekorelované? Resp. neexistuje autokorelácia? Durbin-Watsonov test

𝐻0:𝜌= 0 (nekorelované)

𝐻1:𝜌!= 0

library(car)
durbinWatsonTest(linMod1)
##  lag Autocorrelation D-W Statistic p-value
##    1     -0.09571848      2.179938    0.34
##  Alternative hypothesis: rho != 0

Rezíduá sú nekorelované

Porovnanie modelov

Na to aby sme vybrali najlepší model popisujúci závislosť medzi znakmi, bežne odhadneme viacero modelov a porovnáme ich. Na základe vyššie popísaných charakteristík vyberáme najvhodnejší model (najčastejšie berieme do uváhy kombináciu viacerých).

Keby sme skúsili vzťah medzi cenou a inými premennými popísať kvadratickou funkciou:

linMod2 <- lm(house_data$Price ~ house_data$Size + I(house_data$Size^2) + house_data$Bedrooms + I(house_data$Bedrooms^2) + house_data$Bathrooms + I(house_data$Bathrooms^2) + house_data$Age + + I(house_data$Age^2))  
print(linMod2)
## 
## Call:
## lm(formula = house_data$Price ~ house_data$Size + I(house_data$Size^2) + 
##     house_data$Bedrooms + I(house_data$Bedrooms^2) + house_data$Bathrooms + 
##     I(house_data$Bathrooms^2) + house_data$Age + +I(house_data$Age^2))
## 
## Coefficients:
##               (Intercept)            house_data$Size  
##                -3.959e+04                  1.547e+02  
##      I(house_data$Size^2)        house_data$Bedrooms  
##                -1.509e-03                  6.803e+04  
##  I(house_data$Bedrooms^2)       house_data$Bathrooms  
##                -8.362e+03                  2.409e+03  
## I(house_data$Bathrooms^2)             house_data$Age  
##                 1.307e+02                 -2.675e+03  
##       I(house_data$Age^2)  
##                 9.188e+00
summary(linMod2)
## 
## Call:
## lm(formula = house_data$Price ~ house_data$Size + I(house_data$Size^2) + 
##     house_data$Bedrooms + I(house_data$Bedrooms^2) + house_data$Bathrooms + 
##     I(house_data$Bathrooms^2) + house_data$Age + +I(house_data$Age^2))
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -77693 -17044    259  18741  76995 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -3.959e+04  5.347e+04  -0.740  0.46095    
## house_data$Size            1.547e+02  3.773e+01   4.100 8.98e-05 ***
## I(house_data$Size^2)      -1.509e-03  9.545e-03  -0.158  0.87474    
## house_data$Bedrooms        6.803e+04  2.054e+04   3.313  0.00133 ** 
## I(house_data$Bedrooms^2)  -8.362e+03  2.949e+03  -2.836  0.00563 ** 
## house_data$Bathrooms       2.409e+03  1.471e+04   0.164  0.87032    
## I(house_data$Bathrooms^2)  1.307e+02  3.038e+03   0.043  0.96579    
## house_data$Age            -2.675e+03  3.365e+03  -0.795  0.42881    
## I(house_data$Age^2)        9.188e+00  1.435e+02   0.064  0.94910    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 28310 on 91 degrees of freedom
## Multiple R-squared:  0.8693, Adjusted R-squared:  0.8578 
## F-statistic: 75.65 on 8 and 91 DF,  p-value: < 2.2e-16
lin2AIC<-AIC(linMod2)
lin2BIC<-BIC(linMod2)

linAIC1
## [1] 2345.097
lin2AIC
## [1] 2344.541
linBIC1
## [1] 2360.728
lin2BIC
## [1] 2370.593

Keď porovnáme hodnoty týchto dvoch modelov, zdá sa, že rozdiely v informačných kritériách nie sú príliš veľké. Teda zvolime linearný model.

Cvičenie 11.3 Násobná regresia - dataset

Zadanie: Urobte lineárny regresný model pre premennú mpg z datasetu mtcars, pričom prediktormi budú premenné horsepower, weight a cylinders. Overte predpoklady pre tento model a otestujte jeho štatistickú významnosť. Vizualizujte vzťah reálnych a predikovaných hodnôt pre target. Dávajú iné premenné z datasetu ako prediktory lepší model a ak áno, tak ktoré?

data3 <- datasets::mtcars
data3$lp100km <- 100/(data3$mpg*(1.609344/3.78541)) #prelozim na l/100km
data3$wt <- data3$wt * 1000 / 2.2046 #prelozim na kg
plot(lp100km ~ wt+hp+cyl, data = data3)

Model:

linMod3 <- lm(data3$lp100km ~ data3$wt + data3$cyl + data3$hp)
print(linMod3)
## 
## Call:
## lm(formula = data3$lp100km ~ data3$wt + data3$cyl + data3$hp)
## 
## Coefficients:
## (Intercept)     data3$wt    data3$cyl     data3$hp  
##     1.17574      0.00564      0.18981      0.01482
summary(linMod3)
## 
## Call:
## lm(formula = data3$lp100km ~ data3$wt + data3$cyl + data3$hp)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.7763 -1.1286  0.2647  0.9224  3.4531 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.175742   1.124898   1.045   0.3049    
## data3$wt    0.005640   0.001028   5.487 7.34e-06 ***
## data3$cyl   0.189814   0.346822   0.547   0.5885    
## data3$hp    0.014824   0.007477   1.983   0.0573 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.581 on 28 degrees of freedom
## Multiple R-squared:  0.8487, Adjusted R-squared:  0.8325 
## F-statistic: 52.36 on 3 and 28 DF,  p-value: 1.322e-11

Prvú časť výstupu tvoria charakteristiky týkajúce sa rezíduí. Z týchto hodnôt sa zdá, že by mohli byť normálne rozdelené.

Ďalšia časť je odhad koeficientov. Hodnoty koefiecientov pre X vyjadruju, ak ktorúkoľvek premennú zvýšime o 1, tak sa spotreba paliva zvýši o hodnotu zodpovedajúcej premennej.

V stĺpci Pr(>|t|) sú príslušné p-hodnoty k testom nulovosti regresného koeficienta. H0: ai=0 H1: ai!=0. A teda vhodný je model, kde sú koeficienty štatisticky významne odlišné od 0, teda H0 zamietame. (Predpoklad vhodne použitého testu: Chyby musia byť z rozdelenia N(0,σ^2).) Nie všetky koeficienty sú štatisticky významné. Vidíme, že počet valcov a výkon motora nie su štatisticky významne odlišné od 0.

Ďalej RSE, výberový reziduálny rozptyl s^2- keďže ide o rozptyl bodov okolo regresnej priamky (teda vlastne o rozptyl chýb), vyberáme model, kde je táto charakteristika minimálna. Tomuto hodnota 1,581 vyhovuje. Spotreba paliva sa odchyľuje od odhadnutých hodnôt ležiacich na regresnej priamke približne o ± 1,58.

Ďalšími dôležitými výstupmi sú (adjusted) R-squared. Ide o (korigovaný) koeficient determinácie a popisujú, koľko % variability dát model popisuje. Tu je vhodný model, ktorý popisuje najviac variability, ideálne blízky 1. Na základe R^2 môžeme povedať tiež, že model vhodne popisuje vzťah medzi X a Y. Popisuje takmer 84,9% celkovej variability dát.

Posledná časť výstupu je test vhodnosti celého modelu. Je to vlastne analýza rozptylu (ANOVA)

H0: Všetky koeficienty sú nulové, resp. model nemá zmysel. H1: Model má zmysel.

(Opäť musí platiť predpoklad o normalite chýb.) Na základe p-hodnoty, ktorá je menšia ako hladina významnosti 0.05 H0 zamietame a model je významný.

Informačné kritériá AIC a BIC, majú zmysel hlavne, keď porovnávame viacero modelov. Volíme model s nižšou hodnotou.

linAIC3<-AIC(linMod3)
linBIC3<-BIC(linMod3)
linAIC3
## [1] 125.8594
linBIC3
## [1] 133.188

Keďže počet valcov a výkon motora vyšli štatisticky nevýznamné, skúsime ešte model bez tychto premennych.

linMod4 <- lm(data3$lp100km ~ data3$wt)
print(linMod4)
## 
## Call:
## lm(formula = data3$lp100km ~ data3$wt)
## 
## Coefficients:
## (Intercept)     data3$wt  
##    1.451025     0.007746
summary(linMod4)
## 
## Call:
## lm(formula = data3$lp100km ~ data3$wt)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -4.2300 -1.0420  0.2214  1.0668  2.7422 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.4510246  1.1043204   1.314    0.199    
## data3$wt    0.0077460  0.0007249  10.685 9.57e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.791 on 30 degrees of freedom
## Multiple R-squared:  0.7919, Adjusted R-squared:  0.785 
## F-statistic: 114.2 on 1 and 30 DF,  p-value: 9.566e-12

Porovnanie modelov

c(AIC(linMod3),AIC(linMod4))
## [1] 125.8594 132.0596
c(BIC(linMod3),BIC(linMod4))
## [1] 133.1880 136.4568

Volíme druhý model - má trochu vačšie hodnoty informačných kritérií ale má menej premenných.

Analýza rezíduí

linMod4$residuals 
##           1           2           3           4           5           6 
##  0.54413245 -0.35182582  0.71393607 -1.75579553 -0.95936340 -0.61267396 
##           7           8           9          10          11          12 
##  2.45412403 -3.01934958 -2.20232027 -1.28692367 -0.32338123 -1.40887553 
##          13          14          15          16          17          18 
## -0.96039810  0.74234787  2.71955179  2.10819203 -4.23002826 -1.92115479 
##          19          20          21          22          23          24 
##  0.61189113 -0.95993064  0.82825340  1.35636487  1.95452672  2.74219650 
##          25          26          27          28          29          30 
## -2.70991622  0.36614790  0.07666346  0.97027444  2.29797619  0.75622994 
##          31          32 
##  1.68652435 -0.22739612

Sú normálne rozdelené? Shapiro-Wilk test

shapiro.test(linMod4$residuals) 
## 
##  Shapiro-Wilk normality test
## 
## data:  linMod4$residuals
## W = 0.9741, p-value = 0.6192

Rezíduá je normalne rozdelené na hladine významnosti 0.05.

Je stredná hodnota nula? t-test

t.test(linMod4$residuals, mu=0)
## 
##  One Sample t-test
## 
## data:  linMod4$residuals
## t = -1.5196e-16, df = 31, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.635375  0.635375
## sample estimates:
##     mean of x 
## -4.733898e-17

Áno, stredná hodnota je nula.

Sú náhodné? Runs test a turning point test

library(randtests)
runs.test(linMod3$residuals)
## 
##  Runs Test
## 
## data:  linMod3$residuals
## statistic = 0, runs = 17, n1 = 16, n2 = 16, n = 32, p-value = 1
## alternative hypothesis: nonrandomness
turning.point.test(linMod4$residuals)
## 
##  Turning Point Test
## 
## data:  linMod4$residuals
## statistic = -0.86333, n = 32, p-value = 0.388
## alternative hypothesis: non randomness

Rezíduá sú náhodné.

Platí homoskedasticita? Golfeld-Quandtov test

𝐻0: homoskedasticita

𝐻1: heteroskedasticita

library(lmtest)
gqtest(linMod4)
## 
##  Goldfeld-Quandt test
## 
## data:  linMod4
## GQ = 1.5468, df1 = 14, df2 = 14, p-value = 0.2123
## alternative hypothesis: variance increases from segment 1 to 2

Rezíduá maju konštantne rozptyly.

Sú nekorelované? Resp. neexistuje autokorelácia? Durbin-Watsonov test

𝐻0:𝜌= 0 (nekorelované)

𝐻1:𝜌!= 0

library(car)
durbinWatsonTest(linMod4)
##  lag Autocorrelation D-W Statistic p-value
##    1       0.1322653      1.731857   0.406
##  Alternative hypothesis: rho != 0

Rezíduá sú nekorelované

Vizualizacia:

a=0.00746
b=1.451025
plot(data3$lp100km ~ data3$wt, xlim = c(0,2000), ylim = c(0,25))
abline(a=b, b=a)

Pokúsime sa porovnať náš model s inými modelmi obsahujúcimi iné premenné z datasetu ako prediktory.

linMod5 <- lm(data3$lp100km ~ data3$disp + data3$drat + data3$qsec)
linMod6 <- lm(data3$lp100km ~ data3$am + data3$gear)
linMod7 <- lm(data3$lp100km ~ data3$hp + data3$qsec + data3$vs)
linMod8 <- lm(data3$lp100km ~ data3$disp + data3$qsec + data3$wt)
c(AIC(linMod4),AIC(linMod5),AIC(linMod6),AIC(linMod7),AIC(linMod8))
## [1] 132.0596 138.6195 172.9527 147.8050 125.4084
c(BIC(linMod4),BIC(linMod5),BIC(linMod6),BIC(linMod7),BIC(linMod8))
## [1] 136.4568 145.9482 178.8156 155.1336 132.7371

Vídime, že najlepši 8. model - má nižšie hodnoty informačných kritérií ale ma viac premenných.