Nie zawracałbym stl()
sobie tym głowy - szerokość pasma dla bezużytecznej wygładzarki wykorzystywanej do wydobywania trendu jest daleko, daleko, aż do małych, co powoduje fluktuacje na małą skalę, które widzisz. Użyłbym modelu addytywnego. Oto przykład wykorzystujący dane i kod modelu z książki Simona Wooda na temat GAM:
require(mgcv)
require(gamair)
data(cairo)
cairo2 <- within(cairo, Date <- as.Date(paste(year, month, day.of.month,
sep = "-")))
plot(temp ~ Date, data = cairo2, type = "l")
![dane temperatury Cairo](https://i.stack.imgur.com/ZxNHU.png)
Dopasuj model do trendu i komponentów sezonowych --- ostrzeżenie, że jest to powolne:
mod <- gamm(temp ~ s(day.of.year, bs = "cc") + s(time, bs = "cr"),
data = cairo2, method = "REML",
correlation = corAR1(form = ~ 1 | year),
knots = list(day.of.year = c(0, 366)))
Dopasowany model wygląda następująco:
> summary(mod$gam)
Family: gaussian
Link function: identity
Formula:
temp ~ s(day.of.year, bs = "cc") + s(time, bs = "cr")
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 71.6603 0.1523 470.7 <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(day.of.year) 7.092 7.092 555.407 < 2e-16 ***
s(time) 1.383 1.383 7.035 0.00345 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.848 Scale est. = 16.572 n = 3780
i możemy wizualizować trendy i warunki sezonowe za pośrednictwem
plot(mod$gam, pages = 1)
![Kair dopasował trend i sezonowość](https://i.stack.imgur.com/7B6iy.png)
a jeśli chcemy wykreślić trend na obserwowanych danych, możemy to zrobić z prognozą poprzez:
pred <- predict(mod$gam, newdata = cairo2, type = "terms")
ptemp <- attr(pred, "constant") + pred[,2]
plot(temp ~ Date, data = cairo2, type = "l",
xlab = "year",
ylab = expression(Temperature ~ (degree*F)))
lines(ptemp ~ Date, data = cairo2, col = "red", lwd = 2)
![Trend dopasowany do Kairu](https://i.stack.imgur.com/hH4EN.png)
Lub to samo dla rzeczywistego modelu:
pred2 <- predict(mod$gam, newdata = cairo2)
plot(temp ~ Date, data = cairo2, type = "l",
xlab = "year",
ylab = expression(Temperature ~ (degree*F)))
lines(pred2 ~ Date, data = cairo2, col = "red", lwd = 2)
![Model dopasowany do Kairu](https://i.stack.imgur.com/15BYc.png)
To tylko przykład, a dogłębniejsza analiza może wymagać uwzględnienia kilku brakujących danych, ale powyższe powinno być dobrym punktem wyjścia.
Jeśli chodzi o twoje zdanie na temat kwantyfikacji trendu - to jest problem, ponieważ trend nie jest liniowy, ani w twojej stl()
wersji, ani w wersji GAM, którą pokazuję. Jeśli tak, możesz podać tempo zmian (nachylenie). Jeśli chcesz wiedzieć, o ile zmienił się szacowany trend w okresie próbkowania, możemy wykorzystać zawarte w nim dane pred
i obliczyć różnicę między początkiem i końcem serii tylko w komponencie trendu :
> tail(pred[,2], 1) - head(pred[,2], 1)
3794
1.756163
więc temperatury są średnio o 1,76 stopnia wyższe niż na początku zapisu.
gls()
w pakiecie nlme). Ale jak pokazano powyżej dla Kairu i STL sugeruje dla twoich danych, trend nie jest liniowy. Jako taki trend liniowy nie byłby odpowiedni, ponieważ nie opisuje poprawnie danych. Musisz wypróbować to na swoich danych, ale AM, taki jak ja, pokazałby się jako trend liniowy, gdyby najlepiej pasował do danych.Gavin pod warunkiem bardzo dokładną odpowiedź, ale na łatwiejszą i szybszą rozwiązania, zalecam ustawienie STL funkcji t.window parametr do wartości, która jest wielokrotnością częstotliwości z ts danych. Użyłbym wywnioskowanej okresowości odsetek (np. Wartość 3660 dla trendów dekadalnych z dobowymi danymi rozdzielczości). Może Cię również zainteresować pakiet stl2 opisany w rozprawie autora . Zastosowałem metodę Gavina do własnych danych i jest ona również bardzo skuteczna.
źródło