# 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-log-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-log-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-log-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-log-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.35192961
## 2 Cat I 0.55268946
## 3 Cat II 0.04729054
## 4 Cat III 0.04809038
##
## [[1]]$tick
## [1] 1
##
##
## [[2]]
## [[2]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3519
## 2 Cat I 0.5527
## 3 Cat II 0.0473
## 4 Cat III 0.0481
##
## [[2]]$tick
## [1] 11
##
##
## [[3]]
## [[3]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3520
## 2 Cat I 0.5526
## 3 Cat II 0.0473
## 4 Cat III 0.0481
##
## [[3]]$tick
## [1] 21
##
##
## [[4]]
## [[4]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3519
## 2 Cat I 0.5526
## 3 Cat II 0.0474
## 4 Cat III 0.0481
##
## [[4]]$tick
## [1] 31
##
##
## [[5]]
## [[5]]$proportions
## Category Proportion
## 1 Non-Drinking 0.35186481
## 2 Cat I 0.55264474
## 3 Cat II 0.04739526
## 4 Cat III 0.04809519
##
## [[5]]$tick
## [1] 41
##
##
## [[6]]
## [[6]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3520
## 2 Cat I 0.5525
## 3 Cat II 0.0474
## 4 Cat III 0.0481
##
## [[6]]$tick
## [1] 51
head(alcohol_proportions_df)
## Category Proportion Tick
## 1 Non-Drinking 0.35192961 1
## 2 Cat I 0.55268946 1
## 3 Cat II 0.04729054 1
## 4 Cat III 0.04809038 1
## 5 Non-Drinking 0.35190000 11
## 6 Cat I 0.55270000 11
tail(alcohol_summaries)
## [[1]]
## [[1]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3522
## 2 Cat I 0.5538
## 3 Cat II 0.0472
## 4 Cat III 0.0468
##
## [[1]]$tick
## [1] 10891
##
##
## [[2]]
## [[2]]$proportions
## Category Proportion
## 1 Non-Drinking 0.35226477
## 2 Cat I 0.55374463
## 3 Cat II 0.04719528
## 4 Cat III 0.04679532
##
## [[2]]$tick
## [1] 10901
##
##
## [[3]]
## [[3]]$proportions
## Category Proportion
## 1 Non-Drinking 0.3522
## 2 Cat I 0.5537
## 3 Cat II 0.0473
## 4 Cat III 0.0468
##
## [[3]]$tick
## [1] 10911
##
##
## [[4]]
## [[4]]$proportions
## Category Proportion
## 1 Non-Drinking 0.35226477
## 2 Cat I 0.55354465
## 3 Cat II 0.04729527
## 4 Cat III 0.04689531
##
## [[4]]$tick
## [1] 10921
##
##
## [[5]]
## [[5]]$proportions
## Category Proportion
## 1 Non-Drinking 0.35249425
## 2 Cat I 0.55303409
## 3 Cat II 0.04748575
## 4 Cat III 0.04698590
##
## [[5]]$tick
## [1] 10931
##
##
## [[6]]
## [[6]]$proportions
## Category Proportion
## 1 Non-Drinking 0.35252949
## 2 Cat I 0.55288942
## 3 Cat II 0.04749050
## 4 Cat III 0.04709058
##
## [[6]]$tick
## [1] 10941
tail(alcohol_proportions_df)
## Category Proportion Tick
## 4375 Cat II 0.04748575 10931
## 4376 Cat III 0.04698590 10931
## 4377 Non-Drinking 0.35252949 10941
## 4378 Cat I 0.55288942 10941
## 4379 Cat II 0.04749050 10941
## 4380 Cat III 0.04709058 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.31212373
## 2 Cat I 0.36251285
## 3 Cat II 0.04145661
## 4 Cat III 0.28390681
# 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.30612245
## 2 Cat I 0.37414966
## 3 Cat II 0.02040816
## 4 Cat III 0.29931973
##
## $`7 days after release`
## Category Proportion
## 1 Non-Drinking 0.28585559
## 2 Cat I 0.33432245
## 3 Cat II 0.03956479
## 4 Cat III 0.34025717
##
## $`14 days after release`
## Category Proportion
## 1 Non-Drinking 0.29025745
## 2 Cat I 0.34174659
## 3 Cat II 0.03785967
## 4 Cat III 0.33013629
##
## $`30 days after release`
## Category Proportion
## 1 Non-Drinking 0.29195348
## 2 Cat I 0.33728934
## 3 Cat II 0.03868977
## 4 Cat III 0.33206741
##
## $`90 days after release`
## Category Proportion
## 1 Non-Drinking 0.29886569
## 2 Cat I 0.34226533
## 3 Cat II 0.03912543
## 4 Cat III 0.31974355
##
## $`180 days after release`
## Category Proportion
## 1 Non-Drinking 0.30383184
## 2 Cat I 0.34898183
## 3 Cat II 0.03998248
## 4 Cat III 0.30720385
##
## $`365 days after release`
## Category Proportion
## 1 Non-Drinking 0.31212373
## 2 Cat I 0.36251285
## 3 Cat II 0.04145661
## 4 Cat III 0.28390681
# 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.2993197 0.3402572 0.3301363
## 30 days after release 90 days after release 180 days after release
## 0.3320674 0.3197435 0.3072039
## 365 days after release
## 0.2839068
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.04809038
## 2 11 0.04810000
## 3 21 0.04810000
## 4 31 0.04810000
## 5 41 0.04809519
## 6 51 0.04810000
# 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 = 12, face = "bold"))
text = element_text(size = 12))
p
