[R-sig-Geo] how to fit corSpher with gls?

[Andy Bunn]

Hello all, I'm playing around with using gls and appropriate corr
structure to model spatial autocorrelation in the residuals of a
regression. Below I have an example of a regression of y~x where the
residuals are spatially autocorrelated. I then try to fit a gls() model
with a spherical correlation structure using parameters from a fitted
variogram to the residuals of y~x. My expectation is that I would find the
std error on my estimate of x to decrease and AIC to decrease with the
spatial structure modeled correctly. But I think I am not specifying the
terms correctly on the call to corSpher().

Advice appreciated -Andy

# ~~~~~~~~~~~~~~~
# goal:
# create a variable y that is a function of x and spatially
# autocorrelated noise (espsilon)
# model y~x and show that the residuals are not iid
# apply a gls model with appropriate corr structure

rm(list=ls())
library(gstat)
library(sp)
library(nlme)
library(spdep)
library(ncf)

set.seed(234)
# generate the spatially autocorrelated error term, epsilon
n <- 200
epsilon <- rnorm(n,0,1)

easting <- runif(n,0,100)
northing <- runif(n,0,100)
points <- cbind(easting,northing)
dnb <- dnearneigh(points, 0, 150)
dsts <- nbdists(dnb, points)
idw <- lapply(dsts, function(x) 1/x^1.5) # p=1.5
lw <- nb2listw(dnb, glist=idw, style="W")

rho <- 0.95
inv <- invIrW(lw, rho)

epsilon <- inv %*% epsilon
epsilon <- scale(epsilon[,1])


# generate x as noise.
x <- rnorm(n)
# generate y as a function of x and epsilon (which is spatially
autocorrelated)
B0 <- 10
B1 <- 2
y <- B0 + B1*x + epsilon

# Create a spatial object
dat <- data.frame(easting,northing,y,x)
coordinates(dat)<-c('easting','northing')

# Fit a model (we don't know about spatial component)
lm1 <- lm(y~x, dat)
coefficients(summary(lm1))

# Add residuals to the SPDF dat
dat$lmResids <- residuals(lm1)

# Map the rediduals
bubble(dat, zcol = 'lmResids')

# Test them for spatial autocorrelation using moran.test,
# a correlogram and a variogram.

# Moran's I test
w <- knn2nb(knearneigh(dat, k=8))
moran.test(dat$lmResids, nb2listw(w))

# Moran's I correlogram
residsI <- spline.correlog(x=coordinates(dat)[,1], y=coordinates(dat)[,2],
                         z=dat$lmResids, resamp=20)
plot(residsI)

# Variogram
lmResidsVar <- variogram(lmResids~1, data=dat)
plot(lmResidsVar)
sph.model <- vgm(psill=1, model="Sph", range=20, nugget=0.3)
sph.fit <- fit.variogram(object = lmResidsVar, model = sph.model)
plot(lmResidsVar,model=sph.fit)

# Refit a model with an approparite correlation structure.
gls1 <- gls(y~x,data=dat) # same as ols above (lm1)

# How to get the right corSpher structure?
# I think I should give it values from the variogram fit above
# but that appears to be incorrect!
cs1 <- corSpher(value=c(17,0.3),nugget=TRUE)
gls2 <- gls(y~x,data=dat,correlation=cs1)

summary(gls1)
summary(gls2)





>