GLM_building.Rmd

```{r}
library(tidyverse)
setwd("~/Documents/Python/Datathon/MIT_COVID-19_Datathon")
```

nyt_data: county-level COVID statistics
```{r}
nyt_data = read_csv('~/Documents/Python/Datathon/covid-19-data/us-counties.csv') 
# https://github.com/nytimes/covid-19-data
```


Prison analysis
prison_df: state, prison name, and county
```{r}
prisons <- read_csv('prison_boundaries.csv')
prisons$VAL_DATE <- as.Date(prisons$VAL_DATE)
str_list <- paste(c('CORRECTION','PRISON'), collapse="|")
prison_df <- prisons %>% filter(grepl(str_list,NAME)) %>%
  select(NAME,STATE,COUNTY,POPULATION,CAPACITY,COUNTYFIPS) %>% 
  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) %>%
  group_by(COUNTYFIPS) %>% 
  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,COUNTYFIPS) %>% unique()
  
prison_df
```

Get the states and their abbreviations, add State column to prison_df
```{r}
statenames <- read_csv('state_abbreviations.csv') %>% select(State,Code) %>% rename(STATE=Code)
prison_df <- merge(prison_df,statenames,all.x=TRUE)
prison_df
```

pr: state, fraction of prison population released
prison_counties_frac_released: counties with prisons, fraction of prisoners released statewide
frac_released = releases/prior_pop
```{r}
pr <- read_csv('Data/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') %>%              # not considering national-level releases
  mutate(releases = as.numeric(releases)) %>%   
  mutate(frac_released = releases/prior_pop) %>% # fraction released = # of releases / prior prison population
  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)) %>% # if a state has no reported releases, put 0 for frac_released
  select(COUNTYFIPS,frac_released,prison_pop,prison_capacity,prison_fullness,State) %>%unique()
prison_counties_frac_released
```


data for GLMs

```{r}
county_pop <- read_csv('new_county_demographics_df.csv') %>% rename(COUNTYFIPS = fips)
county_pop

nyt_data %>% arrange(fips,date)

nytdf <- nyt_data %>%
  mutate(prison_county = ifelse(fips %in% prison_counties_frac_released$COUNTYFIPS,1,0)) %>% # add column of 0/1: 1 if county has prison, 0 else
  rename(COUNTYFIPS = fips) %>%
  unique()
nytdf

nyt_with_prisons <- merge(nytdf, prison_counties_frac_released, all.x=TRUE) %>% group_by(COUNTYFIPS) %>%
  arrange(COUNTYFIPS,date) %>%
  #mutate(diff_cases = cases-lag(cases,7)) %>%
  mutate(daily_cases = cases-lag(cases,1)) %>%
  #mutate(frac_change_cases = (cases-lag(cases,7))/lag(cases,7)) %>%
  filter(date >= as.Date('2020-03-01')) %>%
  select(COUNTYFIPS,county,state,date,cases,daily_cases,deaths,prison_county,frac_released,prison_pop,prison_capacity,prison_fullness) %>%
  unique() %>% drop_na(COUNTYFIPS)

prison_pop_df <- merge(nyt_with_prisons,county_pop,all.x=TRUE) %>% 
  group_by(COUNTYFIPS) %>%
  mutate(cases_per_capita = cases/tot_pop) %>%
  mutate(deaths_per_capita = deaths/tot_pop) %>%
  arrange(COUNTYFIPS,date) %>%
  #mutate(daily_cases_7_days_ago = lag(daily_cases,7)) %>%
  #mutate(cases_per_capita_7_days_ago = lag(cases_per_capita,7)) %>% 
  filter(!is.na(daily_cases)) %>%
  mutate(prison_pop = ifelse(is.na(prison_pop),0,prison_pop)) %>%
  filter(daily_cases >= 0)

prison_pop_df

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


#write.csv(prison_pop_df,'Data/full_population_and_prison_data.csv')
```


