[R-sig-eco] Create new corStruct class for lme model

[Shawn O'Neil]

Thanks for the helpful resources.  I have a related question regarding the
use of the available correlation structures in nlme.  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 have tried this approach in
our models 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!

I will also post this to the mixed models list as it may be more relevant
to that forum...

Thank you,
Shawn O.

On Sat, Jan 21, 2012 at 11:40 AM, Ben Bolker <bbolker at gmail.com> wrote:

> Shawn O'Neil <oneil.shawnt at ...> writes:
>
> > I am running mixed effects models using the nlme package.  I am working
> > with a dataset where spatial autocorrelation is inherent in the response
> > variable.  I have tested each of the five correlation structures that are
> > available for mixed models in nlme, but none of these functions fit my
> > spatial correlation pattern very well.  To quote Pinheiro and Bates
> (2000),
> > "New corStruct classes, representing user-defined correlation structures,
> > can be added to the set of standard classes and used with the modeling
> > functions in the nlme library. For this, one must specify a constructor
> > function, generally with the same name as the class, and, at a minimum,
> > methods for the functions coef, corMatrix, and initialize."  I am fairly
> > new to R and could benefit from an example showing how this is done.
>  I've
> > spent a couple of days searching the net and haven't come up with
> > anything.  Can anyone refer me to a good source where a custom corStruct
> > class is created and run with gls or lme?
>
>   The best examples I know of are in the ape package, for doing
> phylogenetic generalized least-squares problems (PGLS) -- I believe
> corBrownian() and corMartins() are in there, there may be others.
>
>  I also started to build some code along these lines for
> 'antecedent' models, as requested by a poster on the r-sig-mixed-models
> list.  I never finished, but what I did get done might be useful:
> http://www.math.mcmaster.ca/~bolker/R/misc/newcorstruct.R
>
>
> http://www.math.mcmaster.ca/~bolker/R/misc/newcorstruct.R
>
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>



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