Daripada mencoba menguraikan deret waktu secara eksplisit, saya lebih suka menyarankan Anda memodelkan data spatio-temporal karena, seperti yang akan Anda lihat di bawah ini, tren jangka panjang kemungkinan bervariasi secara spasial, tren musiman bervariasi dengan tren jangka panjang dan secara spasial.
Saya telah menemukan bahwa model aditif umum (GAM) adalah model yang baik untuk pemasangan seri waktu tidak beraturan seperti yang Anda gambarkan.
Di bawah ini saya menggambarkan model cepat yang saya siapkan untuk data lengkap dari formulir berikut
E (ysaya)= α +f1(Rubahsaya) +f2(Doysaya) +f3(Tahunsaya) +f4(xsaya,ysaya) +f5(Doysaya,Tahunsaya) +f6(xsaya,ysaya,Rubahsaya) +f7(xsaya,ysaya,Doysaya) +f8(xsaya,ysaya,Tahunsaya)
dimana
- α adalah model mencegat,
- f1(Rubahsaya) adalah fungsi waktu yang lancar,
- f2(Doysaya) adalah fungsi hari yang lancar,
- f3(Tahunsaya) adalah fungsi tahun yang lancar,
- f4(xsaya,ysaya) adalah kelancaran 2D garis bujur dan lintang,
- f5(Doysaya,Tahunsaya) adalah kelancaran produk tensor dari tahun ke tahun,
- f6(xsaya,ysaya,Rubahsaya) produk tensor kelancaran lokasi & waktu dalam sehari
- f7(xsaya,ysaya,Doysaya) kelancaran produk tensor lokasi hari tahun &
- f8(xsaya,ysaya,Tahunsaya produk tensor kelancaran lokasi & tahun
Secara efektif, empat smooths pertama adalah efek utama
- waktu hari,
- musim,
- tren jangka panjang,
- variasi spasial
sementara empat produk tensor yang tersisa memodelkan interaksi yang halus antara kovariat yang dinyatakan, model mana
- bagaimana pola musiman suhu bervariasi dari waktu ke waktu,
- bagaimana efek waktu hari bervariasi secara spasial,
- bagaimana efek musiman bervariasi secara spasial, dan
- bagaimana tren jangka panjang bervariasi secara spasial
Data dimuat ke R dan dipijat sedikit dengan kode berikut
library('mgcv')
library('ggplot2')
library('viridis')
theme_set(theme_bw())
library('gganimate')
galveston <- read.csv('gbtemp.csv')
galveston <- transform(galveston,
datetime = as.POSIXct(paste(DATE, TIME),
format = '%m/%d/%y %H:%M', tz = "CDT"))
galveston <- transform(galveston,
STATION_ID = factor(STATION_ID),
DoY = as.numeric(format(datetime, format = '%j')),
ToD = as.numeric(format(datetime, format = '%H')) +
(as.numeric(format(datetime, format = '%M')) / 60))
Model itu sendiri dilengkapi menggunakan bam()
fungsi yang dirancang untuk menyesuaikan GAM ke set data yang lebih besar seperti ini. Anda dapat menggunakan gam()
untuk model ini juga, tetapi akan membutuhkan waktu agak lebih lama.
knots <- list(DoY = c(0.5, 366.5))
M <- list(c(1, 0.5), NA)
m <- bam(MEASUREMENT ~
s(ToD, k = 10) +
s(DoY, k = 30, bs = 'cc') +
s(YEAR, k = 30) +
s(LONGITUDE, LATITUDE, k = 100, bs = 'ds', m = c(1, 0.5)) +
ti(DoY, YEAR, bs = c('cc', 'tp'), k = c(15, 15)) +
ti(LONGITUDE, LATITUDE, ToD, d = c(2,1), bs = c('ds','tp'),
m = M, k = c(20, 10)) +
ti(LONGITUDE, LATITUDE, DoY, d = c(2,1), bs = c('ds','cc'),
m = M, k = c(25, 15)) +
ti(LONGITUDE, LATITUDE, YEAR, d = c(2,1), bs = c('ds','tp'),
m = M), k = c(25, 15)),
data = galveston, method = 'fREML', knots = knots,
nthreads = 4, discrete = TRUE)
The s()
istilah adalah efek utama, sementara ti()
istilah yang tensor produk interaksi menghaluskan mana efek utama dari kovariat bernama telah dihapus dari dasar. Ini ti()
menghaluskan adalah cara untuk memasukkan interaksi dari variabel yang dinyatakan dengan cara numerik yang stabil.
The knots
objek hanya pengaturan titik akhir dari siklus yang halus saya gunakan untuk hari berlaku tahun - kami ingin 23:59 pada 31 Desember untuk bergabung dengan lancar dengan 00:01 1 Jan. Ini menyumbang sampai batas tertentu untuk tahun kabisat.
Ringkasan model menunjukkan semua efek ini signifikan;
> summary(m)
Family: gaussian
Link function: identity
Formula:
MEASUREMENT ~ s(ToD, k = 10) + s(DoY, k = 12, bs = "cc") + s(YEAR,
k = 30) + s(LONGITUDE, LATITUDE, k = 100, bs = "ds", m = c(1,
0.5)) + ti(DoY, YEAR, bs = c("cc", "tp"), k = c(12, 15)) +
ti(LONGITUDE, LATITUDE, ToD, d = c(2, 1), bs = c("ds", "tp"),
m = list(c(1, 0.5), NA), k = c(20, 10)) + ti(LONGITUDE,
LATITUDE, DoY, d = c(2, 1), bs = c("ds", "cc"), m = list(c(1,
0.5), NA), k = c(25, 12)) + ti(LONGITUDE, LATITUDE, YEAR,
d = c(2, 1), bs = c("ds", "tp"), m = list(c(1, 0.5), NA),
k = c(25, 15))
Parametric coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 21.75561 0.07508 289.8 <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(ToD) 3.036 3.696 5.956 0.000189 ***
s(DoY) 9.580 10.000 3520.098 < 2e-16 ***
s(YEAR) 27.979 28.736 59.282 < 2e-16 ***
s(LONGITUDE,LATITUDE) 54.555 99.000 4.765 < 2e-16 ***
ti(DoY,YEAR) 131.317 140.000 34.592 < 2e-16 ***
ti(ToD,LONGITUDE,LATITUDE) 42.805 171.000 0.880 < 2e-16 ***
ti(DoY,LONGITUDE,LATITUDE) 83.277 240.000 1.225 < 2e-16 ***
ti(YEAR,LONGITUDE,LATITUDE) 84.862 329.000 1.101 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
R-sq.(adj) = 0.94 Deviance explained = 94.2%
fREML = 29807 Scale est. = 2.6318 n = 15276
Analisis yang lebih hati-hati ingin memeriksa apakah kita memerlukan semua interaksi ini; beberapa ti()
istilah spasial hanya menjelaskan sejumlah kecil variasi dalam data, seperti yang ditunjukkan olehFstatistik; ada banyak data di sini sehingga bahkan ukuran efek kecil secara statistik signifikan tetapi tidak menarik.
Namun, sebagai pemeriksaan cepat, menghilangkan tiga ti()
smooth spasial ( m.sub
), menghasilkan kecocokan yang secara signifikan lebih buruk seperti yang dinilai oleh AIC:
> AIC(m, m.sub)
df AIC
m 447.5680 58583.81
m.sub 239.7336 59197.05
Kita dapat memplot efek parsial dari lima smooths pertama menggunakan plot()
metode - smooths produk tensor 3D tidak dapat diplot dengan mudah dan tidak secara default.
plot(m, pages = 1, scheme = 2, shade = TRUE, scale = 0)
The scale = 0
argumen di sana menempatkan semua plot pada skala mereka sendiri, untuk membandingkan besaran dampak, kita dapat mematikan ini:
plot(m, pages = 1, scheme = 2, shade = TRUE)
Sekarang kita dapat melihat bahwa efek musiman mendominasi. Tren jangka panjang (rata-rata) ditunjukkan di plot kanan atas. Untuk benar-benar melihat tren jangka panjang, Anda harus memilih stasiun dan kemudian memprediksi dari model untuk stasiun itu, menetapkan waktu hari dan hari dalam setahun ke beberapa nilai representatif (tengah hari, untuk hari dalam setahun di musim panas) mengatakan). Ada awal atau dua tahun dari seri memiliki beberapa nilai suhu rendah relatif terhadap sisa catatan, yang kemungkinan diambil di semua smooths yang terlibat YEAR
. Data-data ini harus dilihat lebih dekat.
Ini bukan tempat yang tepat untuk membahas hal itu, tetapi di sini ada beberapa visualisasi model yang pas. Pertama saya melihat pola suhu spasial dan bagaimana itu bervariasi selama bertahun-tahun seri. Untuk melakukan itu saya memprediksi dari model untuk grid 100x100 di atas domain spasial, pada tengah hari pada hari 180 setiap tahun:
pdata <- with(galveston,
expand.grid(ToD = 12,
DoY = 180,
YEAR = seq(min(YEAR), max(YEAR), by = 1),
LONGITUDE = seq(min(LONGITUDE), max(LONGITUDE), length = 100),
LATITUDE = seq(min(LATITUDE), max(LATITUDE), length = 100)))
fit <- predict(m, pdata)
kemudian saya set ke missing,, NA
nilai prediksi fit
untuk semua titik data yang terletak agak jauh dari pengamatan (proporsional; dist
)
ind <- exclude.too.far(pdata$LONGITUDE, pdata$LATITUDE,
galveston$LONGITUDE, galveston$LATITUDE, dist = 0.1)
fit[ind] <- NA
dan gabungkan prediksi dengan data prediksi
pred <- cbind(pdata, Fitted = fit)
Menetapkan nilai yang diprediksi NA
seperti ini akan menghentikan kami melakukan ekstrapolasi di luar dukungan data.
Menggunakan ggplot2
ggplot(pred, aes(x = LONGITUDE, y = LATITUDE)) +
geom_raster(aes(fill = Fitted)) + facet_wrap(~ YEAR, ncol = 12) +
scale_fill_viridis(name = expression(degree*C), option = 'plasma',
na.value = 'transparent') +
coord_quickmap() +
theme(legend.position = 'top', legend.key.width = unit(2, 'cm'))
kami mendapatkan yang berikut ini
Kita dapat melihat variasi suhu tahun-ke-tahun dengan sedikit lebih detail jika kita menghidupkan daripada melihat sisi plotnya
p <- ggplot(pred, aes(x = LONGITUDE, y = LATITUDE, frame = YEAR)) +
geom_raster(aes(fill = Fitted)) +
scale_fill_viridis(name = expression(degree*C), option = 'plasma',
na.value = 'transparent') +
coord_quickmap() +
theme(legend.position = 'top', legend.key.width = unit(2, 'cm'))+
labs(x = 'Longitude', y = 'Latitude')
gganimate(p, 'galveston.gif', interval = .2, ani.width = 500, ani.height = 800)
Untuk melihat tren jangka panjang secara lebih rinci, kami dapat memprediksi stasiun tertentu. Misalnya, untuk STATION_ID
13364 dan prediksi hari di empat kuartal, kita dapat menggunakan yang berikut untuk menyiapkan nilai kovariat yang ingin kita prediksi pada (tengah hari, pada hari tahun 1, 90, 180, dan 270, di stasiun yang dipilih , dan mengevaluasi tren jangka panjang pada 500 nilai yang berjarak sama)
pdata <- with(galveston,
expand.grid(ToD = 12,
DoY = c(1, 90, 180, 270),
YEAR = seq(min(YEAR), max(YEAR), length = 500),
LONGITUDE = -94.8751,
LATITUDE = 29.50866))
Kemudian kami memperkirakan dan meminta kesalahan standar, untuk membentuk perkiraan interval kepercayaan 95%
fit <- data.frame(predict(m, newdata = pdata, se.fit = TRUE))
fit <- transform(fit, upper = fit + (2 * se.fit), lower = fit - (2 * se.fit))
pred <- cbind(pdata, fit)
yang kami plot
ggplot(pred, aes(x = YEAR, y = fit, group = factor(DoY))) +
geom_ribbon(aes(ymin = lower, ymax = upper), fill = 'grey', alpha = 0.5) +
geom_line() + facet_wrap(~ DoY, scales = 'free_y') +
labs(x = NULL, y = expression(Temperature ~ (degree * C)))
memproduksi
Jelas, ada lebih banyak untuk memodelkan data ini daripada yang saya tunjukkan di sini, dan kami ingin memeriksa autokorelasi residu dan overfitting dari splines, tetapi mendekati masalah sebagai salah satu pemodelan fitur data memungkinkan untuk lebih detail pemeriksaan tren.
Anda tentu saja bisa hanya memodelkan masing STATION_ID
- masing secara terpisah, tetapi itu akan membuang data, dan banyak stasiun memiliki beberapa pengamatan. Di sini model meminjam dari semua informasi stasiun untuk mengisi kekosongan dan membantu memperkirakan tren bunga.
Beberapa catatan bam()
The bam()
model menggunakan semua mgcv trik 's untuk memperkirakan model cepat (beberapa benang 4 ), REML cepat pemilihan kelancaran ( method = 'fREML'
), dan diskritisasi kovariat. Dengan opsi ini diaktifkan, model ini dapat digunakan dalam waktu kurang dari satu menit pada workstation Xeon dual-core Xeon dual-core 2013 saya dengan 64Gb RAM.
Salah satu opsi adalah menggunakan strategi regresi yang memperlakukan data Anda sebagai menarik dari fungsi mendasar yang berkelanjutan (lihat diskusi di sini: Apakah ada standar emas untuk pemodelan deret waktu yang tidak beraturan? ). Dari sana Anda dapat menggunakan metode seperti Singular Spectrum Analysis untuk menguraikan sinyal (untuk R: https://cran.r-project.org/web/packages/Rssa/index.html ).
sumber