# Spatial Regression tutorial

library(sp)
library(gstat)
library(nlme)
library(RColorBrewer)

par(ask=T)

scatter.color <- function(X,Y,Z,brewer.name,breaks,legend=F,...){
# scatter.color(X,Y,Z,brewer.name,breaks,legend,…)
# Input:
#   X: X coordinate
#   Y: Y coordinate
#   Z: value for coloring
#   brewer.name: string for color palette (e.g. “Blues”); see below
#   breaks: cutpoint for colors (e.g. quantile(Z,c(0,.33,.66,1)))
#   legend: add legend to plot? default is False.
#    if true, then click on plot to add a legend using locator()
#   …: extra parameters to be passed to plot(X,Y,…)
# Output:
#  $colors: hexadecimal colors for every element in Z
#  $levels: numeric ranges of levels
#  $palette: hexadecimal colors of palette
# Palette Names:
# The sequential palettes names are
# Blues BuGn BuPu GnBu Greens Greys Oranges OrRd PuBu PuBuGn PuRd Purples RdPu Reds YlGn YlGnBu YlOrBr YlOrRd
# The diverging palettes are
# BrBG PiYG PRGn PuOr RdBu RdGy RdYlBu RdYlGn Spectral
num.cuts <- length(breaks)-1
palette <- brewer.pal(num.cuts,brewer.name)
if (num.cuts==2){palette <- c(palette[1],palette[3])}
Z.level <- cut(Z,breaks)
# Should probably base palette on actual number of breaks!!!
colors <- palette[as.numeric(Z.level)]
levels <-levels(Z.level)
plot(X,Y,pch=16,type="p",col=colors,asp=1,...)
points(X,Y)
if (legend) {
legend(locator(),levels,fill=palette)
}
return(list(colors=colors,levels=levels,palette=palette))
}


data(meuse)
data(meuse.grid)
gridded(meuse.grid) <- ~x+y
summary(meuse)
gr
# Plot lead concentration
scatter.color(meuse$x,meuse$y,log(meuse$lead),'YlGnBu',
	quantile(log(meuse$lead),c(0,.2,.4,.6,.8,1)),legend=T,
	main='Lead Concentration')
#Plot distance
spplot(meuse.grid,'dist', main = 'Distance to River')


# Regress lead on distance to the river

lead.OLS <- lm(sqrt(lead)~sqrt(dist),data=meuse)
summary(lead.OLS)
#distance explains about 40% of the variation. 
# Essentially, this means that we are asking 40% less of 
#  kriging than we were before

# Add the OLS residuals to the dataframe
meuse$OLS.resid <- residuals(lead.OLS)

# Plot the OLS residuals
scatter.color(meuse$x,meuse$y,meuse$OLS.resid,'YlGnBu',quantile(meuse$OLS.resid,c(0,.2,.4,.6,.8,1)),legend=T)

# Calculate and plot the semivariogram of residuals
OLS.vgram.anis <- variogram(OLS.resid~1,~x+y,data=meuse,alpha=c(0,45,90,135))
plot(OLS.vgram.anis)
# Note, the covariate has not solved the anisotropy problems
# Also, it is difficult to determine the direction of max continuity


# It appears that it might be 90 rather than 45, but I have trouble
# justifying this scientifically
# For simplicity, I will fit an isotropic model
fit.eye <- vgm(7,'Sph',900,3)
plot(OLS.vgram.anis,fit.eye)

# Fit spherical model by weighted least squares (WLS)
OLS.vgram.iso <- variogram(OLS.resid~1,~x+y,data=meuse)

OLS.fit.sd <- fit.variogram(OLS.vgram.iso,fit.eye,fit.method=2)
plot(OLS.vgram.iso,OLS.fit.sd)

# We will use the gls function in the nlme library.  
# This function requires a new object that specifies the corr structure
OLS.fit.sd
OLS.cor <- corSpher(c(982,.375),form=~x+y,nugget=T,fixed=T)
# This is a spherical model with range = 982 and nugget = 37.5%
lead.GLS1 <- gls(sqrt(lead)~sqrt(dist),data=meuse,OLS.cor)

meuse$GLS1.resid <- residuals(lead.GLS1)
GLS1.vgram.iso <- variogram(GLS1.resid~1,~x+y,data=meuse)
GLS1.fit.sd <- fit.variogram(GLS1.vgram.iso,OLS.fit.sd,fit.method=2)
plot(GLS1.vgram.iso,GLS1.fit.sd)
GLS1.fit.sd
GLS1.cor <- corSpher(c(800,.318),form=~x+y,nugget=T,fixed=T)
lead.GLS2 <- gls(sqrt(lead)~sqrt(dist),data=meuse,GLS1.cor)
meuse$GLS2.resid <- residuals(lead.GLS2)

# compare the coefficient values
summary(lead.GLS2)
summary(lead.GLS1)
summary(lead.OLS)

#The coefficients are still changing (but relatively little compared to their standard error),
#Also compare the t-statistics of GLS with OLS
#OLS t-statistics are almost always too large


#############################################3
# Regression kriging

# load the prediction raster
data(meuse.grid)

# Plot the trend component
beta <- coefficients(lead.GLS2)
trend <- beta[1] + beta[2]*meuse.grid$dist
pred <- data.frame(cbind(x=as.data.frame(meuse.grid)$x,
				y=as.data.frame(meuse.grid)$y,
				trend=trend))
gridded(pred) <- ~x+y
spplot(pred,'trend',main='Trend Component')

# Calculate SK of errors
SK <- krige(formula=GLS2.resid~1, locations=~x+y,
		model=GLS1.fit.sd, beta=0,
		data=as.data.frame(meuse),
		newdata=as.data.frame(meuse.grid)
)

pred$error <- SK$var1.pred
spplot(pred,'error',main = 'SK of errors')

pred$RK <- pred$trend + pred$error

spplot(pred,'RK',main = 'Regression Kriging Prediction')
############