Naprawdę nie rozumiem heteroscedastyczności. Chciałbym wiedzieć, czy mój model jest odpowiedni, czy nie zgodnie z tym wątkiem.
Naprawdę nie rozumiem heteroscedastyczności. Chciałbym wiedzieć, czy mój model jest odpowiedni, czy nie zgodnie z tym wątkiem.
Odpowiedzi:
Jak skomentował @IrishStat, musisz sprawdzić zaobserwowane wartości pod kątem błędów, aby sprawdzić, czy występują problemy ze zmiennością. Wrócę do tego pod koniec.
Właśnie dlatego masz pojęcie o tym, co rozumiemy przez heteroskedastyczność: kiedy dopasowujesz model liniowy do zmiennej , zasadniczo mówisz, że przyjmujesz założenie, że twoje y ∼ N ( X β , σ 2 ) lub w kategoriach laika, że twoje y ma się równać X β plus niektóre błędy, które mają wariancję σ 2 . To praktycznie twój model liniowy y = X β + ϵ , gdzie błędy ϵ ∼ N ( 0 , σ 2 ). OK, fajne jak dotąd zobaczymy to w kodzie:
set.seed(1); #set the seed for reproducability
N = 100; #Sample size
x = runif(N) #Independant variable
beta = 4; #Regression coefficient
epsilon = rnorm(N); #Error with variance 1 and mean 0
y = x * beta + epsilon #Your generative model
lin_mod <- lm(y ~x) #Your linear model
więc dobrze, jak zachowuje się mój model:
x11(); par(mfrow=c(1,3)); #Make a new 1-by-3 plot
plot(residuals(lin_mod));
title("Simple Residual Plot - OK model")
acf(residuals(lin_mod), main = "");
title("Residual Autocorrelation Plot - OK model");
plot(fitted(lin_mod), residuals(lin_mod));
title("Residual vs Fit. value - OK model");
co powinno dać ci coś takiego: co oznacza, że twoje reszty nie wydają się mieć oczywistego trendu opartego na twoim arbitralnym indeksie (pierwszy wykres - naprawdę najmniej informacyjny), wydają się nie mieć prawdziwej korelacji między nimi (drugi wykres - dość ważny i prawdopodobnie ważniejsze niż homoskedastyczność) i że dopasowane wartości nie mają oczywistej tendencji do niepowodzenia, tj. twoje dopasowane wartości względem twoich reszt wydają się dość losowe. Na tej podstawie powiedzielibyśmy, że nie mamy problemów z heteroskedastycznością, ponieważ nasze pozostałości wydają się mieć wszędzie taką samą wariancję.
OK, chcesz jednak heteroskedastyczności. Biorąc pod uwagę te same założenia liniowości i addytywności, zdefiniujmy inny model generatywny z „oczywistymi” problemami heteroskedastyczności. Mianowicie po niektórych wartościach nasza obserwacja będzie znacznie głośniejsza.
epsilon_HS = epsilon;
epsilon_HS[ x>.55 ] = epsilon_HS[x>.55 ] * 9 #Heteroskedastic errors
y2 = x * beta + epsilon_HS #Your generative model
lin_mod2 <- lm(y2 ~x) #Your unfortunate LM
gdzie proste wykresy diagnostyczne modelu:
par(mfrow=c(1,3)); #Make a new 1-by-3 plot
plot(residuals(lin_mod2));
title("Simple Residual Plot - Fishy model")
acf(residuals(lin_mod2), main = "");
title("Residual Autocorrelation Plot - Fishy model");
plot(fitted(lin_mod2), residuals(lin_mod2));
title("Residual vs Fit. value - Fishy model");
powinien dać coś takiego: Tutaj pierwszy wątek wydaje się nieco „dziwny”; wygląda na to, że mamy kilka reszt, które skupiają się w małych wielkościach, ale nie zawsze jest to problem ... Drugi wykres jest w porządku, oznacza to, że nie mamy korelacji między twoimi resztami w różnych opóźnieniach, więc możemy oddychać przez chwilę. Trzeci spisek rozlewa fasolę: nie wiadomo, czy gdy osiągnęliśmy wyższe wartości, nasze pozostałości eksplodują. Zdecydowanie mamy heteroskedastyczność w pozostałościach tego modelu i musimy coś z tym zrobić (np. IRLS , regresja Theil – Sen itp.)
Tutaj problem był naprawdę oczywisty, ale w innych przypadkach moglibyśmy przeoczyć; Aby zmniejszyć nasze szanse na przeoczenie tego wydarzenia, innym wnikliwym spiskiem był wspomniany przez IrishStat: Wartości rezydualne a wartości zaobserwowane lub dla naszego problemu z zabawkami:
par(mfrow=c(1,2))
plot(y, residuals(lin_mod) );
title( "Residual vs Obs. value - OK model")
plot(y2, residuals(lin_mod2) );
title( "Residual vs Obs. value - Fishy model")
co powinno dać coś takiego:
W uczciwości twojej sytuacji wykres wartości resztkowych w porównaniu z dopasowanymi wartościami wydaje się względnie OK. Sprawdzanie wartości resztkowych w stosunku do zaobserwowanych wartości prawdopodobnie byłoby pomocne, aby upewnić się, że jesteś po bezpiecznej stronie. (Nie wspominałem o wykresach QQ ani nic podobnego, aby nie wprawiać w zakłopotanie bardziej, ale możesz też to krótko sprawdzić). Mam nadzieję, że pomoże to w zrozumieniu heteroskedastyczności i tego, na co powinieneś zwrócić uwagę.
Twoje pytanie wydaje się dotyczyć heteroscedastyczności (ponieważ wspomniałeś je po nazwie i dodałeś tag), ale twoje wyraźne pytanie (np. W tytule i) kończące post jest bardziej ogólne: „czy mój model jest odpowiedni, czy nie zgodnie z tym wątek". Ustalenie, czy model jest nieodpowiedni, polega nie tylko na ocenie heteroscedastyczności.
Zeskrobałem twoje dane za pomocą tej strony (ht @Alexis). Pamiętaj, że dane są sortowane w porządku rosnącym od fitted
. W oparciu o regresję i lewy górny wykres wydaje się wystarczająco wierny:
mod = lm(residuals~fitted)
summary(mod)
# ...
# Residuals:
# Min 1Q Median 3Q Max
# -0.78374 -0.13559 0.00928 0.19525 0.48107
#
# Coefficients:
# Estimate Std. Error t value Pr(>|t|)
# (Intercept) 0.06406 0.35123 0.182 0.856
# fitted -0.01178 0.05675 -0.208 0.836
#
# Residual standard error: 0.2349 on 53 degrees of freedom
# Multiple R-squared: 0.0008118, Adjusted R-squared: -0.01804
# F-statistic: 0.04306 on 1 and 53 DF, p-value: 0.8364
Nie widzę tu żadnych dowodów na heteroscedastyczność. W prawym górnym rogu (wykres qq) wydaje się, że nie ma również żadnych problemów z założeniem normalności.
Z drugiej strony, krzywa „S” w czerwonym dopasowaniu lowess (w lewym górnym wykresie) oraz wykresy acf i pacf (u dołu) wydają się problematyczne. Po lewej stronie większość reszt znajduje się powyżej szarej linii 0. Gdy przesuniesz się w prawo, większość reszt spadnie poniżej 0, następnie powyżej, a następnie ponownie poniżej. Powoduje to, że gdybym ci powiedział, że patrzę na konkretną resztkę i że ma ona wartość ujemną (ale nie powiedziałem ci, na którą patrzę), możesz zgadnąć z dobrą dokładnością, że resztki w pobliżu zostały również negatywnie ocenione. Innymi słowy, reszty nie są niezależne - wiedza o jednym daje informacje o innych.
Oprócz wykresów można to przetestować. Prostym podejściem jest użycie testu przebiegu :
library(randtests)
runs.test(residuals)
# Runs Test
#
# data: residuals
# statistic = -3.2972, runs = 16, n1 = 27, n2 = 27, n = 54, p-value = 0.0009764
# alternative hypothesis: nonrandomness
Aby odpowiedzieć na Twoje wyraźne pytania: Twój wykres pokazuje szereg autokorelacji / nie-niezależności twoich reszt. Oznacza to, że Twój model nie jest odpowiedni w obecnej formie.