SSA Poverty - Spatial Regression

This notebook is to replicate Sara’s district-wise poverty regressions in R and explore alternate approaches to defining spatial weights across districts (testing for spatial autocorrelation in the series of predictors), and to map GWR coefficients. Sara’s original code is taken from 3 Regression analysis district rural poverty - shock reg.do.

library(raster)
library(foreign)
library(data.table)
library(spdep)
library(splm)
library(spgwr)
library(tmap)
library(stringr)

# Presentation options
library(pander)
panderOptions("big.mark", ",")
panderOptions("round", 2)
panderOptions("missing", ".")
panderOptions("table.split.table", 100)

# Load R workspace on BUSTER
setwd("~/Projects/hc-shiny/geospat")
load("./tmp/poverty_r16.05.RData")
#setwd("~/Dropbox (IFPRI)/SDA/Data/analysis/_global_codes")
#load("./temp/2016.05/poverty_r16.05.RData")

# Helper
"%||%" <- function(a, b) if (!is.null(a)) a else b
# Helper - Collect important coefficients from `splm` results
# Could also check if there's an existing broom method for `splm` and `grw` models...
splm.combine <- function(x) data.table(
  model=names(x), 
  do.call(rbind, lapply(x, coef)), 
  do.call(rbind, lapply(x, `[[`, "arcoef")),
  phi=as.numeric(do.call(rbind, lapply(x, function(x) try(x$errcomp[["phi"]])))),
  psi=as.numeric(do.call(rbind, lapply(x, function(x) try(x$errcomp[["psi"]])))),
  rho=as.numeric(do.call(rbind, lapply(x, function(x) try(x$errcomp[["rho"]])))))
# Helper - AIC function for `spml`, show goodness of fit measure
# Credits to https://stat.ethz.ch/pipermail/r-sig-geo/2016-February/024077.html
spml.godf<-function(object, k=2, criterion=c("AIC", "BIC"),  ...) {
  s <- summary(object)
  l <- s$logLik[1,1]
  np <- length(coef(s))
  N <- nrow(s$model)
  if (criterion=="AIC") {
    aic <- -2*l+k*np
    names(aic) <-"AIC"
    return(aic)
  }
  if (criterion=="BIC") {
    bic <- -2*l+log(N)*np
    names(bic) <-"BIC"
    if (k!=2) {
      warning("parameter <k> not used for BIC")
    }
    return(bic)
  }
}

Data

# Load districts and attributes
g2 <- shapefile("./maps/svyMaps_2016.06.22_sara.shp")
dt2 <- read.dta("./tmp/SSApoverty_Dist_forGWR.12.dta")
#g2 <- shapefile("./out/r16.05/svyMaps_2016.06.22_sara.shp")
#dt2 <- read.dta("./temp/2016.05/SSApoverty_Dist_forGWR.12.dta")
# Keep STATA labels for re-use
dt2.lbl <- data.table(varCode=names(dt2), varLabel=attr(dt2, "var.labels"))
setkey(dt2.lbl, varCode)
dt2.lbl[is.na(varLabel), varLabel := varCode]
dt2 <- data.table(dt2)

# Verify that we have all explanatory variables and shocks already constructed, e.g.
rbind(
  dt2.lbl[varCode %like% "ndvi"],
  dt2.lbl[varCode %like% "shock"],
  dt2.lbl[varCode %like% "AEZ"],
  dt2.lbl[varCode %like% "land"])

# Make unique shape IDs explicit
g2$rn <- row.names(g2)
g2.dt <- data.table(g2@data)
# Recode Ethiopia woredas
dt2[svyCode=="eth2010", svyL2Cd := svyL1Cd * 10000 + svyL2Cd]
# Merge shapes and attributes
setkey(g2.dt, ISO3, svyCode, svyL1Cd, svyL2Cd)
setkey(dt2, ISO3, svyCode, svyL1Cd, svyL2Cd)
# Look for possible duplicates
g2.dt[duplicated(g2.dt), .(ISO3, svyCode, svyL1Cd, svyL2Cd)]

dt2[duplicated(dt2), .(ISO3, svyCode, svyL1Cd, svyL2Cd)]

# Drop duplicated vars from Sara's file before merging
dt2[, `:=`(
  rn=NULL, svyL1Nm=NULL, svyL2Nm=NULL, prttyNm=NULL, areakm=NULL, X=NULL, Y=NULL)]
