index.Rmd

---
title: "cl-spatial-epi"
output: 
  flexdashboard::flex_dashboard:
    navbar:
      - { title: "hesscl.com", href: "http://hesscl.com", align: right }
    orientation: columns
    theme: bootstrap
    logo: logo.svg
    source_code: "https://github.com/hesscl/cl-spatial-epi"
---

```{r setup, include=FALSE}

#packages
library(tidyverse)
library(lubridate)
library(sqldf)
library(ggthemes)
library(gridExtra)
library(RColorBrewer)
library(viridis)
library(sp)
library(spdplyr)
library(rgdal)
library(rgeos)
library(spdep)
library(INLA)
library(flexdashboard)
library(DT)
library(plotly)
library(coefINLA)

#load 2012-2016 ACS data (2016 tracts same as 2010 for KC)
census <- read_csv(file = "./input/nhgis0052_ds225_20165_2016_tract.csv")

#read in tract shapefile for king county
kc_shp <- readOGR(dsn = "./input/KCTract2010/KCTract2010.shp",
                  layer = "KCTract2010",
                  GDAL1_integer64_policy = TRUE,
                  stringsAsFactors = F)

#determine config (dev or github)
if(file.exists("../data/cl/craigslistDB.sqlite")){
  DB <- dbConnect(SQLite(), dbname="../data/cl/craigslistDB.sqlite")
  cl <- tbl(DB, "clean") #clean listing table
  
  #compute tract aggregates for CL listing count
  tractCl <- cl %>%
    filter(!is.na(GISJOIN), !is.na(cleanBeds), !is.na(cleanRent), !is.na(cleanSqft)) %>% #only listings with valid Bed/Rent
    filter(GISJOIN %in% kc_shp@data$GISJOIN) %>% #filter to KC only (db has metro area)
    dplyr::select(listingMoYr, GISJOIN, seattle, matchAddress, matchAddress2, matchType, cleanBeds, cleanRent) %>% #SELECT these columns
    collect %>% #bring db query into memory
    mutate(listingMoYr = as.Date(listingMoYr)) %>%
    filter(listingMoYr < "2018-03-12") %>% #everything up till march 12, 2018  (i.e cut off new data)
    filter(!grepl("Google", matchType)) %>% #no Google geocodes, only Smartystreets (precise to Zip9)
    distinct(matchAddress, matchAddress2, cleanBeds, cleanRent, cleanSqft, .keep_all = T) %>% #dedupe unique address-bed-rent combos
    group_by(GISJOIN) %>% #group listings by tract
    summarize(nListings = n(),
              seattle = ifelse(is.na(max(seattle)),0, max(seattle)))
  dbDisconnect(DB)
  
} else{
  stop("Input data needed")
}

#merge tables together, allow missing values for tracts with no listings (i.e. no value in tract)
test <- left_join(census, tractCl)

#compute some new vars/mutate sensibly named ACS vars
test <- test %>%
  mutate(nListings = ifelse(is.na(nListings), 0, nListings), #if NA then 0, since no listings were observed
         nHU = AF7NE001,
         GISJOIN = as.character(GISJOIN)) %>%
  mutate(tpop = AF2LE001,
         tnhw = AF2UE003,
         tnhb = AF2UE004,
         tnha = AF2UE006,
         thsp = AF2UE012,
         tnho =  AF2UE005+AF2UE007+AF2UE008+AF2UE009,
         nforborn = AF95E005,
         nvachu = AF7OE003,
         nforrent = AF7ZM002,
         nocc =  AF7OE002,
         ppov = (AF43E002+AF43E003)/AF43E001,
         nownocc = AF7PE002,
         medGRentACS = AF89M001,
         medStrYr = AF8IE001,
         medVal = AF9LE001,
         medHHInc = AF49M001) %>%
  mutate(pnhw = tnhw/tpop,
         pnhb = tnhb/tpop,
         pnho = tnho/tpop,
         pnha = tnha/tpop,
         phsp = thsp/tpop,
         pforborn = nforborn/tpop,
         pvac = nvachu/nHU,
         pforrent = ifelse(nHU == 0, 0, AF7ZM002/nHU),
         pownocc = nownocc/nocc,
         pblt14lat = AF8HE002/nHU) %>%
  filter(tpop > 0) %>%
  select(-starts_with("AF"))

kc_shp <- readOGR(dsn = "./input/KCTract2010/KCTract2010.shp",
                  layer = "KCTract2010",
                  GDAL1_integer64_policy = TRUE,
                  stringsAsFactors = F, verbose = F)

test <- test[match(kc_shp$GISJOIN, test$GISJOIN),]
test$GISJOIN == kc_shp$GISJOIN #dbl check order
rownames(test) <- row.names(kc_shp)

#spCbind the data frame onto shapefile
kc_shp <- maptools::spCbind(kc_shp, test) #both are complete tables of KC tracts
kc_df <- as.tbl(kc_shp@data) #save table (now in the order of the shapefile)
kc_adj <- poly2nb(kc_shp)
nb2INLA("./output/graphINLA/kctract.graph", kc_adj)

#create a couple numeric ID columns to use later for INLA 
idx <- 1:nrow(kc_df) 
idxx <- idx

#INLA seems to be finicky and want either option?
kc_df$idx <- 1:nrow(kc_df) 
kc_df$idxx <- kc_df$idx

#extract the dataframe to hand to INLA
kc_df <- as.data.frame(kc_shp@data)

#for race/eth viz
blk <- kc_df %>% 
  mutate(race = "Black",
         percent = pnhb) %>%
  select(race, percent, nListings, GISJOIN)

hsp <- kc_df %>% 
  mutate(race = "Hispanic",
         percent = phsp) %>%
  select(race, percent, nListings, GISJOIN)

asi <- kc_df %>% 
  mutate(race = "Asian",
         percent = pnha) %>%
  select(race, percent, nListings, GISJOIN)

wht <- kc_df %>% 
  mutate(race = "White",
         percent = pnhw) %>%
  select(race, percent, nListings, GISJOIN)

oth <- kc_df %>%
  mutate(race = "Other",
         percent = pnho) %>%
  select(race, percent, nListings, GISJOIN)

race_dens <- bind_rows(asi, blk, hsp, oth, wht)
race_dens$race <- as.factor(race_dens$race)

### Model 0: Negative Binomial (no random effects)
form0 <- nListings ~ 1 + log(nHU) + pforrent + pownocc + pblt14lat +
  pnhb + phsp + pnha + pnho + medHHInc + ppov + pforborn + seattle

m0 <- inla(form0, 
           family = "nbinomial", 
           data = kc_df,
           control.predictor = list(compute = TRUE),
           control.compute = list(dic = TRUE, waic = TRUE))

#predicted values
m0median <- exp(m0$summary.linear.predictor["0.5quant"])
m0lower <- exp(m0$summary.linear.predictor["0.025quant"])
m0upper <- exp(m0$summary.linear.predictor["0.975quant"])
kc_shp@data$m0median <- m0median[,1]
kc_shp@data$m0lower <- m0lower[,1]
kc_shp@data$m0upper <- m0upper[,1]
kc_shp@data$m0post95wid <- kc_shp@data$m0upper - kc_shp@data$m0lower


### Model 1: Lognormal Non-Spatial Negative Binomial Model
form1 <- nListings ~ 1 + log(nHU) + pforrent + pownocc + pblt14lat +
  pnhb + phsp + pnha + pnho + medHHInc + ppov + pforborn + seattle +
  f(idx, model = "iid")

m1 <- inla(form1, 
           family = "nbinomial", 
           data = kc_df,
           control.predictor = list(compute = TRUE),
           control.compute = list(dic = TRUE, waic = TRUE))

#lognormal random effect medians
m1RE <- exp(m1$summary.random$idx[5])
kc_shp@data$m1RE <- m1RE[,1]

#predicted values
m1median <- exp(m1$summary.linear.predictor["0.5quant"])
m1lower <- exp(m1$summary.linear.predictor["0.025quant"])
m1upper <- exp(m1$summary.linear.predictor["0.975quant"])
kc_shp@data$m1median <- m1median[,1]
kc_shp@data$m1lower <- m1lower[,1]
kc_shp@data$m1upper <- m1upper[,1]
kc_shp@data$m1post95wid <- kc_shp@data$m1upper - kc_shp@data$m1lower


### Model 2: Lognormal Spatial Negative Binomial Model
form2 <- nListings ~ 1 + log(nHU) + pforrent + pownocc + pblt14lat +
  pnhb + phsp + pnha + pnho + medHHInc + ppov + pforborn + seattle +
  f(idx, model = "iid") +
  f(idxx, model = "besag", graph = "./output/graphINLA/kctract.graph")

m2 <- inla(form2, 
           family = "nbinomial", 
           data = kc_df,
           control.predictor = list(compute = TRUE),
           control.compute = list(dic = TRUE, waic = TRUE))

#lognormal random effect medians
m2RE <- exp(m2$summary.random$idx[5])
kc_shp@data$m2RE <- m2RE[,1]

#spatial effect median
m2SE <- exp(m2$summary.random$idxx[5])
kc_shp@data$m2SE <- m2SE[,1]

#predicted values
m2median <- exp(m2$summary.linear.predictor["0.5quant"])
m2lower <- exp(m2$summary.linear.predictor["0.025quant"])
m2upper <- exp(m2$summary.linear.predictor["0.975quant"])
kc_shp@data$m2median <- m2median[,1]
kc_shp@data$m2lower <- m2lower[,1]
kc_shp@data$m2upper <- m2upper[,1]
kc_shp@data$m2post95wid <- kc_shp@data$m2upper - kc_shp@data$m2lower

### Model 3: GLM Poisson Model
form3 <- nListings ~ 1 + log(nHU) + pforrent + pownocc + pblt14lat +
  pnhb + phsp + pnha + pnho + medHHInc + ppov + pforborn + seattle

m3 <- inla(form3, 
           family = "poisson", 
           data = kc_df,
           control.predictor = list(compute = TRUE),
           control.compute = list(dic = TRUE, waic = TRUE))
summary(m3)
plot(m3)


#predicted values
m3median <- exp(m3$summary.linear.predictor["0.5quant"])
m3lower <- exp(m3$summary.linear.predictor["0.025quant"])
m3upper <- exp(m3$summary.linear.predictor["0.975quant"])
kc_shp@data$m3median <- m3median[,1]
kc_shp@data$m3lower <- m3lower[,1]
kc_shp@data$m3upper <- m3upper[,1]
kc_shp@data$m3post95wid <- kc_shp@data$m3upper - kc_shp@data$m3lower


### Model 4: Lognormal Non-Spatial Poisson Model
form4 <- nListings ~ 1 + log(nHU) + pforrent + pownocc + pblt14lat +
  pnhb + phsp + pnha + pnho + medHHInc + ppov + pforborn + seattle +
  f(idx, model = "iid")

m4 <- inla(form4, 
           family = "poisson", 
           data = kc_df,
           control.predictor = list(compute = TRUE),
           control.compute = list(dic = TRUE, waic = TRUE))
summary(m4)
plot(m4)

#lognormal random effect medians
m4RE <- exp(m4$summary.random$idx[5])
kc_shp@data$m4RE <- m4RE[,1]

#predicted values
m4median <- exp(m4$summary.linear.predictor["0.5quant"])
m4lower <- exp(m4$summary.linear.predictor["0.025quant"])
m4upper <- exp(m4$summary.linear.predictor["0.975quant"])
kc_shp@data$m4median <- m4median[,1]
kc_shp@data$m4lower <- m4lower[,1]
kc_shp@data$m4upper <- m4upper[,1]
kc_shp@data$m4post95wid <- kc_shp@data$m4upper - kc_shp@data$m4lower


### Model 5: Lognormal Spatial Poisson Model
form5 <- nListings ~ 1 + log(nHU) + pforrent + pownocc + pblt14lat +
  pnhb + phsp + pnha + pnho + medHHInc + ppov + pforborn + seattle +
  f(idx, model = "iid") +
  f(idxx, model = "besag", graph = "./output/graphINLA/kctract.graph")

m5 <- inla(form5, 
           family = "poisson", 
           data = kc_df,
           control.predictor = list(compute = TRUE),
           control.compute = list(dic = TRUE, waic = TRUE))
summary(m5)
plot(m5)

#lognormal random effect medians
m5RE <- exp(m5$summary.random$idx[5])
kc_shp@data$m5RE <- m5RE[,1]

#spatial effect median
m5SE <- exp(m5$summary.random$idxx[5])
kc_shp@data$m5SE <- m5SE[,1]

#predicted values
m5median <- exp(m5$summary.linear.predictor["0.5quant"])
m5lower <- exp(m5$summary.linear.predictor["0.025quant"])
m5upper <- exp(m5$summary.linear.predictor["0.975quant"])
kc_shp@data$m5median <- m5median[,1]
kc_shp@data$m5lower <- m5lower[,1]
kc_shp@data$m5upper <- m5upper[,1]
kc_shp@data$m5post95wid <- kc_shp@data$m5upper - kc_shp@data$m5lower

#gof
rmse <- function(error)
{
  sqrt(mean(error^2))
}

#marginal likelihood
m0mlik <- m0$mlik[1]
m1mlik <- m1$mlik[1]
m2mlik <- m2$mlik[1]
m3mlik <- m3$mlik[1]
m4mlik <- m4$mlik[1]
m5mlik <- m5$mlik[1]

#DIC
m0DIC <- m0$dic$dic
m1DIC <- m1$dic$dic
m2DIC <- m2$dic$dic
m3DIC <- m3$dic$dic
m4DIC <- m4$dic$dic
m5DIC <- m5$dic$dic

#WAIC
m0WAIC <- m0$waic$waic
m1WAIC <- m1$waic$waic
m2WAIC <- m2$waic$waic
m3WAIC <- m3$waic$waic
m4WAIC <- m4$waic$waic
m5WAIC <- m5$waic$waic

#RMSE
m0RMSE <- rmse(kc_shp@data$nListings-kc_shp@data$m0median)
m1RMSE <- rmse(kc_shp@data$nListings-kc_shp@data$m1median)
m2RMSE <- rmse(kc_shp@data$nListings-kc_shp@data$m2median)
m3RMSE <- rmse(kc_shp@data$nListings-kc_shp@data$m3median)
m4RMSE <- rmse(kc_shp@data$nListings-kc_shp@data$m4median)
m5RMSE <- rmse(kc_shp@data$nListings-kc_shp@data$m5median)

#compile df
mfit <- data.frame(
  Distribution = c("Neg. Binomial", "Neg. Binomial", "Neg. Binomial", "Poisson", "Poisson", "Poisson"),
  Model = c("GLM", "Non-Spatial RE", "Spatial RE", "GLM", "Non-Spatial RE", "Spatial RE"),
  MLik = c(m0mlik, m1mlik, m2mlik, m3mlik, m4mlik, m5mlik),
  RMSE = c(m0RMSE, m1RMSE, m2RMSE, m3RMSE, m4RMSE, m5RMSE),
  DIC = c(m0DIC, m1DIC, m2DIC, m3DIC, m4DIC, m5DIC),
  WAIC = c(m0WAIC, m1WAIC, m2WAIC, m3WAIC, m4WAIC, m5WAIC)
)


mfit[,3:6] <- round(mfit[,3:6], 2)


pal <- brewer.pal(7, "Purples")

#maps
kc_shp@data$id <- rownames(kc_shp@data)
kc_f <- fortify(kc_shp)
kc_f <- inner_join(kc_f, kc_shp@data, "id")

#for coefINLA
var_labs <- c(
  `(Intercept)`="Intercept",
  `log(nHU)`="log(N HU)",
  pforrent="Prp HU For Rent",
  pblt14lat="Prp HU Blt >2014",
  pownocc="Prp HU Own-Occ",
  pnhb="Prp Black",
  phsp="Prp Hispanic",
  pnha="Prp Asian",
  pnho="Prp Other",
  medHHInc="Med. HH Income",
  ppov="Poverty Rate",
  pforborn="Prp Foreign Born",
  seattle="Seattle"
)

var_labeller <- labeller(
     var = var_labs
)

```

