Vícenásobná regrese a mixed-effects modely

Data ke stažení

data.zip

Použité a doporučené zdroje

Winter 2012: PSEUDOREPLICATION IN PHONETIC RESEARCH

Winter: Linear models and linear mixed effects models in R: Tutorial 1

Winter: A very basic tutorial for performing linear mixed effects analyses (Tutorial 2)

Použity byly také příklady z přednášky Regression Models v rámci specializace Data Science na Coursera (prof. Brian Caffo z Johns Hopkins Bloomberg School of Public Health).

Faktorové proměnné v lineárních modelech

Insect Sprays

library(readxl)
data <- read_excel("sprays.xlsx")
boxplot(count ~ spray, data)

fit <- lm(count ~ spray, data)
summary(fit)
## 
## Call:
## lm(formula = count ~ spray, data = data)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -8.333 -1.958 -0.500  1.667  9.333 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  14.5000     1.1322  12.807  < 2e-16 ***
## sprayB        0.8333     1.6011   0.520    0.604    
## sprayC      -12.4167     1.6011  -7.755 7.27e-11 ***
## sprayD       -9.5833     1.6011  -5.985 9.82e-08 ***
## sprayE      -11.0000     1.6011  -6.870 2.75e-09 ***
## sprayF        2.1667     1.6011   1.353    0.181    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.922 on 66 degrees of freedom
## Multiple R-squared:  0.7244, Adjusted R-squared:  0.7036 
## F-statistic:  34.7 on 5 and 66 DF,  p-value: < 2.2e-16

Protože lm() vidí, že spray je faktorová proměnná se 6 hodnotami (A-F), vyrobí si uvnitř pět dummy variables (binární ano/ne). Pro první hodnotu není nutné dummy variable tvořit, bude přímo součástí interceptu, zatímco těch dalších 5 bude vyhodnoceno jako rozdíl oproti interceptu.

Při vyhodnocování výsledné regresní rovnice u každé dummy variable odpovídáme na otázku “Je to ono?”. Pokud není, násobíme koeficient nulou 0, pokud je, násobíme jedničkou 1. A všechny dílčí členy sčítáme.

Obecný postup sestavování rovnice: věci v řádcích se násobí, řádky mezi sebou sčítají. Za (Intercept) vždy dosazujeme 1. Dummy proměnné (poznáme podle slepeniny názvu původní proměnné a konkrétní hodnoty) = 1, pokud je to ono, jinak 0.

Příklady výsledných rovnic pro spray A, B a E:

# A
count <- 14.5000 + 0.8333*0 - 12.4167*0 - 9.5833*0 - 11.0000*0 + 2.1667*0
count <- 14.5000

# B
count <- 14.5000 + 0.8333*1 - 12.4167*0 - 9.5833*0 - 11.0000*0 + 2.1667*0
count <- 15.3333

# E
count <- 14.5000 + 0.8333*0 - 12.4167*0 - 9.5833*0 - 11.0000*1 + 2.1667*0
count <- 3.5

Interpretace? Lineární model s referencí skupiny A (R si volí samo vždy první skupinu podle pořadí levelů). Koeficienty říkají, jak se změní počet ve srovnání s A (což je intercept), když se použije ten daný sprej. Je vidět, že C až E jsou níž, což odpovídá i obrázku.

Mimochodem: v tomto jednoduchém případě se jedná přímo o průměry daných skupin, takže se použití lm nejeví výhodné. Ale jakmile model zkomplikujeme, přidáme interakce apod., pochopíme smysl v plné síle.

Post-hoc testy: párové testy s korekcí pro multiple-testing

library(emmeans)
em <- emmeans(fit, pairwise ~ spray) # pozn. pokud nepotřebujeme em[[2]], stačí jen
                                     # em <- emmeans(fit, ~spray)
