Loading Data

tx_death<-read.dbf("C:/Users/munta/Google Drive/### Spatial Demography Dem 7263/Spatial Demography Corey Materials/data/2008_Texas_Deaths.dbf")
tx_death$cofips<-substr(tx_death$D_TR2000, 3,5)

#get only bexar county, fips==029
bexar08<-subset(tx_death, tx_death$cofips=="029")

Forming the tract ID

bexar08$tract<-substr(bexar08$D_TR2000,6,11)

#aggregate by tract
bextr08<-data.frame(xtabs(~tract, bexar08, drop.unused.levels=T))

trs<-bextr08$tract

Making the Dataframe

#make a dataframe with an id and deaths data
trmort<-data.frame(tract=bextr08$tract, tx_death=bextr08$Freq)
head(trmort)
##    tract tx_death
## 1 110100       41
## 2 110200        3
## 3 110300       32
## 4 110400       12
## 5 110500       15
## 6 110600       36

Reading Shape File

sa<-readShapePoly("C:/Users/munta/Google Drive/### Spatial Demography Dem 7263/Spatial Demography Corey Materials/data/SA_classdata.shp")
## Warning: use rgdal::readOGR or sf::st_read
mdat<-sa@data

Merging the death data to the shapefile

#Creting a mortality rate
mdat2<-merge(mdat, trmort, by.x="TractID", by.y="tract", all.x=T, sort=F)

sa@data<-mdat2

writePolyShape(sa, fn="C:/Users/munta/Google Drive/### Spatial Demography Dem 7263/Spatial Demography Corey Materials/data/SA_2008merged.shp")
## Warning: use rgdal::readOGR or sf::st_read
mydata <- sa
mydata$mortrate <- 1000*mydata$tx_death / mydata$acs_poptot 

mydata$mortrate<-ifelse(is.na(mydata$mortrate)==T,0, mydata$mortrate)

Performing Exploratory Analysis

Histogram

qplot(mydata$mortrate, geom="histogram", binwidth=.5)+ ggtitle("Distribution of Mortality Rate in San Antonio, 2008")+ xlab("Mortality Rate")

The histogram shows that the mortality rate is around 7 per 1000 population.

Boxplot

mydata$wht<-cut(mydata$PWHITE, quantile(mydata$PWHITE))
mydata$blk<-cut(mydata$PBLACK, quantile(mydata$PBLACK))
mydata$hisp<-cut(mydata$PHISPANIC, quantile(mydata$PHISPANIC))

ggplot(mydata@data, aes(x=wht, y=mydata$mortrate)) + geom_boxplot()+ ggtitle("Mortality rate among White population in San Antonio, 2008")+xlab("White population (quantiles)") +ylab("Mortality Rate")

ggplot(mydata@data, aes(x=blk, y=mydata$mortrate)) + geom_boxplot()+ ggtitle("Mortality rate among Black population in San Antonio, 2008")+xlab("Black population (quantiles)") +ylab("Mortality Rate")

ggplot(mydata@data, aes(x=hisp, y=mydata$mortrate)) + geom_boxplot()+ ggtitle("Mortality rate among Hispanic population in San Antonio, 2008")+xlab("Hispanic population (quantiles)") +ylab("Mortality Rate")

Mortality Rate in San Antonio (Map)

spplot(mydata, "mortrate",
       main="Mortality Rate in San Antonio, 2008",
       at=quantile(mydata$mortrate), 
       col.regions=brewer.pal(n=5, "Reds"))

Moran’s I Analysis of the Mortality Rate

moran.test(mydata$mortrate, listw=saq2)
## 
##  Moran I test under randomisation
## 
## data:  mydata$mortrate  
## weights: saq2    
## 
## 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(mydata$mortrate, listw=saq2, nsim=999)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mydata$mortrate 
## weights: saq2  
## number of simulations + 1: 1000 
## 
## statistic = 0.10625, observed rank = 994, p-value = 0.006
## alternative hypothesis: greater

Moran Scatterplot (Queen)

moran.plot(mydata$mortrate, listw=saq2, main="Moran Scatterplot (Queen)", xlab = "Mortality Rate", ylab="Spatially lagged Mortality Rate")

Moran’s I Analysis of the Mortality Rate

moran.test(mydata$mortrate, listw=sar2)
## 
##  Moran I test under randomisation
## 
## data:  mydata$mortrate  
## weights: sar2    
## 
## 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(mydata$mortrate, listw=sar2, nsim=999)
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mydata$mortrate 
## weights: sar2  
## number of simulations + 1: 1000 
## 
## statistic = 0.12266, observed rank = 997, p-value = 0.003
## alternative hypothesis: greater

Moran Scatterplot (Rook)

moran.plot(mydata$mortrate, listw=sar2, main="Moran Scatterplot (Rook)", xlab = "Mortality Rate", ylab="Spatially lagged Mortality Rate")

Both Queen and Rook analysis show that there’s not much of a visible difference between the two plots; however, the Rook adjancey matrix give us a greater value than the Queen adjacency matrix.

Local Moran’s I (Queen)

locali<-localmoran(mydata$mortrate, saq2, p.adjust.method="fdr")
mydata$locali<-locali[,1]
mydata$localp<-locali[,5]

#Making our own cluster identifiers
mydata$cl<-as.factor(ifelse(mydata$localp<=.05,"Clustered","NotClustered"))


#Plotting the results
spplot(mydata, "locali", main="Local Moran's I", 
       at=quantile(mydata$locali), 
       col.regions=brewer.pal(n=4, "RdBu"))

#Local clusters
spplot(mydata, "cl", main="Local Moran's Clusters (Queen)", 
       col.regions=c(2,0))

Local Moran’s I (Rook)

localir<-localmoran(mydata$mortrate, sar2, p.adjust.method="fdr")
mydata$localir<-localir[,1]
mydata$localpr<-localir[,5]
mydata$clr<-as.factor(ifelse(mydata$localpr<=.05,"Clustered","Not Clustered"))
spplot(mydata, "localir", main="Local Moran's I (Rook)", at=quantile(mydata$localir), col.regions=brewer.pal(n=4, "RdBu"))

spplot(mydata, "clr", main="Local Moran Clusters", col.regions=c(2,0))

The local moran analysis show us that mortality rates are spatially clustered among the census tracts with not much of a difference between the two methods - queen and rook.