Czy autokorelowane wzorce resztkowe pozostają nawet w modelach z odpowiednimi strukturami korelacji i jak wybrać najlepsze modele?

Kontekst

To pytanie używa R, ale dotyczy ogólnych problemów statystycznych.

Analizuję wpływ czynników umieralności (% umieralności z powodu chorób i pasożytnictwa) na tempo wzrostu populacji ćmy w czasie, gdy populacje larw pobierano z 12 miejsc raz w roku przez 8 lat. Dane dotyczące tempa wzrostu populacji pokazują wyraźny, ale nieregularny trend cykliczny w czasie.

Resztki z prostego uogólnionego modelu liniowego (tempo wzrostu ~% choroby +% pasożytnictwo + rok) wykazywały podobnie wyraźny, ale nieregularny trend cykliczny w czasie. Dlatego uogólnione modele najmniejszych kwadratów tej samej postaci zostały również dopasowane do danych z odpowiednimi strukturami korelacji, aby poradzić sobie z czasową autokorelacją, np. Złożoną symetrią, autoregresyjnym porządkiem procesu 1 i autoregresyjnymi strukturami średniej ruchomej.

Wszystkie modele zawierały te same ustalone efekty, zostały porównane za pomocą AIC i zostały dopasowane przez REML (aby umożliwić porównanie różnych struktur korelacji przez AIC). Korzystam z pakietu R nlme i funkcji gls.

Pytanie 1

Resztki modeli GLS nadal wyświetlają prawie identyczne wzory cykliczne, gdy są drukowane w funkcji czasu. Czy takie wzorce zawsze pozostaną, nawet w modelach, które dokładnie uwzględniają strukturę autokorelacji?

Symulowałem niektóre uproszczone, ale podobne dane w R poniżej mojego drugiego pytania, które pokazuje problem w oparciu o moje obecne zrozumienie metod potrzebnych do oceny czasowo autokorelowanych wzorców w resztkach modelu , które teraz wiem, że są błędne (patrz odpowiedź).

pytanie 2

Do moich danych dopasowałem wszystkie modele GLS ze wszystkimi możliwymi prawdopodobnymi strukturami korelacji, ale w rzeczywistości żadne z nich nie są znacznie lepsze niż GLM bez żadnej struktury korelacji: tylko jeden model GLS jest nieznacznie lepszy (wynik AIC = 1,8 niższy), podczas gdy wszystkie pozostałe mają wyższe wartości AIC. Jest tak jednak tylko wtedy, gdy wszystkie modele są dopasowane przez REML, a nie ML, gdzie modele GLS są wyraźnie znacznie lepsze, ale rozumiem z książek statystycznych, że musisz używać REML tylko do porównywania modeli o różnych strukturach korelacji i tych samych stałych efektach z powodów Nie będę tutaj szczegółowo.

Biorąc pod uwagę wyraźnie czasowo autokorelowany charakter danych, jeśli żaden model nie jest nawet umiarkowanie lepszy od zwykłego GLM, jaki jest najbardziej odpowiedni sposób podjęcia decyzji, który model zastosować do wnioskowania, zakładając, że stosuję odpowiednią metodę (ostatecznie chcę użyć AIC, aby porównać różne kombinacje zmiennych)?

Q1 „symulacja” badająca wzorce resztkowe w modelach z odpowiednimi strukturami korelacji i bez nich

Wygeneruj symulowaną zmienną odpowiedzi z efektem cyklicznym „czasu” i dodatnim efektem liniowym „x”:

time <- 1:50
x <- sample(rep(1:25,each=2),50)
y <- rnorm(50,5,5) + (5 + 15*sin(2*pi*time/25)) + (x/1)

y powinien wyświetlać trend cykliczny w „czasie” z losową zmiennością:

plot(time,y)

I dodatni związek liniowy z „x” z losową odmianą:

plot(x,y)

Utwórz prosty liniowy model addytywny „y ~ time + x”:

require(nlme)
m1 <- gls(y ~ time + x, method="REML")

Model wyświetla wyraźne cykliczne wzory w resztach, gdy są wykreślane względem „czasu”, jak można się spodziewać:

plot(time, m1$residuals)

A jaki powinien być ładny, wyraźny brak jakiegokolwiek wzorca lub trendu w pozostałościach, gdy wykreślone w stosunku do „x”:

plot(x, m1$residuals)

Prosty model „y ~ time + x”, który obejmuje autoregresyjną strukturę korelacji rzędu 1, powinien pasować do danych znacznie lepiej niż w poprzednim modelu ze względu na strukturę autokorelacji, przy ocenie za pomocą AIC:

m2 <- gls(y ~ time + x, correlation = corAR1(form=~time), method="REML")
AIC(m1,m2)

Jednak model powinien nadal wyświetlać prawie identycznie „czasowo” autokorelowane reszty:

plot(time, m2$residuals)

Dziękuję bardzo za wszelkie porady.

Pytanie 1

Robisz tutaj dwie rzeczy źle. Pierwsza jest ogólnie rzeczą złą; nie robić w ogóle zagłębić modelowych obiektów i wyrwać komponentów. Naucz się w tym przypadku korzystać z funkcji ekstraktora resid(). W takim przypadku otrzymujesz coś przydatnego, ale jeśli miałbyś inny typ obiektu modelu, na przykład GLM glm(), mod$residualszawierałby resztki robocze z ostatniej iteracji IRLS i jest czymś, czego na ogół nie chcesz!

