Diagnostické nástroje a celá veda okolo prediktorov nám napovedajú o vhodnosti fitu a o tom, či náš model vhodne popisuje (sedí na) dáta. Čo urobíme, ak vizualizácia odhalí značné problémy.

Odpovedať na tieto otázky je ťažké. Teraz predpokladajme vhodnú štruktúru modelu, ale diagnostické nástroje stále napriek tomu vykazujú heteroskedastickú varianciu reziduálov. V takom prípade môžeme do modelu zakomponovať extra informáciu o variancii. \[var(\epsilon) = \sigma^2I\Rightarrow \sigma^2\Sigma\]

Aký má dopad heteroskedastická variancia reziduálov na signifikanciu (p-hodnotu) odhadu?

Pri nezvyčajných pozorovaniach je to už iné, niektoré vedia byť poriadne nebezpečné a odhad parametra môže byť výrazne vychýlený. Sú aj ukážkové prípady, kde ostránenie jedného bodu zmeni znamienka niekoľkých odhadov! Uvažujme situáciu, kde sa nevieme rozhodnúť, či bod je alebo nie je potrebné vyhodiť z dátovej vzorky. Niekedy je rozhodovanie veľmi ťažké a slúžia nám na to napríklad robustná regresia. Vyhadzovanie bodov z dátovej vzorky je veľmi nebezpečné, o čom nás už presvedčila minulosť, preto v tejto úlohe musíme byť vždy opatrný a preto aj robustnú regresiu budeme používať len v prípade potreby.

Čo môže zakryť využitie robustnej regresie?

Heteroskedastická variancia

Pre jednoduchosť použime dáta pipeline z knižnice faraway. Najprv urobte, čo je potrebné pred každým lineárnym modelom.

V našom jednoduchom prípade si môžeme vykresliť odozvu voči prediktoru pomocou ofarbenia skupín. Pozor: Prerušovaná čiara na obrázku nie je fit.

Na mieste sú otázky:

Teraz fitnite linearny regeresny model s Lab ako odozvou a Field ako prediktorom. Urobte diagnostiku modelu.

##              Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.967500   1.574787 -1.2494   0.2143
## Field        1.222968   0.041069 29.7781   <2e-16
## 
## n = 107, p = 2, Residual SE = 7.86476, R-Squared = 0.89

Chceli by sme použit metódu váženia reziduálov na vyriešenie problému s varianciou. Niekedy heteroskedasticitu vytvára jeden prediktor, ako je to v našom prípade. Ale nie vždy je to pravda, hlavne ak prediktorov je viacej. Teraz môžeme urobiť napríklad nasledovné.

Vyskúšame rozdeliť prediktor Field na 12 rovnakých skupín po 9 až na poslednú, ktorá bude obsahovať 8 pozorovaní. V každej skupine spočítame varianciu Lab (varlab) a strednú hodnotu (meanfield). Nakoniec si to nakreslíme, toto sa dá znova urobiť jedným ťahom za pomoci knižnice dplyr a ggplot2. Pokúste sa interpretovať kód:

pipeline %>% arrange(Field) %>%
  mutate(ff = gl(12,9)[-108]) %>% 
  group_by(ff) %>% summarise(meanfield = mean(Field),varlab = var(Lab)) %>%
  ggplot(aes(x=meanfield,y=varlab))+ geom_point() + theme_minimal()

Predpokladajme (alebo vidime z obrázka), že variancia našej odozvy k prediktoru (z obrázku) sa správa nasledovne: \[var(Lab) = a_0*Field^{a_1}.\] Fitnite regresný model log(varlab) na log(meanfield) na odhadnutie \(a_0\) a \(a_1\). Pre lepší fit bude vhodné posledný bod zamlčať. Nakoniec fit znázornite na obrázku

Lepšie vzťah neodmodelujeme, lebo odhadovanie exponencionálnej funkcie je veľmi zložitá optimalizačná úloha (na úlohu podobného charakteru existujú špecializované knižnice). Viac menej pre náš prípad to postačí. Použite tieto váhy v modeli a fitnite Lab ~ Field s váhami \(a_0*Field^{a_1}\) a vysvetlite summary.

## 
## Call:
## lm(formula = Lab ~ Field, data = pipeline, weights = my_exp(Field))
## 
## Weighted Residuals:
##      Min       1Q   Median       3Q      Max 
## -290.075  -33.881  -14.547    4.379  224.886 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.24617    2.94196   0.763    0.447    
## Field        1.13431    0.05689  19.938   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 76.21 on 105 degrees of freedom
## Multiple R-squared:  0.7911, Adjusted R-squared:  0.7891 
## F-statistic: 397.5 on 1 and 105 DF,  p-value: < 2.2e-16

Urobte ešte diagnostiku a sformulute záver, či sme si takýmto odmodelovaním pomohli alebo nie. Ak je nutné, pokúste sa navrhnuť ďalší postup.

Nezvyčajné pozorovania

Častokrát sa stane, že dáta obsahujú “veľa” nezvyčajných pozorovaní. A rovnako často je “veľmi” ťažké usúdiť, či vyhodiť pozorovanie/ia alebo nie. V prípade pákových bodov nemáme informáciu na okraji, vplyvné pozorovania môžu niečo naznačovať (napr. nastupujúci vývoj) a nakoniec outlier nesie extra informáciu.

Naš model obsahuje aj nezvyčajné pozorovanie/a, teraz fitnite rovnaký model pomocou funkcie rg() z knižnice quantreg, pozrite sa na výstup a zakrestlite tento nový fit do obrázku inou farbou. Poznámka: prístupov na fitnutie robustnej regresie je mnoho. pre jednoduchosť požijeme jeden. Doma si môžete vyskúšať aj ostatné.

## 
## Call: rq(formula = Lab ~ Field, data = pipeline)
## 
## tau: [1] 0.5
## 
## Coefficients:
##             coefficients lower bd upper bd
## (Intercept) -1.02059     -2.89777  0.47932
## Field        1.13824      1.09648  1.19711

Teraz zakreslite do dát aj fit lineárneho modelu, kde sme nastavili váhy na varianciu.

Samostatná práca

Vezmite dáta stackloss a urobte, čo je treba na začiatok. Fitnite plný regresný model s stack.loss ako odozvou.

  1. Sústreďte sa na varianciu reziduálov.
  2. Overte nezvyčajné pozorovania.
  3. Navrhnite prípadné kroky a realizujte.
  4. Svoje zmeny potvrďte diagnostikou modelu.