[R] correlation range estimates with nlme::gls

Ben Bolker bolker at zoo.ufl.edu
Mon Apr 11 23:52:50 CEST 2005


   I'm trying to do a simple (?) analysis of a 1D spatial data set, 
allowing for spatial autocorrelation.  (Actually, I'm comparing expected 
vs. observed for a spatial model of a 1D spatial data set.)  I'm using 
models like

gls(obs~exp,correlation=corExp(form=~pos),data=data)

or

gls(obs~exp,correlation=corLin(form=~pos),data=data)

  This form is supposed to fit a linear model of obs=a*exp+b using an 
autocorrelation model based on the position variable "pos".

   The problem: I get reasonable answers for the slope & intercept, but the 
estimated ranges for the autocorrelation functions are huge and seemingly 
unconnected to the pictures I get when I plot acf(resid(M0)) [where M0 is 
the OLS fit to the data].  I tried simulating a bunch of "data" with 
different data sizes: the disturbing thing is that the range results get 
(much) worse, not better, as the number of data points goes up [from n=30, 
around the real size, to n=50, to n=100].  When I am able to get 
confidence intervals on the range, they often don't make sense either. 
On the other hand, the estimates of slope and intercept are reasonable in 
reality and in the simulations.

   I've skipped the gory details at this point.  A more complete write-up 
is at http://www.zoo.ufl.edu/bolker/tiwari-new-reg.pdf if anyone wants 
more information ...

   Should I be concerned about this?  Are there any obvious diagnostics of 
non-convergence etc.?  (Reported AIC values suggest that the models with 
autocorrelation *are* preferable, by quite a bit -- I could pursue this 
further.)  Or should I just not worry about it and move on?

   Ben Bolker


--------------------------
simulation code:

library(MASS)
simdata <- function(sd=200,range=2,n=NULL,mmin=0) {
   if (is.null(n)) mile <- seq(mmin,17.5,by=0.5) else {
     mile <- seq(mmin,17.5,length=n)
   }
   mean <- 3000-50*(mile-10)^2
   v <- sd^2
   dist <- abs(outer(mile,mile,"-"))
   Sigma <- v*exp(-dist/range)
   X <- mvrnorm(1,mu=mean,Sigma=Sigma)
   data.frame(mile=mile,X=X,mean=mean)
}

-- 
620B Bartram Hall                            bolker at zoo.ufl.edu
Zoology Department, University of Florida    http://www.zoo.ufl.edu/bolker
Box 118525                                   (ph)  352-392-5697
Gainesville, FL 32611-8525                   (fax) 352-392-3704




More information about the R-help mailing list