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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
