## ----getData------------------------------------------------------------------

if (!require("rspat")) remotes::install_github("rspatial/rspat")

library(rspat)
h <- spat_data('houses2000')


## -----------------------------------------------------------------------------
dim(h)
names(h)


## ----echo=FALSE---------------------------------------------------------------
dv <- data.frame(matrix(c("nhousingUn", "number of housing units", "recHouses", "number of houses for recreational use", "nMobileHom", "number of mobile homes", "nBadPlumbi", "number of houses with incomplete plumbing", "nBadKitche", "number of houses with incomplete kitchens", "Population", "total population", "Males", "number of males", "Females", "number of females", "Under5", "number of persons under five", "White", "number of persons identifying themselves as white (only)", "Black", "number of persons identifying themselves African-american (only)", "AmericanIn", "number of persons identifying themselves American Indian (only)", "Asian", "number of persons identifying themselves as American Indian (only)", "Hispanic", "number of persons identifying themselves as hispanic (only)", "PopInHouse", "number of persons living in households", "nHousehold", "number of households", "Families", "number of families", "houseValue", "value of the house", "yearBuilt", "year house was built", "nRooms", "median number of rooms per house", "nBedrooms", "median number of bedrooms per house", "medHHinc", "median household income", "MedianAge", "median age of population", "householdS", "median household size", "familySize", "median family size"), ncol=2, byrow=TRUE))
colnames(dv) <- c("variable", "description")
knitr::kable(dv)


## -----------------------------------------------------------------------------
# using a tiny buffer to get a cleaner aggregation 
hb <- buffer(h, 1)
values(hb) <- values(h)
hha <- aggregate(hb, "County")


## -----------------------------------------------------------------------------
d1 <- as.data.frame(h)[, c("nhousingUn", "recHouses", "nMobileHom", "nBadPlumbi", 
 "nBadKitche", "Population", "Males", "Females", "Under5", "White", 
 "Black", "AmericanIn", "Asian", "Hispanic", "PopInHouse", "nHousehold", "Families")]

 d1a <- aggregate(d1, list(County=h$County), sum, na.rm=TRUE)


## -----------------------------------------------------------------------------
d2 <- as.data.frame(h)[, c("houseValue", "yearBuilt", "nRooms", "nBedrooms", 
       "medHHinc", "MedianAge", "householdS",  "familySize")]
d2 <- cbind(d2 * h$nHousehold, hh=h$nHousehold)

d2a <- aggregate(d2, list(County=h$County), sum, na.rm=TRUE)
d2a[, 2:ncol(d2a)] <- d2a[, 2:ncol(d2a)] / d2a$hh


## -----------------------------------------------------------------------------
d12 <- merge(d1a, d2a, by='County')


## -----------------------------------------------------------------------------
hh <- merge(hha[, "County"], d12, by='County')


## ----spreg2-------------------------------------------------------------------
library(RColorBrewer)
grps <- 10
brks <- quantile(h$houseValue, 0:(grps-1)/(grps-1), na.rm=TRUE)
plot(h, "houseValue", breaks=brks, col=rev(brewer.pal(grps, "RdBu")), border=NA)
lines(hh, col="white")


## ----spreg4-------------------------------------------------------------------
brks <- quantile(h$medHHinc, 0:(grps-1)/(grps-1), na.rm=TRUE)
plot(h, "medHHinc", breaks=brks, col=rev(brewer.pal(grps, "RdBu")), border=NA)
lines(hh, col="white")


## -----------------------------------------------------------------------------
hh$fBadP <- pmax(hh$nBadPlumbi, hh$nBadKitche) / hh$nhousingUn
hh$fWhite <- hh$White / hh$Population
hh$age <- 2000 - hh$yearBuilt

f1 <- houseValue ~ age +  nBedrooms 
m1 <- lm(f1, data=as.data.frame(hh))
summary(m1)


## -----------------------------------------------------------------------------
y <- matrix(hh$houseValue)
X <- cbind(1, hh$age, hh$nBedrooms)


## -----------------------------------------------------------------------------
ols <- solve(t(X) %*% X) %*% t(X) %*% y
rownames(ols) <- c('intercept', 'age', 'nBedroom')
ols


## ----spreg6-------------------------------------------------------------------
hh$residuals <- residuals(m1)
brks <- quantile(hh$residuals, 0:(grps-1)/(grps-1), na.rm=TRUE)
plot(hh, "residuals", breaks=brks, col=rev(brewer.pal(grps, "RdBu")))


## ----spreg8, message=FALSE----------------------------------------------------
library(spdep)

sfhh <- sf::st_as_sf(hh)
nb <- poly2nb(sfhh, snap=1/120)
nb[[21]] <- sort(as.integer(c(nb[[21]], 38)))
nb[[38]] <- sort(as.integer(c(21, nb[[38]])))
nb

par(mai=c(0,0,0,0))
plot(hh)
plot(nb, crds(centroids(hh)), col='red', lwd=2, add=TRUE)


## ----spreg10------------------------------------------------------------------
resnb <- sapply(nb, function(x) mean(hh$residuals[x]))
cor(hh$residuals, resnb)
plot(hh$residuals, resnb, xlab="Residuals", ylab="Mean adjacent residuals", pch=20)
abline(lm(resnb ~ hh$residuals), lwd=2, lty=2)


## -----------------------------------------------------------------------------
lw <- nb2listw(nb)
moran.mc(hh$residuals, lw, 999)


## ----spreg, message=FALSE-----------------------------------------------------
library(spatialreg )


## ----spregplot1---------------------------------------------------------------
m1s <- lagsarlm(f1, data=as.data.frame(hh), lw, tol.solve=1.0e-30)

summary(m1s)

hh$residuals <- residuals(m1s)
moran.mc(hh$residuals, lw, 999)

brks <- quantile(hh$residuals, 0:(grps-1)/(grps-1), na.rm=TRUE)
plot(hh, "residuals", breaks=brks, col=rev(brewer.pal(grps, "RdBu")))



## ----spregplotx---------------------------------------------------------------
m1e <- errorsarlm(f1, data=as.data.frame(hh), lw, tol.solve=1.0e-30)
summary(m1e)

hh$residuals <- residuals(m1e)
moran.mc(hh$residuals, lw, 999)

brks <- quantile(hh$residuals, 0:(grps-1)/(grps-1), na.rm=TRUE)
plot(hh, "residuals", breaks=brks, col=rev(brewer.pal(grps, "RdBu")))


## ----spregplot3---------------------------------------------------------------
brks <- quantile(hh$residuals, 0:(grps-1)/(grps-1), na.rm=TRUE)
plot(hh, "residuals", breaks=brks, col=rev(brewer.pal(grps, "RdBu")))