em
## $emmeans
##  spray emmean   SE df lower.CL upper.CL
##  A      14.50 1.13 66   12.240    16.76
##  B      15.33 1.13 66   13.073    17.59
##  C       2.08 1.13 66   -0.177     4.34
##  D       4.92 1.13 66    2.656     7.18
##  E       3.50 1.13 66    1.240     5.76
##  F      16.67 1.13 66   14.406    18.93
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast estimate  SE df t.ratio p.value
##  A - B      -0.833 1.6 66  -0.520  0.9952
##  A - C      12.417 1.6 66   7.755  <.0001
##  A - D       9.583 1.6 66   5.985  <.0001
##  A - E      11.000 1.6 66   6.870  <.0001
##  A - F      -2.167 1.6 66  -1.353  0.7542
##  B - C      13.250 1.6 66   8.276  <.0001
##  B - D      10.417 1.6 66   6.506  <.0001
##  B - E      11.833 1.6 66   7.391  <.0001
##  B - F      -1.333 1.6 66  -0.833  0.9603
##  C - D      -2.833 1.6 66  -1.770  0.4921
##  C - E      -1.417 1.6 66  -0.885  0.9489
##  C - F     -14.583 1.6 66  -9.108  <.0001
##  D - E       1.417 1.6 66   0.885  0.9489
##  D - F     -11.750 1.6 66  -7.339  <.0001
##  E - F     -13.167 1.6 66  -8.223  <.0001
## 
## P value adjustment: tukey method for comparing a family of 6 estimates
plot(em[[1]], comparisons = TRUE)

plot(em[[2]], comparisons = TRUE)
## Warning: Comparison discrepancy in group "1", (B - D) - (D - E):
##     Target overlap = 0.0788, overlap on graph = -0.0534
## Warning: Comparison discrepancy in group "1", (C - D) - (D - F):
##     Target overlap = 0.0873, overlap on graph = -0.0879

# pozn. pro složitější situace lze v plot přidat parametr by = "název_další_proměnné",
# pak to rozdělí graf na více panelů, ovšem srovnání lze dělat vždy jen v rámci 1 panelu

Důležitý je překryv šipek, nikoliv modrých pruhů (to jsou konfidenční intervaly EMM - Estimated marginal means, které jsou pro párová porovnávání velmi zavádějící).

Více viz

https://cran.r-project.org/web/packages/emmeans/vignettes/comparisons.html

https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html

https://cran.r-project.org/web/packages/emmeans/vignettes/sophisticated.html

Vícenásobná regrese

„Existuje vztah mezi užíváním žvýkaček a zhoršenou funkcí plic.“ => skepse

There are known knowns. These are things we know that we know. There are known unknowns. That is to say, there are things that we know we don’t know. But there are also unknown unknowns. There are things we don’t know we don’t know. Donald Rumsfeld

(Known knowns) Regressors that we know we should check to include in the model and have. (Known Unknowns) Regressors that we would like to include in the model, but don’t have. (Unknown Unknowns) Regressors that we don’t even know about that we should have included in the model.

Spíše: kuřáci asi používají častěji žvýkačky než nekuřáci a kouření je propojeno se zhoršenou funkcí plic. Jak to vyřešit? Rozdělit na skupinu kuřáci a nekuřáci, v každé pak hledat vztah.

Švýcarská data porodnosti (1888) ze 47 francouzsky mluvících “provincií”, popis viz ?swiss.

data <- read_excel("swiss.xlsx")

fit <- lm(Fertility ~ . , data)  # model, tečka znamená "vše"
summary(fit)
## 
## Call:
## lm(formula = Fertility ~ ., data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.2743  -5.2617   0.5032   4.1198  15.3213 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      66.91518   10.70604   6.250 1.91e-07 ***
## Agriculture      -0.17211    0.07030  -2.448  0.01873 *  
## Examination      -0.25801    0.25388  -1.016  0.31546    
## Education        -0.87094    0.18303  -4.758 2.43e-05 ***
## Catholic          0.10412    0.03526   2.953  0.00519 ** 
## Infant.Mortality  1.07705    0.38172   2.822  0.00734 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.165 on 41 degrees of freedom
## Multiple R-squared:  0.7067, Adjusted R-squared:  0.671 
## F-statistic: 19.76 on 5 and 41 DF,  p-value: 5.594e-10

Příklad interpretace: agriculture je v procentech (0-100), znamená procento mužů zaměstnaných v zemědělství.

Estimate je -0.1721. Odhadujeme pokles porodnosti o 0.17 na každé 1 % nárůstu počtu mužů v zemědělství, zatímco ostatní hodnoty necháme stejné.