# Any unmatched obs?
dt2[!g2.dt, .N, by=svyCode]

# Seems okay, so let's merge
g2.dt <- dt2[g2.dt]
# Re-attach Sara's attributes to shapes
g2 <- SpatialPolygonsDataFrame(g2, data.frame(g2.dt), match.ID="rn")

# Visually check a few vars
data(World)
p0 <- tm_shape(World) + tm_fill("grey90") + tm_shape(g2, is.master=T)
p1 <- tm_shape(World) + tm_borders("grey50", lwd=.2)
p0 + tm_fill("ndvi_ave", title="NDVI\nLong-term mean",
  style="pretty", n=9, palette=rev(pal.nrwc), auto=F) + p1

p0 + tm_fill("pcexp_ppp_m", title="Per Capita\nExpenditure\n(PPP USD/month)",
  style="jenks", n=9, auto=F) + p1

p0 + tm_fill("foodexp_ppp_m", title="Food Expenditure\n(PPP USD/month)",
  style="jenks", n=9, auto=F) + p1

p0 + tm_fill("pre_lt", title="Precipitation\nlong-term mean\n(mm)",
  style="jenks", n=9, palette=pal.pre, auto=F) + p1

p0 + tm_fill("spei_lt", title="SPEI\nlong-term mean",
  style="jenks", n=9, palette="RdYlGn", auto=F) + p1

p0 + tm_fill("nighlight_yearly", title="Night Light\nReflectance",
  style="kmeans", n=9, palette=pal.earth[1:230], auto=F) + p1

p0 + tm_fill("tt10_20k", title="Travel Time\nto nearest 20K market\n(hours)",
  style="fixed", breaks=c(.5,1,2,4,6,8,10,12,24), palette="-YlOrRd", auto=F) + p1

Model Specifications

Following Sara’s setup, 2 outcome variables Y, 9 alternate groups of biophysical and climatic shock variables and 4 cross-combinations of predictors are considered, for a total of 72 models. Sara chose a SAC form in STATA (equivalent to sacsarlm() in R).

# Define outcome vars
Y <- c("pcexp_ppp_m", "foodexp_ppp_m")
# Define groups of explanatory vars
X <- list(
  socio    = c("hh_female", "agehead", "marriedhead", "meaneduc", "depratio2", "children"),
  agri     = c("LGP_AVG", "ELEVATION", "cropland_pcap", "TLU_pcap"),
  access   = c("nighlight_yearly", "tt10_20k", "PD12_TOT"),
  weather1 = c("spei_lt", "L1_speidif", "L1_speidif2"),
  weather2 = c("spei_lt", "Lcum_speidif", "Lcum_speidif2"), 
  weather3 = c("spei_lt", "L1_speihishock", "L1_speiloshock"),
  weather4 = c("pre_lt", "L1_predif", "L1_predif2", "temp_lt", "L1_tempdif", "L1_tempdif2"),
  weather5 = c("pre_lt", "Lcum_predif", "Lcum_predif2", "temp_lt", "Lcum_tempdif", "Lcum_tempdif2"),
  weather6 = c("pre_lt", "temp_lt", "L1_prehishock", "L1_preloshock", "L1_temphishock"),
  weather7 = c("ndvi_lt", "L1_ndvidif", "L1_ndvidif2"), 
  weather8 = c("ndvi_lt", "Lcum_ndvidif", "Lcum_ndvidif2"), 
  weather9 = c("ndvi_lt", "L1_ndvihishock", "L1_ndviloshock"),
  other    = c("Malaria_yearly"),
  AEZ      = c("arid_warm", "humid_cool", "humid_warm"),
  country  = c("ISO3"),
  yr       = c("year"))
# Create a list of 72 model formulas
models <- list()
for (i in Y) for (j in 1:9) for (k in LETTERS[1:4]) {
  models[[paste0(i, " ~ clim-", k, j)]] <- as.formula(switch(k,
    A = paste0(i, "~", paste(c(X[[paste0("weather", j)]]), collapse="+")),
    B = paste0(i, "~", paste(c(X[[paste0("weather", j)]], X$socio), collapse="+")),
    C = paste0(i, "~", paste(c(X[[paste0("weather", j)]], X$socio, X$access, X$agri), collapse="+")),
    D = paste0(i, "~", paste(c(X[[paste0("weather", j)]], X$socio, X$access, X$agri, X$other), collapse="+"))
  ))
}
# Print first/last 2 formulas
pander(models[c(1:2, 71:72)])