```{r}
single_county = prison_pop_df %>% filter(COUNTYFIPS == '01001')
hist(single_county$daily_cases,20)

sny <- prison_pop_df %>% filter(COUNTYFIPS=='36103')
hist(sny$daily_cases,30)
```

Need to use a Poisson model: Target variable is the daily incidence of cases.

Incorporate history dependence: AIC analysis (probably) to figure out how many days in history to include.

Goodness-of-fit analysis: see how the residuals are distributed over time, by prison county indicatior, location, etc. See if evenly dispersed

This shows that as we increase the number of history terms, the AIC decreases. However, it's not completely fair as we're losing one data point per county each time we add a history term.
```{r}
aic_vec = c()
J = 50
for(j in 1:J){
  print(j)
  model_df <- prison_pop_df %>% select(COUNTYFIPS, date, daily_cases, tot_pop, prison_county, population_density, lowrisk_agegroup_perc)%>%
    group_by(COUNTYFIPS) %>% arrange(COUNTYFIPS, date)
  
  for(i in 1:j){
    varname <- paste("cases",i,sep="")
    model_df <- model_df %>% mutate(!!varname := lag(daily_cases,i))
  }
  
  model_df <- model_df %>% drop_na() %>% ungroup() %>% select(-COUNTYFIPS,-date)
  poiss_model = glm(daily_cases ~ ., family=poisson, data = model_df)
  
  aic_vec <- c(aic_vec,poiss_model$aic)
  
}

plot(1:J,aic_vec,type='l')
```

Set the number of history terms (here j = 7)
```{r}
j = 7
model_df <- prison_pop_df %>% 
  select(COUNTYFIPS, date, daily_cases, tot_pop, prison_county, population_density, lowrisk_agegroup_perc) %>%
  group_by(COUNTYFIPS) %>% arrange(COUNTYFIPS, date)
  
for(i in 1:j){
  varname <- paste("cases",i,sep="")
  model_df <- model_df %>% mutate(!!varname := lag(daily_cases,i))
}

model_df <- model_df %>% drop_na() %>% ungroup() %>% select(-COUNTYFIPS,-date)
poiss_model = glm(daily_cases ~ ., family=poisson, data = model_df)
summary(poiss_model)
```

It's possible that if the above models work we can interpret it this way:
a 0.84 coefficient for prison_county means we multiply the daily incidence of cases by 1 when no prison in the county or by e^0.84~=2 when there is a prison county --> daily incidence of cases doubles!

Now consider just counties with prisons, examine the effect of statewide fraction released:
interpret -1.37 coefficient on frac_released as "if we reduce the prison population by 25%, we can expect a 30% reduction in the surrounding county's COVID-19 cases"

```{r}
j = 7
model_df <- prison_pop_df %>%
  filter(prison_county==1) %>%
  select(COUNTYFIPS, date, daily_cases, tot_pop, frac_released, population_density, lowrisk_agegroup_perc) %>%
  group_by(COUNTYFIPS) %>% arrange(COUNTYFIPS, date)
  
for(i in 1:j){
  varname <- paste("cases",i,sep="")
  model_df <- model_df %>% mutate(!!varname := lag(daily_cases,i))
}

model_df <- model_df %>% drop_na() %>% ungroup() %>% select(-COUNTYFIPS,-date)
poiss_model = glm(daily_cases ~ ., family=poisson, data = model_df)
summary(poiss_model)

```



##########################################################################################################
OLD MODELS ASSUMING NORMAL DISTRIBUTION

Models assessing the impact of having a prison in a given county
Model 1:
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 2:
difference in cases (cases from today - cases from 7 days ago) ~ total prison population (0 if no prison in county), total county population, cases per capita from 7 days ago, population density

```{r}
model1 <- lm(diff_cases ~ tot_pop + prison_county + population_density, data = prison_pop_df)
summary(model1)

model2 <- lm(diff_cases ~  prison_pop + tot_pop + cases_per_capita_7_days_ago + population_density, data = prison_pop_df)
summary(model2)

model2 <- lm(cases_per_capita ~  prison_county + tot_pop + cases_per_capita_7_days_ago + population_density, data = prison_pop_df)
summary(model2)

model2 <- lm(cases_per_capita ~  prison_county + tot_pop + population_density, data = prison_pop_df)
summary(model2)
```

