[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:
[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))
More information about the R-sig-mixed-models
mailing list