• pcexp_ppp_m ~ clim-A1: pcexp_ppp_m ~ spei_lt + L1_speidif + L1_speidif2
• pcexp_ppp_m ~ clim-B1: pcexp_ppp_m ~ spei_lt + L1_speidif + L1_speidif2 + hh_female + agehead + marriedhead + meaneduc + depratio2 + children
• foodexp_ppp_m ~ clim-C9: foodexp_ppp_m ~ ndvi_lt + L1_ndvihishock + L1_ndviloshock + hh_female + agehead + marriedhead + meaneduc + depratio2 + children + nighlight_yearly + tt10_20k + PD12_TOT + LGP_AVG + ELEVATION + cropland_pcap + TLU_pcap
• foodexp_ppp_m ~ clim-D9: foodexp_ppp_m ~ ndvi_lt + L1_ndvihishock + L1_ndviloshock + hh_female + agehead + marriedhead + meaneduc + depratio2 + children + nighlight_yearly + tt10_20k + PD12_TOT + LGP_AVG + ELEVATION + cropland_pcap + TLU_pcap + Malaria_yearly

# Verify that all vars included in these models exist in the dataset
vars.missing <- setdiff(unlist(X,Y), dt2.lbl$varCode)
vars.missing

[1] "Lcum_speidif"  "Lcum_speidif2" "Lcum_predif"   "Lcum_predif2"  "Lcum_tempdif" 
[6] "Lcum_tempdif2" "ndvi_lt"       "Lcum_ndvidif"  "Lcum_ndvidif2"

Nine vars are missing so we need to reconstruct them or else drop the corresponding models from the analysis (done for now).

# Drop models with missing vars
models.drop <- sapply(models, str_locate, vars.missing)
models.drop <- colSums(models.drop, na.rm=T)
models <- models[models.drop==0]

This leaves 56 models to estimate.

Imputations

In R spatial regression commands require all predictor variables to be non-missing. Here we use simple region median value (or country in case regional median is missing). Other interpolation methods for missing values could be used (e.g. triangulation or kriging using district centroids). Since we are interpolating a biophysical phenomenon (SPEI) these techniques might indeed be more appropriate.

X <- unique(unlist(sapply(models, '[[', "X")))
X

g2.dt <- data.table(g2@data)

# Impute missing X values with regional median
g2.dt[, spei_lt_imp := median(spei_lt, na.rm=T), by=.(svyCode, svyL1Cd)]
g2.dt[is.na(spei_lt), spei_lt := spei_lt_imp]
g2.dt[is.na(spei_lt), .N, by=svyCode]

# Impute still missing X values with national median
g2.dt[, spei_lt_imp := median(spei_lt, na.rm=T), by=svyCode]
g2.dt[is.na(spei_lt), spei_lt := spei_lt_imp]

# Verify
g2.dt[is.na(spei_lt), .N, by=svyCode]

# Tabulate


# Re-attach imputed attributes to shapes
g2 <- SpatialPolygonsDataFrame(g2, data.frame(g2.dt), match.ID="rn")

Spatial Effects

This section is to test for the presence of spatial autocorrelation in one or more of the predictors, and to choose a regression approach such that model residuals are not spatially autocorrelated (thus ensuring a better fit).

Moran’s I Statistic

moran.plot(g2.nb@data[, Y], w, zero.policy=T,
  xlim=c(0, 200), ylim=c(0, 200),
  xlab=dt2.lbl[Y, varLabel], 
  ylab=paste("Spatially Lagged", Y))

#plot(
#  variogram(as.formula(paste(Y, "~1")), 
#    locations=coordinates(g2.nb), data=g2.nb, cloud=F), 
#  type="b", pch=16, main=paste("Variogram of", dt2.lbl[Y, varLabel]))

moran.mc(g2.nb@data[, Y], w, zero.policy=T, nsim=999)

moran.plot(g2.nb@data[, X[1]], w, zero.policy=T,
  xlab=dt2.lbl[X[1], varLabel], 
  ylab=paste("Spatially Lagged", X[1]))

moran.mc(g2.nb@data[, X[1]], w, zero.policy=T, nsim=999)

Spatial Weights

Need to choose between QUEEN or ROOK contiguity or else we can experiment with k-nearest points/shapes as neighbors using knn2nb() instead of poly2nb(). Can also assign neighbors based on a given distance threshold using dnearneigh().