Models assessing the impact of release rates in counties with prisons


Model A:
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 B:
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, fraction of population in high risk age group

```{r}
modelA <- lm(diff_cases ~ tot_pop + frac_released + population_density, data = prison_counties)
summary(modelA)

modelA <- lm(cases_per_capita ~ tot_pop + frac_released + population_density + cases_per_capita_7_days_ago, data = prison_counties)
summary(modelA)

modelB <- lm(diff_cases ~ tot_pop + frac_released +  population_density + lowrisk_agegroup_perc,data = prison_counties)
summary(modelB)

```

Repeat these analyses using frac_change_cases as the target variable instead of diff_cases

```{r}
model1 <- lm(frac_change_cases ~ tot_pop + prison_county + population_density, data = prison_pop_df)
summary(model1)

model2 <- lm(frac_change_cases ~  prison_pop + tot_pop + cases_per_capita_7_days_ago + population_density, data = prison_pop_df)
summary(model2)
```

```{r}
modelA <- lm(frac_change_cases ~ tot_pop + frac_released + population_density, data = prison_counties)
summary(modelA)

modelB <- glm(frac_change_cases ~ tot_pop + frac_released +  population_density + lowrisk_agegroup_perc,data = prison_counties)
summary(modelB)

```


```{r}
modelx <- lm(cases_per_capita ~ tot_pop + prison_county + population_density, data = prison_pop_df)
summary(modelx)

# this one has an R^2 value of 0.9!

modely <- lm(cases_per_capita ~  prison_pop + tot_pop + cases_per_capita_7_days_ago + population_density, data = prison_pop_df)
summary(modely)

modelz <- lm(deaths_per_capita ~ tot_pop + prison_county + population_density, data = prison_pop_df)
summary(modelz)

modelw <- lm(deaths_per_capita ~  prison_pop + tot_pop + cases_per_capita_7_days_ago + population_density, data = prison_pop_df)
summary(modelw)
```

Working with ACLU Massachusetts dataset: ACLU.csv gives for each considered DOC facility on each date the total population, detainees tested/positive, staff tested/positive, # released, and fips county code.

ACLU_cases dataframe gives date, county fips code, county-wide prisoner cases and cumulative sum of county-wide prisoner cases

mass_prison_counties adds information from prison_pop_df to ACLU_cases: we now have county-level cases, deaths, and population demographics

```{r}
# excluding hospitals, rehabilitation centers
DOC_Facility = c('MCI-Norfolk','NCCI-Gardn','OCCC','MCI-CJ','MCI-C','Pondville','MCI-F',
                       'MCI-Shirley','MCI-Shirley Min','MTC','SMCC','SBCC')
fips = c(25021,25027,25023,25021,25017,25021,25017,25017,25017,25023,25017,25017)
facilities_fips <- data.frame(DOC_Facility,fips)
facilities_fips

ACLU <- read_csv('ACLU_data.csv')
names(ACLU) <- str_replace_all(names(ACLU)," ","_")
ACLU <- merge(ACLU,facilities_fips,all.y=TRUE) %>% mutate(Date = as.Date(Date,"%m/%d/%y"))
ACLU <- ACLU[,c(1:11,27)] # remove empty columns
ACLU
write.csv(ACLU,'Data/ACLU_prison_data.csv')

ACLU_cases <- ACLU  %>% 
  select(Date,DOC_Facility,fips,`N_Positive_-_Detainees/Inmates`) %>%
  mutate(`N_Positive_-_Detainees/Inmates` = replace_na(`N_Positive_-_Detainees/Inmates`, 0))%>% # replace NA with 0 to get total counts
  group_by(fips,Date) %>%
  mutate(county_wide_prisoner_cases = sum(`N_Positive_-_Detainees/Inmates`)) %>% # sum over all cases in a county
  arrange(fips,Date) %>% 
  select(Date,fips,county_wide_prisoner_cases) %>%
  rename(COUNTYFIPS = fips,date = Date) %>% unique() %>% ungroup %>% group_by(COUNTYFIPS) %>%
  mutate(cumulative_county_wide_prisoner_cases = cumsum(county_wide_prisoner_cases)) # cumuylative cases in a county
ACLU_cases

mass_prison_counties <- merge(prison_pop_df,ACLU_cases,all.y=TRUE)
mass_prison_counties
```

