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

[Andy Bunn]

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