Zajímavé je, že samotná lineární regrese Fertility vs. Agriculture by vyšla přesně naopak.

fit2 <- lm(Fertility ~ Agriculture, data)
summary(fit2)
## 
## Call:
## lm(formula = Fertility ~ Agriculture, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.5374  -7.8685  -0.6362   9.0464  24.4858 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 60.30438    4.25126  14.185   <2e-16 ***
## Agriculture  0.19420    0.07671   2.532   0.0149 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.82 on 45 degrees of freedom
## Multiple R-squared:  0.1247, Adjusted R-squared:  0.1052 
## F-statistic: 6.409 on 1 and 45 DF,  p-value: 0.01492

Jak to? Toto není analýza kauzality! Typický příklad Simpsonova paradoxu.

Umělý příklad

pocetMesicu <- 100
cas <- 1 : pocetMesicu;

penize <- .01*cas + runif(pocetMesicu, -.1, .1); # S rostoucím časem přibývají peníze, i když je tam nějaký šum

stesti <- -penize + cas + rnorm(pocetMesicu, sd = .01)
# S rostoucím časem narůstá štěstí (děti atd.), ale více peněz je více starostí, takže
# rostoucí peníze ovlivňují štěstí negativně.

# peníze jsou zhruba 1/100x čas. Když modelujeme stesti na penezich (bez casu), tak model
# nabyde pocitu, že pozorovaná věc je výsledkem zhruba 100x peníze. I když je ve skutečnosti
# hlavně důsledkem velkého času a peníze jsou jen důsledek, nikoliv příčinou.

# A pozor, to vše je podle modelu velice statisticky významné (extrémně nízké p-hodnoty).

fit1 <- lm(stesti ~ penize)
fit2 <- lm(stesti ~ penize + cas)

summary(fit1)
## 
## Call:
## lm(formula = stesti ~ penize)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -10.1931  -4.6313   0.0673   4.9298  11.6316 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    3.468      1.113   3.115  0.00241 ** 
## penize        92.624      1.899  48.774  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 5.74 on 98 degrees of freedom
## Multiple R-squared:  0.9604, Adjusted R-squared:   0.96 
## F-statistic:  2379 on 1 and 98 DF,  p-value: < 2.2e-16
summary(fit2)
## 
## Call:
## lm(formula = stesti ~ penize + cas)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -0.0270848 -0.0056340 -0.0006313  0.0059300  0.0222275 
## 
## Coefficients:
##               Estimate Std. Error  t value Pr(>|t|)    
## (Intercept) -0.0005095  0.0020863   -0.244    0.808    
## penize      -1.0159708  0.0172453  -58.913   <2e-16 ***
## cas          1.0001368  0.0001806 5538.231   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.01026 on 97 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 3.876e+08 on 2 and 97 DF,  p-value: < 2.2e-16

Interakce

Zpátky k datům “swiss”. Jsou regiony, kde evidentně převládají katolíci, a naopak regiony, kde převládají protestanti. Proto vytváříme binární proměnnou “katolíci”.

data <- read_excel("swiss.xlsx")
data$CatholicBin <- factor(data$Catholic > 50)

Bez faktoru religion

fit1 <- lm(Fertility ~ Agriculture, data)
summary(fit1)
## 
## Call:
## lm(formula = Fertility ~ Agriculture, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -25.5374  -7.8685  -0.6362   9.0464  24.4858 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 60.30438    4.25126  14.185   <2e-16 ***
## Agriculture  0.19420    0.07671   2.532   0.0149 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.82 on 45 degrees of freedom
## Multiple R-squared:  0.1247, Adjusted R-squared:  0.1052 
## F-statistic: 6.409 on 1 and 45 DF,  p-value: 0.01492

Výsledná rovnice je následující:

fertility <- 60.30438 + 0.19420*Agriculture

S faktorem religion

