I'm trying to specify the covariance structure parameters in a linear mixed model (using the correlation structure facilities in nlme). The parameters are estimated from a gstat or geoR variogram model fit to the empirical semivariogram.
My justification for specifying the gstat-derived covariance model in nlme is because nlme uses only the classical semivariance estimator fit to the full data set. gstat and geoR permit the fit of variogram models to the robust (Cressie) semivariogram and the setting of a maximum distance of influence. For example, see my question seeking guidance on the selection between the classical vs. robust semivariogram.
I can produce nearly identical empirical semivariograms to the residuals from an intercept-only model (first figure below). Likewise, the exponential variogram model fits from gstat and geoR are essentially identical. However, fitting the exponential variogram model in nlme produces bizarre parameter estimates.
nlme uses a slightly different specification of the variogram models and correlation structures compared to gstat and geoR: namely, nlme standardizes the within-group errors to unit variance and uses a multiplicative rather than additive nugget effect (see, e.g., p. 230-232 in Pinheiro and Bates, 2000, Mixed-Effects Models in S and S-Plus). Nonetheless, in the absence of a nugget effect (the example below) the exponential fits should be easily translatable among the programs.
Here is a simple illustration of what is perplexing me. This example uses the data set of heavy metal concentrations along the river Meuse in The Netherlands (?meuse in gstat package); the autocorrelation is spatial (heavy metal concentrations closer to each other in space are more similar on average than concentrations farther apart).
The empirical variograms of the model residuals using nlme, gstat and geoR are essentially identical.
# Load necessary packageslibrary(gstat) # variogram model fittinglibrary(sp) # classes and methods for spatial datalibrary(geoR) # variogram model fittinglibrary(nlme) # linear mixed models with autocorrelationlibrary(ggplot2) #graphics# Use heavy metal data from `gstat`# Focus on zinc concentrationsdata(meuse)meuse$dGroup
My justification for specifying the gstat-derived covariance model in nlme is because nlme uses only the classical semivariance estimator fit to the full data set. gstat and geoR permit the fit of variogram models to the robust (Cressie) semivariogram and the setting of a maximum distance of influence. For example, see my question seeking guidance on the selection between the classical vs. robust semivariogram.
I can produce nearly identical empirical semivariograms to the residuals from an intercept-only model (first figure below). Likewise, the exponential variogram model fits from gstat and geoR are essentially identical. However, fitting the exponential variogram model in nlme produces bizarre parameter estimates.
nlme uses a slightly different specification of the variogram models and correlation structures compared to gstat and geoR: namely, nlme standardizes the within-group errors to unit variance and uses a multiplicative rather than additive nugget effect (see, e.g., p. 230-232 in Pinheiro and Bates, 2000, Mixed-Effects Models in S and S-Plus). Nonetheless, in the absence of a nugget effect (the example below) the exponential fits should be easily translatable among the programs.
Here is a simple illustration of what is perplexing me. This example uses the data set of heavy metal concentrations along the river Meuse in The Netherlands (?meuse in gstat package); the autocorrelation is spatial (heavy metal concentrations closer to each other in space are more similar on average than concentrations farther apart).
The empirical variograms of the model residuals using nlme, gstat and geoR are essentially identical.
# Load necessary packageslibrary(gstat) # variogram model fittinglibrary(sp) # classes and methods for spatial datalibrary(geoR) # variogram model fittinglibrary(nlme) # linear mixed models with autocorrelationlibrary(ggplot2) #graphics# Use heavy metal data from `gstat`# Focus on zinc concentrationsdata(meuse)meuse$dGroup