# 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.35256474
## 2 Cat I 0.54864514
## 3 Cat II 0.05119488
## 4 Cat III 0.04759524
##
## [[1]]$tick
## [1] 1
##
##
## [[2]]
## [[2]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3525
## 2 Cat I 0.5486
## 3 Cat II 0.0513
## 4 Cat III 0.0476
##
## [[2]]$tick
## [1] 11
##
##
## [[3]]
## [[3]]$proportions
## Category Proportion
## 1 Non-Drinking 0.35252949
## 2 Cat I 0.54839032
## 3 Cat II 0.05128974
## 4 Cat III 0.04779044
##
## [[3]]$tick
## [1] 21
##
##
## [[4]]
## [[4]]$proportions
## Category Proportion
## 1 Non-Drinking 0.35246475
## 2 Cat I 0.54854515
## 3 Cat II 0.05129487
## 4 Cat III 0.04769523
##
## [[4]]$tick
## [1] 31
##
##
## [[5]]
## [[5]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3525
## 2 Cat I 0.5484
## 3 Cat II 0.0513
## 4 Cat III 0.0478
##
## [[5]]$tick
## [1] 41
##
##
## [[6]]
## [[6]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3525
## 2 Cat I 0.5484
## 3 Cat II 0.0513
## 4 Cat III 0.0478
##
## [[6]]$tick
## [1] 51
head(alcohol_proportions_df)
## Category Proportion Tick
## 1 Non-Drinking 0.35256474 1
## 2 Cat I 0.54864514 1
## 3 Cat II 0.05119488 1
## 4 Cat III 0.04759524 1
## 5 Non-Drinking 0.35250000 11
## 6 Cat I 0.54860000 11
tail(alcohol_summaries)
## [[1]]
## [[1]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3424
## 2 Cat I 0.5655
## 3 Cat II 0.0494
## 4 Cat III 0.0427
##
## [[1]]$tick
## [1] 10891
##
##
## [[2]]
## [[2]]$proportions
## Category Proportion
## 1 Non-Drinking 0.34266573
## 2 Cat I 0.56554345
## 3 Cat II 0.04929507
## 4 Cat III 0.04249575
##
## [[2]]$tick
## [1] 10901
##
##
## [[3]]
## [[3]]$proportions
## Category Proportion
## 1 Non-Drinking 0.34293141
## 2 Cat I 0.56548690
## 3 Cat II 0.04929014
## 4 Cat III 0.04229154
##
## [[3]]$tick
## [1] 10911
##
##
## [[4]]
## [[4]]$proportions
## Category Proportion
## 1 Non-Drinking 0.34289713
## 2 Cat I 0.56553034
## 3 Cat II 0.04918524
## 4 Cat III 0.04238728
##
## [[4]]$tick
## [1] 10921
##
##
## [[5]]
## [[5]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3428
## 2 Cat I 0.5656
## 3 Cat II 0.0492
## 4 Cat III 0.0424
##
## [[5]]$tick
## [1] 10931
##
##
## [[6]]
## [[6]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3430
## 2 Cat I 0.5656
## 3 Cat II 0.0491
## 4 Cat III 0.0423
##
## [[6]]$tick
## [1] 10941
tail(alcohol_proportions_df)
## Category Proportion Tick
## 4375 Cat II 0.0492 10931
## 4376 Cat III 0.0424 10931
## 4377 Non-Drinking 0.3430 10941
## 4378 Cat I 0.5656 10941
## 4379 Cat II 0.0491 10941
## 4380 Cat III 0.0423 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.34336544
## 2 Cat I 0.53384163
## 3 Cat II 0.01596219
## 4 Cat III 0.10683075
# 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.5
## 2 Cat I 0.4
## 3 Cat III 0.1
##
## $`7 days after release`
## Category Proportion
## 1 Non-Drinking 0.42
## 2 Cat I 0.46
## 3 Cat III 0.12
##
## $`14 days after release`
## Category Proportion
## 1 Non-Drinking 0.4047619
## 2 Cat I 0.4761905
## 3 Cat III 0.1190476
##
## $`30 days after release`
## Category Proportion
## 1 Non-Drinking 0.4090909
## 2 Cat I 0.4886364
## 3 Cat III 0.1022727
##
## $`90 days after release`
## Category Proportion
## 1 Non-Drinking 0.4146868
## 2 Cat I 0.4794816
## 3 Cat III 0.1058315
##
## $`180 days after release`
## Category Proportion
## 1 Non-Drinking 0.42216687
## 2 Cat I 0.48069738
## 3 Cat II 0.00249066
## 4 Cat III 0.09464508
##
## $`365 days after release`
## Category Proportion
## 1 Non-Drinking 0.38068966
## 2 Cat I 0.52827586
## 3 Cat II 0.01379310
## 4 Cat III 0.07724138
# 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.10000000 0.12000000 0.11904762
## 30 days after release 90 days after release 180 days after release
## 0.10227273 0.10583153 0.09464508
## 365 days after release
## 0.07724138
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.04759524
## 2 11 0.04760000
## 3 21 0.04779044
## 4 31 0.04769523
## 5 41 0.04780000
## 6 51 0.04780000
# 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
