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

[Ben Bolker]

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:
 

[snip]
 
> # 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,
    trancano=asin(sqrt(canopy/100)),
    replicat=factor(replicat))


> > library(nlme)

[snip]

> 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:

[snip]

Here's my reproducible example:

library(nlme)
set.seed(101)
d <- expand.grid(landposi=c("bottom","mid","top"),
                 replicat=factor(1:5),
                 distanc=1:100)
d$z <- runif(nrow(d))
d <- transform(d,posrep=interaction(landposi,replicat))
m1 <- lme(z~landposi,random=~1|posrep,data=d)
plot(Variogram(m1,form=~distanc|posrep))