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).
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.
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
„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
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)
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
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.
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'
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
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 :-)
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í.