fit2 <- lm(Fertility ~ Agriculture + CatholicBin, data)
summary(fit2)
## 
## Call:
## lm(formula = Fertility ~ Agriculture + CatholicBin, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -29.803  -6.701   1.382   6.855  20.435 
## 
## Coefficients:
##                 Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      60.8322     4.1059  14.816   <2e-16 ***
## Agriculture       0.1242     0.0811   1.531   0.1329    
## CatholicBinTRUE   7.8843     3.7484   2.103   0.0412 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.39 on 44 degrees of freedom
## Multiple R-squared:  0.2046, Adjusted R-squared:  0.1685 
## F-statistic:  5.66 on 2 and 44 DF,  p-value: 0.006492

Výsledné rovnice jsou následující:

# protestanti
fertility <- 60.8322 + 0.1242*Agriculture + 7.8843*0
fertility <- 60.8322 + 0.1242*Agriculture

# katolíci
fertility <- 60.8322 + 0.1242*Agriculture + 7.8843*1
fertility <- 68.7165 + 0.1242*Agriculture

Faktor katolíci pouze jednorázově zvedl intercept.

Interakce: různé intercepts i slopes

fit3 <- lm(Fertility ~ Agriculture * CatholicBin, data)
summary(fit3)
## 
## Call:
## lm(formula = Fertility ~ Agriculture * CatholicBin, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -28.840  -6.668   1.016   7.092  20.242 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                 62.04993    4.78916  12.956   <2e-16 ***
## Agriculture                  0.09612    0.09881   0.973    0.336    
## CatholicBinTRUE              2.85770   10.62644   0.269    0.789    
## Agriculture:CatholicBinTRUE  0.08914    0.17611   0.506    0.615    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.49 on 43 degrees of freedom
## Multiple R-squared:  0.2094, Adjusted R-squared:  0.1542 
## F-statistic: 3.795 on 3 and 43 DF,  p-value: 0.01683

Tabulka koeficientů nám začíná houstnout, připomeňme ještě jednou obecné pravidlo:

Obecný postup sestavování rovnice: věci v řádcích se násobí, řádky mezi sebou sčítají. Za (Intercept) vždy dosazujeme 1. Dummy proměnné (poznáme podle slepeniny názvu původní proměnné a konkrétní hodnoty) = 1, pokud je to ono, jinak 0.

Výsledné rovnice jsou následující:

# protestanti
fertility <- 62.04993 + 0.09612*Agriculture + 2.85770*0 + 0.08914*Agriculture*0
fertility <- 62.04993 + 0.09612*Agriculture

# katolíci
fertility <- 62.04993 + 0.09612*Agriculture + 2.85770*1 + 0.08914*Agriculture*1
fertility <- 64.90763 + 0.18526*Agriculture

Faktor katolíci jednak zvedne intercept = 62.04993 + 2.85770, ale interakce změní také slope = 0.09612 + 0.08914.

Interakce nám umožňují modelovat různé slopes, ale kdybychom toto zapomněli, z odvozených rovnic na to vždy automaticky přijdeme.

Tím dostáváme stejnou situaci, jako kdybychom každou skupinu prokládali regresní přímkou samostatně.

ggplot(data, aes(x = Agriculture, y = Fertility, color = CatholicBin)) +
    geom_point() + geom_smooth(method = "lm")
## `geom_smooth()` using formula 'y ~ x'

Jak porovnávat modely?

fit1 <- lm(Fertility ~ Agriculture, data)
fit3 <- lm(Fertility ~ Agriculture + Examination + Education, data)
fit5 <- lm(Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality, data)
anova(fit1, fit3, fit5)
## Analysis of Variance Table
## 
## Model 1: Fertility ~ Agriculture
## Model 2: Fertility ~ Agriculture + Examination + Education
## Model 3: Fertility ~ Agriculture + Examination + Education + Catholic + 
##     Infant.Mortality
##   Res.Df    RSS Df Sum of Sq      F    Pr(>F)    
## 1     45 6283.1                                  
## 2     43 3180.9  2    3102.2 30.211 8.638e-09 ***
## 3     41 2105.0  2    1075.9 10.477 0.0002111 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Stepwise regression: Automatické nalezení vhodných prediktorů

Má své zastánce i odpůrce. Lepší je datům rozumět a prediktory a interakce vybírat podle smyslu. Jenže, neuteče nám tak něco? Dobré je porovnat oba přístupy.