```{r}
modelx <- lm(cases ~ cumulative_county_wide_prisoner_cases+ tot_pop, data = mass_prison_counties)
summary(modelx)
```


This chunk shows that across all prisons, we only have data on 10 releases --> not powerful enough for our model. We'll have to either A) consider case data that we get from each prison or B) consider jail data.

A) is likely not extremely useful given large gaps in testing
```{r}
tmp <- ACLU %>% select(fips,Date,N_Released) %>% mutate(N_Released = replace_na(N_Released,0))
sum(tmp$N_Released)
tmp <- read_csv('ACLU_data.csv')
names(tmp) <- str_replace_all(names(tmp)," ","_")
tmp <- tmp %>% select(Date,DOC_Facility,N_Released) %>% drop_na()
sum(tmp$N_Released)
```

Now considering ACLU data for county jails: MA_aclu_data has date, county info, population demographics, county-level cases and deaths, release statistics

```{r}
aclu_data <- read_csv('ACLU_data_Counties_All_DOC.csv') %>% rename(county = County,date = Date) %>%
  mutate(date = as.Date(date,"%m/%d/%y"))
aclu_data <- aclu_data[,1:19]
names(aclu_data) <- str_replace_all(names(aclu_data)," ","_")

aclu_data <- aclu_data %>% filter(county != 'DOC') %>% # just look at county jails, not DOC
  select(date,county,`N_Released_Pre-Trial`,N_Released_Sentenced,Total_Population)
aclu_data[is.na(aclu_data)] = 0


MA_prison_pop_df <- prison_pop_df %>% filter(state=='Massachusetts') %>%
  filter(county %in% aclu_data$county)

MA_aclu_data <- merge(MA_prison_pop_df,aclu_data)%>% group_by(county) %>% arrange(county,date)%>%
  mutate(cumulative_pretrial_release = cumsum(`N_Released_Pre-Trial`))%>%
  mutate(cumulative_sentenced_release = cumsum(N_Released_Sentenced)) %>%
  mutate(cumulative_total_release = cumulative_sentenced_release + cumulative_pretrial_release) %>%
  mutate(total_release = `N_Released_Pre-Trial` + N_Released_Sentenced) %>%
  select(date,county,COUNTYFIPS,tot_pop,cases,cases_per_capita,diff_cases,population_density,deaths,cumulative_sentenced_release,cumulative_pretrial_release,total_release,Total_Population,frac_change_cases,cases_per_capita_7_days_ago) %>%
  mutate(frac_released = total_release/Total_Population) %>% filter(Total_Population>0)

prison_pop_df
MA_aclu_data
```

Tough to infer causality --> do we increase frac_released as we increase cases_per_capita because we're experiencing an overall increase? How could we show what the cases_per_capita would be if we released no detainees? include days since first case, something like that?

```{r}
modelx <- lm(frac_change_cases ~ total_release + tot_pop, data = MA_aclu_data)
summary(modelx)

modely0 <- lm(cases_per_capita ~ cases_per_capita_7_days_ago, data = MA_aclu_data)
summary(modely0)
modely <- lm(cases_per_capita ~ frac_released + cases_per_capita_7_days_ago, data = MA_aclu_data)
summary(modely)

AIC(modely0)
AIC(modely)
```

Notes: I think the way I calculate R^2 is wrong