[R-sig-ME] setting corStruct range in a mixed effects model

[Shawn O'Neil]

Hello all,

I have a question regarding the use of the available correlation structures
in package nlme.  We are using mixed effects models to relate information
from lesser scaup kernel density distributions to underlying habitat
covariates.  Based on semivariogram and correlogram examinations, we find
spatial autocorrelation between the residuals of our full linear mixed
effects model.  We have tried to follow an example given by Pinheiro &
Bates (2000, pp 262-264) where a spherical correlation structure is fitted
to a sample semivariogram by imposing a user-defined range and nugget.  We
are trying this approach in an attempt to ignore some potentially spurious
autocorrelation effects occurring at large lag distances.  However, we
cannot achieve the result that the authors do in forcing a specific range.

For example,  we run a random intercept model using all of our covariates
including a surface trend:


FormXY <- formula(logUDval ~ DEP + I(DEP^2) + log(EDGE+0.5) + H2O:SUB +
log(CATTR+0.5) + NHOT +
    I(NHOT^2) + H2O*NDENS + I(NDENS^2) + Year_ + BCIndex + Age + X_ST*Y_ST
+ I(X_ST^2) + I(Y_ST^2) + I(X_ST^3)
    + I(Y_ST^3))

m11.lme<-lme(FormXY, random = ~1 | BirdID, method = "REML", data=hrdata)

Next, we examine the sample semivariogram (Figure 1) and determine that the
range should be approximately 900 and we do not see a significant nugget
effect.  At distances > 3500, there are some odd things occurring that for
now we would like to ignore.  We are primarily concerned with modeling the
spatial autocorrelation that is evident between 0 and 900 so we try to
specify this using the spherical correlation structure (as an example).  I
had to use lmeControl to get around false convergence errors:

m11.lme.spher<-lme(FormXY, random = ~1 | BirdID, correlation =
corSpher(c(900), form = ~X_COORD + Y_COORD | BirdID), method =
"REML",       data=hrdata,
control = lmeControl(opt = c("optim")))

The model finishes but we can see by the sample semivariogram (Figure 2)
and the sample semivariogram of normalized residuals (Figure 3) that the
correlation structure did not fit the way we wanted it to.  Also, the range
in the model summary does not match the value that we put in:

> summary(m11.lme.spher)
Linear mixed-effects model fit by REML
 Data: hrdata
        AIC       BIC   logLik
  -26755.58 -26562.75 13402.79

Random effects:
 Formula: ~1 | BirdID
        (Intercept)  Residual
StdDev:   0.5524643 0.3513228

Correlation Structure: Spherical spatial correlation
 Formula: ~X_COORD + Y_COORD | BirdID
 Parameter estimate(s):
   range
1950.937

Am I making a mistake somewhere?  We have tried all the other corStructs
and also tried setting the nugget in addition to the range but we never get
to a result that models the spatial autocorrelation trend accurately at
distances 0 - 900.  Any advice is greatly appreciated!

Thank you,
Shawn

-- 
Shawn O'Neil
Research Assistant, Applied Ecology
Forest Resources and Environmental Science
Michigan Technological University
701-741-4361