Dear R Help,
I am having difficulty using Variogram within GLS to examine spatial structure
of nested data. 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:
site landposi replicat distanc canopy litdepth
soilmoist shrubcov soildepth
PionB bottom 1 0.5 2 14
0 20 1.82
PionB bottom 1 1 2 20
0 20 4.17
PionB bottom 1 1.5 3 26
0 20 2.6
PionB bottom 1 2 3 23
0 18 2.86
PionB bottom 1 2.5 1 14
0 16.2 2.34
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)
> attach(smith)
> trancano=(asin(sqrt(canopy/100)))
> replicat=factor(replicat)
> library(nlme)
Thus, the data are entered and transformed, and the appropriate library
summoned.
> 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)
> summary(model)
Generalized least squares fit by REML
Model: trancano ~ landposi
Data: structure
AIC BIC logLik
-4968.469 -4947.224 2488.235
Coefficients:
Value Std.Error
t-value p-value
(Intercept) 0.17966306 0.002040385 88.05352
0.0000
landposiridge 0.00095435 0.002885540 0.33073
0.7409
landposislope 0.01226755 0.002885540 4.25139
0.0000
Correlation:
(Intr) lndpsr
Landposiridge -0.707
landposislope -0.707 0.500
Standardized residuals:
Min Q1 Med
Q3 Max
-3.9587908 -0.7265938 -0.0819489 0.7360712
3.8154715
Residual standard error: 0.04562439
Degrees of freedom: 1500 total; 1497 residual
> plot(Variogram(model, form=~distanc) )
Error in as.data.frame.default(data, optional = TRUE) :
cannot coerce class '"function"' into a
data.frame
And here it stops. For some reason, R seems to think there is a function in the
grouped data object. Despite a considerable amount of time spent in debugging,
I cannot discover the reason. I notice that Christina Silva described the same
problem in the R-Help log a few years ago (11/19/2002), but her question
doesn't seem to have been answered.
I greatly appreciate any advice or suggestions anyone can give me on this!
Many thanks,
Glenn Matlack
Athens, Ohio
[[alternative HTML version deleted]]
