HM1 Exploratory Analysis

# San Antonio in Bexar county
bexar<-subset(sa08, sa08$cofips=="029")
# tract ID
bexar$tract<-substr(bexar$D_TR2000,6,11)
# aggregate deaths by tract in Bexar
bextr<-data.frame(xtabs(~tract, bexar, drop.unused.levels=T))
trs<-bextr$tract
tractmort<-data.frame(tract=bextr$tract, sa08=bextr$Freq)

library(spdep)

## Loading required package: sp

## Loading required package: Matrix

## Loading required package: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source'))`

library(maptools)

## Checking rgeos availability: FALSE
##      Note: when rgeos is not available, polygon geometry     computations in maptools depend on gpclib,
##      which has a restricted licence. It is disabled by default;
##      to enable gpclib, type gpclibPermit()

library(RColorBrewer)
library(classInt)
library(ggplot2)

sashp<-readShapePoly("SA_classdata.shp")

## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read

mdat<-sashp@data

# We merge deaths to the shapefile
mdat2<-merge(mdat, tractmort, by.x="TractID", by.y="tract", all.x=T, sort=F)

sashp@data<-mdat2
writePolyShape(sashp, fn="SA_2008merged.shp")

## Warning: writePolyShape is deprecated; use rgdal::writeOGR or sf::st_write

data<-sashp

# Crude Mortality Rate, SA 2008
data$mortrate<-1000*data$sa08/data$acs_poptot
data$mortrate<-ifelse(is.na(data$mortrate)==T,0, data$mortrate)

Exploratory Analysis

# Mortality Rate 2008
qplot(data$mortrate, geom="histogram", binwidth=.15)+ggtitle("Histogram of CMR in SA, 2008") + xlab("CMR")

# Boxplots
# boxplot of mortality rate 2008
data$hispq<-cut(data$PHISPANIC, quantile(data$PHISPANIC))
data$whiteq<-cut(data$PWHITE, quantile(data$PWHITE))
data$blackq<-cut(data$PBLACK, quantile(data$PBLACK))

ggplot(data@data, aes(x=hispq, y=data$mortrate)) + geom_boxplot()+ ggtitle("Hispanic population and CMR", "San Antonio, 2008")+xlab("Hispanic population (quantiles)") +ylab("CMR")

ggplot(data@data, aes(x=whiteq, y=data$mortrate)) + geom_boxplot()+ ggtitle("Non-Hispanic White population and CMR", "San Antonio, 2008")+xlab("Non-Hispanic White population (quantiles)") +ylab("CMR")

ggplot(data@data, aes(x=blackq, y=data$mortrate)) + geom_boxplot()+ ggtitle("Non-Hispanic Black population and CMR", "San Antonio, 2008")+xlab("Non-Hispanic Black population (quantiles)") +ylab("CMR")

# Quantile maps of SA
cmrmap2

cmrmap

Global Moran’s I

Choosing a Rook adjacency matrix gives us a greater Moran’s I than the selection of a Queen adjacency matrix. However the difference is barely visible bewteen the scatterplots.

# Moran's I for the CMR
moran.test(data$mortrate, listw=saq) # Queen 

## 
##  Moran I test under randomisation
## 
## data:  data$mortrate  
## weights: saq    
## 
## Moran I statistic standard deviate = 2.8086, p-value = 0.002488
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.106251726      -0.004273504       0.001548625

moran.mc(data$mortrate, listw=saq, nsim=999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  data$mortrate 
## weights: saq  
## number of simulations + 1: 1000 
## 
## statistic = 0.10625, observed rank = 994, p-value = 0.006
## alternative hypothesis: greater

moran.plot(data$mortrate, listw=saq, main="Moran Scatterplot (Queen)", xlab = "CMR", ylab="Spatially lagged CMR")

moran.test(data$mortrate, listw=sar) # Rook 

## 
##  Moran I test under randomisation
## 
## data:  data$mortrate  
## weights: sar    
## 
## Moran I statistic standard deviate = 2.9408, p-value = 0.001637
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.122655508      -0.004273504       0.001862871

moran.mc(data$mortrate, listw=sar, nsim=999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  data$mortrate 
## weights: sar  
## number of simulations + 1: 1000 
## 
## statistic = 0.12266, observed rank = 997, p-value = 0.003
## alternative hypothesis: greater

moran.plot(data$mortrate, listw=sar, main="Moran Scatterplot (Rook)", xlab = "CMR", ylab="Spatially lagged CMR")

Local Moran’s I (Queen)

Local Moran’s I analysis showed a significant cluster of CMR by tracts. There is no big difference between the decision of select Queen or Rook adjacency matrices.

locali<-localmoran(data$mortrate, saq, p.adjust.method="fdr")
data$locali<-locali[,1]
data$localp<-locali[,5]
data$cl<-as.factor(ifelse(data$localp<=.05,"Clustered","Not Clustered"))
# Local Moran's I
lmq<-spplot(data, "locali", main="Local Moran's I", at=quantile(data$locali), col.regions=brewer.pal(n=4, "RdBu"))
lmqcl<-spplot(data, "cl", main="Local Moran Clusters", col.regions=c(2,0))
lmq

lmqcl

Local Moran’s I (Rook)

localir<-localmoran(data$mortrate, sar, p.adjust.method="fdr")
data$localir<-localir[,1]
data$localpr<-localir[,5]
data$clr<-as.factor(ifelse(data$localpr<=.05,"Clustered","Not Clustered"))
# Local Moran's I
lmr<-spplot(data, "localir", main="Local Moran's I", at=quantile(data$localir), col.regions=brewer.pal(n=4, "RdBu"))
lmrcl<-spplot(data, "clr", main="Local Moran Clusters", col.regions=c(2,0))
lmr

lmrcl