R Notebook
27 April, 2023
We confirm below various attributes of the agents that are logged.
We will do this in data.table for speed.
rm(list=ls())
renv::load()
library(data.table)
library(ggplot2)
library(yaml)Read input parameters:¬
input_params <- read_yaml("~/code/cadre/python/myparams/model_params.yaml")Read the agent data into data-table:
agent_dt <- fread("~/code/cadre/python/output/agent_log_9.csv")
str(agent_dt)## Classes 'data.table' and 'data.frame': 10965743 obs. of 16 variables:
## $ tick : int 1 1 1 1 1 1 1 1 1 1 ...
## $ id : int 0 1 2 3 4 5 6 7 8 9 ...
## $ age : int 72 41 72 56 26 28 72 37 54 61 ...
## $ race : chr "White" "Asian" "Hispanic" "Black" ...
## $ female : int 0 1 1 1 0 0 1 0 1 1 ...
## $ alc_use_status : int 1 0 1 1 1 1 1 2 1 0 ...
## $ smoking_status : chr "Former" "Never" "Never" "Never" ...
## $ last_incarceration_tick : int -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
## $ last_release_tick : int -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
## $ current_incarceration_status: int 0 0 0 0 0 0 0 0 0 0 ...
## $ entry_at_tick : int 0 0 0 0 0 0 0 0 0 0 ...
## $ exit_at_tick : int -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 ...
## $ n_incarcerations : int 0 0 0 0 0 0 0 0 0 0 ...
## $ n_releases : int 0 0 0 0 0 0 0 0 0 0 ...
## $ n_smkg_stat_trans : int 0 0 0 0 2 0 0 0 0 0 ...
## $ n_alc_use_stat_trans : int 0 0 0 0 0 0 0 0 0 0 ...
## - attr(*, ".internal.selfref")=<externalptr>last_tick <- max(agent_dt$tick)
print(last_tick)## [1] 10950How many unique agents are in the dataset?
length(unique(agent_dt$id))## [1] 15743range(unique(agent_dt$id))## [1] 0 15742There are 15743 unique agents in the dataset. Their IDs range from 0 to 15742.
Demographics
agent_dt_by_race <- agent_dt[tick == last_tick, .N, by=race][,"per_cent":=round(N/sum(N)*100)][order(race)]
target_values <- data.frame(race = names(input_params$RACE_DISTRIBUTION),
target_pct = unlist(input_params$RACE_DISTRIBUTION))# Calculate the percentages
agent_dt_by_race[, .(count = .N), by = race][, percentage := round(count/sum(count) * 100, 2)]
# Order the data by race
agent_dt_by_race <- agent_dt_by_race[order(race)]
# Create the bar chart
ggplot(agent_dt_by_race, aes(x = race, y = per_cent, fill = race)) +
geom_bar(stat = "identity", color = "black") +
scale_fill_manual(values = c("#7bccc4", "#43a2ca", "#0868ac", "#f0f9e8")) +
labs(title = "Race Distribution",
x = "Race",
y = "Percentage",
fill = "Race")+
geom_segment(data = target_values,
aes(x = race, xend = race, y = 0, yend = target_pct),
color = "red", size = 1, linetype = "dashed") +
geom_text(data = target_values,
aes(x = race, y = target_pct*97, label = paste0(round(target_pct*100, 1), "%")),
vjust = -0.5, size = 4, color = "red")## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
agent_dt[tick==last_tick, .N, by=female][,"%":=N/sum(N)][order(female)]female_target <- input_params$FEMALE_PROP
gender_pct <- agent_dt[tick == last_tick, .N, by = female][, "%":= N/sum(N)][order(female)]# calculate percentages
gender_pct <- agent_dt[tick == last_tick, .N, by = female][, "%":= N/sum(N)][order(female)]
# read in YAML file with target value
female_target <- input_params$FEMALE_PROP
# calculate target percentages for each gender
female_target_pct <- female_target
male_target_pct <- 1 - female_target
# calculate actual percentages for each gender
female_actual_pct <- gender_pct[female == 1, `%`]
male_actual_pct <- gender_pct[female == 0, `%`]
# create plot
ggplot(gender_pct, aes(x = ifelse(female == 0, "Male", "Female"), y = `%`, fill = ifelse(female == 0, "Male", "Female"))) +
geom_bar(stat = "identity", color = "black") +
scale_fill_manual(values = c("#1b9e77", "#d95f02")) +
annotate("text", x = 1, y = female_actual_pct,
label = paste0("Female Target = ", round(female_actual_pct*100, 1), "%"),
color = "#1b9e77", size = 4, vjust = 0) +
annotate("text", x = 2, y = male_actual_pct,
label = paste0("Male Target = ", round(male_actual_pct*100, 1), "%"),
color = "#d95f02", size = 4, vjust = 0) +
labs(title = "Agent Gender Distribution",
x = "",
y = "Percentage",
fill = "") +
theme_minimal()
agebreaks <- c(18, 25, 35, 45, 55, 65)
agelabels <- c("18-24", "25-34", "35-44", "45-54", "55-64")agent_dt[tick == last_tick, .N, ]## [1] 10000setDT(agent_dt)[ , age_groups := cut(age,
breaks = agebreaks,
include.lowest = TRUE,
right = FALSE,
labels = agelabels)]
nrow(agent_dt[tick == last_tick])## [1] 10000agent_dt[tick == last_tick, .N, by=c("age_groups", "race", "female")][order(age_groups, race, female)]agent_dt[tick == last_tick, .N, by=c("race", "age_groups", "female")][order(race, age_groups)]agent_dt[tick==last_tick, .N, by=c("race", "female")][,"%":=N/sum(N)*100][order(race)]agent_dt[tick==last_tick, .N, by=c("age_groups")][,"%":=round(N/sum(N)*100)][order(age_groups)]Median age at the start:
agent_dt[tick==1, median(age)]## [1] 51Median age at the end:
nrow(agent_dt[tick==last_tick])## [1] 10000agent_dt[tick==last_tick, median(age)]## [1] 63# Create age groups
agent_dt_by_age_group <- agent_dt[tick==last_tick, .N, by = .(age_group = cut(age, breaks = c(18, 25, 35, 45, 55, 65, 75, 80, 84)))]
# Create bar chart
ggplot(agent_dt_by_age_group, aes(x = age_group, y = N, fill = age_group)) +
geom_bar(stat = "identity", color = "black") +
scale_fill_manual(values = c("#e41a1c", "#377eb8", "#4daf4a", "#984ea3", "#ff7f00", "#ffff33", "#a65628", "#f781bf")) +
labs(title = "Agent Age Distribution",
x = "Age Group",
y = "Count",
fill = "Age Group") +
theme_minimal()
Baseline Smoking Behaviors
At the first tick:
agent_dt[tick==1 & race == "Black" & female == 0, .N, by=c("smoking_status")][]agent_dt[tick==1 & race == "Black" & female == 0, .N, by=c("smoking_status")][,"sim_prop":= round(N/sum(N), 3)][]# targets
BLACK_MALE_CURRENT = input_params$SMOKING_PREV$BLACK_MALE_CURRENT
BLACK_MALE_FORMER = input_params$SMOKING_PREV$BLACK_MALE_FORMER
BLACK_MALE_NEVER = input_params$SMOKING_PREV$BLACK_MALE_NEVER
smkg_prev_black_men_target <- c(BLACK_MALE_CURRENT, BLACK_MALE_FORMER, BLACK_MALE_NEVER)
smkg_prev_black_men_target## [1] 0.205 0.473 0.322## Add difference between the two target and the simAt the last tick:
smkg_blackmen_last_tick <- agent_dt[tick==last_tick & race == "Black" & female == 0, .N, by=c("smoking_status")][,"prop":=round(N/sum(N), 3)][]
smkg_blackmen_last_ticksmkg_prev_black_men_target## [1] 0.205 0.473 0.322# Count the number and percentage of each smoking status
smoking_dt <- agent_dt[tick == last_tick, .N, by = smoking_status][, `:=` (percentage = round(N/sum(N)*100, 2))]
# Create the bar chart
ggplot(smoking_dt, aes(x = smoking_status, y = percentage, fill = smoking_status)) +
geom_bar(stat = "identity", color = "black") +
scale_fill_manual(values = c("#fec44f", "#43a2ca", "#7bccc4")) +
labs(title = "Agent Smoking Status Distribution",
x = "Smoking Status",
y = "Percentage",
fill = "Smoking Status") +
theme_minimal()
Baseline Alcohol Use
alc_use_last_tick <-
agent_dt[tick == last_tick, .N, by = c("alc_use_status")][,
"%" := round(N /sum(N), 3)][
order(alc_use_status)][
]
# targets
alc_use_props <- input_params$ALC_USE_PROPS
alc_use_targets <- c(alc_use_props$A, alc_use_props$O, alc_use_props$R, alc_use_props$D)
cbind(alc_use_last_tick[["%"]], alc_use_targets)## alc_use_targets
## [1,] 0.083 0.083
## [2,] 0.730 0.729
## [3,] 0.131 0.132
## [4,] 0.055 0.056# Calculate the percentages for alcohol use status
alc_use_dt <- agent_dt[tick == last_tick, .N, by = alc_use_status][, `:=` (percentage = round(N/sum(N)*100, 2))]
# Order the data table by alcohol use status
setorder(alc_use_dt, alc_use_status)
# Rename the alcohol use status codes to descriptive labels
alc_use_dt[, alc_use_status := factor(alc_use_status, levels = c(0, 1, 2, 3),
labels = c("Lifetime Abstainer",
"Ocassional",
"Regular",
"Heavy"))]
ggplot(alc_use_dt, aes(x = alc_use_status, y = percentage, fill = alc_use_status)) +
geom_bar(stat = "identity", color = "black") +
labs(title = "Agent Alcohol Use Status Distribution",
x = "Use Frequency",
y = "Percentage",
fill = "Alcohol Use Status") +
theme_minimal() +
guides(fill = "none") +
geom_text(aes(label = paste0(as.numeric(input_params$ALC_USE_PROPS)*100, "%")),
position = position_stack(vjust = 1.02))+
scale_y_continuous(limits = c(0, 100), breaks = seq(0, 100, 10))
Smoking and Alcohol Co-Use
Conditional Probabilities on Co-Use
Approach 1
Write functions to compute co-use probabilities at any time step:
# Functions for computing probabilities at a specific tick
compute_p_smoker_given_alc_user <- function(tick_data) {
alc_users <- tick_data[tick_data$alc_use_status != 0,]
n_smokers_and_alc_users <-
sum(alc_users$smoking_status == "Current")
p_smoker_given_alc_user <-
n_smokers_and_alc_users / nrow(alc_users)
return(p_smoker_given_alc_user)
}
compute_p_alc_user_given_smoker <- function(tick_data) {
smokers <- tick_data[tick_data$smoking_status == "Current",]
n_smokers_and_alc_users <- sum(smokers$alc_use_status != 0)
p_alc_user_given_smoker <-
n_smokers_and_alc_users / nrow(smokers)
return(p_alc_user_given_smoker)
}
compute_p_abstainer_given_never_smoker <- function(tick_data) {
never_smokers <- tick_data[tick_data$smoking_status == "Never",]
n_abstainers_and_never_smokers <-
sum(never_smokers$alc_use_status == 0)
p_abstainer_given_never_smoker <-
n_abstainers_and_never_smokers / nrow(never_smokers)
return(p_abstainer_given_never_smoker)
}
compute_p_never_smoker_given_abstainer <- function(tick_data) {
abstainers <- tick_data[tick_data$alc_use_status == 0,]
n_abstainers_and_never_smokers <-
sum(abstainers$smoking_status == "Never")
p_never_smoker_given_abstainer <-
n_abstainers_and_never_smokers / nrow(abstainers)
return(p_never_smoker_given_abstainer)
}
compute_p_current_or_former_smoker_given_alc_user <-
function(tick_data) {
alc_users <- tick_data[tick_data$alc_use_status != 0,]
n_current_or_former_smokers_and_alc_users <-
sum(alc_users$smoking_status %in% c("Current", "Former"))
p_current_or_former_smoker_given_alc_user <-
n_current_or_former_smokers_and_alc_users / nrow(alc_users)
return(p_current_or_former_smoker_given_alc_user)
}
compute_p_alc_user_given_current_or_former_smoker <-
function(tick_data) {
current_or_former_smokers <-
tick_data[tick_data$smoking_status %in% c("Current", "Former"),]
n_current_or_former_smokers_and_alc_users <-
sum(current_or_former_smokers$alc_use_status != 0)
p_alc_user_given_current_or_former_smoker <-
n_current_or_former_smokers_and_alc_users / nrow(current_or_former_smokers)
return(p_alc_user_given_current_or_former_smoker)
}Apply these functions to the first and last tick data:
first_tick_data <- agent_dt[agent_dt$tick == 1,]
last_tick_data <- agent_dt[agent_dt$tick == max(tick)]compute_p_smoker_given_alc_user(first_tick_data)## [1] 0.1152335compute_p_alc_user_given_smoker(first_tick_data)## [1] 0.9134948compute_p_abstainer_given_never_smoker(first_tick_data)## [1] 0.08655536compute_p_never_smoker_given_abstainer(first_tick_data)## [1] 0.4712919compute_p_current_or_former_smoker_given_alc_user(first_tick_data)## [1] 0.546268compute_p_alc_user_given_current_or_former_smoker(first_tick_data)## [1] 0.9188693compute_p_smoker_given_alc_user(last_tick_data)## [1] 0.1411431compute_p_alc_user_given_smoker(last_tick_data)## [1] 0.915782compute_p_abstainer_given_never_smoker(last_tick_data)## [1] 0.08287536compute_p_never_smoker_given_abstainer(last_tick_data)## [1] 0.453125compute_p_current_or_former_smoker_given_alc_user(last_tick_data)## [1] 0.5449389compute_p_alc_user_given_current_or_former_smoker(last_tick_data)## [1] 0.9165291# Create a contingency table of smoking and alcohol use status
smoking_alc_use_table <- dcast(agent_dt[tick == last_tick], smoking_status ~ alc_use_status, value.var = "id", fun.aggregate = length)
smoking_alc_use_table <- smoking_alc_use_table[, .(smoking_status, `Lifetime Abstainer` = `0`, Occasional = `1`, Regular = `2`, Heavy = `3`)]
smoking_alc_use_melt <- melt(smoking_alc_use_table, id.vars = "smoking_status",
variable.name = "alcohol_use")Approach II
Create a contingency table at the last tick:
smkg_alc_use_xtab <- dcast(smoking_alc_use_melt, smoking_status ~ alcohol_use, value.var = "value")
smkg_alc_use_xtabsmoking_alc_col_props <- prop.table(as.matrix(smkg_alc_use_xtab[,-1]), margin = 1)
smoking_alc_row_props <- prop.table(as.matrix(smkg_alc_use_xtab[,-1]), margin = 2)
smoking_alc_col_props <- cbind(smkg_alc_use_xtab[,1], round(smoking_alc_col_props, 3))
smoking_alc_row_props <- cbind(smkg_alc_use_xtab[,1], round(smoking_alc_row_props, 3))
print(smoking_alc_col_props)## smoking_status Lifetime Abstainer Occasional Regular Heavy
## 1: Current 0.084 0.713 0.154 0.049
## 2: Former 0.083 0.738 0.123 0.055
## 3: Never 0.083 0.729 0.131 0.058print(smoking_alc_row_props)## smoking_status Lifetime Abstainer Occasional Regular Heavy
## 1: Current 0.143 0.138 0.166 0.125
## 2: Former 0.404 0.408 0.380 0.401
## 3: Never 0.453 0.454 0.454 0.475rowSums(smoking_alc_col_props[,-1])## [1] 1.000 0.999 1.001colSums(smoking_alc_row_props[,-1])## Lifetime Abstainer Occasional Regular Heavy
## 1.000 1.000 1.000 1.001P(Current Smoker |Heavy Drinker) = 0.125
P(Current Smoker | Regular Drinker) = 0.166 P(Current Smoker |
Occasional Drinker) = 0.138
Alcohol Use conditioned on Current Smoking: P(Heavy Drinker | Current
Smoker) = 0.049 P(Regular Drinker| Current Smoker) = 0.154
P(Occasional Drinker | Current Smoker) = 0.713
Marginal Distributions on Alcohol and Tobacco Use
Start by computing the mean distributions of smoking and alcohol use statesa cross the simulation
# Compute the proportion of each smoking status
smoking_prop <- prop.table(table(agent_dt$smoking_status))
# Compute the mean proportion of current, former, and never smokers
mean_current_smokers <- mean(smoking_prop["Current"])
mean_former_smokers <- mean(smoking_prop["Former"])
mean_never_smokers <- mean(smoking_prop["Never"])
# Print the mean smoking status distribution
cat("Mean smoking status distribution:\n")## Mean smoking status distribution:cat("Current Smokers: ", mean_current_smokers, "\n")## Current Smokers: 0.1388275cat("Former Smokers: ", mean_former_smokers, "\n")## Former Smokers: 0.4040381cat("Never Smokers: ", mean_never_smokers, "\n")## Never Smokers: 0.4571344# Compute the proportion of each alcohol use status
alc_prop <- prop.table(table(agent_dt$alc_use_status))
# Compute the mean proportion of each alcohol use status
mean_abstainers <- mean(alc_prop[1])
mean_occasional <- mean(alc_prop[2])
mean_regular <- mean(alc_prop[3])
mean_heavy <- mean(alc_prop[4])
# Print the mean alcohol use status distribution
cat("Mean alcohol use status distribution:\n")## Mean alcohol use status distribution:cat("Abstainers: ", mean_abstainers, "\n")## Abstainers: 0.08269079cat("Occasional Users: ", mean_occasional, "\n")## Occasional Users: 0.729248cat("Regular Users: ", mean_regular, "\n")## Regular Users: 0.131906cat("Heavy Users: ", mean_heavy, "\n")## Heavy Users: 0.05615525Joint Probability Computations of Alcohol and Tobacco Use
Joint probability of current smoking and alcohol use (i.e., not lifetime abstainer):
# Filter the dataset to include only current smokers who are also alcohol users
current_smokers_alcohol_users_dt <- agent_dt[agent_dt$smoking_status == "Current" & agent_dt$alc_use_status %in% c(1, 2, 3),]
# Compute the proportion of individuals who are both current smokers and alcohol users
joint_prob_current_smokers_alcohol_users <- nrow(current_smokers_alcohol_users_dt) / nrow(agent_dt)
# Print the joint probability of current smokers and alcohol users
cat("Joint probability of current smokers and alcohol users: ", joint_prob_current_smokers_alcohol_users, "\n")## Joint probability of current smokers and alcohol users: 0.1273249Joint probability of former or current/former smoking and alcohol use (i.e., not lifetime abstainer):
# Filter the dataset to include only current smokers who are also alcohol users
former_or_current_smokers_alcohol_users_dt <- agent_dt[agent_dt$smoking_status %in% c("Current", "Former") & agent_dt$alc_use_status %in% c(1, 2, 3),]
# Compute the proportion of individuals who are both current smokers and alcohol users
joint_prob_former_or_current_smokers_alcohol_users <- nrow(former_or_current_smokers_alcohol_users_dt) / nrow(agent_dt)
# Print the joint probability of current smokers and alcohol users
cat("Joint probability of former or current smokers and alcohol users: ", joint_prob_former_or_current_smokers_alcohol_users, "\n")## Joint probability of former or current smokers and alcohol users: 0.4978723Joint probability of never smoking and lifetime abstention from alcohol use:
# Filter the dataset to include only current smokers who are also alcohol users
no_smoking_or_alcohol_users_dt <- agent_dt[agent_dt$smoking_status %in% c("Never") & agent_dt$alc_use_status %in% c(0),]
# Compute the proportion of individuals who are both current smokers and alcohol users
joint_prob_no_smoking_or_alcohol_users <- nrow(no_smoking_or_alcohol_users_dt) / nrow(agent_dt)
# Print the joint probability of current smokers and alcohol users
cat("Joint probability of never smoker and lifetime alcohol abstainer: ",
joint_prob_no_smoking_or_alcohol_users, "\n")## Joint probability of never smoker and lifetime alcohol abstainer: 0.03769749Joint probability of never/former smoking and lifetime abstention from alcohol use:
# Filter the dataset to include only current smokers who are also alcohol users
former_or_never_smoking_and_no_alcohol_users_dt <- agent_dt[agent_dt$smoking_status %in%
c("Never", "Former") &
agent_dt$alc_use_status %in% c(0),]
# Compute the proportion of individuals who are both current smokers and alcohol users
joint_prob_former_or_never_smoking_and_no_alcohol_users <- nrow(former_or_never_smoking_and_no_alcohol_users_dt) /
nrow(agent_dt)
# Print the joint probability of current smokers and alcohol users
cat("Joint probability of former/never smoker and lifetime alcohol abstainer: ",
joint_prob_former_or_never_smoking_and_no_alcohol_users, "\n")## Joint probability of former/never smoker and lifetime alcohol abstainer: 0.07118815Joint probability of current/former smoking and lifetime abstention from alcohol use:
# Filter the dataset to include only current/former smokers who are lifetime abstainers from alcohol
current_former_smokers_lifetime_abstainers_dt <- agent_dt[agent_dt$smoking_status %in% c("Current", "Former") & agent_dt$alc_use_status == 0,]
# Compute the proportion of individuals who are both current/former smokers and lifetime abstainers from alcohol
joint_prob_current_former_smokers_lifetime_abstainers <- nrow(current_former_smokers_lifetime_abstainers_dt) / nrow(agent_dt)
# Print the joint probability of current/former smokers and lifetime abstainers from alcohol use
cat("Joint probability of current/former smokers and lifetime abstainers from alcohol use: ", joint_prob_current_former_smokers_lifetime_abstainers, "\n")## Joint probability of current/former smokers and lifetime abstainers from alcohol use: 0.0449933Joint probability of never smoking and using alcohol:
# Filter the dataset to include only alcohol users who have never smoked
alcohol_users_never_smokers_dt <- agent_dt[agent_dt$alc_use_status %in% c(1, 2, 3) & agent_dt$smoking_status == "Never",]
# Compute the proportion of individuals who are alcohol users but never smokers
joint_prob_alcohol_users_never_smokers <- nrow(alcohol_users_never_smokers_dt) / nrow(agent_dt)
# Print the joint probability of alcohol users but never smokers
cat("Joint probability of alcohol users but never smokers: ", joint_prob_alcohol_users_never_smokers, "\n")## Joint probability of alcohol users but never smokers: 0.4194369Compare the conditional, marginal, and joint probability distributions to the Falk et al. (or other empirical) data.
# Create a stacked bar chart of smoking and alcohol use status
ggplot(smoking_alc_use_melt, aes(x = smoking_status, y = value, fill = alcohol_use)) +
geom_bar(stat = "identity", color = "black") +
labs(title = "Agent Smoking and Alcohol Use Status Distribution",
x = "Smoking Status",
y = "Count",
fill = "Alcohol Use Status") +
theme_minimal() +
scale_fill_manual(values = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c")) +
guides(fill = "none")
# Create a contingency table of smoking and alcohol use status
smoking_alc_use_table <- dcast(agent_dt[tick == last_tick], smoking_status ~ alc_use_status, value.var = "id", fun.aggregate = length)
smoking_alc_use_table <- smoking_alc_use_table[, .(smoking_status, `Lifetime Abstainer` = `0`, Occasional = `1`, Regular = `2`, Heavy = `3`)]
smoking_alc_use_melt <- melt(smoking_alc_use_table, id.vars = "smoking_status")
# Create a stacked bar chart of smoking and alcohol use status
ggplot(smoking_alc_use_melt, aes(x = smoking_status, y = value, fill = variable)) +
geom_bar(stat = "identity", color = "black") +
labs(title = "Agent Smoking and Alcohol Use Status Distribution",
x = "Smoking Status",
y = "Count") +
scale_fill_manual(values = c("#a6cee3", "#1f78b4", "#b2df8a", "#33a02c"),
labels = c("Lifetime Abstainer", "Occasional", "Regular", "Heavy"),
name = "Alcohol Use Status") +
theme_minimal()
# Age Assignment at Entry
We check below the ages of new agents, and if their initial age assignments and age increments are working correctly.
How many agents are present at the last tick:
agent_dt[tick == last_tick, .N] ## [1] 10000How many agents present at tick 991 entered after model initialization?
agent_dt[tick == 991 & id > 9999, .N] ## [1] 422The mean age of agents at time by time of entry at tick 1:
agent_dt[tick == 1,
.("mean_age_today"=mean(age), .N),
by=.(entry_at_tick),
]At tick 491:
agent_dt[tick == 491,
.("mean_age_today"=mean(age), .N),
by=.(entry_at_tick),
]And at tick 991:
agent_dt[tick == 991,
.("mean_age_today"=mean(age), .N),
by=.(entry_at_tick),
]These data seem to match expectations of aging of newly entering agents. See that at tick 1, only one agent had entered after model initialization. and this agent was 18 years old.
At tick 491, this agent was 19 years old (491-365 years older than 18, rounded), and at tick 991, this agent was 21 years old (991-365 years older than 18, rounded).
So, the aging process seems to work as expected.
Incarceration Statistics
Number and proportion of incarcerated persons at last tick:
agent_dt[tick==last_tick, .N, by=c("current_incarceration_status")][,"prop":=round(N/sum(N)*100, 0)][]Break down number and proportions of incarcerated persons at last tick by race, sex and age:
agent_dt[tick==last_tick & current_incarceration_status == 1,
.N,
by=c("race")][,"prop":=round(N/sum(N)*100, 0)][]agent_dt[tick==last_tick & current_incarceration_status == 1,
.N,
by=c("female")][,"prop":=round(N/sum(N)*100, 0)][]agent_dt[tick==last_tick & current_incarceration_status == 1,
.N,
by=c("age_groups")][,"prop":=round(N/sum(N)*100, 0)][]agent_dt[tick==last_tick & current_incarceration_status == 1,
.N,
by=c("race", "female", "age_groups")][order(race, female, age_groups)][,"prop":=round(N/sum(N)*100, 0)][]How many persons have been ever incarcerated?
To answer this question, we can group the last incarceration times by agent ID. If last incarceration time is -1, the person has never been incarcerated. If last incarceration time is > -1, the person has been incarcerated at least once.
# Group data by agent ID and compute the last incarceration time
last_incarceration_time <- agent_dt[,
.(last_incarceration_time = max(last_incarceration_tick)),
by = id]
# Count the number of agents who have been incarcerated at least once
n_ever_incarcerated <- sum(last_incarceration_time$last_incarceration_time != -1)
# Compute the proportion of agents who have been incarcerated at least once
prop_ever_incarcerated <- n_ever_incarcerated / nrow(last_incarceration_time)
# Agents who have been incarcerated at least once
ever_incarcerated_times <- last_incarceration_time[last_incarceration_time > 1]
dim(ever_incarcerated_times)## [1] 740 2# Extract demographic information for ever incarcerated agents
ever_incarcerated_ids <-
unique(agent_dt[id %in% ever_incarcerated_times$id, id])
ever_incarcerated_info <- agent_dt[id %in% ever_incarcerated_ids,
.SD[.N, list(
id,
age = last(age),
race = last(race),
female = last(female)
)],
by = id]
# Compute distributions
## mean age
age_summary <- ever_incarcerated_info[, .(mean_age=mean(age),
min_age=min(age),
max_age=max(age)
)
]
sex_dist <-
ever_incarcerated_info[, .N, by = female][, .(
Sex = c("Female", "Male"),
Count = c(N[female == 0], N[female == 1]),
"Proportion of ever incarcerated" = round(N / sum(N), 2)
)]
print(sex_dist)## Sex Count Proportion of ever incarcerated
## 1: Female 382 0.48
## 2: Male 358 0.52race_dist <- ever_incarcerated_info[, .N, by = race][,
.(Race = race,
Count = N,
Proportion = round(N / sum(N), 2)
)
]
race_dist[]Thus, 740 persons have been incarcerated at least once. The proportion of persons who have been ever incarcerated is 0.047005.
Now let us compute the number of females and males whoa have ever been incarcerated, and the proportions by race (% of White, Black, …)
first_tick_race_sex <-
agent_dt[, .SD[1], by = id][, .(
ID = id,
Tick = tick,
Female = female,
Race = race
)]
# Count number of females and males
sex_count <- first_tick_race_sex[, .N, by = Female][,
.(Sex = c("Female", "Male"),
Count = c(N[Female == 1], N[Female == 0])) ]
print(sex_count)## Sex Count
## 1: Female 7925
## 2: Male 7818# Count number of individuals in each race
race_count <- first_tick_race_sex[, .N, by = Race][,
.(Race, Count = N)]
print(race_count)## Race Count
## 1: White 11177
## 2: Asian 583
## 3: Hispanic 2636
## 4: Black 1347Thus, the distribution of ever incarcerated persons by sex is
cbind(sex_count$Sex, round(as.numeric(sex_dist$Count/sex_count$Count), 3))## [,1] [,2]
## [1,] "Female" "0.048"
## [2,] "Male" "0.046"The distribution of ever incarcerated persons by race is:
cbind(race_count$Race, round(as.numeric(race_dist$Count/race_count$Count), 3))## [,1] [,2]
## [1,] "White" "0.048"
## [2,] "Asian" "0.101"
## [3,] "Hispanic" "0.044"
## [4,] "Black" "0.018"What is the duration of incarceration for persons who have been ever incarcerated? Restricting our data to only persons who have been released from their last incarceration episode, we can estimate this as:
# Filter the data to include only the observations for agents who have been ever incarcerated and released
ever_incarcerated_released <-
agent_dt[last_release_tick > last_incarceration_tick,
.SD[1, ], # ensure times for each incarceration, release episode appear only once
by = .(id, last_incarceration_tick)]
# Compute the length of incarceration time for each ever incarcerated and released agent
ever_incarcerated_released[, incarceration_length := last_release_tick - last_incarceration_tick,
by = id]
# Compute the summary statistics of the lengths of incarceration time for the ever incarcerated and released agents
incarceration_summary <-
ever_incarcerated_released[, .(
min_incarceration_length = min(incarceration_length),
median_incarceration_length = median(incarceration_length),
max_incarceration_length = max(incarceration_length)
),
by = id]
# Compute the distribution of incarceration times for all incarceration episodes followed by release
summary(ever_incarcerated_released$incarceration_length)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 7.0 21.0 55.0 114.9 128.0 2166.0Assess the total number of incarceration-release events and the number of unique agents who experienced them:
length(ever_incarcerated_released$id)## [1] 1219length(unique(ever_incarcerated_released$id))## [1] 730Distribution of incarceration events:
hist(ever_incarcerated_released$incarceration_length,
breaks = seq(0, 3000, by = 30),
xlab = "Length of incarceration (days)",
ylab = "Number of individuals",
main = "Distribution of incarceration lengths")
Break down number and proportions of persons ever incarcerated at last tick by race, sex and age:
ever_inc_dt <-
agent_dt[tick==last_incarceration_tick & n_incarcerations > 0,
.N,
by=c("race", "female", "age_groups")][order(race, female, age_groups)][
,"prop":=round(N/sum(N)*100, 0)][]
ever_inc_dtever_inc_dt[, colSums(.SD), .SDcols = 4:5]## N prop
## 99 99Mean number of incarcerations at last tick:
agent_dt[tick == last_tick,
.(
"mean_n_incarcerations" = mean(n_incarcerations),
"min_n_incarcerations" = min(n_incarcerations),
"max_n_incarcerations" = max(n_incarcerations),
"median_n_incarcerations" = median(n_incarcerations),
"n_ever_incarcerated_agents"=sum(n_incarcerations > 0),
.N)] Incarceration Duration
Compute statistics on agents whose last release date is after their last incarceration date.
agent_dt[tick == last_tick & last_release_tick >= last_incarceration_tick,
.(
"mean_duration" = mean(last_release_tick - last_incarceration_tick),
"max_duration" = max(last_release_tick - last_incarceration_tick),
"min_duration" = min(last_release_tick - last_incarceration_tick),
"median_duration" = median(last_release_tick - last_incarceration_tick),
.N)] Release Statistics
Number
agent_dt[tick == last_tick,
.(
"mean_n_releases" = mean(n_releases),
"min_n_releases" = min(n_releases),
"max_n_releases" = max(n_releases),
"median_n_releases" = median(n_releases),
"n_ever_released_agents"=sum(n_releases > 0),
.N)] Duration
Compute statistics on agents whose last release date is *before* their last incarceration date.
agent_dt[tick == last_tick & last_release_tick < last_incarceration_tick,
.(
"mean_rel_dur" = mean(abs(last_release_tick - last_incarceration_tick)),
"max_rel_dur" = max(abs(last_release_tick - last_incarceration_tick)),
"min_rel_dur" = min(abs(last_release_tick - last_incarceration_tick)),
"median_rel_dur" = median(abs(last_release_tick - last_incarceration_tick)),
.N)] Smoking transition behavior
agent_dt[tick == last_tick,
.(
"mean_n_smkg_trans" = mean(n_smkg_stat_trans),
"max_n_smkg_trans" = max(n_smkg_stat_trans),
"min_n_smkg_trans" = min(n_smkg_stat_trans),
"median_n_smkg_trans" = median(n_smkg_stat_trans),
.N)] Smoking State Distributions:
ans <- agent_dt[tick == last_tick, .N, by=c("race", "smoking_status")][order(race)]
ans ans[, prop := N / sum(N), by = race][]ans_num_smokers <- ans[smoking_status %in% c("Former", "Current", "Never"),
.(num_smokers = sum(N)),
by = race]
ans_num_smokerssum(ans_num_smokers$num_smokers)## [1] 10000sim_prop_smokers_by_race <- ans_num_smokers$num_smokers/sum(ans_num_smokers$num_smokers)Alcohol transition behavior
agent_dt[tick == last_tick,
.(
"mean_n_alc_trans" = mean(n_alc_use_stat_trans),
"max_n_alc_trans" = max(n_alc_use_stat_trans),
"min_n_alc_trans" = min(n_alc_use_stat_trans),
"median_n_alc_trans" = median(n_alc_use_stat_trans),
.N)] Alcohol Use State Distributions:
ans <- agent_dt[tick == last_tick, .N, by=alc_use_status]
prop <- ans$N/sum(ans$N)
cbind(ans[,1], prop)Smoking and Alcohol Use among Incarcerated Persons vs Not-
Smoking
Compare the proportion of the current smokers (to all persons) at the last tick in the n_releases = 0 group vs the n_releases = 1 group.
agent_dt[tick == last_tick, .N, by = n_releases][]agent_dt[tick == last_tick, .N, by = smoking_status][]agent_dt[tick == last_tick, .N, by = c("smoking_status", "n_releases")][]# Filter data at the last_tick
last_tick_data <- agent_dt[tick == last_tick]
# Calculate the ratio for each group
ratios <- last_tick_data[, .(
n_current_smokers = sum(smoking_status == "Current"),
n_total_persons = sum(smoking_status %in% c("Current", "Former", "Never"))
), by = .(group = n_releases >= 1)]
# Compute the ratio of current smokers to total smokers in each group
ratios[, ratio := n_current_smokers / n_total_persons]
# Print the results
ratiosThe population size at each time tick is
sum(agent_dt[tick == last_tick, .N, by = c("smoking_status", "n_releases")][["N"]])## [1] 10000Now let us compute this comparison-of-ratios across all time steps, i.e., mean across all time steps of the current smoking proportion in the n_releases = 0 group vs. the n_releases >= 1 group. We add the standard deviation and the interquartile range as measures of variation across both groups.
# Calculate the ratio for each group and each tick
ratios_by_tick <- agent_dt[, .(
n_current_smokers = sum(smoking_status == "Current"),
n_total_persons = sum(smoking_status %in% c("Current", "Former", "Never"))
), by = .(tick, group = n_releases >= 1)]
# Compute the ratio of current smokers to total smokers in each group and each tick
ratios_by_tick[, ratio := n_current_smokers / n_total_persons]
# Calculate the mean ratio, standard deviation, and interquartile range across all time steps for each group
summary_by_group <- ratios_by_tick[, .(
mean_ratio = mean(ratio, na.rm = TRUE),
sd_ratio = sd(ratio, na.rm = TRUE),
iqr_ratio = IQR(ratio, na.rm = TRUE)
), by = group]
# Print the results
summary_by_groupTest if the means for the two groups are statistically different:
# Extract the ratio data for each group
group_0_ratios <- ratios_by_tick[group == FALSE, ratio]
group_1_ratios <- ratios_by_tick[group == TRUE, ratio]
# Perform the two-sample t-test
t_test_result <- t.test(group_0_ratios, group_1_ratios)
# Print the t-test result
t_test_result##
## Welch Two Sample t-test
##
## data: group_0_ratios and group_1_ratios
## t = 0.8956, df = 4387.4, p-value = 0.3705
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.006108857 0.016384181
## sample estimates:
## mean of x mean of y
## 0.1446523 0.1395146Compute change in current smoking rate over time since release
calculate_proportions_smoking <- function(data) {
proportions <- prop.table(table(data$smoking_status))
df <- data.frame(Category = names(proportions),
Proportion = as.numeric(proportions))
df$Category <- factor(df$Category,
levels = c("Never", "Former", "Current"),
labels = c("Never", "Former", "Current"))
return(df)
}
summary_after_release <- function(time_period) {
agents_within_time_period <- agent_dt[last_release_tick > 0 &
tick - last_release_tick <= time_period]
proportions <- calculate_proportions_smoking(agents_within_time_period)
agent_count <- nrow(agents_within_time_period)
return(list("proportions" = proportions, "agent_count" = agent_count))
}
time_periods <- c(1, 7, 14, 30, 90, 180, 365)
summaries_smoking <- lapply(time_periods, function(x) summary_after_release(x))
summaries_smoking## [[1]]
## [[1]]$proportions
## Category Proportion
## 1 Current 0.2283465
## 2 Former 0.2519685
## 3 Never 0.5196850
##
## [[1]]$agent_count
## [1] 127
##
##
## [[2]]
## [[2]]$proportions
## Category Proportion
## 1 Current 0.2240566
## 2 Former 0.3125000
## 3 Never 0.4634434
##
## [[2]]$agent_count
## [1] 848
##
##
## [[3]]
## [[3]]$proportions
## Category Proportion
## 1 Current 0.2004648
## 2 Former 0.3317838
## 3 Never 0.4677513
##
## [[3]]$agent_count
## [1] 1721
##
##
## [[4]]
## [[4]]$proportions
## Category Proportion
## 1 Current 0.1817684
## 2 Former 0.3607993
## 3 Never 0.4574322
##
## [[4]]$agent_count
## [1] 3653
##
##
## [[5]]
## [[5]]$proportions
## Category Proportion
## 1 Current 0.1504653
## 2 Former 0.3901225
## 3 Never 0.4594121
##
## [[5]]$agent_count
## [1] 10853
##
##
## [[6]]
## [[6]]$proportions
## Category Proportion
## 1 Current 0.1400676
## 2 Former 0.3972493
## 3 Never 0.4626831
##
## [[6]]$agent_count
## [1] 21304
##
##
## [[7]]
## [[7]]$proportions
## Category Proportion
## 1 Current 0.1394853
## 2 Former 0.3942471
## 3 Never 0.4662676
##
## [[7]]$agent_count
## [1] 41266Alcohol Use
Compare the alcohol use rates in the: (1) entire agent population; (2) agent population that has been never released; (3) agent population has been released at least once:
# Entire agent population
all_agents <- agent_dt
# Agents who have never been released (n_releases == 0)
never_released_agents <- agent_dt[n_releases == 0]
# Agents who have been released at least once (n_releases >= 1)
released_agents <- agent_dt[n_releases >= 1]
# Function to calculate proportions
calculate_proportions <- function(data) {
proportions <- prop.table(table(data$alc_use_status))
df <- data.frame(Category = names(proportions),
Proportion = as.numeric(proportions))
df$Category <- factor(df$Category,
levels = c(0, 1, 2, 3),
labels = c("Lifetime Abstainer", "Occasional User", "Regular User", "Heavy User"))
return(df)
}
# Entire agent population
all_agents_summary <- calculate_proportions(all_agents)
# Agents who have never been released
never_released_agents_summary <- calculate_proportions(never_released_agents)
# Agents who have been released at least once
released_agents_summary <- calculate_proportions(released_agents)
# Display summaries
print("Summary for Entire Agent Population")## [1] "Summary for Entire Agent Population"print(all_agents_summary)## Category Proportion
## 1 Lifetime Abstainer 0.08269079
## 2 Occasional User 0.72924799
## 3 Regular User 0.13190597
## 4 Heavy User 0.05615525print("Summary for Agents Who Have Never Been Released")## [1] "Summary for Agents Who Have Never Been Released"print(never_released_agents_summary)## Category Proportion
## 1 Lifetime Abstainer 0.08273326
## 2 Occasional User 0.72941171
## 3 Regular User 0.13185987
## 4 Heavy User 0.05599516print("Summary for Agents Who Have Been Released at Least Once")## [1] "Summary for Agents Who Have Been Released at Least Once"print(released_agents_summary)## Category Proportion
## 1 Lifetime Abstainer 0.08133877
## 2 Occasional User 0.72403571
## 3 Regular User 0.13337366
## 4 Heavy User 0.06125187Once agents are released, the baseline state transition probabilities take effect, which force the alcohol use rates to equilibrate at the baseline state distribution of alcohol use. This might be the reason why the heavy alcohol use rates do not differ between the groups that have never been released vs. the group that has been released at least once.
To account for this effect of baseline transition probabilities, we will instead compare the alcohol use rates shortly after release, i.e., we will compute the alcohol use rates for a certain period after the agents are released. This may help uncover the immediate effect of release on alcohol use rates.
# Define a one-year period after release
fixed_time_after_release <- 1
# 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 Lifetime Abstainer 0.11023622
## 2 Occasional User 0.58267717
## 3 Regular User 0.06299213
## 4 Heavy User 0.24409449# Note for The proportion of heavy users one day after release is in the right range: approximately 24.4094488%. But the proportion of lifetime abstainers at 11.023622 seems higher than it should be.
Sensitivity Analysis of time since release:
# Function to filter agents and calculate proportions for the given time period after release
summary_after_release <- function(time_period) {
agents_within_time_period <- agent_dt[last_release_tick > 0 & tick <= last_release_tick + time_period]
proportions <- calculate_proportions(agents_within_time_period)
return(proportions)
}
# 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))
# Display summaries
names(summaries) <- paste(time_periods, "days after release")
summaries## $`1 days after release`
## Category Proportion
## 1 Lifetime Abstainer 0.11023622
## 2 Occasional User 0.58267717
## 3 Regular User 0.06299213
## 4 Heavy User 0.24409449
##
## $`7 days after release`
## Category Proportion
## 1 Lifetime Abstainer 0.08490566
## 2 Occasional User 0.62146226
## 3 Regular User 0.11674528
## 4 Heavy User 0.17688679
##
## $`14 days after release`
## Category Proportion
## 1 Lifetime Abstainer 0.0889018
## 2 Occasional User 0.6275421
## 3 Regular User 0.1051714
## 4 Heavy User 0.1783847
##
## $`30 days after release`
## Category Proportion
## 1 Lifetime Abstainer 0.08020805
## 2 Occasional User 0.64193813
## 3 Regular User 0.11251027
## 4 Heavy User 0.16534355
##
## $`90 days after release`
## Category Proportion
## 1 Lifetime Abstainer 0.08430849
## 2 Occasional User 0.65603980
## 3 Regular User 0.12411315
## 4 Heavy User 0.13553856
##
## $`180 days after release`
## Category Proportion
## 1 Lifetime Abstainer 0.08045437
## 2 Occasional User 0.68221930
## 3 Regular User 0.12739392
## 4 Heavy User 0.10993241
##
## $`365 days after release`
## Category Proportion
## 1 Lifetime Abstainer 0.08009015
## 2 Occasional User 0.70101294
## 3 Regular User 0.13061600
## 4 Heavy User 0.08828091time_points <- c(1, 7, 14, 30, 90, 180, 365)
labels <- c("1D", "1W", "2W", "1M", "3M", "6M", "1Y")
heavy_user_proportions <- sapply(summaries, function(x) x[x$Category == "Heavy User", "Proportion"])
heavy_user_df <- data.frame(Time = time_points,
Labels = labels,
Proportion = heavy_user_proportions)
ggplot(heavy_user_df, aes(x = Time, y = Proportion, group = 1)) +
geom_line() +
geom_point() +
scale_x_continuous(breaks = time_points, labels = labels) +
labs(title = "Heavy User Proportions After Release",
x = "Time After Release",
y = "Proportion of Heavy Users") +
theme_minimal()