R Notebook

rm(list=ls())
library(data.table)
library(here)
library(DiagrammeR)
library(yaml)
library(igraph)
library(ggraph)
library(ggplot2)
here::here("~/code/cadre/data-analysis-plotting/Simulated-Data-Analysis/r")

[1] "~/code/cadre/data-analysis-plotting/Simulated-Data-Analysis/r"

input_params <- read_yaml("../../../python/myparams/model_params.yaml") #input params
network_dt <- fread("../../../python/output/network_log_11.csv") #simulated network data
agent_dt <- fread("../../../python/output/agent_log_11.csv") #simulated network data
last_tick <- max(network_dt$tick) 

Let’s focus on the last tick for which we have agent and network data:

last_tick_network_dt <- network_dt[tick==last_tick,]
last_tick_agent_dt <- agent_dt[tick == last_tick,]
print(last_tick_network_dt)

print(last_tick_agent_dt)

Visualize using DiagrammeR:

# Identify the IDs of the recently released agents
recently_released_agents <- last_tick_agent_dt$id[last_tick_agent_dt$last_release_tick > (last_tick - 365)]

# Get the network data for the recently released agents
network_recently_released <- last_tick_network_dt[(last_tick_network_dt$p1 %in% recently_released_agents) | 
                                                    (last_tick_network_dt$p2 %in% recently_released_agents),]

# Get the first-degree network (neighbors) for each agent
first_degree_neighbors <- unique(c(network_recently_released$p1, network_recently_released$p2))

# Get agent data for the first-degree neighbors
first_degree_neighbors_agent_data <- last_tick_agent_dt[last_tick_agent_dt$id %in% first_degree_neighbors,]

# Create an edge data frame from your network data
edf <- data.frame(from = network_recently_released$p1, 
                  to = network_recently_released$p2)

# Create a node data frame with your agent data
ndf <- data.frame(id = first_degree_neighbors_agent_data$id, 
                  smoking_status = first_degree_neighbors_agent_data$smoking_status,
                  alc_use_status = first_degree_neighbors_agent_data$alc_use_status,
                  recently_released = ifelse(first_degree_neighbors_agent_data$id %in% recently_released_agents, "yes", "no"))

# Create a graph with DiagrammeR
graph <- create_graph(nodes_df = ndf, 
                      edges_df = edf,
                      directed = FALSE)

# Customize the node fillcolor based on smoking_status, alc_use_status and recently_released
graph <- set_node_attrs(graph, "fillcolor", ifelse(ndf$smoking_status == "Current", "red", 
                                                   ifelse(ndf$alc_use_status == 3, "blue", 
                                                          ifelse(ndf$recently_released == "yes", "black", "white"))))

# Customize the node color based on smoking_status, alc_use_status and recently_released
graph <- set_node_attrs(graph, "color", ifelse(ndf$smoking_status == "Current", "red", 
                                               ifelse(ndf$alc_use_status == 3, "blue", 
                                                      ifelse(ndf$recently_released == "yes", "black", "gray"))))

# Customize the node shape based on recently_released
graph <- set_node_attrs(graph, "shape", ifelse(ndf$recently_released == "yes", "square", "circle"))

# Remove labels from nodes to make them stand out more
graph <- set_node_attrs(graph, "label", "")

# Customize the edge color and line type
graph <- set_edge_attrs(graph, "color", "black")  # Change the edge color to gray
graph <- set_edge_attrs(graph, "style", "solid")  # Change the edge line type to dashed

# Customize the edge thickness
graph <- set_edge_attrs(graph, "penwidth", 2)  # Increase the edge thickness

# Customize the outline color of all nodes
graph <- set_node_attrs(graph, "color", "darkgray") 

# Render the graph
render_graph(graph)

# Create a legend
legend <- rbind(data.frame(group = "Recently Released", shape = "square", color = "black"),
                data.frame(group = "Currently Smoking", shape = "circle", color = "red"),
                data.frame(group = "AUD", shape = "circle", color = "blue"))

Prevalence of Current Smoking and AUD among network members:

# Number of agents with Current Smoking status
num_current_smoking <- sum(first_degree_neighbors_agent_data$smoking_status == "Current")

# Number of agents with AUD
num_aud <- sum(first_degree_neighbors_agent_data$alc_use_status == 3)

# Total number of agents in the network
total_agents <- nrow(first_degree_neighbors_agent_data)

# Calculate prevalence
prevalence_smoking <- num_current_smoking / total_agents
prevalence_aud <- num_aud / total_agents

