Jail_Prison_analysis_R.Rmd

```{r}
library(tidyverse)
```

nyt_data: county-level COVID statistics
jp: jail population, surrounding county population
full data: nyt_data and jp merged; jail population, cases, deaths for each county over time
```{r}
nyt_data = read_csv('NYT_data/us-counties.csv') # https://github.com/nytimes/covid-19-data
jp = read_csv('https://raw.githubusercontent.com/vera-institute/jail-population-data/master/jail_population.csv')
jp <- jp %>% rename(county = county_name) %>% rename(state=state_name) %>% arrange(date)
jp$county<-str_remove(jp$county," County")
jp
nyt_data
full_data <- merge(jp,nyt_data) %>% 
  select(date,county,state,fips,jail_population,resident_population,cases,deaths,jail_incarceration_rate_per_100k,urbanicity)
#full_data <- merge(jp,nyt_data,by = c('date','county','state')) %>% select(date,county,state,fips,jail_population,resident_population,cases,deaths,jail_incarceration_rate_per_100k,urbanicity)
full_data %>% arrange(county,state)

write.csv(full_data,"jail_county_df.csv")
```
Cluster: predict covid rates, average increase? something as a function of incaraceration rate, urbanicity, jail population

recent: most recent statistics
```{r}
full_data_recent <- full_data %>% filter(date == as.Date("2020-05-11"))
ggplot(data = full_data_recent) + 
  geom_point(aes(log(jail_incarceration_rate_per_100k),log(cases),color=log(resident_population)))+
  scale_color_gradient(low="blue", high="red")
```


% difference in jail population and cases from March 1st to now
```{r}
full_data %>% arrange((date)) 
diff_df <- full_data %>% arrange((date)) %>% filter(date >= as.Date("2020-03-01")) %>% group_by(county,state) %>%
  mutate(frac_jail_red = (first(jail_population)-last(jail_population))/first(jail_population)) %>%
  mutate(frac_increase_cases = (last(cases)-first(cases))/first(cases)) %>%
  select(county,state,frac_jail_red,frac_increase_cases) %>% distinct %>% arrange(desc(frac_increase_cases))
diff_df
```

Relationship between % difference in cases and % difference in jail population?
```{r}
diff_df
ggplot(data = diff_df) + geom_point(aes((log(frac_jail_red)),log(frac_increase_cases)))
```


```{r}
jail_counties <- jp%>% select(county,state) %>% unique()
jail_counties$countystate <- paste(jail_counties$county, jail_counties$state, sep=", ")
nyt_data$countystate <- paste(nyt_data$county, nyt_data$state, sep=", ")
```

Look at a particular county
```{r}
countyname = 'Clackamas'
statename = 'Oregon'
county_data <- full_data %>% filter(county==countyname,state==statename)%>% select(county,date,jail_population,cases)
county_data$release_rate <- (lag(county_data$jail_population,3)-county_data$jail_population)/lag(county_data$jail_population,3)*500
county_data
library(reshape2)
ggplot(data=melt(county_data, id=c('date','county'))) + geom_line(aes(date,value,color=variable))
```

Look at an entire state
```{r}
statename = 'Oregon'
state_data <- full_data %>% filter(state==statename) %>% group_by(date) %>% mutate(total_jailed = sum(jail_population)) %>%
  mutate(total_cases = sum(cases)) %>% select(date,total_jailed,total_cases) %>% unique()
state_data
ggplot(data=melt(state_data, id=c('date'))) + geom_line(aes(date,value,color=variable))
```

