Vícenásobná regrese a mixed-effects modely
Data ke stažení
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-16Protož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.5Interpretace? 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 estimatesplot(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 paneluDů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-10Pří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.01492Jak 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-16summary(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-16Interakce
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.01492Vý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.006492Vý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*AgricultureFaktor 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.01683Tabulka 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*AgricultureFaktor 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 ' ' 1Stepwise 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-16fitStep <- 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-16Má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
- 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.
- 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.
- 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).
- Normalita residuálů (Normality of residuals) Nejméně důležitá podmínka, protože metoda nejmenších čtverců je poměrně robustní.
- 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.
- 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í.
Model splňuje všechny podmínky, zkontrolovali jsme residuals. Skutečně?
data <- read_excel("orly_owl_Lin_4p_5_flat.xlsx")
fit <- lm(V1 ~ ., data)
summary(fit)##
## Call:
## lm(formula = V1 ~ ., data = data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0652 -0.9878 0.1419 0.7958 1.9560
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.001194 0.021170 0.056 0.955
## V2 0.986328 0.128630 7.668 2.56e-14 ***
## V3 0.972295 0.127506 7.625 3.54e-14 ***
## V4 0.861429 0.120431 7.153 1.14e-12 ***
## V5 0.927474 0.084432 10.985 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.001 on 2293 degrees of freedom
## Multiple R-squared: 0.05, Adjusted R-squared: 0.04834
## F-statistic: 30.17 on 4 and 2293 DF, p-value: < 2.2e-16Máme koeficienty a p-hodnoty, stačí nám to?
plot(predict(fit), resid(fit))
Mixed-effects modely
Mixed-effects = Fixed and random effects
Zajímá nás vztah frequency ~ politeness + sex + ε, kde politeness je faktor se dvěma hodnotami (pol = formal register, inf = informal register), jehož vliv je hlavní otázkou studie. Do modelu je ještě zahrnut fixní efekt sex. Fixed-effect je to kvůli tomu, že to máme pod kontrolou. Můžeme to ovlivnit a v experimentu máme kontrolovaně zastoupeny všechny možné hodnoty této proměnné, není to tedy malý náhodný výběr z obrovské populace. Navíc se dá očekávat, že efekt pohlaví je systematický a vlastně i predikovatelný. Člen ε je „error term“ zachycující zbylé „náhodné“, nebo spíše nevysvětlené odchylky.
Komplikace: více měření od každé osoby (subject) (porušení podmínky nezávislosti pro lineární regresní modely). Vyřešíme to hladce, subject prohlásíme za random-effect (něco, čehož vlastnosti nemáme pod kontrolou, něco, co jsme náhodně vybrali z populace). Toto umožní modelu uvažovat rozdílnou frequency pro každou osobu zvlášť. Např. subject 1 může mít průměr 232 Hz napříč promluvami, zatímco subject 2 třeba 258 Hz.
library(readxl)
data <- read_excel("politeness_data.xlsx")
which(!complete.cases(data)) # v 39. řádku sice chybí hodnota, ale následující metody si s tím umí poradit## [1] 39library(tidyverse)
data %>% group_by(subject) %>% summarise(mean(frequency, na.rm = TRUE))## # A tibble: 6 x 2
## subject `mean(frequency, na.rm = TRUE)`
## <chr> <dbl>
## 1 F1 232.
## 2 F2 258.
## 3 F3 251.
## 4 M3 169.
## 5 M4 146.
## 6 M7 102.boxplot(frequency ~ subject, data)
Je sice patrné, že muži jsou níže, ale zřejmá je také vysoká individuální variabilita. Tu můžeme modelovat různými random intercepts. Tím tedy ubíráme něco z ε, něco, co je sice náhodné, ale v rámci konkrétní osoby již charakteristické. Nový vzorce vypadá frequency ~ politeness + sex + (1|subject) + ε. Jednička znamená, že se bavíme o intercept, a ten má být různý pro každý subject.
V designu studie je ještě jeden problém – různé položky, scénáře. Například požádejte profesora o laskavost (formální situace) vs. požádejte kamaráda o laskavost (neformální situace). Nebo omluvte se, že jdete pozdě – jiný scénář. Očekáváme tedy variaci v rámci těchto scénářů.
boxplot(frequency ~ scenario, data)
To tedy také přidáme do předpisu modelu frequency ~ politeness + sex + (1|subject) + (1|scenario) + ε. Takový model je mnohem efektivnější než dříve používané průměrování. Průměrování přes více proměnných je problém. Tady vidí model vše, jak skutečně je, a má i představu o variabilitě.
Vztah mezi frequency a politness v závislosti na pohlaví.
boxplot(frequency ~ politeness*sex, col=c("white", "lightgray"), data)
library(lme4)
# Toto by nešlo: lmer(frequency ~ politeness, data), jelikož není vložen žádný "random-effect"
politeness.model <- lmer(frequency ~ politeness + (1|subject) + (1|scenario), data)
summary(politeness.model)## Linear mixed model fit by REML ['lmerMod']
## Formula: frequency ~ politeness + (1 | subject) + (1 | scenario)
## Data: data
##
## REML criterion at convergence: 793.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2006 -0.5817 -0.0639 0.5625 3.4385
##
## Random effects:
## Groups Name Variance Std.Dev.
## scenario (Intercept) 219 14.80
## subject (Intercept) 4015 63.36
## Residual 646 25.42
## Number of obs: 83, groups: scenario, 7; subject, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 202.588 26.754 7.572
## politenesspol -19.695 5.585 -3.527
##
## Correlation of Fixed Effects:
## (Intr)
## politensspl -0.103coef(politeness.model)## $scenario
## (Intercept) politenesspol
## 1 189.0772 -19.69454
## 2 209.1573 -19.69454
## 3 213.9799 -19.69454
## 4 223.1972 -19.69454
## 5 200.6423 -19.69454
## 6 190.4838 -19.69454
## 7 191.5789 -19.69454
##
## $subject
## (Intercept) politenesspol
## F1 241.4365 -19.69454
## F2 267.2893 -19.69454
## F3 259.9240 -19.69454
## M3 179.0959 -19.69454
## M4 154.7280 -19.69454
## M7 113.0550 -19.69454
##
## attr(,"class")
## [1] "coef.mer"Random effects Zajímavý sloupec Variance: míra variability (rozptylu) zachycené danými random-effects. Residual je naše ε.
Fixed effects To už známe z klasických lineárních regresních modelů. (Intercept) přísluší první variantě, tedy inf. Varianta pol ho sníží o 19.695. Co znamená hodnota 202.588? Je to průměr mužů a žen dohromady při inf. To není moc chytré a chtělo by to modelu přidat informaci o pohlaví (farma s tuctem slepic a tuctem krav má v průměru tři nohy, což příliš neinformuje o tom, co se skutečně na farmě děje).
politeness.model <- lmer(frequency ~ politeness + sex + (1|subject) + (1|scenario), data)
summary(politeness.model)## Linear mixed model fit by REML ['lmerMod']
## Formula: frequency ~ politeness + sex + (1 | subject) + (1 | scenario)
## Data: data
##
## REML criterion at convergence: 775.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2591 -0.6236 -0.0772 0.5388 3.4795
##
## Random effects:
## Groups Name Variance Std.Dev.
## scenario (Intercept) 219.5 14.81
## subject (Intercept) 615.6 24.81
## Residual 645.9 25.41
## Number of obs: 83, groups: scenario, 7; subject, 6
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 256.846 16.116 15.938
## politenesspol -19.721 5.584 -3.532
## sexM -108.516 21.013 -5.164
##
## Correlation of Fixed Effects:
## (Intr) pltnss
## politensspl -0.173
## sexM -0.652 0.004coef(politeness.model)## $scenario
## (Intercept) politenesspol sexM
## 1 243.3398 -19.72111 -108.5163
## 2 263.4292 -19.72111 -108.5163
## 3 268.2541 -19.72111 -108.5163
## 4 277.4757 -19.72111 -108.5163
## 5 254.9102 -19.72111 -108.5163
## 6 244.6724 -19.72111 -108.5163
## 7 245.8426 -19.72111 -108.5163
##
## $subject
## (Intercept) politenesspol sexM
## F1 242.9386 -19.72111 -108.5163
## F2 267.2654 -19.72111 -108.5163
## F3 260.3348 -19.72111 -108.5163
## M3 285.2283 -19.72111 -108.5163
## M4 262.2248 -19.72111 -108.5163
## M7 223.0858 -19.72111 -108.5163
##
## attr(,"class")
## [1] "coef.mer"Random effects Velice se snížil rozptyl (Intercept) subject, což je dobré, protože to má vzhledem k odděleným pohlaví logiku. Jednoduché rozdělení do dvou skupin značně snížilo variabilitu. Nicméně rozptyl odchylek zůstává u obou modelů prakticky stejný. Jen se rozptyl subject (vysvětlený díky modelování jejich jedinečných interceptů random faktorem) nyní přesunul do vysvětlení fixním faktorem sex. Takže pro model to není nic zásadního, ale pro nás, pokud chceme na první pohled vidět rozdíl mezi pohlavími a neluštit to z jednotlivých koeficientů random effects, to přináší zajímavou a přehlednou informaci.
Fixed effects (Intercept) znamená evidentně F + inf. Muži M jsou o 108.5 Hz níže a situace pol má podobný efekt, jako předtím, opět snížení o cca 20 Hz.
Signifikance výsledků
Bohužel zde není výpočet p-hodnot tak přímočarý, jako byl u lineárních modelů. Je několik různých přístupů, o kterých se bouřlivě diskutuje a zatím neexistuje odpověď, co je ideální řešení. Připravme se tedy na to, že ať zvolíme cokoliv, vždy se to někomu nebude líbit. Zkusíme Likelihood Ratio Test. Porovnáme model se všemi faktory a redukovaný model bez faktoru, jehož vliv chceme vyhodnotit.
politeness.model <- lmer(frequency ~ politeness + sex + (1|subject) + (1|scenario), data, REML=FALSE)
politeness.0politeness <- lmer(frequency ~ sex + (1|subject) + (1|scenario), data, REML=FALSE)
# REML=FALSE je potřeba uvádět, když chceme používat Likelihood Ratio Test
anova(politeness.model, politeness.0politeness)## Data: data
## Models:
## politeness.0politeness: frequency ~ sex + (1 | subject) + (1 | scenario)
## politeness.model: frequency ~ politeness + sex + (1 | subject) + (1 | scenario)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## politeness.0politeness 5 816.72 828.81 -403.36 806.72
## politeness.model 6 807.10 821.61 -397.55 795.10 11.618 1 0.0006532
##
## politeness.0politeness
## politeness.model ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1A už můžeme psát “politeness affected frequency (χ2(1)=11.62, p=0.00065), lowering it by about 19.7 Hz ± 5.6 (standard errors)”
A jak by to bylo s interakcemi?
Pokud máme podezření, že mluvení zdvořilé vs. neformální má opačný (nebo jinak velký) efekt u mužů než u žen, je dobré zkusit model s interakcí. Pro test signifikance by se pak porovnávaly modely
- full model: frequency ~ politeness * sex
- reduced model: frequency ~ politeness + sex
Pokud vyjde signifikantní rozdíl, interakce tam je.
Random slopes modely
Toto je skoro to nejdůležitější. Když se podíváme na koeficienty modelu, tak každý scénář a každý subjekt má jiný Intercept. Proto se nazývá random intercept model. Problém je, že předpokládáme, že efekt zdvořilosti je ve všech případech stejný, což nemusí být pravda. Různí lidé budou určitě různě zdvořilí, různé situace asi také povedou na různé rozdíly.
coef(politeness.model)## $scenario
## (Intercept) politenesspol sexM
## 1 243.4860 -19.72207 -108.5173
## 2 263.3592 -19.72207 -108.5173
## 3 268.1322 -19.72207 -108.5173
## 4 277.2546 -19.72207 -108.5173
## 5 254.9319 -19.72207 -108.5173
## 6 244.8015 -19.72207 -108.5173
## 7 245.9618 -19.72207 -108.5173
##
## $subject
## (Intercept) politenesspol sexM
## F1 243.3684 -19.72207 -108.5173
## F2 266.9442 -19.72207 -108.5173
## F3 260.2276 -19.72207 -108.5173
## M3 284.3535 -19.72207 -108.5173
## M4 262.0575 -19.72207 -108.5173
## M7 224.1293 -19.72207 -108.5173
##
## attr(,"class")
## [1] "coef.mer"Pokud je dost dat, je dobré použít random slope model.
politeness.model <- lmer(frequency ~ politeness + sex +
(1+politeness|subject) + (1+politeness|scenario), data, REML=FALSE)## boundary (singular) fit: see help('isSingular')Povolujeme nejen různý intercept napříč osobami a situacemi, ale také různý vliv politeness, což je hlavní proměnná, jejíž vliv zkoumáme.
Toto je sice hezký nápad, ale model nám zahlásí varování
boundary (singular) fit: see ?isSingular, což znamená
problém s konvergencí, modelu tedy nelze věřit. Důvodem bude malé
množství dat na to, abychom si mohli dovolit modelovat různý slope
politeness pro každou osobu. Chtělo by to od každé osoby více nahrávek,
na což by bývalo bylo dobré myslet už od začátku při designu
experimentu.
Model tedy zjednodušíme a budeme modelovat jen různé slope politeness pro jednotlivé scenario, kde už máme pro každou situaci dat dostatečně.
politeness.model <- lmer(frequency ~ politeness + sex +
(1|subject) + (1+politeness|scenario), data, REML=FALSE)coef(politeness.model)## $scenario
## (Intercept) politenesspol sexM
## 1 244.0806 -20.35867 -108.5499
## 2 262.1605 -15.97479 -108.5499
## 3 268.0862 -20.72032 -108.5499
## 4 275.7626 -16.42730 -108.5499
## 5 254.9207 -19.39768 -108.5499
## 6 245.6879 -21.89365 -108.5499
## 7 247.3428 -23.51028 -108.5499
##
## $subject
## (Intercept) politenesspol sexM
## F1 243.3651 -19.75467 -108.5499
## F2 266.9752 -19.75467 -108.5499
## F3 260.2488 -19.75467 -108.5499
## M3 284.4393 -19.75467 -108.5499
## M4 262.0225 -19.75467 -108.5499
## M7 224.1273 -19.75467 -108.5499
##
## attr(,"class")
## [1] "coef.mer"Takže rozdíly sice jsou, ale ne zase až tak velké. A příjemné pro náš výzkum je, že všechny změny mají stejné znaménko.
Pojďme otestovat signifikanci.
politeness.0politeness <- lmer(frequency ~ sex +
(1|subject) + (1+politeness|scenario), data, REML=FALSE)
# vyhodili jsme akorát politeness z fixních. ale struktura random effects zůstává!
anova(politeness.model, politeness.0politeness)## Data: data
## Models:
## politeness.0politeness: frequency ~ sex + (1 | subject) + (1 + politeness | scenario)
## politeness.model: frequency ~ politeness + sex + (1 | subject) + (1 + politeness | scenario)
## npar AIC BIC logLik deviance Chisq Df Pr(>Chisq)
## politeness.0politeness 7 815.67 832.60 -400.84 801.67
## politeness.model 8 810.93 830.28 -397.47 794.93 6.739 1 0.009433
##
## politeness.0politeness
## politeness.model **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Výsledek je signifikantní.
Přestože mnoho lidí používá pouze intercept-only models, uvažování random slopes je velice důležité. Lidé přeci určitě reagují odlišně na dané situace. A lze očekávat, že v různých položkách to také dopadně rozdílně.
Simulacemi se ukazuje, že intercept-only models jsou antikonzervativní, mají vyšší riziko chyby 1. druhu, takže z nich lezou častěji nesprávně signifikantní výsledky. Je tedy lákavé pokušení najít „signifikantní objev,“ nicméně s větší pravděpodobností to bude jen omyl.
Vkládat by se tedy měly random-slopes, které dávají vzhledem ke studii smysl, a to pro všechny fixed effects, které jsou důležité pro celkovou interpretaci výsledků studie.
Ve výše uvedené studii šlo jen o politeness, autory nezajímaly rozdíly mezi pohlavími, přestože se vyplatilo mít je pod kontrolou. Proto použili random slopes jen pro efekt politeness, nikoliv pro sex.
Podmínky
Dobrá zpráva je, že zůstávají stejné, jako u lineárních regresních modelů. Stačí tedy provádět stejné testy.
Nezávislost není potřeba, pokud ji máme pod kontrolou. To je výhoda mixed-effects modelů. Problém by nastal, kdybychom nevložili důležité fixed nebo random efekty. Např. bez efektu subject by model nevěděl, že je více položek od stejné osoby. Pak je to porušení podmínek a model nefunguje správně.
Poznámka: funkce dfbeta pro tyto modely nefunguje, v “bw_LME_tutorial2.pdf” na str. 18 a 19 je alternativní funkce.
A jak to citovat?
citation()##
## To cite R in publications use:
##
## R Core Team (2022). R: A language and environment for statistical
## computing. R Foundation for Statistical Computing, Vienna, Austria.
## URL https://www.R-project.org/.
##
## A BibTeX entry for LaTeX users is
##
## @Manual{,
## title = {R: A Language and Environment for Statistical Computing},
## author = {{R Core Team}},
## organization = {R Foundation for Statistical Computing},
## address = {Vienna, Austria},
## year = {2022},
## url = {https://www.R-project.org/},
## }
##
## We have invested a lot of time and effort in creating R, please cite it
## when using it for data analysis. See also 'citation("pkgname")' for
## citing R packages.citation("lme4")##
## To cite lme4 in publications use:
##
## Douglas Bates, Martin Maechler, Ben Bolker, Steve Walker (2015).
## Fitting Linear Mixed-Effects Models Using lme4. Journal of
## Statistical Software, 67(1), 1-48. doi:10.18637/jss.v067.i01.
##
## A BibTeX entry for LaTeX users is
##
## @Article{,
## title = {Fitting Linear Mixed-Effects Models Using {lme4}},
## author = {Douglas Bates and Martin M{\"a}chler and Ben Bolker and Steve Walker},
## journal = {Journal of Statistical Software},
## year = {2015},
## volume = {67},
## number = {1},
## pages = {1--48},
## doi = {10.18637/jss.v067.i01},
## }“We used R (R Core Team, 2020) and lme4 (Bates, Maechler & Bolker, 2015) to perform a linear mixed effects analysis of the relationship between frequency and politeness. As fixed effects, we entered politeness and gender (without interaction term) into the model. As random effects, we had intercepts for subjects and items, as well as by-subject and by-item random slopes for the effect of politeness. Visual inspection of residual plots did not reveal any obvious deviations from homoscedasticity or normality. P-values were obtained by likelihood ratio tests of the full model with the effect in question against the model without the effect in question.”
Post-hoc testy
Pokud má fixed-effect tři a více možných hodnot, může nás kromě celkové statistické významnosti zajímat, které skupiny se konkrétně významně liší od jiných. K tomu slouží tzv. post-hoc testy, obdoba párových testů, které ale zároveň berou v potaz celkový model (nehrozí tedy problém zvýšené chyby 1. druhu důsledkem tzv. multiple-testing).
K tomu slouží např. balíček emmeans, který umí pracovat s modely lm() i lmer(). Na této stránce najdeme sekci Vignettes s návody na konkrétních příkladech, pro smíšené modely konkrétně https://cran.r-project.org/web/packages/emmeans/vignettes/sophisticated.html
Závěrem
Když už v textu napíšeme, že jsme zkontrolovali residuály, aby odpovídaly podmínkám modelu, tak je to dobré také skutečně udělat :-) Ať už se jedná o jednoduchý lineární model s jednou nezávislou proměnnou nebo o obrovskou vícenásobnou regresi s interakcemi či smíšený model včetně důležitého random-slopes rozšíření, vizuální kontrolu residuálů provedeme jednoduše takto:
plot(predict(fit), resid(fit))© 10. 10. 2018 a 1. 12. 2020 Tomáš Bořil, borilt@gmail.com