# 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.34606539
## 2 Cat I 0.55444456
## 3 Cat II 0.04919508
## 4 Cat III 0.05029497
##
## [[1]]$tick
## [1] 1
##
##
## [[2]]
## [[2]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3461
## 2 Cat I 0.5545
## 3 Cat II 0.0491
## 4 Cat III 0.0503
##
## [[2]]$tick
## [1] 11
##
##
## [[3]]
## [[3]]$proportions
## Category Proportion
## 1 Non-Drinking 0.34596540
## 2 Cat I 0.55474453
## 3 Cat II 0.04899510
## 4 Cat III 0.05029497
##
## [[3]]$tick
## [1] 21
##
##
## [[4]]
## [[4]]$proportions
## Category Proportion
## 1 Non-Drinking 0.34589623
## 2 Cat I 0.55483355
## 3 Cat II 0.04898530
## 4 Cat III 0.05028491
##
## [[4]]$tick
## [1] 31
##
##
## [[5]]
## [[5]]$proportions
## Category Proportion
## 1 Non-Drinking 0.34603079
## 2 Cat I 0.55478904
## 3 Cat II 0.04889022
## 4 Cat III 0.05028994
##
## [[5]]$tick
## [1] 41
##
##
## [[6]]
## [[6]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3458
## 2 Cat I 0.5549
## 3 Cat II 0.0490
## 4 Cat III 0.0503
##
## [[6]]$tick
## [1] 51
head(alcohol_proportions_df)
## Category Proportion Tick
## 1 Non-Drinking 0.34606539 1
## 2 Cat I 0.55444456 1
## 3 Cat II 0.04919508 1
## 4 Cat III 0.05029497 1
## 5 Non-Drinking 0.34610000 11
## 6 Cat I 0.55450000 11
tail(alcohol_summaries)
## [[1]]
## [[1]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3457
## 2 Cat I 0.5477
## 3 Cat II 0.0552
## 4 Cat III 0.0514
##
## [[1]]$tick
## [1] 7541
##
##
## [[2]]
## [[2]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3458
## 2 Cat I 0.5477
## 3 Cat II 0.0552
## 4 Cat III 0.0513
##
## [[2]]$tick
## [1] 7551
##
##
## [[3]]
## [[3]]$proportions
## Category Proportion
## 1 Non-Drinking 0.34559632
## 2 Cat I 0.54793562
## 3 Cat II 0.05518344
## 4 Cat III 0.05128461
##
## [[3]]$tick
## [1] 7561
##
##
## [[4]]
## [[4]]$proportions
## Category Proportion
## 1 Non-Drinking 0.34536546
## 2 Cat I 0.54794521
## 3 Cat II 0.05529447
## 4 Cat III 0.05139486
##
## [[4]]$tick
## [1] 7571
##
##
## [[5]]
## [[5]]$proportions
## Category Proportion
## 1 Non-Drinking 0.34516548
## 2 Cat I 0.54814519
## 3 Cat II 0.05529447
## 4 Cat III 0.05139486
##
## [[5]]$tick
## [1] 7581
##
##
## [[6]]
## [[6]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3450
## 2 Cat I 0.5483
## 3 Cat II 0.0553
## 4 Cat III 0.0514
##
## [[6]]$tick
## [1] 7591
tail(alcohol_proportions_df)
## Category Proportion Tick
## 3035 Cat II 0.05529447 7581
## 3036 Cat III 0.05139486 7581
## 3037 Non-Drinking 0.34500000 7591
## 3038 Cat I 0.54830000 7591
## 3039 Cat II 0.05530000 7591
## 3040 Cat III 0.05140000 7591
# 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.26449348
## 2 Cat I 0.40819785
## 3 Cat II 0.05135757
## 4 Cat III 0.27595110
# 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.30081301
## 2 Cat I 0.34146341
## 3 Cat II 0.04065041
## 4 Cat III 0.31707317
##
## $`7 days after release`
## Category Proportion
## 1 Non-Drinking 0.24058920
## 2 Cat I 0.39934534
## 3 Cat II 0.04746318
## 4 Cat III 0.31260229
##
## $`14 days after release`
## Category Proportion
## 1 Non-Drinking 0.25767918
## 2 Cat I 0.37883959
## 3 Cat II 0.05119454
## 4 Cat III 0.31228669
##
## $`30 days after release`
## Category Proportion
## 1 Non-Drinking 0.25857200
## 2 Cat I 0.38725292
## 3 Cat II 0.05082695
## 4 Cat III 0.30334812
##
## $`90 days after release`
## Category Proportion
## 1 Non-Drinking 0.25720683
## 2 Cat I 0.39378673
## 3 Cat II 0.05191716
## 4 Cat III 0.29708928
##
## $`180 days after release`
## Category Proportion
## 1 Non-Drinking 0.2612413
## 2 Cat I 0.4005602
## 3 Cat II 0.0516733
## 4 Cat III 0.2865251
##
## $`365 days after release`
## Category Proportion
## 1 Non-Drinking 0.26449348
## 2 Cat I 0.40819785
## 3 Cat II 0.05135757
## 4 Cat III 0.27595110
# 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.3170732 0.3126023 0.3122867
## 30 days after release 90 days after release 180 days after release
## 0.3033481 0.2970893 0.2865251
## 365 days after release
## 0.2759511
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.05029497
## 2 11 0.05030000
## 3 21 0.05029497
## 4 31 0.05028491
## 5 41 0.05028994
## 6 51 0.05030000
# 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