# Print the results
cat("Prevalence of Current Smoking in networks of recently released agents: ", prevalence_smoking, "\n")

Prevalence of Current Smoking in networks of recently released agents:  0.1899441 

cat("Prevalence of AUD in networks of recently released agents: ", prevalence_aud, "\n")

Prevalence of AUD in networks of recently released agents:  0.0726257 

Prevalence of Cyrrent Smoking and AUD among this set of egos:

# Define selected_agents_df, which contains the data for the selected agents
recently_released_agents_df <- last_tick_agent_dt[last_tick_agent_dt$id %in% recently_released_agents,]

# Number of selected agents who are current smokers
num_ego_current_smoking <- sum(recently_released_agents_df$smoking_status == "Current")

# Number of selected agents with AUD
num_ego_aud <- sum(recently_released_agents_df$alc_use_status == 3)

# Total number of selected agents
total_ego_agents <- length(recently_released_agents)

# Calculate prevalence
prevalence_ego_smoking <- num_ego_current_smoking / total_ego_agents
prevalence_ego_aud <- num_ego_aud / total_ego_agents

# Print the results
cat("Prevalence of Current Smoking among recently released agents: ", prevalence_ego_smoking, "\n")

Prevalence of Current Smoking among recently released agents:  0.1944444 

cat("Prevalence of AUD among recently released agents: ", prevalence_ego_aud, "\n")

Prevalence of AUD among recently released agents:  0.1111111 

# Create a legend
legend <- rbind(
  data.frame(group = "Recently Released", shape = "square", color = "black"),
  data.frame(group = "Recently Released & Currently Smoking", shape = "square", color = "red"),
  data.frame(group = "Recently Released & AUD", shape = "square", color = "blue"),
  data.frame(group = "Currently Smoking", shape = "circle", color = "red"),
  data.frame(group = "AUD", shape = "circle", color = "blue")
)


# The ggplot call
p <- ggplot() +
  geom_point(data = legend, 
             aes(x = 1, y = group, color = color, fill = color, shape = shape), 
             size = 5) +
  geom_text(data = legend, 
            aes(x = 1.1, y = group, label = group), 
            hjust = 0) +
  scale_color_identity() +
  scale_fill_identity() +
  scale_shape_manual(values=c("square"=22, "circle"=21)) +
  theme_void() +
  theme(legend.position = "none") +
  coord_cartesian(xlim = c(1, 2))  # Adjust as necessary to ensure labels fit

# Print the plot
p

Alternate legend:

# Create a legend
legend <- data.frame(
  group = c("Recently Released", "Network Member"),
  shape = c("square", "circle"),
  fillcolor = c("black", "white"),
  color = c("black", "black")
)

# Map the attributes for smoking and alcohol use to the legend
smoking_legend <- data.frame(
  group = c("Current Smoking"),
  shape = c("circle"),
  fillcolor = c("blue"),
  color = c("blue")
)

alc_legend <- data.frame(
  group = c("Alcohol Use Disorder"),
  shape = c("circle"),
  fillcolor = c("red"),
  color = c("red")
)

# Combine the legend data frames
legend <- rbind(legend, smoking_legend, alc_legend)

# The ggplot call
legend_p_2 <- ggplot() +
  geom_point(data = legend, 
             aes(x = 1, y = group, fill = fillcolor, color = color, shape = shape), 
             size = 5) +
  geom_text(data = legend, 
            aes(x = 1.1, y = group, label = group), 
            hjust = 0) +
  scale_fill_identity() +
  scale_color_identity() +
  scale_shape_manual(values = c("square" = 22, "circle" = 21)) +
  theme_void() +
  theme(
    legend.position = "none") +
  coord_cartesian(xlim = c(1, 2))  # Adjust as necessary to ensure labels fit

# Print the plot
legend_p_2

legend_p_2

Horizontal legend alignment:

# Create a legend
legend <- data.frame(
  group = c("Recently Released", "Network Member", "Current Smoking", "Alcohol Use Disorder"),
  shape = c("square", "circle", "circle", "circle"),
  fillcolor = c("black", "white", "blue", "red"),
  color = c("black", "black", "blue", "red")
)

# The ggplot call
horizontal_legend_p <- ggplot() +
  geom_point(data = legend, 
             aes(x = group, y = 1, fill = fillcolor, color = color, shape = shape), 
             size = 5) +
  geom_text(data = legend, 
            aes(x = group, y = 1.1, label = group), 
            vjust = -0.5, angle = 45, hjust = 1) +
  scale_fill_identity() +
  scale_color_identity() +
  scale_shape_manual(values = c("square" = 22, "circle" = 21)) +
  theme_void() +
  theme(
    legend.position = "none"
  ) +
  coord_cartesian(ylim = c(0.8, 2))  # Adjust as necessary to ensure labels fit

