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

[Andy Bunn]

One more typo. Here is the complete script:

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.
# Corr structure with range and nugget from above and specifying
# the form
cs1 <- corSpher(c(17,0.3),form=~easting+northing,nugget=TRUE)
gls1 <- gls(y~x,data=dat,correlation=cs1)

summary(gls1)

# add the residuals to the spatial dat object
dat$glsResids <- residuals(gls1,type="normalized")

# Map the rediduals from gls
bubble(dat, zcol = 'glsResids')

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






On 5/18/15, 3:03 PM, "Andy Bunn" <Andy.Bunn at wwu.edu> wrote:

>Ugh! Forgot to specify the form of the object. Fixed below.
>
>On 5/18/15, 12:57 PM, "Andy Bunn" <Andy.Bunn at wwu.edu> wrote:
>
>>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.
># Corr structure with range and nugget from above and specifying
># the form
>cs1 <- corSpher(c(17,0.3),form=~easting+northing,nugget=TRUE)
>gls1 <- gls(y~x,data=dat,correlation=cs1)
>summary(gls1)
># add the residuals to the spatial dat object
>dat$glsResids <- residuals(gls2,type="normalized")
># Moran's I test comes back clean
>moran.test(dat$glsResids, nb2listw(w))
>
>
>
>
>
>
>
>
>
>
>
>>
>># 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)
>>
>>
>>
>>
>>
>>>
>>
>>_______________________________________________
>>R-sig-Geo mailing list
>>R-sig-Geo at r-project.org
>>https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>