19
$\begingroup$

How do I fit a linear model with autocorrelated errors in R? In stata I would use the prais command, but I can't find an R equivalent...

$\endgroup$

closed as off-topic by gung Apr 14 at 14:13

This question appears to be off-topic. The users who voted to close gave this specific reason:

  • "This question appears to be off-topic because EITHER it is not about statistics, machine learning, data analysis, data mining, or data visualization, OR it focuses on programming, debugging, or performing routine operations within a statistical computing platform. If the latter, you could try the support links we maintain." – gung
If this question can be reworded to fit the rules in the help center, please edit the question.

23
$\begingroup$

Have a look at gls (generalized least squares) from the package nlme

You can set a correlation profile for the errors in the regression, e.g. ARMA, etc:

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

for ARMA(1,1) errors.

$\endgroup$
  • 1
    $\begingroup$ Can I use the "predict" function on a new dataset, and retain the same error structure? How does the gls command know what order the observations are in? $\endgroup$ – Zach Jan 24 '11 at 0:23
27
$\begingroup$

In addition to the gls() function from nlme, you can also use the arima() function in the stats package using MLE. Here is an example with both functions.

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 

The advantage of the arima() function is that you can fit a much larger variety of ARMA error processes. If you use the auto.arima() function from the forecast package, you can automatically identify the ARMA error:

require(forecast)    
fit3 <- auto.arima(y, xreg=x)
$\endgroup$
  • 1
    $\begingroup$ What is the 0.5 for in "corr=corAR1(0.5,form=~1)?" $\endgroup$ – Zach Jan 24 '11 at 0:30
  • 1
    $\begingroup$ It gives a starting value for the optimization. It makes almost no difference if it is omitted. $\endgroup$ – Rob Hyndman Jan 24 '11 at 2:11
  • 1
    $\begingroup$ +1 The arima option looks more different from Stata's prais at first glance, but it's more flexible and you can also use tsdiag to get a nice visual of how well your AR(1) assumption actually fits. $\endgroup$ – Wayne Aug 24 '13 at 14:37
  • 1
    $\begingroup$ How do I use arima() if I regress y on a constant only? $\endgroup$ – user43790 Nov 9 '15 at 20:06
7
$\begingroup$

Use function gls from package nlme. Here is the example.

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

Since model is fitted using maximum likelihood you need to supply starting values. The default starting value is 0, but as always it is good to try several values to ensure the convergence.

As Dr. G pointed out you can also use other correlation structures, namely ARMA.

Note that in general least squares estimates are consistent if covariance matrix of regression errors is not multiple of identity matrix, so if you fit model with specific covariance structure, first you need to test whether it is appropriate.

$\endgroup$
  • $\begingroup$ What is the 0.5 for in "corr=corAR1(0.5,form=~1)?" $\endgroup$ – Zach Jan 24 '11 at 0:29
3
$\begingroup$

You can use predict on gls output. See ?predict.gls. Also you can specify the order of the observation by the "form" term in the correlation structure. For example:
corr=corAR1(form=~1) indicates that order of the data is the one they are in the table. corr=corAR1(form=~Year) indicates that the order is the one of factor Year.. Finally the "0.5" value in corr=corAR1(0.5,form=~1)? is generally set up to the value of the parameter estimated to represent the variance structure (phi, in case of AR, theta in case of MA...). It is optional to set it up and use for optimization as Rob Hyndman mentioned.

$\endgroup$
  • $\begingroup$ Predicting or fitted values will be far different though correct (gls() versus Arima())? As gls will only use fixed effects while Arima will include the arima errors to the fitted values. $\endgroup$ – B_Miner Sep 4 '15 at 17:40

Not the answer you're looking for? Browse other questions tagged or ask your own question.