# Print the plot
horizontal_legend_p

Plot AUD and Current Smking in the networks of a same number of agents who have never been incarcerated:

# Set seed for reproducibility
set.seed(Sys.time())

# Identify the IDs of the agents who have never been incarcerated
never_incarcerated_agents <- last_tick_agent_dt$id[last_tick_agent_dt$n_incarcerations == 0]

# Randomly select 36 agents who have never been incarcerated
selected_agents <- sample(never_incarcerated_agents, 
                          size = length(recently_released_agents))

# Get the network data for the selected agents
network_selected_agents <- last_tick_network_dt[(last_tick_network_dt$p1 %in% selected_agents) | 
                                                (last_tick_network_dt$p2 %in% selected_agents),]

# Get the first-degree network (neighbors) for each agent
first_degree_neighbors <- unique(c(network_selected_agents$p1, network_selected_agents$p2))

# Get agent data for the first-degree neighbors
first_degree_neighbors_agent_data <- last_tick_agent_dt[last_tick_agent_dt$id %in% first_degree_neighbors,]

# Create an edge data frame from your network data
edf <- data.frame(from = network_selected_agents$p1, 
                  to = network_selected_agents$p2)

# Create a node data frame with your agent data
ndf <- data.frame(id = first_degree_neighbors_agent_data$id, 
                  smoking_status = first_degree_neighbors_agent_data$smoking_status,
                  alc_use_status = first_degree_neighbors_agent_data$alc_use_status)

# Create a graph with DiagrammeR
graph <- create_graph(nodes_df = ndf, 
                      edges_df = edf,
                      directed = FALSE)

# Customize the node fillcolor based on smoking_status and alc_use_status
graph <- set_node_attrs(graph, "fillcolor", ifelse(ndf$smoking_status == "Current", "red", 
                                                   ifelse(ndf$alc_use_status == 3, "blue", "white")))

# Customize the node color based on smoking_status and alc_use_status
graph <- set_node_attrs(graph, "color", ifelse(ndf$smoking_status == "Current", "red", 
                                               ifelse(ndf$alc_use_status == 3, "blue", "gray")))

# Remove labels from nodes to make them stand out more
graph <- set_node_attrs(graph, "label", "")

# Render the graph
render_graph(graph)

NA

Numerical prevalence of AUD and Current Smoking in this network:

# Number of agents with Current Smoking status
num_current_smoking <- sum(first_degree_neighbors_agent_data$smoking_status == "Current")

# Number of agents with AUD
num_aud <- sum(first_degree_neighbors_agent_data$alc_use_status == 3)

# Total number of agents in the network
total_agents <- nrow(first_degree_neighbors_agent_data)

# Calculate prevalence
prevalence_smoking <- num_current_smoking / total_agents
prevalence_aud <- num_aud / total_agents

# Print the results
cat("Prevalence of Current Smoking in the networks of randomly 
    selected never incarcerated persons: ", prevalence_smoking, "\n")

Prevalence of Current Smoking in the networks of randomly 
    selected never incarcerated persons:  0.1414141 

cat("Prevalence of AUD in the networks of randomly 
    selected never incarcerated persons: ", prevalence_aud, "\n")

Prevalence of AUD in the networks of randomly 
    selected never incarcerated persons:  0.07575758 

Smoking and AUD among the egos:

# Define selected_agents_df, which contains the data for the selected agents
selected_agents_df <- last_tick_agent_dt[last_tick_agent_dt$id %in% selected_agents,]

# Number of selected agents who are current smokers
num_ego_current_smoking <- sum(selected_agents_df$smoking_status == "Current")

# Number of selected agents with AUD
num_ego_aud <- sum(selected_agents_df$alc_use_status == 3)

# Total number of selected agents
total_ego_agents <- length(selected_agents)

# Calculate prevalence
prevalence_ego_smoking <- num_ego_current_smoking / total_ego_agents
prevalence_ego_aud <- num_ego_aud / total_ego_agents

# Print the results
cat("Prevalence of Current Smoking among selected agents: ", prevalence_ego_smoking, "\n")

Prevalence of Current Smoking among selected agents:  0.1111111 

cat("Prevalence of AUD among selected agents: ", prevalence_ego_aud, "\n")

Prevalence of AUD among selected agents:  0.08333333