Každý další prediktor zlepší R2. Zároveň tím ale narůstá složitost modelu a ne vždy je tak přidání další proměnné výhodou. Účelem je najít co nejmenší počet prediktorů, aby R2 bylo hodně dobré. Jako míra kvality může sloužit např. AIC kritérium, nebo raději BIC.

data <- read_excel("Housing.xlsx")  # hodnota nemovitostí podle mnoha kritérií
fit1 <- lm(Medv ~ ., data = data)  # základní model zahrnující všechny proměnné bez interakcí
summary(fit1) 
## 
## Call:
## lm(formula = Medv ~ ., data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.8948  -2.7585  -0.4663   1.7963  26.0911 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      36.461352   5.100994   7.148 3.21e-12 ***
## CensusTract      -0.002526   0.002080  -1.215 0.225046    
## Crime            -0.108762   0.032855  -3.310 0.001000 ** 
## Zoning            0.048031   0.013785   3.484 0.000538 ***
## Industrial        0.019932   0.061468   0.324 0.745871    
## CharlesRiverYes   2.705245   0.861298   3.141 0.001786 ** 
## NitricOxide     -17.541602   3.822390  -4.589 5.66e-06 ***
## Rooms             3.839225   0.418422   9.175  < 2e-16 ***
## HouseAge         -0.001938   0.013380  -0.145 0.884866    
## Distance         -1.493304   0.199892  -7.471 3.68e-13 ***
## Access            0.324925   0.068111   4.771 2.43e-06 ***
## TaxRate          -0.011598   0.003807  -3.046 0.002443 ** 
## PupilTeacher     -0.947985   0.130822  -7.246 1.67e-12 ***
## BBT               0.009357   0.002685   3.485 0.000536 ***
## LowerStatus      -0.526184   0.050704 -10.377  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.743 on 491 degrees of freedom
## Multiple R-squared:  0.7414, Adjusted R-squared:  0.734 
## F-statistic: 100.6 on 14 and 491 DF,  p-value: < 2.2e-16
fitStep <- step(fit1, scope = ~ .^2, k = log(nrow(data)), direction = "both", trace = 0) # kritérium BIC
# ^2 povolí přidat i interakce vždy 2 proměnných
summary(fitStep)
## 
## Call:
## lm(formula = Medv ~ Crime + Industrial + CharlesRiver + NitricOxide + 
##     Rooms + HouseAge + Distance + Access + TaxRate + PupilTeacher + 
##     BBT + LowerStatus + Rooms:LowerStatus + Access:LowerStatus + 
##     Rooms:Access + Distance:Access + BBT:LowerStatus + Distance:PupilTeacher + 
##     Crime:CharlesRiver + CharlesRiver:NitricOxide + CharlesRiver:Rooms + 
##     CharlesRiver:PupilTeacher + Rooms:PupilTeacher + HouseAge:BBT + 
##     Industrial:Distance + Industrial:LowerStatus + Crime:Rooms + 
##     Crime:LowerStatus, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -8.5845 -1.6797 -0.3157  1.5433 19.4311 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  -9.673e+01  1.350e+01  -7.167 2.93e-12 ***
## Crime                        -1.454e+00  3.147e-01  -4.620 4.95e-06 ***
## Industrial                    7.647e-01  1.237e-01   6.182 1.36e-09 ***
## CharlesRiverYes               6.341e+01  1.115e+01   5.687 2.26e-08 ***
## NitricOxide                  -1.691e+01  3.020e+00  -5.598 3.67e-08 ***
## Rooms                         1.946e+01  1.730e+00  11.250  < 2e-16 ***
## HouseAge                      2.233e-01  5.898e-02   3.786 0.000172 ***
## Distance                     -2.462e+00  6.776e-01  -3.634 0.000309 ***
## Access                        3.461e+00  3.109e-01  11.132  < 2e-16 ***
## TaxRate                      -1.401e-02  2.536e-03  -5.522 5.52e-08 ***
## PupilTeacher                  1.207e+00  7.085e-01   1.704 0.089111 .  
## BBT                           7.946e-02  1.262e-02   6.298 6.87e-10 ***
## LowerStatus                   2.939e+00  2.707e-01  10.857  < 2e-16 ***
## Rooms:LowerStatus            -3.793e-01  3.592e-02 -10.559  < 2e-16 ***
## Access:LowerStatus           -4.804e-02  4.465e-03 -10.760  < 2e-16 ***
## Rooms:Access                 -3.490e-01  4.370e-02  -7.986 1.05e-14 ***
## Distance:Access              -9.236e-02  2.603e-02  -3.548 0.000427 ***
## BBT:LowerStatus              -8.337e-04  3.355e-04  -2.485 0.013292 *  
## Distance:PupilTeacher         1.371e-01  3.719e-02   3.686 0.000254 ***
## Crime:CharlesRiverYes         2.544e+00  3.813e-01   6.672 7.01e-11 ***
## CharlesRiverYes:NitricOxide  -3.706e+01  6.202e+00  -5.976 4.48e-09 ***
## CharlesRiverYes:Rooms        -3.774e+00  7.402e-01  -5.099 4.94e-07 ***
## CharlesRiverYes:PupilTeacher -1.185e+00  3.701e-01  -3.203 0.001451 ** 
## Rooms:PupilTeacher           -3.792e-01  1.067e-01  -3.555 0.000415 ***
## HouseAge:BBT                 -7.107e-04  1.552e-04  -4.578 5.99e-06 ***
## Industrial:Distance          -1.316e-01  2.533e-02  -5.197 3.00e-07 ***
## Industrial:LowerStatus       -2.580e-02  5.204e-03  -4.959 9.88e-07 ***
## Crime:Rooms                   1.605e-01  4.001e-02   4.011 7.00e-05 ***
## Crime:LowerStatus             1.511e-02  4.954e-03   3.051 0.002408 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.045 on 477 degrees of freedom
## Multiple R-squared:  0.8964, Adjusted R-squared:  0.8904 
## F-statistic: 147.5 on 28 and 477 DF,  p-value: < 2.2e-16

