R help - Mar 2004 - corARMA and ACF in nlme

Hi R-sters,

Just wondering what I might be doing wrong.  I'm trying to fit a multiple 
linear regression model, and being ever mindful about the possibilities of 
autocorrelation in the errors (it's a time series), the errors appear to 
follow an AR1 process (ar(ts(glsfit$residuals)) selected order 1).  So, 
when I go back and try to do the simultaneous regression and error fit with 
gls, the acf and pacf plots of residuals from the old model (glsfit) and 
those plots of the new model (glsAR1fit, below)  look exactly the same (a 
significant autocorrelation at lag of 1).

Any ideas out there as to what I may be doing wrong?  Is there an error in 
my code?

Here's my R code for the simultaneous model fit (taking a phi estimate=0.6 
from a previous step <ACF(glsfit)>):

glsAR1fit<-gls(y~x1+x2+x3+x4, na.action = na.omit, subset=12:54,
                 correlation = corARMA(0.6, p=1, q=0, fixed = FALSE))

Thanks much,

Jeff

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Jeff Jorgensen

Center for Limnology
University of Wisconsin Madison           ph (608) 263-2304
680 North Park Street                           fx (608) 265-2340
Madison, Wisconsin 53706                http://limnology.wisc.edu

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

On Tue, 9 Mar 2004, Jeff Jorgensen wrote:
> Just wondering what I might be doing wrong.  I'm trying to fit a
multiple
> linear regression model, and being ever mindful about the possibilities of 
> autocorrelation in the errors (it's a time series), the errors appear
to
> follow an AR1 process (ar(ts(glsfit$residuals)) selected order 1).  So, 
> when I go back and try to do the simultaneous regression and error fit with
> gls, 
That's not really what you did: gls fits a multivariate normal 
distribution with covariance matrix like that of observations from an AR1 
rather than fit an AR model.  A fine distinction, but as the rest of the 
para shows, it does matter.
> the acf and pacf plots of residuals from the old model (glsfit) and 
> those plots of the new model (glsAR1fit, below)  look exactly the same (a 
> significant autocorrelation at lag of 1).
They should.  Those are the residuals from the regression, not the 
innovations of the autoregression fitted to the residuals.  So the problem 
lies in how you interpreted what you did, I believe.

You may find this clearer if you use arima which does fit an ARIMA model.
> Any ideas out there as to what I may be doing wrong?  Is there an error in 
> my code?
> 
> Here's my R code for the simultaneous model fit (taking a phi
estimate=0.6
> from a previous step <ACF(glsfit)>):
> 
> glsAR1fit<-gls(y~x1+x2+x3+x4, na.action = na.omit, subset=12:54,
>                  correlation = corARMA(0.6, p=1, q=0, fixed = FALSE))
-- 
Brian D. Ripley,                  ripley at stats.ox.ac.uk
Professor of Applied Statistics,  http://www.stats.ox.ac.uk/~ripley/
University of Oxford,             Tel:  +44 1865 272861 (self)
1 South Parks Road,                     +44 1865 272866 (PA)
Oxford OX1 3TG, UK                Fax:  +44 1865 272595

Jeff Jorgensen <jcjorgensen at wisc.edu> writes:
> Hi R-sters,
> 
> Just wondering what I might be doing wrong.  I'm trying to fit a
> multiple linear regression model, and being ever mindful about the
> possibilities of autocorrelation in the errors (it's a time series),
> the errors appear to follow an AR1 process (ar(ts(glsfit$residuals))
> selected order 1).  So, when I go back and try to do the simultaneous
> regression and error fit with gls, the acf and pacf plots of residuals
> from the old model (glsfit) and those plots of the new model
> (glsAR1fit, below)  look exactly the same (a significant
> autocorrelation at lag of 1).
> 
> Any ideas out there as to what I may be doing wrong?  Is there an
> error in my code?
This is one of the dangers of accessing model fit structures directly:
There is more than one way to define residuals for correlated data. Try
looking at help(residuals.gls).

(It's been a while, but as far as I remember, you can plot the
Variogram of the raw residuals and overlay the theoretical Variogram
of the fitted model, so raw residuals are useful too. Also,
standardized residuals are not uniquely defined since matrix square
roots aren't.)

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard at biostat.ku.dk)             FAX: (+45) 35327907