Konteks
Pertanyaan ini menggunakan R, tetapi tentang masalah statistik umum.
Saya sedang menganalisis efek dari faktor kematian (% kematian karena penyakit dan parasitisme) pada tingkat pertumbuhan populasi ngengat dari waktu ke waktu, di mana populasi larva diambil sampelnya dari 12 lokasi setahun sekali selama 8 tahun. Data tingkat pertumbuhan populasi menampilkan tren siklus yang jelas namun tidak teratur dari waktu ke waktu.
Sisa dari model linear umum sederhana (tingkat pertumbuhan ~% penyakit +% parasitisme + tahun) menunjukkan tren siklus yang jelas tetapi tidak teratur dari waktu ke waktu. Oleh karena itu, model kuadrat terkecil umum dari bentuk yang sama juga dipasang pada data dengan struktur korelasi yang sesuai untuk menangani autokorelasi temporal, misalnya simetri senyawa, urutan proses autoregresif 1 dan struktur korelasi rata-rata bergerak autoregresif.
Semua model berisi efek tetap yang sama, dibandingkan dengan menggunakan AIC, dan dipasang oleh REML (untuk memungkinkan perbandingan struktur korelasi yang berbeda dengan AIC). Saya menggunakan paket R nlme dan fungsi gls.
pertanyaan 1
Residu model GLS masih menampilkan pola siklus yang hampir identik ketika diplot dengan waktu. Apakah pola seperti itu akan selalu ada, bahkan dalam model yang secara akurat menjelaskan struktur autokorelasi?
Saya telah mensimulasikan beberapa data yang disederhanakan tetapi serupa dalam R di bawah pertanyaan kedua saya, yang menunjukkan masalah berdasarkan pemahaman saya saat ini tentang metode yang diperlukan untuk menilai pola autokorelasi sementara dalam residual model , yang sekarang saya tahu salah (lihat jawaban).
Pertanyaan 2
Saya telah memasang model GLS dengan semua struktur korelasi yang mungkin masuk akal ke data saya, tetapi tidak ada yang secara substansial lebih pas daripada GLM tanpa struktur korelasi apa pun: hanya satu model GLS yang sedikit lebih baik (skor AIC = 1,8 lebih rendah), sementara sisanya memiliki nilai AIC yang lebih tinggi. Namun, ini hanya terjadi ketika semua model dipasang oleh REML, bukan ML di mana model GLS jelas jauh lebih baik, tetapi saya mengerti dari buku statistik Anda hanya harus menggunakan REML untuk membandingkan model dengan struktur korelasi yang berbeda dan efek tetap yang sama karena alasan Saya tidak akan merinci di sini.
Mengingat sifat data yang secara autokorelasi jelas sementara, jika tidak ada model yang bahkan lebih baik daripada GLM sederhana apa cara yang paling tepat untuk memutuskan model mana yang akan digunakan untuk inferensi, dengan asumsi saya menggunakan metode yang sesuai (akhirnya saya ingin menggunakan AIC untuk membandingkan kombinasi variabel yang berbeda)?
Q1 'simulasi' mengeksplorasi pola residual dalam model dengan dan tanpa struktur korelasi yang sesuai
Hasilkan variabel respons simulasi dengan efek siklus 'waktu', dan efek linear positif '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 harus menampilkan tren siklus selama 'waktu' dengan variasi acak:
plot(time,y)
Dan hubungan linier positif dengan 'x' dengan variasi acak:
plot(x,y)
Buat model aditif linier sederhana "y ~ time + x":
require(nlme)
m1 <- gls(y ~ time + x, method="REML")
Model ini menampilkan pola siklus yang jelas dalam residu ketika diplot terhadap 'waktu', seperti yang diharapkan:
plot(time, m1$residuals)
Dan apa yang seharusnya menjadi kekurangan, baik jelas dari pola atau tren dalam residu ketika diplot terhadap 'x':
plot(x, m1$residuals)
Model sederhana "y ~ waktu + x" yang mencakup struktur korelasi autoregresif urutan 1 harus sesuai dengan data yang jauh lebih baik daripada model sebelumnya karena struktur autokorelasi, ketika dinilai menggunakan AIC:
m2 <- gls(y ~ time + x, correlation = corAR1(form=~time), method="REML")
AIC(m1,m2)
Namun, model tersebut harus tetap menampilkan residu autokorelasi yang 'identik sementara':
plot(time, m2$residuals)
Terima kasih banyak atas sarannya.
Jawaban:
Q1
Anda melakukan dua hal yang salah di sini. Yang pertama adalah hal yang umumnya buruk; jangan secara umum mempelajari objek model dan merobek komponen. Pelajari cara menggunakan fungsi ekstraktor, dalam hal ini
resid()
. Dalam hal ini Anda mendapatkan sesuatu yang bermanfaat tetapi jika Anda memiliki jenis objek model yang berbeda, seperti GLM dariglm()
, makamod$residuals
akan berisi residu yang berfungsi dari iterasi IRLS terakhir dan merupakan sesuatu yang umumnya tidak Anda inginkan!Hal kedua yang Anda lakukan salah adalah sesuatu yang membuat saya keluar juga. Residu yang Anda ekstrak (dan juga akan diekstraksi jika Anda akan gunakan
resid()
) adalah residu mentah atau respons. Pada dasarnya ini adalah perbedaan antara nilai yang dipasang dan nilai yang diamati dari respons, dengan mempertimbangkan ketentuan efek tetap saja . Nilai-nilai ini akan mengandung autokorelasi residu yang sama denganm1
karena efek tetap (atau jika Anda lebih suka, prediktor linier) adalah sama pada kedua model (~ time + x
).Untuk mendapatkan residu yang menyertakan istilah korelasi yang Anda tentukan, Anda perlu residu yang dinormalisasi . Anda mendapatkan ini dengan melakukan:
Ini (dan jenis residu lain yang tersedia) dijelaskan dalam
?residuals.gls
:Dengan perbandingan, berikut adalah ACFs dari raw (response) dan residual yang dinormalisasi
Untuk melihat mengapa ini terjadi, dan di mana residu mentah tidak menyertakan istilah korelasi, pertimbangkan model yang Anda paskan
dimana
Residual mentah, default dikembalikan oleh
resid(m2)
hanya dari bagian prediktor linier, maka dari ini sedikitQ2
Tampaknya Anda mencoba menyesuaikan tren non-linear dengan fungsi linier
time
dan memperhitungkan kurangnya kecocokan dengan "tren" dengan AR (1) (atau struktur lain). Jika data Anda sama dengan contoh data yang Anda berikan di sini, saya akan cocok dengan GAM untuk memungkinkan fungsi kovariat yang lancar. Model ini akan menjadidi mana
select = TRUE
berlaku beberapa penyusutan tambahan untuk memungkinkan model untuk menghapus salah satu syarat dari model.Model ini memberi
dan memiliki istilah halus yang terlihat seperti ini:
Residu dari model ini juga berperilaku lebih baik (residu mentah)
Sekarang kata hati-hati; ada masalah dengan smoothing time series dalam hal metode yang memutuskan seberapa halus atau goyahnya fungsi diasumsikan bahwa data independen. Apa artinya ini secara praktis adalah bahwa fungsi waktu yang mulus (
s(time)
) dapat cocok dengan informasi yang benar-benar kesalahan otokorelasi acak dan tidak hanya tren yang mendasarinya. Oleh karena itu Anda harus sangat berhati-hati ketika menyesuaikan data deret waktu.Ada beberapa cara di sekitar ini, tetapi satu cara adalah beralih untuk menyesuaikan model melalui
gamm()
mana panggilanlme()
internal dan yang memungkinkan Anda untuk menggunakancorrelation
argumen yang Anda gunakan untukgls()
model. Berikut ini sebuah contohs(time)
s(time)
s(time)
Model dengan AR (1) tidak mewakili peningkatan yang signifikan atas model tanpa AR (1):
Jika kita melihat perkiraan untuk $ \ hat {\ rho}} kita lihat
Phi
sumber