Home {data-icon="fa-home"}
=====================================

Column
-------------------------------------

### Welcome!

This website provides supplementary information for the paper _Sociodemographic and Spatial Variations in Craigslist Rental Housing Activity_

You can download the R source code and data for this project from Github [here](https://github.com/hesscl/cl-spatial-epi)

Contact Chris Hess at [hesscl@uw.edu](mailto://hesscl@uw.edu) for more information about this report.

_This page was last updated: `r Sys.Date()`_

Column
-------------------------------------

### Histogram of rental listings scraped per neighborhood (census tracts)

```{r nlistings, message=F}
nListings.gg <- ggplot(kc_shp@data, aes(x = nListings)) +
  geom_histogram(bins = 50, fill = pal[3]) +
  theme_minimal() +
  xlab("N Listings in Neighborhood") +
  ylab("Frequency")

ggplotly(nListings.gg, tooltip = NULL)
```


Total Population and Socioeconomics {data-navmenu="Description"}
=====================================

Column
-------------------------------------

### Total Population (ACS 2012-2016)

```{r pop.dens, message=F}
pop.gg <- ggplot(kc_shp@data, aes(x = log(tpop))) +
  geom_density(fill = pal[7]) +
  theme_minimal() +
  xlab("\nLog(Total Population)") +
  ylab("Density\n")
ggplotly(pop.gg)
```

### Prop. Persons below Federal Poverty Line (ACS 2012-2016)

```{r pov.dens, message=F}
pov.gg <- ggplot(kc_shp@data, aes(x = ppov)) +
  geom_density(fill = pal[7]) +
  theme_minimal() +
  xlab("\nProp. Persons below Federal Poverty Line") +
  ylab("Density\n")
ggplotly(pov.gg)
```

Column
-------------------------------------

### Median HH Income (ACS 2012-2016)

```{r medhhinc.dens, message=F}
inc.gg <- ggplot(kc_shp@data, aes(x = medHHInc)) +
  geom_density(fill = pal[7]) +
  theme_minimal() +
  xlab("\nMedian HH Income") +
  ylab("Density\n")
ggplotly(inc.gg)
```

### Prop. Foreign-Born (ACS 2012-2016)

```{r pforborn.dens, message=F}
pforborn.gg <- ggplot(kc_shp@data, aes(x = pforborn)) +
  geom_density(fill = pal[7]) +
  theme_minimal() +
  xlab("\nProp. Persons Foreign-born") +
  ylab("Density\n")
ggplotly(pforborn.gg)
```


Housing Stock {data-navmenu="Description"}
=====================================

Column
-------------------------------------

### N Housing Units in Tract (ACS 2012-2016)

```{r hu.dens, message=F}
hu.gg <- ggplot(kc_shp@data, aes(x = nHU)) +
  geom_density(fill = pal[7]) +
  theme_minimal() +
  xlab("\nNumber of Housing Units") +
  ylab("Density\n")
ggplotly(hu.gg)
```

### Prop. HU For Rent (ACS 2012-2016)

```{r forrent.dens, message=F}
forrent.gg <- ggplot(kc_shp@data, aes(x = pforrent)) +
  geom_density(fill = pal[7]) +
  theme_minimal() +
  xlab("\nProp. HU For Rent") +
  ylab("Density\n")
ggplotly(forrent.gg)
```

Column
-------------------------------------

### Prop. HU Built 2014 or Later (ACS 2012-2016)

```{r pblt.dens, message=F}
#nHU ACS density
pblt.gg <- ggplot(kc_shp@data, aes(x = pblt14lat)) +
  geom_density(fill = pal[7]) +
  theme_minimal() +
  xlab("\nProp. HU Built 2014 or Later") +
  ylab("Density\n")
ggplotly(pblt.gg)
```

### Prop. HU Owner-Occupied (ACS 2012-2016)

```{r pownocc.dens, message=F}
pownocc.gg <- ggplot(kc_shp@data, aes(x = pownocc)) +
  geom_density(fill = pal[7]) +
  theme_minimal() +
  xlab("\nProp. HU Owner-occupied") +
  ylab("Density\n")
ggplotly(pownocc.gg)
```


Racial/Ethnic Composition {data-navmenu="Description"}
=====================================

Column
-------------------------------------

### Boxplot of Racial/Ethnic Compositions

```{r raceeth1}
raceeth1 <- ggplot(race_dens, aes(x = race, y = percent, group = race, color = race)) + 
  geom_boxplot() +
  scale_fill_brewer(palette = "Set1") +
  theme_minimal() +
  theme(legend.margin=margin(r = 1, unit='cm'))

ggplotly(raceeth1) %>%
  highlight(on = "plotly_hover")
```

Column
-------------------------------------

### Smoothed Fit of N Listings ~ Racial/Ethnic Composition

```{r raceeth2}
raceeth2 <- ggplot(race_dens, aes(x = percent, y = nListings, group = race, color = race)) + 
  geom_smooth(se = F) + 
  scale_x_continuous(labels = scales::percent) +
  scale_fill_brewer(palette = "Set1") +
  theme_minimal() +
  theme(legend.margin=margin(r = 1, unit='cm')) +
  xlab("\n% of Neighborhood Population") +
  ylab("N Listings\n")

ggplotly(raceeth2) %>%
  highlight(on = "plotly_hover")
```


Goodness of Fit {data-navmenu="Models"}
=====================================

Column
-------------------------------------

### Actual versus Predicted

```{r avp, message=F}
m0avp <- kc_shp@data %>%
  mutate(model = "GLM",
         dist = "NB",
         actual = nListings,
         predicted = m0median,
         GISJOIN = GISJOIN) 

m1avp <- kc_shp@data %>%
  mutate(model = "Non-Spatial RE",
         dist = "NB",
         actual = nListings,
         predicted = m1median,
         GISJOIN = GISJOIN) 

m2avp <- kc_shp@data %>%
  mutate(model = "Spatial RE",
         dist = "NB",
         actual = nListings,
         predicted = m2median,
         GISJOIN = GISJOIN) 

m3avp <- kc_shp@data %>%
  mutate(model = "GLM",
         dist = "Poisson",
         actual = nListings,
         predicted = m3median,
         GISJOIN = GISJOIN) 

m4avp <- kc_shp@data %>%
  mutate(model = "Non-Spatial RE",
         dist = "Poisson",
         actual = nListings,
         predicted = m4median,
         GISJOIN = GISJOIN) 

m5avp <- kc_shp@data %>%
  mutate(model = "Spatial RE",
         dist = "Poisson",
         actual = nListings,
         predicted = m5median,
         GISJOIN = GISJOIN) 
avp <- bind_rows(m0avp, m1avp, m2avp, m3avp, m4avp, m5avp)

avp.gg <- ggplot(avp, aes(x = actual, y = predicted, color = model, shape = dist, group = GISJOIN)) +
  facet_grid(dist ~ .) +
  geom_abline(slope = 1, intercept = 0, color = "grey80") +
  geom_point() +
  coord_cartesian(xlim = c(0, 8500), ylim = c(0, 8500)) +
  xlab("Actual") +
  ylab("Predicted") +
  scale_color_brewer(palette = "Purples") +
  theme_minimal()

ggplotly(avp.gg) %>%
  highlight()
```

Column
-------------------------------------

### Goodness-of-Fit Statistics

```{r mfit}
datatable(mfit,
          rownames = F,
          style = "bootstrap",
          autoHideNavigation = T,
          options = list(paging = F, searching = F, autoWidth = T))

```


Fixed Effects {data-navmenu="Models"}
=====================================

Column
-------------------------------------

### Poisson Non-Spatial RE Model posterior distributions (with 95% intervals)

```{r m4coef, fig.dpi = 300, fig.height=6, fig.width=6}
#like a ropeladder, but bayes-y 
c4 <- coefINLA(m4, labeller = var_labeller) + 
  scale_x_continuous(limits = c(-20, 20)) +
  xlab("Exponentiated Coefficients")

ggplotly(c4, tooltip = c("var", "med")) %>%
  layout(margin = list(r = 125, b = 100))
```


Column
-------------------------------------

### Poisson Spatial RE Model posterior distributions (with 95% intervals)

```{r m5coef, fig.dpi = 300, fig.height=6, fig.width=6}
c5 <- coefINLA(m5, labeller = var_labeller) + 
  scale_x_continuous(limits = c(-20, 20)) +
  xlab("Exponentiated Coefficients")

ggplotly(c5, tooltip = c("var", "med")) %>%
  layout(margin = list(r = 125, b = 100))
```


Random Effects {data-navmenu="Models"}
=====================================
Column
-------------------------------------

### Poisson Non-Spatial RE Model --- I.I.D. Random Effect

```{r m4re}
ggplot(kc_f, aes(x = long, y = lat, group = group, fill = m4RE, label = GISJOIN)) +
  geom_polygon(color = "white", lwd = .15) +
  scale_fill_distiller(name = "Median I.I.D. Random Effect", palette = "Purples", 
                       direction = 1, na.value = "grey80") + 
  coord_map() +
  theme_minimal() +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        panel.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(angle = 45)) +
  xlab("") +
  ylab("")
```

### 

Column
-------------------------------------

### Poisson Spatial RE Model --- I.I.D. Random Effect

```{r m5re}
ggplot(kc_f, aes(x = long, y = lat, group = group, fill = m5RE, label = GISJOIN)) +
  geom_polygon(color = "white", lwd = .15) +
  scale_fill_distiller(name = "Median I.I.D. Random Effect", palette = "Purples", 
                       direction = 1, na.value = "grey80") + 
  coord_map() +
  theme_minimal() +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        panel.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(angle = 45)) +
  xlab("") +
  ylab("")
```

### Poisson Spatial RE Model --- ICAR Spatial Effect

```{r m5se}
ggplot(kc_f, aes(x = long, y = lat, group = group, fill = m5SE, label = GISJOIN)) +
  geom_polygon(color = "white", lwd = .15) +
  scale_fill_distiller(name = "Median Spatial Random Effect", palette = "Purples", 
                       direction = 1, na.value = "grey80") + 
  coord_map() +
  theme_minimal() +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        panel.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(angle = 45)) +
  xlab("") +
  ylab("")
```


Fitted Values {data-navmenu="Models"}
=====================================

Column
-------------------------------------

### Poisson Non-Spatial RE Model --- Median Fitted Values

```{r m4fit}
#median fitted value (for faceted plot)
ggplot(kc_f, aes(x = long, y = lat, group = group, fill = log10(m4median), label = GISJOIN)) +
  geom_polygon(color = "white", lwd = .15) +
  scale_fill_distiller(name = "Median Posterior Prediction\nN Listings (Non-Spatial)", palette = "Purples", 
                       limits = c(0, 4),
                       breaks = c(1, 2, 3, 4),
                       labels = c("10", "100", "1000", "10000"),
                       direction = 1, na.value = "grey80") + 
  coord_map() +
  theme_minimal() +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        panel.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(angle = 45)) +
  xlab("") +
  ylab("")
```

### Poisson Non-Spatial RE Model --- 95% Interval

```{r m4cred}
#95% credible interval width
ggplot(kc_f, aes(x = long, y = lat, group = group, fill = log10(m4post95wid), label = GISJOIN)) +
  geom_polygon(color = "white", lwd = .15) +
  scale_fill_distiller(name = "95% Credible Int. Width\n(Non-Spatial)", palette = "Purples", 
                       limits = c(0, 4),
                       breaks = c(1, 2, 3, 4),
                       labels = c("10", "100", "1000", "10000"),
                       direction = 1, na.value = "grey80") +  
  coord_map() +
  theme_minimal() +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        panel.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(angle = 45)) +
  xlab("") +
  ylab("") 
```


Column
-------------------------------------

### Poisson Spatial RE Model --- Median Fitted Values

```{r m5fit}
#median fitted value (for faceted plot)
ggplot(kc_f, aes(x = long, y = lat, group = group, fill = log10(m2median), label = GISJOIN)) +
  geom_polygon(color = "white", lwd = .15) +
  scale_fill_distiller(name = "Median Posterior Prediction\nN Listings (Spatial)", palette = "Purples", 
                       limits = c(0, 4),
                       breaks = c(1, 2, 3, 4),
                       labels = c("10", "100", "1000", "10000"),
                       direction = 1, na.value = "grey80") + 
  coord_map() +
  theme_minimal() +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        panel.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(angle = 45)) +
  xlab("") +
  ylab("")
```

### Poisson Spatial RE Model --- 95% Interval

```{r m5cred}
#95% credible interval width
ggplot(kc_f, aes(x = long, y = lat, group = group, fill = log10(m5post95wid), label = GISJOIN)) +
  geom_polygon(color = "white", lwd = .15) +
  scale_fill_distiller(name = "95% Credible Int. Width\n(Spatial)", palette = "Purples", 
                       limits = c(0, 4),
                       breaks = c(1, 2, 3, 4),
                       labels = c("10", "100", "1000", "10000"),
                       direction = 1, na.value = "grey80") +  
  coord_map() +
  theme_minimal() +
  theme(axis.ticks.y = element_blank(),axis.text.y = element_blank(), # get rid of x ticks/text
        axis.ticks.x = element_blank(),axis.text.x = element_blank(), # get rid of y ticks/text
        panel.background = element_blank(),
        panel.grid = element_blank(),
        legend.position = "bottom",
        legend.text = element_text(angle = 45)) +
  xlab("") +
  ylab("") 
```