Prison analysis
prison_df: state, prison name, and county
```{r}
prisons <- read_csv('prison_boundaries.csv')
prisons$VAL_DATE <- as.Date(prisons$VAL_DATE)
prisons$COUNTY <- tolower(prisons$COUNTY)
prisons <- prisons %>% rename(county = COUNTY,Code=STATE)
str_list <- paste(c('CORRECTION','PRISON'), collapse="|")
prisons <- prisons %>% filter(grepl(str_list,NAME)) %>%
  select(NAME,Code,county,POPULATION,CAPACITY) %>% 
  mutate(POPULATION = ifelse(POPULATION<0, NA, POPULATION)) %>%  # convert negative values to NA
  mutate(CAPACITY = ifelse(CAPACITY<0, NA, CAPACITY)) %>%  # convert negative values to NA
  mutate(fullness = POPULATION/CAPACITY)

statenames <- read_csv('state_abbreviations.csv') %>% select(State,Code)
prison_df <- merge(statenames,prisons) %>% select(State,NAME,county) %>% unique()
prison_df <- merge(statenames,prisons) %>% 
  select(State,NAME,county, POPULATION, CAPACITY,fullness) %>% 
  unique() %>% group_by(State,county) %>%
  mutate(prison_pop = sum(POPULATION)) %>%     # sum over all prisons in a county
  mutate(prison_capacity = sum(CAPACITY)) %>%  # sum over all prisons in a county
  mutate(prison_fullness = mean(fullness)) %>%  # calculate average prison fullness
  select(State,county,prison_pop,prison_capacity,prison_fullness) %>% unique()
prison_df
```

pr: state, fraction of prison population released
prison_counties_frac_released: counties with prisons, fraction of prisoners released statewide
```{r}
pr <- read_csv('covid_prison_releases.csv') %>%
  rename(prior_pop = 'Population Prior to Releases',releases = "Overall Pop. Reduction / \nTotal Number of Releases") %>% 
  select(State,Facility,prior_pop,releases) %>% 
  filter(Facility=='Statewide') %>%
  mutate(releases = as.numeric(releases)) %>%
  mutate(frac_released = releases/prior_pop) %>%
  select(State,frac_released) %>% drop_na()
pr

prison_counties_frac_released <- merge(prison_df, pr,by='State',all.x = TRUE) %>%
  mutate(frac_released = ifelse(is.na(frac_released),0,frac_released)) %>%
  mutate(countystate = paste(county,State,sep=", ")) %>%
  #select(countystate,frac_released) %>% unique()
  select(countystate,frac_released,prison_pop,prison_capacity,prison_fullness) %>%unique()
prison_counties_frac_released
```


GLMs

Model 1:
Restricted to counties containing a prison
difference in cases (cases from today - cases from 7 days ago) ~ total county population, fraction of prisoners released statewide, county population density

