# Analyze incarceration data
rm(list=ls())
renv::activate()
# Load packages ---------
library(here)
## here() starts at /users/akhann16/code/cadre/data-analysis-plotting/Simulated-Data-Analysis/r/agent-analysis
library(data.table)
library(yaml)
library(ggplot2)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:data.table':
##
## between, first, last
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
# Utility functions ------------
#source(here("agent-analysis", "utils", "post-release-alcohol.R"))
#print(here("agent-analysis", "utils", "post-release-alcohol.R"))
source("/users/akhann16/code/cadre/data-analysis-plotting/Simulated-Data-Analysis/r/agent-analysis/utils/post-release-alcohol.R")
# Read RDS file ------------
#agent_log_env <- readRDS(here("agent-analysis", "rds-outs", "agent_log_env.RDS"))
print(here("agent-analysis", "rds-outs", "agent_log_env.RDS"))
## [1] "/users/akhann16/code/cadre/data-analysis-plotting/Simulated-Data-Analysis/r/agent-analysis/agent-analysis/rds-outs/agent_log_env.RDS"
agent_log_env <- readRDS("/users/akhann16/code/cadre/data-analysis-plotting/Simulated-Data-Analysis/r/agent-analysis/rds-outs/agent_log_env.RDS")
# Load data ------------
agent_dt <- agent_log_env[["agent_dt"]]
input_params <- agent_log_env[["input_params"]]
# Ticks of interest ----------
selected_ticks <- seq(1, max(agent_dt$tick), by = 10)
last_tick <- max(agent_dt$tick)
# Population subsets by release status --------
desired_tick=last_tick
all_agents <- agent_dt[tick==desired_tick,] #everyone at all times
never_released_agents <- agent_dt[tick==desired_tick & n_releases == 0] #never release
released_agents <- agent_dt[tick == desired_tick & n_releases >= 1] #released at least once
# Alcohol use state distributions in above subsets --------
all_agents_summary <- calculate_proportions(all_agents)
never_released_agents_summary <- calculate_proportions(never_released_agents)
released_agents_summary <- calculate_proportions(released_agents)
# Summarize alcohol use rates in general population --------
alcohol_summaries <-
lapply(selected_ticks, function(x) {
agents_at_tick <- agent_dt[tick == x]
proportions <- calculate_proportions(agents_at_tick)
return(list("proportions" = proportions, "tick" = x))
})
alcohol_proportions_df <-
do.call(rbind, lapply(alcohol_summaries, function(x) {
df <- x$proportions
df$Tick <- x$tick
return(df)
}
)
)
head(alcohol_summaries)
## [[1]]
## [[1]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3437
## 2 Cat I 0.5565
## 3 Cat II 0.0511
## 4 Cat III 0.0487
##
## [[1]]$tick
## [1] 1
##
##
## [[2]]
## [[2]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3405
## 2 Cat I 0.5573
## 3 Cat II 0.0524
## 4 Cat III 0.0498
##
## [[2]]$tick
## [1] 11
##
##
## [[3]]
## [[3]]$proportions
## Category Proportion
## 1 Non-Drinking 0.33803239
## 2 Cat I 0.55628874
## 3 Cat II 0.05438912
## 4 Cat III 0.05128974
##
## [[3]]$tick
## [1] 21
##
##
## [[4]]
## [[4]]$proportions
## Category Proportion
## 1 Non-Drinking 0.33516648
## 2 Cat I 0.55664434
## 3 Cat II 0.05619438
## 4 Cat III 0.05199480
##
## [[4]]$tick
## [1] 31
##
##
## [[5]]
## [[5]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3316
## 2 Cat I 0.5569
## 3 Cat II 0.0588
## 4 Cat III 0.0527
##
## [[5]]$tick
## [1] 41
##
##
## [[6]]
## [[6]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3285
## 2 Cat I 0.5576
## 3 Cat II 0.0608
## 4 Cat III 0.0531
##
## [[6]]$tick
## [1] 51
head(alcohol_proportions_df)
## Category Proportion Tick
## 1 Non-Drinking 0.3437 1
## 2 Cat I 0.5565 1
## 3 Cat II 0.0511 1
## 4 Cat III 0.0487 1
## 5 Non-Drinking 0.3405 11
## 6 Cat I 0.5573 11
tail(alcohol_summaries)
## [[1]]
## [[1]]$proportions
## Category Proportion
## 1 Non-Drinking 0.11388861
## 2 Cat I 0.66513349
## 3 Cat II 0.09549045
## 4 Cat III 0.12548745
##
## [[1]]$tick
## [1] 10891
##
##
## [[2]]
## [[2]]$proportions
## Category Proportion
## 1 Non-Drinking 0.1136
## 2 Cat I 0.6654
## 3 Cat II 0.0951
## 4 Cat III 0.1259
##
## [[2]]$tick
## [1] 10901
##
##
## [[3]]
## [[3]]$proportions
## Category Proportion
## 1 Non-Drinking 0.1135
## 2 Cat I 0.6666
## 3 Cat II 0.0953
## 4 Cat III 0.1246
##
## [[3]]$tick
## [1] 10911
##
##
## [[4]]
## [[4]]$proportions
## Category Proportion
## 1 Non-Drinking 0.1130
## 2 Cat I 0.6670
## 3 Cat II 0.0956
## 4 Cat III 0.1244
##
## [[4]]$tick
## [1] 10921
##
##
## [[5]]
## [[5]]$proportions
## Category Proportion
## 1 Non-Drinking 0.11278872
## 2 Cat I 0.66773323
## 3 Cat II 0.09529047
## 4 Cat III 0.12418758
##
## [[5]]$tick
## [1] 10931
##
##
## [[6]]
## [[6]]$proportions
## Category Proportion
## 1 Non-Drinking 0.1127
## 2 Cat I 0.6697
## 3 Cat II 0.0945
## 4 Cat III 0.1231
##
## [[6]]$tick
## [1] 10941
tail(alcohol_proportions_df)
## Category Proportion Tick
## 4375 Cat II 0.09529047 10931
## 4376 Cat III 0.12418758 10931
## 4377 Non-Drinking 0.11270000 10941
## 4378 Cat I 0.66970000 10941
## 4379 Cat II 0.09450000 10941
## 4380 Cat III 0.12310000 10941
# Compare alcohol use rates in short time frame after release (1 year) --------
fixed_time_after_release <- 365
# filter agents who were released at least once and within one year after release
within_fixed_time_after_release_agents <- agent_dt[last_release_tick > 0 & tick <= last_release_tick + fixed_time_after_release]
# calculate proportions for agents within one year after release
within_fixed_time_after_release_summary <- calculate_proportions(within_fixed_time_after_release_agents)
# display the summary
print("Summary for Agents Within One Year After Release")
## [1] "Summary for Agents Within One Year After Release"
print(within_fixed_time_after_release_summary)
## Category Proportion
## 1 Non-Drinking 0.2672188
## 2 Cat I 0.5347901
## 3 Cat II 0.0605046
## 4 Cat III 0.1374864
# Sensitivity Analysis on parameter of time after release ---------
# time periods in ticks (days): 1 day, 1 week, 2 weeks, 1 month, 3 months, 6 months, 1 year
time_periods <- c(1, 7, 14, 30, 90, 180, 365)
# calculate summaries for each time period
summaries <- lapply(time_periods, function(x)
summary_after_release(x, agent_dt = agent_dt))
# display summaries
names(summaries) <- paste(time_periods, "days after release")
summaries
## $`1 days after release`
## Category Proportion
## 1 Non-Drinking 0.4
## 2 Cat I 0.6
##
## $`7 days after release`
## Category Proportion
## 1 Non-Drinking 0.33333333
## 2 Cat I 0.58333333
## 3 Cat III 0.08333333
##
## $`14 days after release`
## Category Proportion
## 1 Non-Drinking 0.2933333
## 2 Cat I 0.5866667
## 3 Cat III 0.1200000
##
## $`30 days after release`
## Category Proportion
## 1 Non-Drinking 0.3235294
## 2 Cat I 0.5470588
## 3 Cat III 0.1294118
##
## $`90 days after release`
## Category Proportion
## 1 Non-Drinking 0.34279919
## 2 Cat I 0.51724138
## 3 Cat II 0.01014199
## 4 Cat III 0.12981744
##
## $`180 days after release`
## Category Proportion
## 1 Non-Drinking 0.3430736
## 2 Cat I 0.4859307
## 3 Cat II 0.0465368
## 4 Cat III 0.1244589
##
## $`365 days after release`
## Category Proportion
## 1 Non-Drinking 0.33452169
## 2 Cat I 0.46880570
## 3 Cat II 0.06714201
## 4 Cat III 0.12953060
# Plot heavy drinking/AUD after release ---------
labels <- c("1D", "1W", "2W", "1M", "3M", "6M", "1Y")
heavy_use_or_AUD_proportions <- sapply(summaries, function(x) sum(x[x$Category == "Cat III", "Proportion"]))
heavy_use_or_AUD_proportions
## 1 days after release 7 days after release 14 days after release
## 0.00000000 0.08333333 0.12000000
## 30 days after release 90 days after release 180 days after release
## 0.12941176 0.12981744 0.12445887
## 365 days after release
## 0.12953060
heavy_use_AUD_df <- data.frame(Time = time_periods,
Labels = labels,
Proportion = heavy_use_or_AUD_proportions
)
ggplot(heavy_use_AUD_df, aes(x = Time, y = Proportion, group = 1)) +
geom_line() +
geom_point() +
scale_x_continuous(breaks = time_periods, labels = labels) +
scale_y_continuous(breaks = seq(0, 0.5, 0.1), limits = c(0, 0.5)) +
labs(title = "",
x = "Time After Release",
y = "Proportion of Agents with AUD") +
theme_minimal()+
theme(axis.text.x = element_text(angle = 90, hjust = 1))
# Compare Heavy Use/AUD in general population with the post-release group ---------
# extract rows for heavy use/AUD only
general_heavy_AUD_df <-
alcohol_proportions_df[alcohol_proportions_df$Category == 'Cat III', ]
# aggregate the heavy/AUD data
combined_proportions <- aggregate(Proportion ~ Tick, data = general_heavy_AUD_df, sum)
head(combined_proportions)
## Tick Proportion
## 1 1 0.04870000
## 2 11 0.04980000
## 3 21 0.05128974
## 4 31 0.05199480
## 5 41 0.05270000
## 6 51 0.05310000
# map the 'Tick' column to the same scale as the original 'heavy_use_AUD_df' dataframe
# compute average of AUD rates for the last 10 ticks
last_combined_rate <- mean(combined_proportions$Proportion[(nrow(combined_proportions) - 9):nrow(combined_proportions)])
# Create a new dataframe with these rates for the time points
last_combined_rate_df <- data.frame(Time = max(heavy_use_AUD_df$Time):365,
Proportion = rep(last_combined_rate, 365 - max(heavy_use_AUD_df$Time) + 1))
general_population <- data.frame(
Time = heavy_use_AUD_df$Time,
Proportion = last_combined_rate_df$Proportion)
# Plotting both datasets
p <-
ggplot() +
geom_line(data = heavy_use_AUD_df, aes(x = Time, y = Proportion, color = "Previously Incarcerated"), linewidth=1.5) +
geom_line(data = general_population, aes(x = Time, y = Proportion, color = "Overall"), linewidth=1.5)+
ylim(c(0,1))+
scale_color_manual(values = c("Previously Incarcerated" = "blue", "Overall" = "red")) +
theme_minimal() +
labs(title = "",
x = "Time After Release",
y = "Proportion",
color = "Population") +
scale_x_continuous(breaks = heavy_use_AUD_df$Time, labels = heavy_use_AUD_df$Labels) +
theme(axis.text.x = element_text(angle = 90, hjust = 1),
text = element_text(size = 24, face = "bold"))
p
