Prosty model liniowy z błędami autokorelowanymi w R [zamknięty]

19

Jak dopasować model liniowy z błędami autokorelowanymi w R? W stacie użyłbym praispolecenia, ale nie mogę znaleźć odpowiednika R ...

Zach
źródło

Odpowiedzi:

23

Spójrz na gls(uogólnione najmniejsze kwadraty) z pakietu nlme

Możesz ustawić profil korelacji dla błędów w regresji, np. ARMA, itp:

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

dla błędów ARiMR (1,1).

Dr G.
źródło
1
Czy mogę użyć funkcji „przewidywania” w nowym zbiorze danych i zachować tę samą strukturę błędów? Skąd polecenie gls wie, w jakiej kolejności są obserwacje?
Zach
27

Oprócz gls()funkcji from nlmemożesz także użyć arima()funkcji w statspakiecie za pomocą MLE. Oto przykład z obiema funkcjami.

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 

Zaletą funkcji arima () jest to, że można dopasować znacznie większą różnorodność procesów błędów ARMA. Jeśli użyjesz funkcji auto.arima () z pakietu prognozy, możesz automatycznie zidentyfikować błąd ARMA:

require(forecast)    
fit3 <- auto.arima(y, xreg=x)
Rob Hyndman
źródło
1
Do czego służy 0,5 w „corr = corAR1 (0,5, forma = ~ 1)?”
Zach
1
Daje wartość początkową optymalizacji. Nie ma prawie żadnej różnicy, jeśli zostanie pominięty.
Rob Hyndman
1
+1 arimaOpcja praisna pierwszy rzut oka wygląda inaczej niż Stata , ale jest bardziej elastyczna i można jej również użyć, tsdiagaby uzyskać ładny obraz tego, jak dobrze pasuje twoje założenie AR (1).
Wayne
1
Jak korzystać z arima (), jeśli regresuję y tylko na stałej?
user43790,
7

Użyj funkcji gls z pakietu nlme . Oto przykład.

##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 

Ponieważ model jest montowany z maksymalnym prawdopodobieństwem, musisz podać wartości początkowe. Domyślna wartość początkowa to 0, ale jak zawsze dobrze jest wypróbować kilka wartości, aby zapewnić zbieżność.

Jak zauważył dr G, możesz również użyć innych struktur korelacji, a mianowicie ARMA.

Zauważ, że ogólnie szacunki najmniejszych kwadratów są spójne, jeśli macierz kowariancji błędów regresji nie jest wielokrotnością macierzy tożsamości, więc jeśli dopasujesz model do określonej struktury kowariancji, najpierw musisz sprawdzić, czy jest on odpowiedni.

mpiktas
źródło
Do czego służy 0,5 w „corr = corAR1 (0,5, forma = ~ 1)?”
Zach
3

Możesz użyć przewidywania na wyjściu gls. Widzisz? Prognozuj.gls. Możesz także określić kolejność obserwacji według terminu „forma” w strukturze korelacji. Na przykład:
corr=corAR1(form=~1)wskazuje, że kolejność danych odpowiada kolejności w tabeli. corr=corAR1(form=~Year)wskazuje, że kolejność jest zgodna z współczynnikiem Rok. Wreszcie wartość „0,5” corr=corAR1(0.5,form=~1)?jest na ogół ustawiana na wartość parametru oszacowanego jako reprezentujący strukturę wariancji (phi, w przypadku AR, theta w przypadku MA. .). Opcjonalnie można go skonfigurować i używać do optymalizacji, jak wspomniał Rob Hyndman.

Xochitl C.
źródło
Przewidywane lub dopasowane wartości będą znacznie różne, choć prawidłowe (gls () w porównaniu z Arima ())? Ponieważ gls będą używać tylko stałych efektów, podczas gdy Arima będzie zawierać błędy arima do dopasowanych wartości.
B_Miner