Jak interpretować ten wykres dopasowania i resztek?

Naprawdę nie rozumiem heteroscedastyczności. Chciałbym wiedzieć, czy mój model jest odpowiedni, czy nie zgodnie z tym wątkiem.

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$y$$y \sim N(X\beta,\sigma^2)$$y$$X\beta$$\sigma^2$$y = X\beta + \epsilon$$\epsilon \sim N(0,\sigma^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:

$R^2$$R^2$$0.5989$$0.03919$

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

$X^2$$X^3$

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.