Drugą rzeczą, którą robisz źle, jest coś, co mnie przyłapało. Resztki, które wyodrębniłeś (a także byś je wyodrębnił, gdybyś użył resid()), to resztki surowe lub resztkowe z odpowiedzi. Zasadniczo jest to różnica między dopasowanymi wartościami i zaobserwowanych wartości odpowiedzi, biorąc pod uwagę efekty ustalone terminy tylko . Wartości te będą zawierały tę samą resztkową autokorelację, co m1dlatego, że ustalone efekty (lub, jeśli wolisz, predyktor liniowy) są takie same w dwóch modelach ( ~ time + x).

Aby uzyskać resztki zawierające podany termin korelacji, potrzebujesz znormalizowanych reszt. Dostajesz je wykonując:

resid(m1, type = "normalized")

Ten (i inne dostępne rodzaje pozostałości) opisano w ?residuals.gls:

type: an optional character string specifying the type of residuals
      to be used. If ‘"response"’, the "raw" residuals (observed -
      fitted) are used; else, if ‘"pearson"’, the standardized
      residuals (raw residuals divided by the corresponding
      standard errors) are used; else, if ‘"normalized"’, the
      normalized residuals (standardized residuals pre-multiplied
      by the inverse square-root factor of the estimated error
      correlation matrix) are used. Partial matching of arguments
      is used, so only the first character needs to be provided.
      Defaults to ‘"response"’.

Dla porównania, tutaj są ACF surowego (odpowiedź) i znormalizowanych reszt

layout(matrix(1:2))
acf(resid(m2))
acf(resid(m2, type = "normalized"))
layout(1)

Aby zobaczyć, dlaczego tak się dzieje i gdzie surowe reszty nie zawierają terminu korelacji, rozważ model, który pasowałeś

gdzie

$\Lambda$$\hat{\rho}$$\rho^{|d|}$$d$

Surowe reszty, domyślnie zwrócone, resid(m2)pochodzą tylko z liniowej części predyktora, a więc z tego bitu

$\Lambda$

Q2

Wygląda na to, że próbujesz dopasować trend nieliniowy za pomocą funkcji liniowej timei wyjaśnić brak dopasowania do „trendu” za pomocą AR (1) (lub innych struktur). Jeśli twoje dane są jak dane przykładowe, które podajesz tutaj, dopasowałbym GAM, aby umożliwić płynne działanie zmiennych towarzyszących. Ten model byłby

$\Lambda = \mathbf{I}$

library("mgcv")
m3 <- gam(y ~ s(time) + s(x), select = TRUE, method = "REML")

gdzie select = TRUEstosuje pewien dodatkowy skurcz, aby umożliwić modelowi usunięcie jednego z terminów z modelu.

Ten model daje

> summary(m3)

Family: gaussian 
Link function: identity 

Formula:
y ~ s(time) + s(x)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  23.1532     0.7104   32.59   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
          edf Ref.df      F  p-value    
s(time) 8.041      9 26.364  < 2e-16 ***
s(x)    1.922      9  9.749 1.09e-14 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

i ma gładkie terminy, które wyglądają tak:

Reszty z tego modelu są również lepiej zachowane (surowe reszty)

acf(resid(m3))

Teraz słowo ostrzeżenia; istnieje problem z wygładzaniem szeregów czasowych polegający na tym, że metody decydujące o tym, jak gładkie lub faliste są funkcje, zakładają, że dane są niezależne. W praktyce oznacza to, że płynna funkcja time ( s(time)) może pasować do informacji, które są naprawdę przypadkowym błędem autokorelowanym, a nie tylko leżącym u ich podstaw trendem. Dlatego należy zachować szczególną ostrożność przy dopasowywaniu wygładzaczy do danych szeregów czasowych.

Można to obejść na wiele sposobów, ale jednym ze sposobów jest przejście do dopasowania modelu, za pomocą gamm()którego wywołuje się lme()wewnętrznie i który pozwala użyć correlationargumentu użytego dla gls()modelu. Oto przykład

mm1 <- gamm(y ~ s(time, k = 6, fx = TRUE) + s(x), select = TRUE,
            method = "REML")
mm2 <- gamm(y ~ s(time, k = 6, fx = TRUE) + s(x), select = TRUE,
            method = "REML", correlation = corAR1(form = ~ time))

s(time)s(time)$\rho = 0$s(time)$\rho >> .5$

Model z AR (1) nie stanowi znaczącej poprawy w stosunku do modelu bez AR (1):

> anova(mm1$lme, mm2$lme)
        Model df      AIC      BIC    logLik   Test   L.Ratio p-value
mm1$lme     1  9 301.5986 317.4494 -141.7993                         
mm2$lme     2 10 303.4168 321.0288 -141.7084 1 vs 2 0.1817652  0.6699

Jeśli spojrzymy na szacunkową wartość $ \ hat {\ rho}}, zobaczymy

> intervals(mm2$lme)
....

 Correlation structure:
         lower      est.     upper
Phi -0.2696671 0.0756494 0.4037265
attr(,"label")
[1] "Correlation structure:"

gdzie Phi$\rho$$\rho$