Máme solidní model a můžeme tak dělat kvalifikovaný odhad cen nemovitostí. Můžeme stanovit běžnou cenu, ale také poznat nadsazené či výrazně výhodné nabídky :-)

Podmínky

  1. Linearita (Linearity) Residual plot nesmí vykazovat žádné vzory. Měly by to být náhodné body plus mínus kolem nuly. Řešením je provést vhodnou transformaci (např. log), nebo možná stačí jen zahrnout ještě nějaký další prediktor.
  2. Absence kolinearity (Absence of collinearity) Prediktory nesmějí být korelované. Jakmile jsou, je to vážný problém. Zvýhodňování jednoho znevýhodňuje druhý a model je ve slepé uličce. Řešením je vzdát se jednoho z páru korelovaných prediktorů i za cenu snížení R2. A nebo použití PCA, tím se ale drasticky sníží intuitivní interpretovatelnost modelu.
  3. Homoskedasticita (absence of heteroskedasticity) Rozptyl sledované proměnné musí být stejný napříč celým rozsahem hodnot prediktorů. Řešením může být vhodná transformace (např. log).
  4. Normalita residuálů (Normality of residuals) Nejméně důležitá podmínka, protože metoda nejmenších čtverců je poměrně robustní.
  5. Zbavení se vlivných odlehlých hodnot (Absence of influential outliers) Pokud se jedná o chybu měření, je jednoduché takové body identifikovat a odstranit. Horší je, pokud se o chybu nejedná. Je pak otázkou, zda má vůbec smysl používat regresní model. A nebo se rozhodnout, že model zkrátka postihne jen část populace a outliers odstranit. Regresní modely mají smysl pro popis běžných věcí, nahraditelných bodů „z davu“. Odlišné body mohou model pěkně rozhodit.
  6. Nezávislost (Independence) Extrémně důležitá podmínka, kterou mnozí omylem přehlížejí. Každý bod musí být produktem nezávislého výběru. Jakmile bude více bodů například od stejné osoby, výsledky se ovlivní a p-hodnoty jsou naprosto nevěrohodné. Řešením problému jsou mixed-effects modely :-)

A hlavně: data by měla vždy pocházet z náhodného výběru. Pokud toto porušíme, žádná statistická metoda nás nezachrání.