Model 2: 
Considering all counties
difference in cases (cases from today - cases from 7 days ago) ~ total county population, indicator function (contains a prison = 1, doesn't contain a prison = 0), county population density

Model 3:
Considering all counties 
difference in cases (cases from today - cases from 7 days ago) ~ total county population, prison indicator function * fraction released, prison indicator function, county population density
```{r}
county_pop <- read_csv('county_demographics_df.csv')%>% mutate(countystate = paste(tolower(county),state,sep = ", ")) #%>%
  #select(countystate,tot_pop,population_density)

nytdf <- nyt_data %>% mutate(countystate = paste(tolower(county),state,sep=", ")) %>%
  mutate(prison_county = ifelse(countystate %in% prison_counties_frac_released$countystate,1,0)) %>% unique()

nyt_with_prisons <- merge(nytdf, prison_counties_frac_released, all.x=TRUE) %>% group_by(countystate) %>%
  arrange(countystate,date) %>%
  mutate(diff_cases = cases-lag(cases,7)) %>%
  filter(date >= as.Date('2020-03-01')) %>%
  filter(!is.na(diff_cases)) %>%
  #select(countystate,date,cases,prison_county,frac_released,diff_cases)
  select(countystate,fips,date,cases,prison_county,frac_released,diff_cases,prison_pop,prison_capacity,prison_fullness)

#nyt_with_prisons

prison_pop_df <- merge(nyt_with_prisons,county_pop,all.x=TRUE) #%>% mutate(cases_per_capita = cases/tot_pop) %>%
  #mutate(cases_per_capita_7_days_ago = lag(cases_per_capita,7))
#nyt_with_prisons
write.csv(prison_pop_df,"prison_county_df.csv")
prison_pop_df

prison_counties <- prison_pop_df %>% filter(prison_county==1)
#prison_counties

length(unique(prison_pop_df$countystate))
length(unique(prison_counties$countystate))
nrow(prison_counties)

model1 <- glm(diff_cases ~ tot_pop + frac_released + population_density, data = prison_counties)
summary(model1)

model2 <- glm(diff_cases ~ tot_pop + prison_county + population_density, data = prison_pop_df)
summary(model2)

prison_pop_df_with_indicator <-prison_pop_df %>% mutate(frac_released_for_indicator = ifelse(is.na(frac_released),0,frac_released)) %>%
  select(countystate,date,prison_county,frac_released_for_indicator,diff_cases,tot_pop,population_density) %>% drop_na() %>%
  mutate(int_prison_county_frac_released = prison_county*frac_released_for_indicator)

model3 <- glm(diff_cases ~ tot_pop + int_prison_county_frac_released + prison_county + population_density ,data = prison_pop_df_with_indicator)
summary(model3)

with(summary(model1), 1 - deviance/null.deviance)
with(summary(model2), 1 - deviance/null.deviance)
with(summary(model3), 1 - deviance/null.deviance)
```

Model 4:
Considering all counties 
difference in cases (cases from today - cases from 7 days ago) ~ total county population, prison indicator function * fraction released, prison indicator function, county population density, high risk age group

Model 5:
Restricted to counties containing a prison
difference in cases (cases from today - cases from 7 days ago) ~ total county population, fraction of prisoners released statewide, county population density, prison fullness

Model 6:
Restricted to counties containing a prison
difference in cases (cases from today - cases from 7 days ago) ~ total county population, fraction of prisoners released statewide, county population density, prison fullness, high risk age group

Model 7:
Restricted to counties containing a prison
difference in cases (cases from today - cases from 7 days ago) ~ total county population, fraction of prisoners released statewide, county population density, high risk age group

```{r}
prison_pop_df_with_indicator <-prison_pop_df %>% mutate(frac_released_for_indicator = ifelse(is.na(frac_released),0,frac_released)) %>%
  mutate(int_prison_county_frac_released = prison_county*frac_released_for_indicator)

model4 <- glm(diff_cases ~ tot_pop + int_prison_county_frac_released + prison_county + population_density + highrisk_agegroup_perc,data = prison_pop_df_with_indicator)
model5 <- glm(diff_cases ~ tot_pop + frac_released + population_density + prison_fullness, data = prison_counties)
model6 <- glm(diff_cases ~ tot_pop + frac_released + population_density + prison_fullness + highrisk_agegroup_perc, data = prison_counties)
model7 <- glm(diff_cases ~ tot_pop + frac_released + population_density + highrisk_agegroup_perc, data = prison_counties)
model8 <- glm(diff_cases ~ tot_pop + frac_released + population_density + prison_capacity, data = prison_counties)

summary(model8)
summary(model7)
summary(model6)
summary(model5)
summary(model4)

extractAIC(model7)
extractAIC(model6)
extractAIC(model5)
extractAIC(model4)

with(summary(model7), 1 - deviance/null.deviance)
with(summary(model6), 1 - deviance/null.deviance)
with(summary(model5), 1 - deviance/null.deviance)
with(summary(model4), 1 - deviance/null.deviance)

```



Some data visualization, consider the mean change in cases as a function of fraction released, population density
```{r}
prison_pop_df_today <- prison_pop_df  %>% group_by(countystate) %>%
  mutate(mean_diff_cases = mean(diff_cases)) %>% filter(date == as.Date('2020-05-11')) 
prison_pop_df_today
ggplot(prison_pop_df_today) + geom_point(aes(log(population_density),log(mean_diff_cases),color=prison_county))

prison_counties_today <- prison_counties  %>% group_by(countystate) %>%
  mutate(mean_diff_cases = mean(diff_cases)) %>% filter(date == as.Date('2020-05-11'))%>%
  mutate(log_frac_released = ifelse(is.infinite(log(frac_released)),-9,log(frac_released)))
prison_counties_today
ggplot(prison_counties_today) + geom_point(aes(log(population_density),log(mean_diff_cases),color=(log_frac_released)))
```




```{r}
str_list <- paste(c('CORRECTION','PRISON'), collapse="|")
prison_counties <- read_csv('prison_boundaries.csv') %>% filter(grepl(str_list,NAME))
prison_county_fips <- prison_counties$COUNTYFIPS

data <- read_csv('county_populations.csv') %>% mutate(fips = paste(STATE,COUNTY,sep="")) %>% select(fips) %>%
  mutate(prison_county = ifelse(fips %in% prison_county_fips,1,0)) %>% mutate(prison_county = as.factor(prison_county))
data

library(usmap)

plot_usmap(data = data,values="prison_county", size = .11) +
  scale_fill_manual(name='Prison County',labels = c("No Prison","Prison"),values=c("white","black"))
 
state_fips <- read_csv('NYT_data/us-states.csv') %>% select(state,fips) %>% rename(State = state)
data2 <- merge(state_fips,pr,by='State',all.x=TRUE) %>% unique() %>% 
  mutate(frac_released = ifelse(is.na(frac_released),0,frac_released)) %>%
  select(fips,frac_released)
data2

plot_usmap(data = data2,values="log_frac_released", size = .11,labels=FALSE) +
    scale_fill_continuous(low = "white", high = "darkgreen", name = "Fraction of prisoners released") + theme(legend.position = "right")

plot_usmap(data = data2,values="frac_released", size = .11,labels=FALSE) +
    scale_fill_gradient(low='white',high='darkgreen',name = "Fraction of prisoners released",trans='sqrt') + theme(legend.position = "right")
```
```{r}
county_populations <- read_csv('county_populations.csv') %>% mutate(fips = paste(STATE,COUNTY,sep="")) %>% select(fips,CENSUS2010POP)

county_covid_data <- nyt_data %>% group_by(fips) %>% filter(date == max(date)) %>% select(fips,cases)
county_covid_data <- merge(county_populations,county_covid_data) %>% mutate(infection_rate = cases/CENSUS2010POP) %>% select(fips,infection_rate)
plot_usmap(data = county_covid_data,values="infection_rate", size = .11) +
    scale_fill_continuous(low = "white", high = "red", name = "infection_rate",trans='sqrt') + theme(legend.position = "right")

```


* MOSTLY UNRELATED WORK HERE AND BELOW, THROWING THINGS AT THE WALL, MUCH DOES NOT STICK *

```{r}
df <- full_data %>% arrange(date) %>% group_by(county,state) %>% 
  mutate(jail_pop_diff = jail_population-lag(jail_population,7)) %>%
  mutate(ir = cases/resident_population) %>%
  mutate(ir_diff = ir-lag(ir,7)) %>%
  arrange(county,state,date) %>% drop_na() %>%
  mutate(density = ifelse(urbanicity=='urban',3,
                           ifelse(urbanicity=='small/mid',2,1))) %>%
  select(date,county,state,jail_pop_diff,resident_population,ir_diff,density)
  #filter(date >= as.Date('2020-04-01')) %>%
df

model <- glm(ir_diff ~ jail_pop_diff + resident_population + density, data = df)
summary(model)
```


(Not 100% sure about this analysis--please help!)
```{r}
# jail population dataframe with fraction released

frac_rel <- full_data %>% group_by(county,state) %>% filter(first(jail_population) >20)%>% 
  mutate(num_released = first(jail_population)-last(jail_population)) %>% 
  mutate(frac_released = num_released/first(jail_population)) %>%
  mutate(infection_rate=cases/resident_population) %>%
  select(county,state,frac_released,num_released,resident_population,infection_rate) %>%
  unique() %>% arrange(desc(frac_released))
  
frac_rel

# 5 quantiles
(quantile(frac_rel$frac_released, c(.2,.4,.6,.8,1)))


frac_rel$quantile <- ifelse(frac_rel$frac_released<.08,5,
                               ifelse(frac_rel$frac_released<.16,4,
                                      ifelse(frac_rel$frac_released<.24,3,
                                             ifelse(frac_rel$frac_released<.34,2,1))))

frac_rel$quantile <- as.factor(frac_rel$quantile)

release_df <- merge(full_data_recent,frac_rel) %>%
  mutate(infection_rate=cases/resident_population) %>% select(county,state,resident_population,infection_rate,frac_released,quantile) %>% unique()
release_df$quantile <- as.factor(release_df$quantile) 

table(release_df$quantile)

release_df
ggplot(release_df) + geom_violin(aes(quantile,log(infection_rate)))
```



```{r}
ggplot(release_df) + geom_point(aes(log(frac_released),log(infection_rate)))
```