Note that contiguity requires valid topology (which is surely not the case in this shapefile). Argument snap may be used to correct for slivers. Else one can use edit(nn) to make manual corrections to the matrix. Another approach to address non-contiguous splatial features is to use the feature centroid, or any weighted centroid (e.g. population weighted centroid of the admin unit).

There’s also the issue of choosing a method for the spatial weights (row-standardized, binary). Typically Row standardization is used to create proportional weights in cases where features have an unequal number of neighbors. Use Binary when you want comparable spatial parameters across different data sets with different connectivity structures.

# Generate spatial neighbour list for SSA
nb2 <- poly2nb(g2.nb, row.names=paste(g2.nb$ISO3, g2.nb$rn, sep="."))
summary(nb2)

# Verify the 5 discontiguous districts
bad <- c("477", "1052", "1343", "1746", "1947")
g2.nb.dt[rn %in% bad, .(rn, ISO3, svyL1Nm, svyL2Nm, pcexp_ppp_m)]

for(i in c("AGO", "ETH", "SEN")) {
  bb <- bbox(g2.nb[g2.nb$ISO3==i & g2.nb$rn %in% bad,])
  bb <- bb + c(-.5, -.5, .5, .5)
  print(
    tm_shape(g2, bbox=bb) + tm_borders() +
      tm_shape(g2.nb) + tm_fill("pcexp_ppp_m") + 
      tm_shape(g2.nb[g2.nb$ISO3==i & g2.nb$rn %in% bad,]) +
      tm_borders(col="red") + tm_text("svyL2Nm", col="red") +
      tm_layout(legend.outside=T)
  )}

They’re not islands, but surrounding districts have no data. Need to check a little more what’s going on, fix in QGIS if needed. Also refer to Bivand:

I did look at this 15 years ago with Boris Portnov, in the context of ESDA:

@incollection{bivand+portnov:04,
author = {R. S. Bivand and B. A. Portnov}, editor = {L. Anselin and R. J. G. M. Florax and S. J. Rey}, title = {Exploring spatial data analysis techniques using {}: the case of observations with no neighbours}, booktitle = {Advances in Spatial Econometrics: Methodology, Tools, Applications},
year = {2004},
publisher = {Springer},
address = {Berlin},
pages = {121–142} }

There are oddities in the Moran scatterplot, and also in mapping the graph-based neighbour representation into matrix form, say with the spatial lag of a no-neighbour observation’s value being zero (for zero.policy=TRUE). That paper was the basis for the zero.policy= framework. There are other consequences that you’ve found with respect to the number of subgraphs, which may or may not break formal assumptions of analysis methods. In addition, we don’t know how far the broken assumptions actually matter. This would probably be a good candidate for proper study including simulation.

Plot Contiguities

# Plot contiguities in a few countries
for (i in c("GHA", "ETH", "AGO", "SEN")) {
  
  tmp <- g2.nb[g2.nb$ISO3==i,]
  coords.tmp <- coordinates(tmp)
  nb2.tmp <- poly2nb(tmp)
  nb2.tmp
  
  plot(g2[g2$ISO3==i,], col="red", lwd=0.1)
  plot(tmp, col="grey90", lwd=0.1, add=T)
  plot(nb2.tmp, coords.tmp, col="blue", add=T)
  title(main=paste("Contiguity -", i), font.main=1)
}

# Save distance matrix for SSA (W=row standardized) to STATA for re-use
w2 <- nb2mat(nb2, style="W", zero.policy=T)
w2 <- as.data.frame(w2)
attr(w2, "var.labels") <- paste(g2$svyCode, g2$svyL1Cd, g2$svyL2Cd, sep=".")
write.dta(w2, "../out/r16.05/poverty_continguity.dta", version=12)

# Check population weights
summary(g2.nb$pop)

# Spatial weights for SSA (check doco for how to include pop weights)
# Note that If zero.policy is set to TRUE, weights vectors of zero length are inserted
# for regions without neighbour in the neighbours list. 
w <- nb2listw(nb2, zero.policy=T)

Spatial Regressions (OLS, LAG, SAC)

Batch run all models using OLS, LAG, and SAC regressions and save results to draw comparisons.

# Drop admin units with missing outcome values (not included in model)
g2.nb <- g2[!is.na(g2[[x]]),]

