---
title: "Error and mapping ACS data"
author: "Jenna"
date: "September 19, 2020"
output: word_document
---

```{r setup, include=FALSE}
library(rmarkdown)
knitr::opts_chunk$set(echo = FALSE, warning = FALSE, message = FALSE)

```

### Objective   

Error for small-area geography data can be very high. For example, data on the poverty rate from 2014-2018 ACS has a 15% coefficient of variation (CV) or below for only 105/1126 census tracts in Louisiana. How can we address this issue in order map the most accurate data possible?  
   
### Background   

In 2012, Ben Horowitz (then a TDC employee) [evaluated methodologies for smoothing ACS data](https://gnocdc.s3.amazonaws.com/reports/GNOCDC_ACS_Mapping_Methodology.pdf) (by block group) within New Orleans. His analysis concluded that averaging 'true' neighboring census tracts' data from the 2007 - 2011 ACS produced estimates that were closer to 2010 Census data for the same geographies than three other methods. 'True' neighbors are only those neighbors that actually touch; i.e., are not separated by a body of water. 

There are some concerns about directly applying this methodology to the project at hand, our statewide, maps:   
1. Identifying true neighbors statewide would likely be labor-intensive.   
2. The error of the final measurements wasn't evaluated in this method, and we won't have current decennial data to compare our smoothed estimates to.   

Tools exist in R and other software to apply different forms of smoothing to geographic data.   

[Empirical Bayes smoothing](https://geodacenter.github.io/workbook/3b_rates/lab3b.html#empirical-bayes-eb-smoothed-rate-map) shrinks estimates towards the overall average. Small-population geographies will be more affected by this method, as they are more likely to be spurious outliers. 

[Spatial average smoothing](https://spatialanalysis.github.io/lab_tutorials/Applications_of_Spatial_Weights.html#spatial-rate-smoothing) is the same method Ben used - each tract's smoothed estimate is an average of its raw data with all of its neighbors' raw data.  

One researcher, Spielman, [developed an interesting algorithm](https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0115626) that perfoms what he calls "regionalization" - the smaller geographies are aggregated one unit at a time until that new region passes a user-defined reliability threshold, such as a CV of <= 15%. Open source code of the algorithm is available, but implementing it is tricky. 

An example of how this problem is sometimes addressed given that accurate data rarely exists for small-area geographies is discussed in another [paper by Spielman](https://www.tandfonline.com/doi/full/10.1080/00045608.2015.1052335): "[O]ne way to overcome variable-level uncertainty in the ACS is to focus on a composite of multiple variables, rather than examining individual variables independently. This “contextual” approach has a long tradition in the social scientific literature, extending to the early twentieth century." I point this out not to suggest that we give up on reporting single small-area statistics but to show that this is a difficult and long-standing problem with no easy solution. 

At the Population Division of NYC planning, maps are colored according to bins (e.g. 0 - 5% poverty; 5 - 10% poverty, etc.) that are tested for their reliability with [an interesting tool](https://atcoordinates.info/2017/12/17/the-map-reliability-calculator-for-classifying-acs-data/). Taking into account the MOE of each measurement within the bin, a reliable map is one in which 10% or fewer of the measurements might actually fall outside of their bin. Unfortunately, when I input our LA tract-level poverty data, Having 2 breaks (so 3 categories) was considered reliable only if I filtered out tracts with a CV <=30%. That is to say, this approach may be better suited the urban data. 

Another approach is one we've discussed before: [use shading/patterns](http://sites.tufts.edu/gis/files/2013/11/Amercian-Community-Survey_Margin-of-error-tutorial.pdf) to show reliability of each geography's measurement with an overlay on the map.   

Let's look at two of the methods discussed above for tract-level poverty data. First, we look at the raw data and its MOE at the state level and in New Orleans (so we can consider the data in the context of our own local knowledge). Then we see the same geographies with the Empirical Bayes' smoothing and the simple average smoothing, as well as the difference between these smoothed estimates and their raw counterparts. Finally,we see the impact that simple average smoothing has on the data's error.

I recommend opening the document in two separate windows so you can keep the raw data maps handy for comparing.

```{r}
library(here)
source(here("sandbox/jenna/libraries jenna.R"))

tracts.la <- sf::st_read(here("inputs/tl_2010_22_tract10/tl_2010_22_tract10.shp")) 
tracts.la.OPclipped <- sf::st_read(here("inputs/clippedriver/clippedriver.shp")) 
```


#### Raw data   

```{r}
povTractRaw <- get_acs(geography = "tract",
                      year = 2018,
                      variables = c('B17001_001E','B17001_002E'), 
                      state = c("LA"),
                      survey = "acs5")


povTract <- povTractRaw %>%
  dplyr::select(-NAME) %>%
  gather(-variable, -GEOID, key = type, value = val) %>%
  mutate(vartype = paste0(variable, str_sub(type,1,1))) %>%
  dplyr::select(-variable, -type) %>%
  spread(key = vartype,value = val) %>%
  mutate(pov = `B17001_002e`/`B17001_001e`,
         pov = ifelse(is.nan(pov), 0, pov)) %>%
  mutate(povmoe = moeprop(`B17001_001e`,`B17001_002m`,`B17001_001m`,pov)) %>%
  rename(pop = `B17001_001e`, popMOE = `B17001_001m`, povPop = `B17001_002e`, povPopMOE = `B17001_002m`) %>%
  filter(pop != 0)
```
```{r}
tracts.la %>% 
  left_join(povTract, by= c("GEOID10" = "GEOID")) %>%
  ggplot() +
  geom_sf(aes(fill = pov), color = "gray70", size = .05)+
  themeDC_map() +
  # scale_color_brewer(palette =  "Spectral") +
  # scale_fill_brewer(palette = "Spectral", direction = 1,  na.value = "gray70") +
  scale_fill_viridis_c(limits = c(0,.85)) +
  coord_sf(datum = NA) +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Raw poverty rate",  
       subtitle  =  "Louisiana, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/Louisiana raw pov.png", type = "cairo",
         width = 5, height = 4, units = "in")

```

```{r}
tracts.la %>% 
  left_join(povTract, by= c("GEOID10" = "GEOID")) %>%
  ggplot() +
  geom_sf(aes(fill = povmoe), color = "gray70", size = .05)+
  themeDC_map() +
  # scale_color_brewer(palette =  "Spectral") +
  # scale_fill_brewer(palette = "Spectral", direction = 1,  na.value = "gray70") +
  scale_fill_viridis_c() +
  coord_sf(datum = NA) +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Poverty rate MOE",  
       subtitle  =  "Louisiana, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/Louisiana pov moe.png", type = "cairo",
         width = 5, height = 4, units = "in")

```

```{r}
tracts.la %>% 
  left_join(povTract, by= c("GEOID10" = "GEOID")) %>%
  filter(str_sub(GEOID10,1,5)=="22071") %>%
  ggplot() +
  geom_sf(aes(fill = pov), color = "gray70", size = .05)+
  themeDC_map() +
  scale_fill_viridis_c(limits = c(0,.85)) +
  coord_sf(datum = NA) +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Raw poverty rate",  
       subtitle  =  "New Orleans, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/New Orleans raw pov.png", type = "cairo",
         width = 5, height = 4, units = "in")
```

```{r}
tracts.la %>% 
  left_join(povTract, by= c("GEOID10" = "GEOID")) %>%
  filter(str_sub(GEOID10,1,5)=="22071") %>%
  ggplot() +
  geom_sf(aes(fill = povmoe), color = "gray70", size = .05)+
  themeDC_map() +
  scale_fill_viridis_c() +
  coord_sf(datum = NA) +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Poverty rate MOE",  
       subtitle  =  "New Orleans, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/New Orleans pov moe.png", type = "cairo",
         width = 5, height = 4, units = "in")
```

#### Empirical Bayes smoothing   

```{r}
library(spdep)
q1weights <- read.gal("povTract_nopop0_q1.gal", povTract$GEOID)

```

```{r}
povEBtest <- EBlocal(ri = povTract$povPop, ni = povTract$pop, nb = q1weights)

povEBmatchtest <- bind_cols(povTract, povEBtest) %>%
  mutate(dif = est-raw)
```

```{r}
tracts.la %>% 
  left_join(povEBmatchtest, by= c("GEOID10" = "GEOID")) %>% 
  ggplot() +
  geom_sf(aes(fill = est), color = "gray70", size = .05)+
  coord_sf(datum = NA) +
  themeDC_map() +
  themeDC_map() +
  scale_fill_viridis_c(limits = c(0,.85)) +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Empirical Bayes smoothed pov data",  
       subtitle  =  "Louisiana, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/Louisiana empirical bayes smoothed pov.png", type = "cairo",
         width = 5, height = 4, units = "in")
```

```{r}
tracts.la %>% 
  left_join(povEBmatchtest, by= c("GEOID10" = "GEOID")) %>% 
  ggplot() +
  geom_sf(aes(fill = dif), color = "gray70", size = .05)+
  coord_sf(datum = NA) +
  themeDC_map() +
  themeDC_map() +
  scale_fill_viridis_c() +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Difference btwn raw and Empirical Bayes smoothed pov data \n(smoothed - raw)",  
       subtitle  =  "Louisiana, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/Louisiana empirical bayes smoothed pov DIFF.png", type = "cairo",
         width = 5, height = 4, units = "in")
```

```{r}
tracts.la %>% 
  left_join(povEBmatchtest, by= c("GEOID10" = "GEOID")) %>% 
  # st_as_sf(.) %>%
  filter(str_sub(GEOID10, 1, 5) == "22071") %>%
  ##Filtering for census tracts that are water in Orleans Parish
  # filter(GEOID != "22071980100") %>% 
  # filter(GEOID != "22071990000") %>% 
  # filter(GEOID != "22071980000") %>% 
  ggplot() +
  geom_sf(aes(fill = est), color = "gray70", size = .05)+
  coord_sf(datum = NA) +
  themeDC_map() +
  themeDC_map() +
  scale_fill_viridis_c(limits = c(0,.85)) +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Empirical Bayes smoothed pov data",  
       subtitle  =  "New Orleans, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/New Orleans empirical bayes smoothed pov.png", type = "cairo",
         width = 5, height = 4, units = "in")

tracts.la %>% 
  left_join(povEBmatchtest, by= c("GEOID10" = "GEOID")) %>% 
  # st_as_sf(.) %>%
  filter(str_sub(GEOID10, 1, 5) == "22071") %>%
  ##Filtering for census tracts that are water in Orleans Parish
  # filter(GEOID10 != "22071980100") %>%
  # filter(GEOID10 != "22071990000") %>%
  # filter(GEOID10 != "22071980000") %>%
  ggplot() +
  geom_sf(aes(fill = dif), color = "gray70", size = .05)+
  coord_sf(datum = NA) +
  themeDC_map() +
  scale_fill_viridis_c() +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Difference btwn raw and Empirical Bayes smoothed pov data \n(smoothed - raw)",  
       subtitle  =  "New Orleans, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/New Orleans empirical bayes smoothed pov DIFF.png", type = "cairo",
         width = 5, height = 4, units = "in")
```   

##### Another way to look at the effect of the smoothing method:   

```{r}
ggplot(povEBmatchtest, aes(x = dif)) +
  geom_histogram() +
  labs(title = "Impact of Empirical Bayes smoothing on statewide, tract-level data:",
       subtitle = "Histogram of smoothed - raw data by tract")
```

#### Simple average smoothing   

```{r}

q1weights_listw <- nb2listw(q1weights)

lag_povPop <- lag.listw(q1weights_listw,povTract$povPop)
lag_pop <- lag.listw(q1weights_listw,povTract$pop)
povTract$R_SPAT_R <- lag_povPop/ lag_pop  
povTract <- povTract %>% mutate(rspatdiff = R_SPAT_R - pov)
```


```{r}
tracts.la %>% 
  left_join(povTract, by= c("GEOID10" = "GEOID")) %>%
  ggplot() +
  geom_sf(aes(fill = R_SPAT_R), color = "gray70", size = .05)+
  themeDC_map() +
  # scale_color_brewer(palette =  "Spectral") +
  # scale_fill_brewer(palette = "Spectral", direction = 1,  na.value = "gray70") +
  scale_fill_viridis_c(limits = c(0,.85)) +
  coord_sf(datum = NA) +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Simple average of neighboring tracts smoothed poverty rate",  
       subtitle  =  "Louisiana, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/Louisiana simple smoothe pov.png", type = "cairo",
         width = 5, height = 4, units = "in")
```

```{r}
tracts.la %>% 
  left_join(povTract, by= c("GEOID10" = "GEOID")) %>%
  ggplot() +
  geom_sf(aes(fill = rspatdiff), color = "gray70", size = .05)+
  themeDC_map() +
  # scale_color_brewer(palette =  "Spectral") +
  # scale_fill_brewer(palette = "Spectral", direction = 1,  na.value = "gray70") +
  scale_fill_viridis_c() +
  coord_sf(datum = NA) +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Difference btwn raw and simple avg smoothed pov rate \n(smoothed - raw)",  
       subtitle  =  "Louisiana, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/Louisiana simple smoothe pov DIFF.png", type = "cairo",
         width = 5, height = 4, units = "in")
```

```{r}
tracts.la %>% 
  left_join(povTract, by= c("GEOID10" = "GEOID")) %>%
  filter(str_sub(GEOID10,1,5)=="22071") %>%
  # st_as_sf(.) %>%
  # filter(str_sub(GEOID, 1, 5) == "22071") %>% 
  ##Filtering for census tracts that are water in Orleans Parish
  # filter(GEOID != "22071980100") %>% 
  # filter(GEOID != "22071990000") %>% 
  # filter(GEOID != "22071980000") %>% 
  ggplot() +
  geom_sf(aes(fill = R_SPAT_R), color = "gray70", size = .05)+
  themeDC_map() +
  # scale_color_brewer(palette =  "Spectral") +
  # scale_fill_brewer(palette = "Spectral", direction = 1,  na.value = "gray70") +
  scale_fill_viridis_c(limits = c(0,.85)) +
  coord_sf(datum = NA) +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Simple average of neighboring tracts smoothed poverty rate \n(same legend as others)",  
       subtitle  =  "New Orleans, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/New Orleans simple smoothe pov.png", type = "cairo",
         width = 5, height = 4, units = "in")
```

```{r}
tracts.la %>% 
  left_join(povTract, by= c("GEOID10" = "GEOID")) %>%
  filter(str_sub(GEOID10,1,5)=="22071") %>%
  # st_as_sf(.) %>%
  # filter(str_sub(GEOID, 1, 5) == "22071") %>% 
  ##Filtering for census tracts that are water in Orleans Parish
  # filter(GEOID != "22071980100") %>% 
  # filter(GEOID != "22071990000") %>% 
  # filter(GEOID != "22071980000") %>% 
  ggplot() +
  geom_sf(aes(fill = R_SPAT_R), color = "gray70", size = .05)+
  themeDC_map() +
  # scale_color_brewer(palette =  "Spectral") +
  # scale_fill_brewer(palette = "Spectral", direction = 1,  na.value = "gray70") +
  scale_fill_viridis_c() +
  coord_sf(datum = NA) +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Simple average of neighboring tracts smoothed poverty rate \n(own legend, to show greater detail)",  
       subtitle  =  "New Orleans, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/New Orleans simple smoothe pov.png", type = "cairo",
         width = 5, height = 4, units = "in")
```


```{r}

tracts.la %>% 
  left_join(povTract, by= c("GEOID10" = "GEOID")) %>%
  filter(str_sub(GEOID10,1,5)=="22071") %>%
  ggplot() +
  geom_sf(aes(fill = rspatdiff), color = "gray70", size = .05)+
  themeDC_map() +
  scale_fill_viridis_c() +
  coord_sf(datum = NA) +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Difference btwn raw and simple avg smoothed pov rate \n(smoothed - raw)",  
       subtitle  =  "New Orleans, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/New Orleans simple smoothe pov DIFF.png", type = "cairo",
         width = 5, height = 4, units = "in")
```

##### Another way to look at the effect of the smoothing method:   

```{r}
ggplot(povTract, aes(x = rspatdiff)) +
  geom_histogram() +
  labs(title = "Impact of simple average smoothing on statewide, tract-level data:",
       subtitle = "Histogram of smoothed - raw data by tract")
```

#### Simple average smoothing and error

```{r}
povShp <- tracts.la %>%
  right_join(povTract, by = c("GEOID10" = "GEOID")) 

pov.sf <- st_as_sf(povShp, coords = c("INTPLON10", "INTPLAT10"))

st_queen <- function(a, b = a) st_relate(a, b, pattern = "F***T****")
as.nb.sgbp <- function(x, ...) {
  attrs <- attributes(x)
  x <- lapply(x, function(i) { if(length(i) == 0L) 0L else i } )
  attributes(x) <- attrs
  class(x) <- "nb"
  x
}
queen.sgbp <- st_queen(pov.sf)
queen.nb <- as.nb.sgbp(queen.sgbp)
queen.weights <- nb2listw(queen.nb)

aggs <- data.frame(popAgg = numeric(), povPopAgg = numeric(), popMOEagg = numeric(), povPopMOEagg = numeric())
for (i in 1:dim(povTract)[1]){
  nbs = queen.weights$neighbours[i]
  m <- data.frame(popAgg = numeric(), povPopAgg = numeric(), popMOEagg = numeric(), povPopMOEagg = numeric())
  for (j in 1:lengths(nbs)) {
    m <- bind_rows(m, data.frame(valPop = povTract$pop[nbs[[1]][j]], 
                                 valpovPop = povTract$povPop[nbs[[1]][j]],
                                 valPopMOE = povTract$popMOE[nbs[[1]][j]], 
                                 valpovPopMOE = povTract$povPopMOE[nbs[[1]][j]]))
  }
  new = as.data.frame(m) %>%
    summarise(popAgg = sum(valPop), povPopAgg = sum(valpovPop), popMOEagg = sqrt(sum(valPopMOE^2)), povPopMOEagg = sqrt(sum(valpovPopMOE^2)))
  aggs <- bind_rows(aggs, new)
}

aggs.err <- aggs %>%
  mutate(aggPctPov = povPopAgg/popAgg) %>%
  mutate(moeProp = moeprop(popAgg, povPopMOEagg, popMOEagg, aggPctPov)) %>%
  mutate(cvAgg = (moeProp/1.645)/aggPctPov)

errcheck <- bind_cols(povTract, aggs.err) %>%
  mutate(cvRaw = (povmoe/1.645)/pov,
         cvDiff = cvRaw - cvAgg,
         cvProp = cvAgg/cvRaw)
```

```{r}

tracts.la %>% 
  left_join(errcheck, by= c("GEOID10" = "GEOID")) %>%
  ggplot() +
  geom_sf(aes(fill = cvRaw), color = "gray70", size = .05)+
  themeDC_map() +
  scale_fill_viridis_c(limits = c(0,.65)) +
  coord_sf(datum = NA) +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Raw poverty rate CV",  
       subtitle  =  "Louisiana, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/Louisiana raw cv.png", type = "cairo",
         width = 5, height = 4, units = "in")
```

```{r}

tracts.la %>% 
  left_join(errcheck, by= c("GEOID10" = "GEOID")) %>%
  ggplot() +
  geom_sf(aes(fill = cvAgg), color = "gray70", size = .05)+
  themeDC_map() +
  scale_fill_viridis_c(limits = c(0,.65)) +
  coord_sf(datum = NA) +
  theme(legend.position = "right",
        plot.title=element_text(size=10, hjust=0.5, vjust=-1, family = "Asap", color ="#002F44"),
        plot.subtitle=element_text(size=9, hjust=0.5, face="bold", color="#002F44", family = "Asap"),
        plot.caption =element_text(size=9, hjust=0, face="italic", color = "#002F44", family = "Asap"),
        legend.text=element_text(size=9, family = "Asap", color = "#002F44"),
        legend.title= element_blank()) +
  labs(title = "Simple average smothed poverty rate CV",  
       subtitle  =  "Louisiana, 2018 ACS 5") +
  ggsave(filename = "outputs/smoothing/Louisiana avg smoothed cv.png", type = "cairo",
         width = 5, height = 4, units = "in")
```
#### For Orleans, water clipped vs. unclipped simple average smoothing and error

```{r}

povTract.op <- povTract %>%
  filter(str_sub(GEOID,1,5) == "22071")
povShp.clp <- tracts.la.OPclipped %>%
  right_join(povTract.op, by = c("GEOID10" = "GEOID"))

pov.sf.clp <- st_as_sf(povShp.clp, coords = c("INTPLON10", "INTPLAT10"))


queen.sgbp.clp <- st_queen(pov.sf.clp)
queen.nb.clp <- as.nb.sgbp(queen.sgbp.clp)
queen.weights.clp <- nb2listw(queen.nb.clp)

aggs.clp <- data.frame(popAgg = numeric(), povPopAgg = numeric(), popMOEagg = numeric(), povPopMOEagg = numeric())
for (i in 1:dim(povTract.op)[1]){
  nbs = queen.weights.clp$neighbours[i]
  m <- data.frame(popAgg = numeric(), povPopAgg = numeric(), popMOEagg = numeric(), povPopMOEagg = numeric())
  for (j in 1:lengths(nbs)) {
    m <- bind_rows(m, data.frame(valPop = povTract.op$pop[nbs[[1]][j]], 
                                 valpovPop = povTract.op$povPop[nbs[[1]][j]],
                                 valPopMOE = povTract.op$popMOE[nbs[[1]][j]], 
                                 valpovPopMOE = povTract.op$povPopMOE[nbs[[1]][j]]))
  }
  new = as.data.frame(m) %>%
    summarise(popAgg = sum(valPop), povPopAgg = sum(valpovPop), popMOEagg = sqrt(sum(valPopMOE^2)), povPopMOEagg = sqrt(sum(valpovPopMOE^2)))
  aggs.clp <- bind_rows(aggs.clp, new)
}

aggs.err.clp <- aggs.clp %>%
  mutate(aggPctPov.clp = povPopAgg/popAgg) %>%
  mutate(moeProp.clp = moeprop(popAgg, povPopMOEagg, popMOEagg, aggPctPov.clp)) %>%
  mutate(cvAgg.clp = (moeProp.clp/1.645)/aggPctPov.clp)

errcheck.clp <- bind_cols(povTract.op, aggs.err.clp) %>%
  mutate(cvRaw = (povmoe/1.645)/pov,
         cvDiff = cvRaw - cvAgg,
         cvProp = cvAgg/cvRaw)
```


```{r}
povShp.op <- tracts.la %>%
  right_join(povTract.op, by = c("GEOID10" = "GEOID"))

pov.sf.op <- st_as_sf(povShp.op, coords = c("INTPLON10", "INTPLAT10"))


queen.sgbp.op <- st_queen(pov.sf.op)
queen.nb.op <- as.nb.sgbp(queen.sgbp.op)
queen.weights.op <- nb2listw(queen.nb.op)

aggs.op <- data.frame(popAgg = numeric(), povPopAgg = numeric(), popMOEagg = numeric(), povPopMOEagg = numeric())
for (i in 1:dim(povTract.op)[1]){
  nbs = queen.weights.op$neighbours[i]
  m <- data.frame(popAgg = numeric(), povPopAgg = numeric(), popMOEagg = numeric(), povPopMOEagg = numeric())
  for (j in 1:lengths(nbs)) {
    m <- bind_rows(m, data.frame(valPop = povTract.op$pop[nbs[[1]][j]], 
                                 valpovPop = povTract.op$povPop[nbs[[1]][j]],
                                 valPopMOE = povTract.op$popMOE[nbs[[1]][j]], 
                                 valpovPopMOE = povTract.op$povPopMOE[nbs[[1]][j]]))
  }
  new = as.data.frame(m) %>%
    summarise(popAgg = sum(valPop), povPopAgg = sum(valpovPop), popMOEagg = sqrt(sum(valPopMOE^2)), povPopMOEagg = sqrt(sum(valpovPopMOE^2)))
  aggs.op <- bind_rows(aggs.op, new)
}

aggs.err.op <- aggs.op %>%
  mutate(aggPctPov = povPopAgg/popAgg) %>%
  mutate(moeProp = moeprop(popAgg, povPopMOEagg, popMOEagg, aggPctPov)) %>%
  mutate(cvAgg = (moeProp/1.645)/aggPctPov)

errcheck.op <- bind_cols(povTract.op, aggs.err.op) %>%
  mutate(cvRaw = (povmoe/1.645)/pov,
         cvDiff = cvRaw - cvAgg,
         cvProp = cvAgg/cvRaw)
```

```{r}
clp.stattest <- bind_cols(aggs.err.clp, aggs.err.op, povTract.op) %>%
  select(GEOID, aggPctPov, moeProp,aggPctPov.clp, moeProp.clp) %>%
  mutate(dif = aggPctPov - aggPctPov.clp)

##stat testing for 2000 vs 201* data
##input: columns of estimates and their MOEs (zeros for Census 2000)
##output: column of yes or no if significant
stattest <- function(x, moex = matrix(0, length(x)), y, moey){
  significant <- matrix(0, length(x))
  v <- abs((x-y)/sqrt((moex/1.645)^2+(moey/1.645)^2))
  significant <- ifelse(v>1.645,"yes","no")
  return(as.list(significant))
}

clp.stattest.tested <- stattest(clp.stattest$aggPctPov, clp.stattest$moeProp, clp.stattest$aggPctPov.clp, clp.stattest$moeProp.clp)
```

### Conclusion

Actually measuring the error of the output of these methods is not a well-documented area of this research, as far as I saw. However, it's easy enough to do ourselves!   

The improvement was huge - 65% of tracts saw a decrease of 10 percentage points or more from their raw CV to their smoothed CV using the simple averaging method. However, about 90 tracts still have CVs >=15%. One possibilty is that we combine methods - we map smoothed values with pattern over tracts that are less reliable. We could also just drop the ~8% of tracts that don't meet the 15% cutoff. 

