Lecture 42 (lab)—Friday, April 27, 2007

What was covered?

R functions and commands demonstrated

R function options

R packages used

Preliminaries

Checking for spatial correlation with a Mantel test

> library(nlme)
> data(Wheat2)
> names(Wheat2)

[1] "Block" "variety" "yield" "latitude" "longitude"

Fig. 1  Spatial arrangement of plots

> plot(latitude~longitude, col=as.numeric(Block), data=Wheat2)

#geographic distances
> dist1<-dist(cbind(Wheat2$latitude, Wheat2$longitude))
#differences in yield
> dist2<-dist(Wheat2$yield)

> library(vegan)
> mantel(dist1,dist2)

Mantel statistic based on Pearson's product-moment correlation

Call:
mantel(xdis = dist1, ydis = dist2)

Mantel statistic r: 0.2005
      Significance: < 0.001

Empirical upper confidence limits of r:
   90%    95%  97.5%    99%
0.0338 0.0416 0.0500 0.0596

Based on 1000 permutations

> model2<-lm(yield~variety+Block, data=Wheat2)
> dist3<-dist(residuals(model2))
> mantel(dist1,dist3)

Mantel statistic based on Pearson's product-moment correlation

Call:
mantel(xdis = dist1, ydis = dist3)

Mantel statistic r: 0.1061
      Significance: < 0.001

Empirical upper confidence limits of r:
   90%    95%  97.5%    99%
0.0302 0.0415 0.0472 0.0560

Based on 1000 permutations

> model4<-lm(yield~variety+Block+latitude+longitude, data=Wheat2)
> dist4<-dist(residuals(model4))
> mantel(dist1,dist4)

Mantel statistic based on Pearson's product-moment correlation

Call:
mantel(xdis = dist1, ydis = dist4)

Mantel statistic r: 0.04557
      Significance: 0.024

Empirical upper confidence limits of r:
   90%    95%  97.5%    99%
0.0313 0.0384 0.0444 0.0571

Based on 1000 permutations

Fitting the empirical semivariogram

> gmodel2<-gls(yield~variety+Block, data=Wheat2, method='ML')

> newdat<-data.frame(Wheat2$latitude, Wheat2$longitude, residuals(gmodel2, type='n'))

> library(geoR)
> geodat<-as.geodata(newdat)

> variog(geodat, uvec=1:30, max.dist=30, option='bin')->geodat.v1

> names(geodat.v1)
 [1] "u"                "v"              "n"
 [4] "sd"               "bins.lim"       "ind.bin"
 [7] "var.mark"         "beta.ols"       "output.type"
[10] "max.dist"         "estimator.type" "n.data"
[13] "lambda"           "trend"          "pairs.min"
[16] "nugget.tolerance" "direction"      "tolerance"
[19] "uvec"             "call"

> geodat.v1$n
 [1] 210  199  766 538 832 462  434 1406 695 633  589 1556 1000 813
[15] 628 1207 1200 837 731 474 1619  668 558 428 1131  547  465 196
[29] 165

> variog(geodat, uvec=1:50, max.dist=50, option='bin')->geodat.v1
variog: computing omnidirectional variogram
> geodat.v1$n
 [1] 210  199  766 538 832 462  434 1406 695 633  589 1556 1000 813
[15] 628 1207 1200 837 731 474 1619  668 558 428 1131  547  465 196
[29] 855  480  322 192 417 449  259  124 117 386  130  122   58 106
[43]  56   30   23  13   9   5

Fig. 2  Empirical semivariogram for residuals of yield model

> variog(geodat, uvec=1:50, max.dist=30, option='bin')-> geodat.v1
> plot(geodat.v1)

Modeling spatial correlation with the semivariogram

> ols2<-variofit(geodat.v1, ini=c(.8,30), nugget=.4, cov.model='exponential')
> lines(ols2,col=2)
> ols3<-variofit(geodat.v1, ini=c(.8,30), nugget=.4, cov.model='gaussian')
> lines(ols3,col=4)
> ols4<-variofit(geodat.v1, ini=c(.8,30), nugget=.4, cov.model='spherical')
> lines(ols4,col=3)
> ols5<-variofit(geodat.v1, ini=c(.8,30), nugget=.4, cov.model='linear')
> lines(ols5,col='gray70')
> legend(18,.5, c('exponential', 'Gaussian', 'spherical', 'linear'), col=c(2:4, 'gray70'), lty=rep(1,4), cex=.8, bty='n')

Model Formula
Exponential
Gaussian
Linear
Rational quadratic
Spherical

Here I( ) is the indicator function such that .

The approach in nlme is different than in geoR (where the nugget is an additive term rather than a multiplicative one).

#original model
> gmodel2<-gls(yield~variety+Block, data=Wheat2, method='ML')
> out1a<-update(gmodel2, corr=corSpher(c(20,.4), form=~latitude+longitude, nugget=T))

> out1b <-update(gmodel2, corr=corExp(c(20,.4), form=~latitude+longitude, nugget=T))
> out1c<-update(gmodel2, corr=corGaus(c(20,.4), form=~latitude+longitude, nugget=T))
> out1d<-update(gmodel2, corr=corLin(c(20,.4), form=~latitude+longitude, nugget=T))
> out1e<-update(gmodel2, corr=corRatio(c(20,.4), form=~latitude+longitude, nugget=T))

> sapply(list(gmodel2,out1a,out1b,out1c,out1e),AIC)
[1] 1561.622 1370.701 1372.923 1369.888 1369.571

> gmodel4<-update(gmodel2, .~.+latitude+longitude)
> AIC(gmodel4)

[1] 1459.622

# plot variogram for the residuals original model
> plot(Variogram(gmodel2, form=~latitude+longitude, maxDist=30, resType='n'))
# plot variogram for residuals of final model
> plot(Variogram(out1e, form=~latitude+longitude, resType='n'))

(a)
(b)
Fig. 4 (a) Semivariogram for residuals from an ordinary regression model in which the residuals are assumed to be uncorrelated (b) Semivariogram for residuals from a regression model in which the residuals are given a spatial correlation structure

> newdat2<-data.frame(Wheat2$latitude, Wheat2$longitude, residuals(out1e, type='n'))
> geodat<-as.geodata(newdat2)
> variog(geodat,uvec=1:30,max.dist=30,option='bin')->geodat.v2

> plot(geodat.v2)

Fig. 5  Semivariogram for residuals for model with rational quadratic correlation structure (geoR)

Cited reference

Stroup, W. W., P. S. Baenziger, and D. K. Mulitze. 1994. Removing spatial variation from wheat yield trials: a comparison of methods. Crop Science 86: 62–66.

Course Home Page


Jack Weiss
Phone: (919) 962-5930
E-Mail: jack_weiss@unc.edu
Address: Curriculum in Ecology, Box 3275, University of North Carolina, Chapel Hill, 27516
Copyright © 2007
Last Revised--May 1, 2007
URL: http://www.unc.edu/courses/2007spring/enst/562/001/docs/lectures/lecture42.htm