# Compare models
fm <- as.formula(paste(Y, "~", paste(X, collapse="+")))
fm

# Model 1: simple OLS
m <- lm(fm, data=g2.nb.dt, weights=1/pop)
summary(m)

# Examine spatial autocorrelation among the residuals
lm.morantest(m, listw=w, zero.policy=T)

# Model 2: LAG model
mlagsar <- lagsarlm(fm, w, zero.policy=T, data=g2.nb.dt)
summary(mlagsar)

# Model 3: SAC model
msacsar <- sacsarlm(fm, w, zero.policy=T, data=g2.nb.dt)
summary(msacsar)

# Also show impact effects of spatial models
# To understand the direct (local), indirect(spill-over), and total effect of a unit 
# change in each of the predictor variables
W <- as(w, "CsparseMatrix")
trMatc <- trW(W, type="mult")

impacts(mlagsar, tr=trMatc, R=2000)
impacts(msacsar, tr=trMatc, R=2000)

summary(impacts(mlagsar, tr=trMatc, R=2000), zstats=T, short=T)
summary(impacts(msacsar, tr=trMatc, R=2000), zstats=T, short=T)

The output from impacts() in the LAG model says that a 1 point increase in long-term SPEI leads to an increase in expenditure of PPP $11/month. A 1 sd increase in drought leads to a reduction in expenditure of PPP $20/month.

# Label 20 districts at random
rnd <- sample(1:nrow(g2.nb.dt), 20)

# Plot OLS
plot(m$model$pcexp_ppp_m, m$fitted.values,
  main=m$call, xlab=fm, cex=.5, pch=16,
  xlim=c(0, 300), ylim=c(0, 300))

text(m$model$pcexp_ppp_m[rnd], m$fitted.values[rnd], 
  labels=g2.nb.dt[rnd, svyL2Nm],
  cex=.6, pos=4)

# And residuals x fitted values
plot(m, which=3)

# Plot LAG
plot(mlagsar$y, mlagsar$fitted.values,
  main=mlagsar$call, xlab=fm, cex=.5, pch=16,
  xlim=c(0, 300), ylim=c(0, 300))

text(mlagsar$y[rnd], mlagsar$fitted.values[rnd], 
  labels=g2.nb.dt[rnd, svyL2Nm],
  cex=.6, pos=4)

# Plot SAC
plot(msacsar$y, msacsar$fitted.values,
  main=msacsar$call, xlab=fm, cex=.5, pch=16,
  xlim=c(0, 300), ylim=c(0, 300))

text(msacsar$y[rnd], msacsar$fitted.values[rnd], 
  labels=g2.nb.dt[rnd, svyL2Nm],
  cex=.6, pos=4)

Geographically Weighted Regression

See Bivand at https://cran.r-project.org/web/packages/spgwr/vignettes/GWR.pdf and Anselin http://www.csiss.org/gispopsci/workshops/2011/PSU/readings/W15_Anselin2007.pdf. Also Brunsdon http://rpubs.com/chrisbrunsdon/101305.

Note that sampling weights are not implemented. Choosing a method to estimate optimal bandwidth is unclear (check doco), also not clear how to choose a kernel function (default Gaussian).

# Load package
library(spgwr)

# Try GWR/LM on same model as above (pass shapes, will return shapes with coeff)
bwG <- gwr.sel(fm, data=g2.nb, gweight=gwr.Gauss, verbose=FALSE)
gwrG <- gwr(fm, data=g2.nb, bandwidth=bwG, gweight=gwr.Gauss, hatmatrix=TRUE)
gwrG

# Map coefficients
data(World)

for (i in X) {print(
    tm_shape(gwrG$SDF, is.master=T) + 
      tm_fill(i, palette="RdYlGn", style="jenks", n=9,
        title=stringr::str_wrap(dt2.lbl[Y, varLabel], 30)) +
      tm_shape(World) + tm_borders(lwd=0.1) +
      tm_layout(
        title=dt2.lbl[i, varLabel],
        title.snap.to.legend=T, legend.outside=T)
  )}

# Could also try a GWR/GLM on same model as above
bwG <- ggwr.sel(fm, data=g2.nb, gweight=gwr.Gauss, verbose=FALSE)
gwrG <- ggwr(fm, data=g2.nb, bandwidth=bwG, gweight=gwr.Gauss, hatmatrix=TRUE,
  family="poisson")

save.image("./tmp/poverty_r16.05.RData")