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:
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]]
[[alternative HTML version deleted]]