Liniowość modelu.•Składnik losowy: heteroskedastyczność, normalność, autokorelacja.•Nietypowe obserwacje: leverage, outliers, influential observations.•Wybór zmiennych na podstawie kryteriów.•Współliniowość (condition number)
Problemy
Generalnie podstawowe problemy z modelem regresji liniowej można podzielić na trzy szerokie grupy. Pierwsza klasa problemów to kwestiapoprawnej specyfikacji części strukturalnej modelu, o której zakładamy, że $Ey = X \beta$.
Druga klasa problemów jest związana z założeniami poczynionymi o składniku losowym. Ściśle, zakładamy, że $\epsilon \sim N(0, \sigma^2 I)$czyli, że błędy są niezależne o rozkładzie normalnym ze stałą wariancją (tzw. sferyczność składnika losowego).
Ostatnia grupa problemów dotyczy kwestii nietypowych obserwacji lub obserwacji, które nie pasują do wybranego modelu (outliers).
Zagadnieniem dodatkowym jest kwestia wyboru zmiennych, które powinny być włączone do modelu. Tutaj będziemy się zajmowali jedynieprocedurami opartymi o pewne kryteria "jakości modelu", np. AIC.
Liniowość modelu
Generalnym założeniem modelu regresji jest liniowość, tj. $Ey = X \beta$. Jeżeli założenie to nie jest spełnione, to uzyskane wynikisą niepoprawnie interpretowane, a uzyskane predykcje, zwłaszcza te poza poza próbą, będą najprawdopodobniej błędne. Dwa typowe testy służące do sprawdzania liniowości modeli to test Rainbow i tzw. test RESET.
Test Rainbow
Test Rainbow jest zaimplementowany w pakiecie lmtest jako funkcja raintest(). Idea tego testu polega na tym, żejeżeli rzeczywista zależność między zmienną wyjaśnianą i regresorami jest liniowa, to regresja przeprowadzona na podpróbie nie powinna być lepsza. Jeżeli model wyestymowany na całej próbie jest znacząco gorszy niż model wyestymowany na podpróbie, to należy odrzucić hipotezę o liniowości modelu.Podpróbę wybiera się na podstawie leverage, tj. wybiera się te obserwacje, które mają najnmniejszą leverage.
Rozważmy następujący sztuczny przykład.
> library(lmtest)> n <- 10^3> x <- 1:n> err <- rnorm(n, 0, 1)> y1 <- x+err # model liniowy> y2 <- x^2 + err # model nieliniowy> plot(x, y1)> plot(x, y2)> model1 <- lm(y1~x)> model2 <- lm(y2~x)> raintest(model1) Rainbow testdata: model1 Rain = 1.0038, df1 = 500, df2 = 498, p-value = 0.4833> raintest(model2) Rainbow testdata: model2 Rain = 30.8765, df1 = 500, df2 = 498, p-value < 2.2e-16>
Jak widać w pierwszym przypadku należy przyjąć hipotezę o liniowości modelu. W drugim przypadku należy odrzucić hipotezę o liniowości modelu.
Test RESET
Drugim popularnym testem stosowanym do sprawdzania liniowości modelu jest tzw. test RESET. Główna idea tego testu opiera się na szacowaniu większegomodelu rozszerzonego o (standardowo) drugie i trzecie potęgi regresorów. Natępnie jest stosowane statystyka F do sprawdzenia czydodane regresory są istotne.Test ten jest zaimplementowany w pakiecie lmtest jako funkcjaresettest().
Rozważmy następujący sztuczny przykład (jak powyżej).
> library(lmtest)> n <- 10^3> x <- 1:n> err <- rnorm(n, 0, 1)> y1 <- x+err # model liniowy> y2 <- x^2 + err # model nieliniowy> plot(x, y1)> plot(x, y2)> model1 <- lm(y1~x)> model2 <- lm(y2~x)> resettest(model1, power=2, type="regressor") RESET testdata: model1 RESET = 1.034, df1 = 1, df2 = 997, p-value = 0.3095> resettest(model2, power=2, type="regressor") RESET testdata: model2 RESET = 5.216049e+12, df1 = 1, df2 = 997, p-value < 2.2e-16>
Podobnie jak poprzednio, w pierwszym przypadku należy przyjąć hipotezę zerową o liniowości modelu, a w drugim ją odrzucić.
Składnik losowy
Kolejna grupa zagadnień, które wymagają sprawdzenia to założenia związane ze składnikiem losowym.
Heteroskedastyczność
Heteroskedastyczność występuje gdy wariancja składnika losowego nie jest stała. Podstawową rzeczą, którą należy wykonać jestdokładne przyjrzenie się wykresowi reszt vs. wartości dopasowane (pierwszy wykres w grupie standardowych wykresów związanych z regresją liniową).W sytuacji optymalnej wykres ten powinien przestawiać symetrycznie rozrzucone punkty na około poziomej lini. Jeżeli jest zachowana symetria, ale zwiększasię rozrzucenie punktów, to oznacza to występowania heteroskedastyczności. Jeżeli punktu są rozłożone zgodnie z jakimś trendem, to najprawdopodbniejwystępują problemy z liniowością modelu.
Typowym testem związanym z heteroskedastycznością jest test Goldfeld-Quandt'a. Idea tego testu polega na wyestymowaniu dwóch podmodeli,powstałych z podzielenia oryginalnego modelu i odrzuceniu hipotezy zerowej jeżeli wariancje w dwóch modelach różnią się istotnie. Test ten jestzaimplementowany w pakiecie lmtest jako funkcja gqtest(). Rozważmy następujący przykład.
> library(lmtest)> x <- rep(c(-1,1), 50)> err1 <- c(rnorm(50, sd=1), rnorm(50, sd=2)) # hetoroskedastycznosc> err2 <- rnorm(100) # hom*oskedastycznosc> y1 <- 1 + x + err1> y2 <- 1 + x + err2> gqtest(y1 ~ x) Goldfeld-Quandt testdata: y1 ~ x GQ = 4.4743, df1 = 48, df2 = 48, p-value = 3.72e-07> gqtest(y2 ~ x) Goldfeld-Quandt testdata: y2 ~ x GQ = 0.4918, df1 = 48, df2 = 48, p-value = 0.9923>
W pierwszym przypadku należy udrzucić hipotezę zerową o równości wariancji (a więc składnik losowy nie ma stałej wariancji). W drugim przypadku należy przyjąć hipotezę zerową o równości wariancji.
Inne typowe testy (wszystkie zaimplementowane w pakiecie lmtest) stosowane do wykrywania heteroskedastyczności to test Harrison-McCabe (funkcja hmctest()) oraz test Breusch-Pagan (funkcja bptest()).
Normalność
Graficzną pomocą w wykrywaniu odchyleń od założenia o normalności jest drugi wykres w zespole standardowych wykresów dla regresji liniowej.Jest to tzw. wykres Q-Q (quantile-quantile plot). Jest to wykres estymowanych kwantyli vs. kwantyli teoretycznych z rozkładu normalnego. Jeżeli resztymają rozkład normalny, to na wykresie tym powinniśmy obserować punkty położone niemal na przekątnej (tj. prostej $y=x$). W przypadku gdy obserwuje się odchylenia, oznacza to, że mamy problemy z tym założeniem. (Uwaga: oczywiście w ogonach rozkładów będziemy obserwowalitakie odchylenia, próba na której pracujemy jest skończona!).
Oczwiście można sprawdzić to założenie formalnie używając jednego z licznych testów. Jednym z typowych testów jest test Shapiro-Wilk. Test ten jest zaimplementowany w pakiecie stats jako funkcja shapiro.test(). Rozważmy następujący przykład.
> shapiro.test(rnorm(1000, mean = 5, sd = 3)) Shapiro-Wilk normality testdata: rnorm(1000, mean = 5, sd = 3) W = 0.9977, p-value = 0.1734> shapiro.test(runif(1000, min = 2, max = 4)) Shapiro-Wilk normality testdata: runif(1000, min = 2, max = 4) W = 0.9539, p-value < 2.2e-16>
W pierwszym przypadku nie ma podstaw do odrzucenia hipotezy o normalności. W drugim przypadku hipotezę o normalności należy odrzucić.
Warto zauważyć, że test Shapiro-Wilk jest stosunkowo mocny a dodatkowo jest mało wrażliwy na heteroskedastyczność i autokorelację składnika losowego.
Inny popularny test to test Jarque-Bera zaimplementowany w pakiecie tseries jako funkcja jarque.bera.test().Rozważmy następujący przykład.
> library(tseries)> x <- rnorm(100) # H0> jarque.bera.test(x) Jarque Bera Testdata: x X-squared = 0.39, df = 2, p-value = 0.8229> x <- runif(100) # Ha> jarque.bera.test(x) Jarque Bera Testdata: x X-squared = 6.9225, df = 2, p-value = 0.03139>
Autokorelacja
Ostatnim założeniem dotyczącym składnika losowego jest założenie o braku autokorelacji. Tutaj typowym testem jest test Durbin-Watsonzaimplementowany w pakiecie lmtest jako funkcja dwtest(). W teście tym hipotezą zerową jest zerowa wartość autokorelacji.Rozważmy następujący przykład.
> library(lmtest)> err1 <- rnorm(100) # bialy szum> err2 <- filter(err1, 0.7, method="recursive") # model AR(1) z parametrem 0.7> x <- rep(c(-1,1), 50)> y1 <- 1 + x + err1> y2 <- 1 + x + err2> dwtest(y1 ~ x) Durbin-Watson testdata: y1 ~ x DW = 2.0422, p-value = 0.6233alternative hypothesis: true autocorrelation is greater than 0 > dwtest(y2 ~ x) Durbin-Watson testdata: y2 ~ x DW = 0.7055, p-value = 6.016e-11alternative hypothesis: true autocorrelation is greater than 0 >
W pierwszym przykłądzie należy przyjąć hipotezę zerową, a w drugim ją odrzucić. Warto zaznaczyć, że test Durbin-Watson ma obszar nierozstrzygalności.
Innym typowym testem związanym ze zjawiskiem autokorelacji jest test Breusch-Godfrey zaimplementowany w pakiecie lmtestjako funkcja bgtest(). Rozważmy poprzedni przykład i dla zdefiniowanych tam modeli przeprowadźmy test Breusch-Godfrey.
> bgtest(y1~x) Breusch-Godfrey test for serial correlation of order 1data: y1 ~ x LM test = 0.0485, df = 1, p-value = 0.8256> bgtest(y1~x, order=2) Breusch-Godfrey test for serial correlation of order 2data: y1 ~ x LM test = 2.6845, df = 2, p-value = 0.2613> bgtest(y2~x) Breusch-Godfrey test for serial correlation of order 1data: y2 ~ x LM test = 41.8322, df = 1, p-value = 9.945e-11>
Ten test nie ma niedogodności związanych z obszarem niekonkuzywności. Ten test występuje czasami pod nazwą test mnożnika Lagrage'aautokorelacji składnika losowego.
Nietypowe obserwacje: leverage, outliers, influential observations
Dana obserwacja może nie pasować do modelu (outliers) lub znacznie zmieniać oszacowania modelu (influential observation). Istnienie takichobserwacji może wynikać z błędów w próbie (np. niepoprawnie wpisana dana) i wtedy należy to poprawić. Inna możliwość to konkretna przyczyna,dla której dana obserwacja jest nietypowa. Wtedy może to prowadzić do ciekawych odkryć lub do wyłączenia danej obserwacji. W każdym przypadkunależy się takim danym dokładniej przyjrzeć.
Przypomnijym, że $H = X (X^T X)^{-1} X^T$ (tzw. hat matrix) oraz, że wielkości $h_i = H_{ii}$ nazywają się leverages i mają bezpośredni wpływna wariancje składnika losowego $var \epsilon_i = \sigma^2 (1-h_i)$. Wiadomo, że $\sum_{i} \; h_i = p$, gdzie $p$ jest liczbą parametrów w modelu.Średnią wartością leverage jest $p/n$, więc przyjmuje się, że leverage większy od $2 p/n$ wskazuje na coś nietypowego (tj. duże wartości obserwacji).
Drugą standardową miarą czegoś nietypowego jest tzw. Cook distance zdefiniowany jako
$D_i = \frac{(\hat{y} - \hat{y}_{(i)})^T (\hat{y} - \hat{y}_{(i)})}{p \hat{\sigma}^2}, $
gdzie subskrypt $\hat{y}_{(i)}$ oznacza wartości teoretyczne wystymowane bez obserwacji $i$-tej. Jasne jest, że im większa wartość tak zdefiniowanejstatystyki, tym większy wpływ na modelu ma wyrzucona obserwacja.
Obie w/w charakterystyki są przedstawione w sposób graficzny na ostatnim wykresie w zespole standardowych wykresów związanych z regresją liniową. Wykres ten pozwala w prosty sposób zidentyfikować obserwacje, które są nietypowe. Dalsze postępowanie zależy od konkretnego zastosowania modelu i konkretnych danych, z którymi pracujemy.
Wybór zmiennych
W sytuacji gdy mamy pewien zbiór regresorów powstaje problem, które z nich należy włączyć do modelu, a które nie. Generlanie do tego problemu sądwa podejścia. Pierwsze, którego nie będziemy omawiali, polega na testowaniu zagnieżdżonych hipotez o istotności kolejnych zmiennych. Najprostsza tego typu metoda polega na wyestymowaniu modelu ze wszystkimi zmiennymi. Z tak wyestymowanego modelu należy usunąć najbardziej nieistotną zmienną, tj. zmienną, którajest nieistotna i ma największą wartość p-value. Następnie model należy powtórnie wyestymować. Ta procedura jest powtarzana tak długo, aż wszystkie regresory, które nadal są w modelu są istotne.
Inna metoda polaga na zdefiniowaniu pewnej miary jakości modelu i takiego doboru zmiennych, które zmaksymalizują (lub zminimalizują) taką miarę. Typowymwyborem takiej miary jest Akaike Information Criterion (AIC) lub Bayes Information Criterion (BIC).
Rozważmy następujący przykład, w którym na podstawie minimalizaji AIC szukamy najlepszego modelu. Wykonujemy to komendą step().
> data(state)> dane <- data.frame(state.x77, row.names=state.abb)> model <- lm(Life.Exp ~ ., data=dane)> step(model)Start: AIC= -22.18 Life.Exp ~ Population + Income + Illiteracy + Murder + HS.Grad + Frost + Area Df Sum of Sq RSS AIC- Area 1 0.001 23.298 -24.182- Income 1 0.004 23.302 -24.175- Illiteracy 1 0.005 23.302 -24.17423.297 -22.185- Population 1 1.747 25.044 -20.569- Frost 1 1.847 25.144 -20.371- HS.Grad 1 2.441 25.738 -19.202- Murder 1 23.141 46.438 10.305Step: AIC= -24.18 Life.Exp ~ Population + Income + Illiteracy + Murder + HS.Grad + Frost Df Sum of Sq RSS AIC- Illiteracy 1 0.004 23.302 -26.174- Income 1 0.006 23.304 -26.170 23.298 -24.182- Population 1 1.760 25.058 -22.541- Frost 1 2.049 25.347 -21.968- HS.Grad 1 2.980 26.279 -20.163- Murder 1 26.272 49.570 11.568Step: AIC= -26.17 Life.Exp ~ Population + Income + Murder + HS.Grad + Frost Df Sum of Sq RSS AIC- Income 1 0.006 23.308 -28.161 23.302 -26.174- Population 1 1.887 25.189 -24.280- Frost 1 3.037 26.339 -22.048- HS.Grad 1 3.495 26.797 -21.187- Murder 1 34.739 58.041 17.457Step: AIC= -28.16 Life.Exp ~ Population + Murder + HS.Grad + Frost Df Sum of Sq RSS AIC 23.308 -28.161- Population 1 2.064 25.372 -25.920- Frost 1 3.122 26.430 -23.876- HS.Grad 1 5.112 28.420 -20.246- Murder 1 34.816 58.124 15.528Call:lm(formula = Life.Exp ~ Population + Murder + HS.Grad + Frost, data = dane)Coefficients:(Intercept) Population Murder HS.Grad Frost 7.103e+01 5.014e-05 -3.001e-01 4.658e-02 -5.943e-03 >
Inne typowe kryteria to adjusted $R^2$ i statystyka Mallowa $C_p$.
Współliniowość
Dokładna współliniowość zachodzi wtedy gdy conajmniej jeden z regresorów jest liniową kombinacją pozostałych. Wtedy macierz $X^T X$ jest singularna, tj. nie istnieje jej odwrotność, a więc nie jesteśmy w stanie znaleźć estymatorów $\beta$. Tego typu problemy mogą powstać, gdy omyłkowo włączymy do modeluten sam regresor dwukrotnie, ten sam regresor, ale mierzony w przeskalowanych jednostakach (np. m, km), itp. Generalnie jednak dokładna współliniowośćjest zawzyczaj wynikiem błędu i jest wykrywana automatycznie. W praktyce znacznie częściej pojawia się tzw. przybliżona współliniowość, którapowoduje problemy z estymacją i interpretacją $\beta$.
Jedną z metod wykrywania współliniowości jest sprawdzenie wartości własnych macierzy $X^T X$. Niech $\lambda_{max}$ będzie największą wartością własną, a $\lambda_{min}$ najmniejszą wartością własną. Definujemy tzw. condition number $\kappa = \sqrt{ \lambda_{max}/\lambda_{min}}$. Przy dokładnej współliniowości$\lambda_{min}=0$ i $\kappa$ jest źle zdefiniowana, ale dla przybliżonej współliniowości $\kappa$ jest poprawnie zdefiniowana i im "mocniejsza" współliniowość, tym jest to liczba większa. Przyjmuje się, że wartości $\kappa>25$ są duże (w niektórych źródłach wartość progowa to $\kappa>30$).
Rozważmy następujący przykład.
> library(faraway)> data(seatpos)> model <- lm(hipcenter~., seatpos)> summary(model)Call:lm(formula = hipcenter ~ ., data = seatpos)Residuals: Min 1Q Median 3Q Max -73.827 -22.833 -3.678 25.017 62.337 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 436.43213 166.57162 2.620 0.0138 *Age 0.77572 0.57033 1.360 0.1843 Weight 0.02631 0.33097 0.080 0.9372 HtShoes -2.69241 9.75304 -0.276 0.7845 Ht 0.60134 10.12987 0.059 0.9531 Seated 0.53375 3.76189 0.142 0.8882 Arm -1.32807 3.90020 -0.341 0.7359 Thigh -1.14312 2.66002 -0.430 0.6706 Leg -6.43905 4.71386 -1.366 0.1824 ---Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Residual standard error: 37.72 on 29 degrees of freedomMultiple R-Squared: 0.6866, Adjusted R-squared: 0.6001 F-statistic: 7.94 on 8 and 29 DF, p-value: 1.306e-05 > X <- model.matrix(model)[,-1]> e <- eigen( t(X) %*% X)> e$val[1] 3.653671e+06 2.147948e+04 9.043225e+03 2.989526e+02 1.483948e+02 8.117397e+01 5.336194e+01 7.298209e+00> sqrt(e$val[1]/e$val)[1] 1.00000 13.04226 20.10032 110.55123 156.91171 212.15650 261.66698 707.54911>
Condition number $\kappa$ ma wartość ok. 707, a więc jest bardzo duży. Warto jednak zwrócić uwagę, że poza nim jest jeszcze kilka innychdużych wartości, co oznacza, że więcej zmiennych jest ze sobą współliniowych.
