Model linier sederhana dengan kesalahan autokorelasi dalam R [ditutup]

19

Bagaimana saya menyesuaikan model linier dengan kesalahan autokorelasi dalam R? Di stata saya akan menggunakan praisperintah, tapi saya tidak dapat menemukan setara R ...

Zach
sumber

Jawaban:

23

Lihat gls(kotak paling umum) dari paket nlme

Anda dapat mengatur profil korelasi untuk kesalahan dalam regresi, misalnya ARMA, dll:

 gls(Y ~ X, correlation=corARMA(p=1,q=1))

untuk kesalahan ARMA (1,1).

Dr G
sumber
1
Bisakah saya menggunakan fungsi "prediksi" pada dataset baru, dan mempertahankan struktur kesalahan yang sama? Bagaimana perintah gls mengetahui urutan pengamatannya?
Zach
27

Selain gls()fungsi dari nlme, Anda juga dapat menggunakan arima()fungsi dalam statspaket menggunakan MLE. Berikut ini adalah contoh dengan kedua fungsi tersebut.

x <- 1:100
e <- 25*arima.sim(model=list(ar=0.3),n=100)
y <- 1 + 2*x + e

###Fit the model using gls()
require(nlme)
(fit1 <- gls(y~x, corr=corAR1(0.5,form=~1)))
Generalized least squares fit by REML
  Model: y ~ x 
  Data: NULL 
  Log-restricted-likelihood: -443.6371

Coefficients:
(Intercept)           x 
   4.379304    1.957357 

Correlation Structure: AR(1)
 Formula: ~1 
 Parameter estimate(s):
      Phi 
0.3637263 
Degrees of freedom: 100 total; 98 residual
Residual standard error: 22.32908 

###Fit the model using arima()
(fit2 <- arima(y, xreg=x, order=c(1,0,0)))

Call:
arima(x = y, order = c(1, 0, 0), xreg = x)

Coefficients:
         ar1  intercept       x
      0.3352     4.5052  1.9548
s.e.  0.0960     6.1743  0.1060

sigma^2 estimated as 423.7:  log likelihood = -444.4,  aic = 896.81 

Keuntungan dari fungsi arima () adalah Anda dapat memuat berbagai proses kesalahan ARMA yang jauh lebih besar. Jika Anda menggunakan fungsi auto.arima () dari paket perkiraan, Anda dapat secara otomatis mengidentifikasi kesalahan ARMA:

require(forecast)    
fit3 <- auto.arima(y, xreg=x)
Rob Hyndman
sumber
1
Apa 0,5 untuk dalam "corr = corAR1 (0,5, form = ~ 1)?"
Zach
1
Ini memberikan nilai awal untuk optimasi. Hampir tidak ada bedanya jika dihilangkan.
Rob Hyndman
1
+1 arimaOpsi ini terlihat lebih berbeda dari Stata praispada pandangan pertama, tetapi lebih fleksibel dan Anda juga dapat menggunakan tsdiaguntuk mendapatkan visual yang bagus tentang seberapa baik asumsi AR (1) Anda benar-benar cocok.
Wayne
1
Bagaimana saya menggunakan arima () jika saya mundur hanya dengan konstanta?
user43790
7

Gunakan fungsi gls dari paket nlme . Inilah contohnya.

##Generate data frame with regressor and AR(1) error. The error term is 
## \eps_t=0.3*\eps_{t-1}+v_t
df <- data.frame(x1=rnorm(100), err=filter(rnorm(100)/5,filter=0.3,method="recursive"))

##Create ther response
df$y <- 1 + 2*df$x + df$err

###Fit the model
gls(y~x, data=df, corr=corAR1(0.5,form=~1))

Generalized least squares fit by REML
  Model: y ~ x 
  Data: df 
  Log-restricted-likelihood: 9.986475

 Coefficients:
 (Intercept)           x 
   1.040129    2.001884 

 Correlation Structure: AR(1)
 Formula: ~1 
 Parameter estimate(s):
     Phi 
 0.2686271 
Degrees of freedom: 100 total; 98 residual
Residual standard error: 0.2172698 

Karena model dipasang menggunakan kemungkinan maksimum, Anda perlu memberikan nilai awal. Nilai awal default adalah 0, tetapi seperti biasa, baik untuk mencoba beberapa nilai untuk memastikan konvergensi.

Seperti yang ditunjukkan oleh Dr. G, Anda juga dapat menggunakan struktur korelasi lainnya, yaitu ARMA.

Perhatikan bahwa secara umum estimasi kuadrat terkecil konsisten jika matriks kovarians dari kesalahan regresi bukan kelipatan dari matriks identitas, jadi jika Anda mencocokkan model dengan struktur kovarians tertentu, pertama-tama Anda perlu menguji apakah itu sesuai.

mpiktas
sumber
Apa 0,5 untuk dalam "corr = corAR1 (0,5, form = ~ 1)?"
Zach
3

Anda dapat menggunakan prediksi pada output gls. Lihat? Memprediksi. Anda juga dapat menentukan urutan pengamatan dengan istilah "bentuk" dalam struktur korelasi. Sebagai contoh:
corr=corAR1(form=~1)menunjukkan bahwa urutan data adalah yang ada di tabel. corr=corAR1(form=~Year)menunjukkan bahwa urutan adalah salah satu faktor Tahun .. Akhirnya nilai "0,5" pada corr=corAR1(0.5,form=~1)?umumnya diatur hingga nilai parameter yang diperkirakan mewakili struktur varians (phi, dalam kasus AR, theta dalam kasus MA .. .). Ini opsional untuk mengaturnya dan digunakan untuk optimasi seperti yang disebutkan oleh Rob Hyndman.

Xochitl C.
sumber
Nilai yang diprediksi atau dipasang akan jauh berbeda meskipun benar (gls () versus Arima ())? Karena gls hanya akan menggunakan efek tetap sementara Arima akan memasukkan kesalahan arima ke nilai yang dipasang.
B_Miner