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.
Odpowiedzi:
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 GLMglm()
,mod$residuals
zawierał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ę, com1
dlatego, ż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:
Ten (i inne dostępne rodzaje pozostałości) opisano w
?residuals.gls
:Dla porównania, tutaj są ACF surowego (odpowiedź) i znormalizowanych reszt
Aby zobaczyć, dlaczego tak się dzieje i gdzie surowe reszty nie zawierają terminu korelacji, rozważ model, który pasowałeś
gdzie
Surowe reszty, domyślnie zwrócone,
resid(m2)
pochodzą tylko z liniowej części predyktora, a więc z tego bituQ2
Wygląda na to, że próbujesz dopasować trend nieliniowy za pomocą funkcji liniowej
time
i 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łbygdzie
select = TRUE
stosuje pewien dodatkowy skurcz, aby umożliwić modelowi usunięcie jednego z terminów z modelu.Ten model daje
i ma gładkie terminy, które wyglądają tak:
Reszty z tego modelu są również lepiej zachowane (surowe reszty)
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ćcorrelation
argumentu użytego dlagls()
modelu. Oto przykłads(time)
s(time)
s(time)
Model z AR (1) nie stanowi znaczącej poprawy w stosunku do modelu bez AR (1):
Jeśli spojrzymy na szacunkową wartość $ \ hat {\ rho}}, zobaczymy
gdzieρ ρ
Phi
źródło