[R-sig-ME] Mixed model including spatially autocorrelated response

Ben Bolker bbolker at gmail.com
Sat Dec 8 17:19:33 CET 2012

Matlack, Glenn <matlack at ...> writes:

> Dear R Sig Mixed,

> I am having difficulty using Variogram within GLS to examine spatial
> structure within a mixed model.  My data frame consists of
> ecological measurements of a forest, in which three landscape
> positions ("landposi") are compared.  Each landscape position is
> replicated five times ("replicat"), and the environment is measured
> ("canopy", "litdepth", etc.) one hundred times at 50-cm intervals
> along a transect ("distanc").  Thus, there are 1,500 rows of data.
> The raw data file begins:

> # The general strategy follows Pinheiro and Bates' (2000;
> pp. 260-266) and Crawley's (2007; pp. 778-785) wheat-trials example:
> 1) groupedData is used to specify the nesting structure, 2) a GLS
> model is fit assuming no spatial autocorrelation, 3) spatial
> covariance is added using the Variogram function, and 4) # # #
> spatial correlation structures are compared.  Before beginning,
> "canopy" (a percentage) is arcsin sqareroot transformed to
> "trancano", and "replicat" is specified as a factor.
> smith<-read.table("C:\\Users\\matlack\\Desktop\\
>  Documents\\Undergrad Research\\Nicole Smith\\smith.txt",header=T)

  Right now this isn't reproducible.  Is there any chance you
can post the data somewhere, or extract a subset?
http://tinyurl.com/reproducible-000 is useful.

> attach(smith)

   I have a strong preference against attach(); it is generally
more coherent to add new variables or transform variables *within*
the data frame, as in

smith2 <- transform(smith,

> > library(nlme)


> structure<-groupedData(trancano~landposi|landposi/replicat)

> "Trancano" is the response variable, "landposi" is the factor of
> interest, and "replicat" is a random effect nested within
> "landposi".
> > model<-gls(trancano~landposi,structure)

I have a couple of comments here.  I generally don't use
groupedData because I don't understand them; I prefer to 
specify random effects structures explicitly in the model.
I'm also a bit confused, because as far as I know gls()
doesn't handle random effects??  But maybe gls() just
ignores groupedData information?

  Also, it doesn't make sense to fit landposi as both
a fixed effect and a random effect.  You may want
landposi:replicat rather than landposi/replicat as your
grouping factor ...

> > summary(model)
> Generalized least squares fit by REML
>                 Model: trancano ~ landposi
>                 Data: structure
>              AIC                  BIC                         logLik
>            -4968.469             -4947.224             2488.235
> Coefficients:


Here's my reproducible example:

d <- expand.grid(landposi=c("bottom","mid","top"),
d$z <- runif(nrow(d))
d <- transform(d,posrep=interaction(landposi,replicat))
m1 <- lme(z~landposi,random=~1|posrep,data=d)

More information about the R-sig-mixed-models mailing list