Processing and Visualizing Simulation Results-Kenya CVD Microsimulation Model
James Oguta
December 1, 2025
- 0.1 Code Description
- 1 Clear the workspace
- 2 Processing results for patients with no intervention
- 3 Plot the state occupancy over cycles
- 4 Plot the distribution of state occupancy across the six states for each cycle
- 5 Plot the distribution of CVD events-To be fixed so as to define the maximum number of events
- 6 Plot the distribution of MI events
- 7 Count only the first instance of MI and stroke events per individual
- 8 Summarize people with each event and both events
- 9 Plot the costs over cycles
- 10 Renaming the dataframes for patients on Empower Health intervention
- 11 Plot the state occupancy over cycles for patients on Empower Health intervention
- 12 Plot the distribution of CVD events for patients on Empower Health intervention
- 13 Plot the distribution of MI events for patients on Empower Health intervention
- 14 Count only the first instance of MI and stroke events per individual for patients on Empower Health intervention
- 15 Summarize people with each event and both events for patients on Empower Health intervention
- 16 Plot the costs over cycles for patients on Empower Health intervention
- 17 Renaming the dataframes for patients on Usual Care intervention
- 18 Plot the state occupancy over cycles for patients on Usual Care intervention
- 19 Plot the distribution of CVD events for patients on Usual Care intervention
- 20 Count the number of CVD events(MI and Stroke) for each individual-separate stroke and MI events for patients on Usual Care intervention
- 21 Plot the distribution of MI events for patients on Usual Care intervention
- 22 Count only the first instance of MI and stroke events per individual for patients on Usual Care intervention
- 23 Summarize people with each event and both events for patients on Usual Care intervention
- 24 Plot the costs over cycles for patients on Usual Care intervention
- 25 Comparing the results between No intervention, Empower Health and Usual Care interventions
- 26 Plot the distribution of state occupancy across the six states for each cycle for none and both interventions-label proportions in each state for the last cycle
- 27 Summarise the proportion state occupancy at the final cycle for all interventions
- 28 Plot stacked area chart for all interventions
- 29 Plot the distribution of CVD events for all interventions
- 30 Draw graph of MI, stroke and deaths averted for Empower Health and Usual Care interventions compared to No intervention
- 31 Plot the costs over cycles for all interventions
- 32 Plot mean costs per cycle for All interventions
- 33 Plot mean DALYs per cycle for both interventions
- 34 Plot a cost effectiveness plane for the three interventions
- 35 Plot the cost-effectiveness plane
- 36 Plot the cost-effectiveness plane comparing Empower Health and Usual Care interventions using DALYs averted and costs-add willingness to pay (WTP) threshold line at $1000
- 37 Compute the incremental cost-effectiveness ratio (ICER) for the two interventions
- 38 Alternatively
- 39 Rename datasets
- 40 Add strategy to the intervention datasets
- 41 Combine into one dataset for plotting
- 42 Reshape data to wide format
- 43 Step 3: Calculate Incremental Cost and Effect
- 44 Cost-Effectiveness Plane (CE Plane) Plot
- 45 Add no intervention options
- 46 CE Plane for all the three interventions
- 47 Cost-Effectiveness Acceptability Curve (All Interventions)
0.1 Code Description
The following R script processes the results of the simulation. It takes in the simulation results(m_States, m_Costs and m_Effs) matrices and summarizes them into a data frame. The script also generates a plot of the simulation results.
1 Clear the workspace
1.1 Load Libraries
The script begins by loading the necessary R libraries:
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 4.0.0 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.4
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors##
## Attaching package: 'reshape2'
##
## The following object is masked from 'package:tidyr':
##
## smiths2 Processing results for patients with no intervention
2.1 Load Data
The script loads the simulation results from the rda file
# Load the simulation results from rda file
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_States.rda') # No intervention
# load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_Costs.rda') # No intervention
# load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_Effs.rda') # No intervention
# # load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_States_UC.rda')2.2 Renaming the dataframes
v_states_names <- c("NE", "ST", "PS", "MI", "PM", "D")
state_occupancy <- as.data.frame(m_States)
# Add individual id variable taken from rownames without the prefix "ind_"
state_occupancy$id <- gsub("ind_", "", rownames(state_occupancy))
# Reshape the data to long format
state_occupancy_long <- melt(state_occupancy, id.vars = "id", variable.name = "cycle", value.name = "state")
# Convert cycle to numeric
state_occupancy_long$cycle <- as.numeric(gsub("cycle_", "", state_occupancy_long$cycle))
# Convert state to factor
state_occupancy_long$state <- factor(state_occupancy_long$state, levels = v_states_names)
head(state_occupancy_long)## id cycle state
## 1 1 0 NE
## 2 2 0 NE
## 3 3 0 NE
## 4 4 0 NE
## 5 5 0 NE
## 6 6 0 NE3 Plot the state occupancy over cycles
ggplot(state_occupancy_long, aes(x = cycle, y = id, fill = state)) +
geom_tile() +
scale_fill_manual(values = c("NE" = "green", "ST" = "yellow", "PS" = "orange", "MI" = "red", "PM" = "purple", "D" = "black")) +
labs(title = "Health State Occupancy Over Cycles",
x = "Cycle",
y = "Individual ID",
fill = "Health State") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Plot distribution of state occupancy across the six states for each
cycle # ———- 2. Compute proportions by state per cycle ———- # Compute
the proportions of individuals in each state for each cycle
state_proportions <- state_occupancy_long %>%
group_by(cycle, state) %>%
summarise(proportion = n() / nrow(state_occupancy_long)) %>%
ungroup()## `summarise()` has grouped output by 'cycle'. You can override using the
## `.groups` argument.state_proportions <- state_occupancy_long %>%
group_by(cycle, state) %>%
summarise(count = n(), .groups = "drop") %>%
group_by(cycle) %>%
mutate(proportion = count / sum(count))
head(state_proportions)## # A tibble: 6 × 4
## # Groups: cycle [3]
## cycle state count proportion
## <dbl> <fct> <int> <dbl>
## 1 0 NE 100000 1
## 2 1 NE 98981 0.990
## 3 1 ST 176 0.00176
## 4 1 MI 159 0.00159
## 5 1 D 684 0.00684
## 6 2 NE 97888 0.9794 Plot the distribution of state occupancy across the six states for each cycle
ggplot(state_proportions, aes(x = cycle, y = proportion, fill = state)) +
geom_bar(stat = "identity", position = "fill") +
scale_fill_manual(values = c("NE" = "green", "ST" = "yellow", "PS" = "orange", "MI" = "red", "PM" = "purple", "D" = "black")) +
labs(title = "Distribution of Health State Occupancy per Cycle",
x = "Cycle",
y = "Proportion",
fill = "Health State") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# ———- Plot stacked area chart ———-
# Plot stacked area chart
ggplot(state_proportions, aes(x = cycle, y = proportion, fill = state)) +
geom_area() +
scale_fill_manual(values = c("NE" = "green", "ST" = "yellow", "PS" = "orange", "MI" = "red", "PM" = "purple", "D" = "black")) +
labs(title = "Health State Distribution Over Time",
x = "Cycle",
y = "Proportion",
fill = "Health State") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Counting number of CVD events
# Count the number of CVD events(MI and Stroke) for each individual
cvd_events <- state_occupancy_long %>%
filter(state %in% c("MI", "ST")) %>%
group_by(id) %>%
summarise(cvd_events = n()) %>%
ungroup()
# Print the first few rows of the cvd_events data frame
head(cvd_events)## # A tibble: 6 × 2
## id cvd_events
## <chr> <int>
## 1 10 1
## 2 100 1
## 3 1000 1
## 4 10000 2
## 5 100000 2
## 6 10001 15 Plot the distribution of CVD events-To be fixed so as to define the maximum number of events
ggplot(cvd_events, aes(x = cvd_events)) +
geom_histogram(binwidth = 1, fill = "blue", color = "black", alpha = 0.7) +
labs(title = "Distribution of CVD Events",
x = "Number of CVD Events",
y = "Frequency") +
theme_minimal()
# Count the number of CVD events(MI and Stroke) for each
individual-separate stroke and MI events
# Count the number of MI events for each individual
mi_events <- state_occupancy_long %>%
filter(state == "MI") %>%
group_by(id) %>%
summarise(mi_events = n()) %>%
ungroup()
# Count the number of Stroke events for each individual
stroke_events <- state_occupancy_long %>%
filter(state == "ST") %>%
group_by(id) %>%
summarise(stroke_events = n()) %>%
ungroup()
# Count the number of deaths
deaths <- state_occupancy_long %>%
filter(state == "D") %>%
group_by(id) %>%
summarise(deaths = n()) %>%
ungroup()
# Print the first few rows of the mi_events and stroke_events data frames
head(mi_events)## # A tibble: 6 × 2
## id mi_events
## <chr> <int>
## 1 10 1
## 2 1000 1
## 3 10003 1
## 4 10004 1
## 5 10006 1
## 6 10008 3## # A tibble: 6 × 2
## id stroke_events
## <chr> <int>
## 1 100 1
## 2 10000 2
## 3 100000 2
## 4 10001 1
## 5 10003 1
## 6 10007 1## # A tibble: 6 × 2
## id deaths
## <chr> <int>
## 1 1 47
## 2 100 1
## 3 1000 18
## 4 10000 35
## 5 100000 9
## 6 10001 76 Plot the distribution of MI events
ggplot(mi_events, aes(x = mi_events)) +
geom_histogram(binwidth = 1, fill = "red", color = "black", alpha = 0.7) +
labs(title = "Distribution of MI Events_No Intervention",
x = "Number of MI Events",
y = "Frequency") +
theme_minimal()
# Plot the distribution of Stroke events
ggplot(stroke_events, aes(x = stroke_events)) +
geom_histogram(binwidth = 1, fill = "green", color = "black", alpha = 0.7) +
labs(title = "Distribution of Stroke Events_No Intervention",
x = "Number of Stroke Events",
y = "Frequency") +
theme_minimal()
# Plot the distribution of MI and Stroke events-No intervention
# Combine the MI and Stroke events data frames
cvd_events_combined <- merge(mi_events, stroke_events, by = "id", all = TRUE)
# Replace NA values with 0
cvd_events_combined[is.na(cvd_events_combined)] <- 0
# Reshape the data to long format
cvd_events_long <- melt(cvd_events_combined, id.vars = "id", variable.name = "event_type", value.name = "event_count")
# Convert event_type to factor
cvd_events_long$event_type <- factor(cvd_events_long$event_type, levels = c("mi_events", "stroke_events"), labels = c("MI Events", "Stroke Events"))
# Plot the distribution of MI and Stroke events
ggplot(cvd_events_long, aes(x = event_count, fill = event_type)) +
geom_histogram(binwidth = 1, position = "identity", alpha = 0.7) +
labs(title = "Distribution of MI and Stroke Events (No Intervention)",
x = "Number of Events",
y = "Frequency") +
scale_fill_manual(values = c("MI Events" = "red", "Stroke Events" = "green")) +
theme_minimal()
# Plot the number of CVD events per person over their lifetime
# Combine the MI and Stroke events data frames
cvd_events_combined <- merge(mi_events, stroke_events, by = "id", all = TRUE)
# Replace NA values with 0
cvd_events_combined[is.na(cvd_events_combined)] <- 0
# Reshape the data to long format
cvd_events_long <- melt(cvd_events_combined, id.vars = "id", variable.name = "event_type", value.name = "event_count")
# Convert event_type to factor
cvd_events_long$event_type <- factor(cvd_events_long$event_type, levels = c("mi_events", "stroke_events"), labels = c("MI Events", "Stroke Events"))
# Summarise the data to get the total number of events per person-separate MI and Stroke events
cvd_events_summary <- cvd_events_long %>%
group_by(id) %>%
summarise(total_mi_events = sum(event_count[event_type == "MI Events"]),
total_stroke_events = sum(event_count[event_type == "Stroke Events"])) %>%
ungroup()
# Plot the number of CVD events per person over their lifetime
ggplot(cvd_events_summary, aes(x = id)) +
geom_bar(aes(y = total_mi_events), stat = "identity", fill = "red", alpha = 0.7) +
geom_bar(aes(y = total_stroke_events), stat = "identity", fill = "green", alpha = 0.7) +
labs(title = "Number of CVD Events per Person Over Their Lifetime",
x = "Individual ID",
y = "Number of Events") +
scale_y_continuous(sec.axis = sec_axis(~., name = "Number of Events")) +
theme_minimal()
hist(cvd_events_summary$total_mi_events, breaks = 20, main = "Distribution of MI Events", xlab = "Number of MI Events", col = "red")
# Count the number of deaths
# Count the total number of deaths for all individuals- Only take the first instance of death per individual
deaths <- state_occupancy_long %>%
filter(state == "D") %>%
distinct(id) %>%
mutate(deaths = 1)
# Compute the total number of deaths across all individuals
total_deaths <- nrow(deaths)
# Print the total number of deaths
cat("Total number of deaths:", total_deaths, "\n")## Total number of deaths: 752627 Count only the first instance of MI and stroke events per individual
# Count only the first instance of MI events per individual
mi_first_instance <- state_occupancy_long %>%
filter(state == "MI") %>%
group_by(id) %>%
summarise(mi_first_instance = n() > 0) %>%
ungroup()
# Count only the first instance of Stroke events per individual
stroke_first_instance <- state_occupancy_long %>%
filter(state == "ST") %>%
group_by(id) %>%
summarise(stroke_first_instance = n() > 0) %>%
ungroup()
# View the first few rows of the mi_first_instance and stroke_first_instance data frames
head(mi_first_instance)## # A tibble: 6 × 2
## id mi_first_instance
## <chr> <lgl>
## 1 10 TRUE
## 2 1000 TRUE
## 3 10003 TRUE
## 4 10004 TRUE
## 5 10006 TRUE
## 6 10008 TRUE## # A tibble: 6 × 2
## id stroke_first_instance
## <chr> <lgl>
## 1 100 TRUE
## 2 10000 TRUE
## 3 100000 TRUE
## 4 10001 TRUE
## 5 10003 TRUE
## 6 10007 TRUE# Count the number of individuals with MI events
mi_first_count <- nrow(mi_first_instance[mi_first_instance$mi_first_instance == TRUE, ])
# Count the number of individuals with Stroke events
stroke_first_count <- nrow(stroke_first_instance[stroke_first_instance$stroke_first_instance == TRUE, ])
# Print the counts
cat("Number of individuals with first instance of MI events:", mi_first_count, "\n")## Number of individuals with first instance of MI events: 37617## Number of individuals with first instance of Stroke events: 437048 Summarize people with each event and both events
# Count the number of individuals with MI events
mi_count <- nrow(cvd_events[cvd_events$cvd_events > 0, ])
# Count the number of individuals with Stroke events
stroke_count <- nrow(stroke_events[stroke_events$stroke_events > 0, ])
# Count the number of individuals with both MI and Stroke events
both_count <- nrow(cvd_events_combined[cvd_events_combined$mi_events > 0 & cvd_events_combined$stroke_events > 0, ])
# Print the counts
cat("Number of individuals with MI events:", mi_count, "\n")## Number of individuals with MI events: 73109## Number of individuals with Stroke events: 43704## Number of individuals with both MI and Stroke events: 8212## Number of individuals with first instance of MI events: 37617## Number of individuals with first instance of Stroke events: 43704## Total number of deaths: 752629 Plot the costs over cycles
# Clear the workspace from previous data frames
rm(m_States)
# Load the simulation results from rda file
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_Costs.rda') # No intervention
costs_df <- as.data.frame(m_Costs)
# Add individual id variable taken from rownames without the prefix "ind_"
costs_df$id <- gsub("ind_", "", rownames(costs_df))
costs_long <- melt(costs_df, id.vars = "id", variable.name = "cycle", value.name = "cost")
# Convert cycle to numeric
costs_long$cycle <- as.numeric(gsub("cycle_", "", costs_long$cycle))
# Plot the costs over cycles
ggplot(costs_long, aes(x = cycle, y = id, fill = cost)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "black") +
labs(title = "Costs Over Cycles",
x = "Cycle",
y = "Individual ID",
fill = "Cost") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Plot mean costs per cycle
# Calculate mean costs per cycle
mean_costs <- costs_long %>%
group_by(cycle) %>%
summarise(mean_cost = mean(cost, na.rm = TRUE)) %>%
ungroup()
# Plot mean costs per cycle
ggplot(mean_costs, aes(x = cycle, y = mean_cost)) +
geom_line() +
labs(title = "Mean Costs per Cycle",
x = "Cycle",
y = "Mean Cost") +
theme_minimal()
# Plot the DALYs over cycles
## Warning in rm(m_States, m_Costs): object 'm_States' not found# Load the simulation results from rda file
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_Effs.rda') # No intervention
effs_df <- as.data.frame(m_Effs)
# Add individual id variable taken from rownames without the prefix "ind_"
effs_df$id <- gsub("ind_", "", rownames(effs_df))
effs_long <- melt(effs_df, id.vars = "id", variable.name = "cycle", value.name = "eff")
# Convert cycle to numeric
effs_long$cycle <- as.numeric(gsub("cycle_", "", effs_long$cycle))
# Plot the DALYs over cycles
ggplot(effs_long, aes(x = cycle, y = id, fill = eff)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "black") +
labs(title = "DALYs Over Cycles",
x = "Cycle",
y = "Individual ID",
fill = "DALY") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Plot mean DALYs per cycle
# Calculate mean DALYs per cycle
mean_effs <- effs_long %>%
group_by(cycle) %>%
summarise(mean_eff = mean(eff, na.rm = TRUE)) %>%
ungroup()
# Plot mean DALYs per cycle
ggplot(mean_effs, aes(x = cycle, y = mean_eff)) +
geom_line() +
labs(title = "Mean DALYs per Cycle",
x = "Cycle",
y = "Mean DALY") +
theme_minimal()
## Repeat process for patients on Empower Health intervention
# Clear the workspace from previous data frames for no intervention
rm(m_Effs)
# Load the simulation results from rda file for patients on Empower Health intervention
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_States_Emp.rda') # On Empower Health intervention
# load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_Costs_Emp.rda') # On Empower Health intervention
# load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_Effs_Emp.rda') # On Empower Health intervention
# load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_MI_Risks_Emp.rda') # On Empower Health intervention10 Renaming the dataframes for patients on Empower Health intervention
v_states_names <- c("NE", "ST", "PS", "MI", "PM", "D")
state_occupancy_emp <- as.data.frame(m_States_Emp)
# Add individual id variable taken from rownames without the prefix "ind_"
state_occupancy_emp$id <- gsub("ind_", "", rownames(state_occupancy_emp))
# Reshape the data to long format
state_occupancy_emp_long <- melt(state_occupancy_emp, id.vars = "id", variable.name = "cycle", value.name = "state")
# Convert cycle to numeric
state_occupancy_emp_long$cycle <- as.numeric(gsub("cycle_", "", state_occupancy_emp_long$cycle))
# Convert state to factor
state_occupancy_emp_long$state <- factor(state_occupancy_emp_long$state, levels = v_states_names)
head(state_occupancy_emp_long)## id cycle state
## 1 1 0 NE
## 2 2 0 NE
## 3 3 0 NE
## 4 4 0 NE
## 5 5 0 NE
## 6 6 0 NE11 Plot the state occupancy over cycles for patients on Empower Health intervention
ggplot(state_occupancy_emp_long, aes(x = cycle, y = id, fill = state)) +
geom_tile() +
scale_fill_manual(values = c("NE" = "green", "ST" = "yellow", "PS" = "orange", "MI" = "red", "PM" = "purple", "D" = "black")) +
labs(title = "Health State Occupancy Over Cycles (Empower Health Intervention)",
x = "Cycle",
y = "Individual ID",
fill = "Health State") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Plot distribution of state occupancy across the six states for each
cycle for patients on Empower Health intervention
# Compute the proportions of individuals in each state for each cycle for patients on Empower Health intervention
state_proportions_emp <- state_occupancy_emp_long %>%
group_by(cycle, state) %>%
summarise(proportion = n() / nrow(state_occupancy_emp_long)) %>%
ungroup()## `summarise()` has grouped output by 'cycle'. You can override using the
## `.groups` argument.state_proportions_emp <- state_occupancy_emp_long %>%
group_by(cycle, state) %>%
summarise(count = n(), .groups = "drop") %>%
group_by(cycle) %>%
mutate(proportion = count / sum(count))
head(state_proportions_emp)## # A tibble: 6 × 4
## # Groups: cycle [3]
## cycle state count proportion
## <dbl> <fct> <int> <dbl>
## 1 0 NE 100000 1
## 2 1 NE 98902 0.989
## 3 1 ST 225 0.00225
## 4 1 MI 173 0.00173
## 5 1 D 700 0.007
## 6 2 NE 97797 0.978# Plot the distribution of state occupancy across the six states for each cycle for patients on Empower Health intervention
ggplot(state_proportions_emp, aes(x = cycle, y = proportion, fill = state)) +
geom_bar(stat = "identity", position = "fill") +
scale_fill_manual(values = c("NE" = "green", "ST" = "yellow", "PS" = "orange", "MI" = "red", "PM" = "purple", "D" = "black")) +
labs(title = "Distribution of Health State Occupancy per Cycle (Empower Health Intervention)",
x = "Cycle",
y = "Proportion",
fill = "Health State") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Plot stacked area chart for patients on Empower Health
intervention
# Plot stacked area chart for patients on Empower Health intervention
ggplot(state_proportions_emp, aes(x = cycle, y = proportion, fill = state)) +
geom_area() +
scale_fill_manual(values = c("NE" = "green", "ST" = "yellow", "PS" = "orange", "MI" = "red", "PM" = "purple", "D" = "black")) +
labs(title = "Health State Distribution Over Time (Empower Health Intervention)",
x = "Cycle",
y = "Proportion",
fill = "Health State") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Counting number of CVD events for patients on Empower Health
intervention
# Count the number of CVD events(MI and Stroke) for each individual for patients on Empower Health intervention
cvd_events_emp <- state_occupancy_emp_long %>%
filter(state %in% c("MI", "ST")) %>%
group_by(id) %>%
summarise(cvd_events = n()) %>%
ungroup()
# Print the first few rows of the cvd_events_emp data frame
head(cvd_events_emp)## # A tibble: 6 × 2
## id cvd_events
## <chr> <int>
## 1 10 1
## 2 100 1
## 3 10000 1
## 4 100000 2
## 5 10002 2
## 6 10003 212 Plot the distribution of CVD events for patients on Empower Health intervention
ggplot(cvd_events_emp, aes(x = cvd_events)) +
geom_histogram(binwidth = 1, fill = "blue", color = "black", alpha = 0.7) +
labs(title = "Distribution of CVD Events (Empower Health Intervention)",
x = "Number of CVD Events",
y = "Frequency") +
theme_minimal()
# Count the number of CVD events(MI and Stroke) for each
individual-separate stroke and MI events for patients on Empower Health
intervention
# Count the number of MI events for each individual for patients on Empower Health intervention
mi_events_emp <- state_occupancy_emp_long %>%
filter(state == "MI") %>%
group_by(id) %>%
summarise(mi_events = n()) %>%
ungroup()
# Count the number of Stroke events for each individual for patients on Empower Health intervention
stroke_events_emp <- state_occupancy_emp_long %>%
filter(state == "ST") %>%
group_by(id) %>%
summarise(stroke_events = n()) %>%
ungroup()
# Print the first few rows of the mi_events_emp and stroke_events_emp data frames
head(mi_events_emp)## # A tibble: 6 × 2
## id mi_events
## <chr> <int>
## 1 10 1
## 2 100000 2
## 3 10002 2
## 4 10004 1
## 5 1001 1
## 6 10010 2## # A tibble: 6 × 2
## id stroke_events
## <chr> <int>
## 1 100 1
## 2 10000 1
## 3 10003 2
## 4 10011 1
## 5 10015 1
## 6 10018 113 Plot the distribution of MI events for patients on Empower Health intervention
ggplot(mi_events_emp, aes(x = mi_events)) +
geom_histogram(binwidth = 1, fill = "red", color = "black", alpha = 0.7) +
labs(title = "Distribution of MI Events (Empower Health Intervention)",
x = "Number of MI Events",
y = "Frequency") +
theme_minimal()
# Plot the distribution of Stroke events for patients on Empower Health
intervention
ggplot(stroke_events_emp, aes(x = stroke_events)) +
geom_histogram(binwidth = 1, fill = "green", color = "black", alpha = 0.7) +
labs(title = "Distribution of Stroke Events (Empower Health Intervention)",
x = "Number of Stroke Events",
y = "Frequency") +
theme_minimal()
# Plot the distribution of MI and Stroke events for patients on Empower
Health intervention
# Combine the MI and Stroke events data frames for patients on Empower Health intervention
cvd_events_combined_emp <- merge(mi_events_emp, stroke_events_emp, by = "id", all = TRUE)
# Replace NA values with 0
cvd_events_combined_emp[is.na(cvd_events_combined_emp)] <- 0
# Reshape the data to long format
cvd_events_long_emp <- melt(cvd_events_combined_emp, id.vars = "id", variable.name = "event_type", value.name = "event_count")
# Convert event_type to factor
cvd_events_long_emp$event_type <- factor(cvd_events_long_emp$event_type, levels = c("mi_events", "stroke_events"), labels = c("MI Events", "Stroke Events"))
# Plot the distribution of MI and Stroke events for patients on Empower Health intervention
ggplot(cvd_events_long_emp, aes(x = event_count, fill = event_type)) +
geom_histogram(binwidth = 1, position = "identity", alpha = 0.7) +
labs(title = "Distribution of MI and Stroke Events (Empower Health Intervention)",
x = "Number of Events",
y = "Frequency") +
scale_fill_manual(values = c("MI Events" = "red", "Stroke Events" = "green")) +
theme_minimal()
# Plot the number of CVD events per person over their lifetime for
patients on Empower Health intervention
# Combine the MI and Stroke events data frames for patients on Empower Health intervention
cvd_events_combined_emp <- merge(mi_events_emp, stroke_events_emp, by = "id", all = TRUE)
# Replace NA values with 0
cvd_events_combined_emp[is.na(cvd_events_combined_emp)] <- 0
# Reshape the data to long format
cvd_events_long_emp <- melt(cvd_events_combined_emp, id.vars = "id", variable.name = "event_type", value.name = "event_count")
# Convert event_type to factor
cvd_events_long_emp$event_type <- factor(cvd_events_long_emp$event_type, levels = c("mi_events", "stroke_events"), labels = c("MI Events", "Stroke Events"))
# Summarise the data to get the total number of events per person for patients on Empower Health intervention
cvd_events_summary_emp <- cvd_events_long_emp %>%
group_by(id) %>%
summarise(total_mi_events = sum(event_count[event_type == "MI Events"]),
total_stroke_events = sum(event_count[event_type == "Stroke Events"])) %>%
ungroup()
# Plot the number of CVD events per person over their lifetime for patients on Empower Health intervention
ggplot(cvd_events_summary_emp, aes(x = id)) +
geom_bar(aes(y = total_mi_events), stat = "identity", fill = "red", alpha = 0.7) +
geom_bar(aes(y = total_stroke_events), stat = "identity", fill = "green", alpha = 0.7) +
labs(title = "Number of CVD Events per Person Over Their Lifetime (Empower Health Intervention)",
x = "Individual ID",
y = "Number of Events") +
scale_y_continuous(sec.axis = sec_axis(~., name = "Number of Events")) +
theme_minimal()
hist(cvd_events_summary_emp$total_mi_events, breaks = 20, main = "Distribution of MI Events (Empower Health Intervention)", xlab = "Number of MI Events", col = "red")
# Count the number of deaths for patients on Empower Health
intervention
# Count the total number of deaths for all individuals for patients on Empower Health intervention- Only take the first instance of death per individual
deaths_emp <- state_occupancy_emp_long %>%
filter(state == "D") %>%
distinct(id) %>%
mutate(deaths = 1)
# Compute the total number of deaths across all individuals for patients on Empower Health intervention
total_deaths_emp <- nrow(deaths_emp)
# Print the total number of deaths for patients on Empower Health intervention
cat("Total number of deaths (Empower Health Intervention):", total_deaths_emp, "\n")## Total number of deaths (Empower Health Intervention): 6569814 Count only the first instance of MI and stroke events per individual for patients on Empower Health intervention
# Count only the first instance of MI events per individual for patients on Empower Health intervention
mi_first_instance_emp <- state_occupancy_emp_long %>%
filter(state == "MI") %>%
group_by(id) %>%
summarise(mi_first_instance = n() > 0) %>%
ungroup()
# Count only the first instance of Stroke events per individual for patients on Empower Health intervention
stroke_first_instance_emp <- state_occupancy_emp_long %>%
filter(state == "ST") %>%
group_by(id) %>%
summarise(stroke_first_instance = n() > 0) %>%
ungroup()
# View the first few rows of the mi_first_instance_emp and stroke_first_instance_emp data frames
head(mi_first_instance_emp)## # A tibble: 6 × 2
## id mi_first_instance
## <chr> <lgl>
## 1 10 TRUE
## 2 100000 TRUE
## 3 10002 TRUE
## 4 10004 TRUE
## 5 1001 TRUE
## 6 10010 TRUE## # A tibble: 6 × 2
## id stroke_first_instance
## <chr> <lgl>
## 1 100 TRUE
## 2 10000 TRUE
## 3 10003 TRUE
## 4 10011 TRUE
## 5 10015 TRUE
## 6 10018 TRUE# Count the number of individuals with MI events for patients on Empower Health intervention
mi_first_count_emp <- nrow(mi_first_instance_emp[mi_first_instance_emp$mi_first_instance == TRUE, ])
# Count the number of individuals with Stroke events for patients on Empower Health intervention
stroke_first_count_emp <- nrow(stroke_first_instance_emp[stroke_first_instance_emp$stroke_first_instance == TRUE, ])
# Print the counts for patients on Empower Health intervention
cat("Number of individuals with first instance of MI events (Empower Health Intervention):", mi_first_count_emp, "\n")## Number of individuals with first instance of MI events (Empower Health Intervention): 29268cat("Number of individuals with first instance of Stroke events (Empower Health Intervention):", stroke_first_count_emp, "\n")## Number of individuals with first instance of Stroke events (Empower Health Intervention): 3563815 Summarize people with each event and both events for patients on Empower Health intervention
# Count the number of individuals with MI events for patients on Empower Health intervention
mi_count_emp <- nrow(cvd_events_emp[cvd_events_emp$cvd_events > 0, ])
# Count the number of individuals with Stroke events for patients on Empower Health intervention
stroke_count_emp <- nrow(stroke_events_emp[stroke_events_emp$stroke_events > 0, ])
# Count the number of individuals with both MI and Stroke events for patients on Empower Health intervention
both_count_emp <- nrow(cvd_events_combined_emp[cvd_events_combined_emp$mi_events > 0 & cvd_events_combined_emp$stroke_events > 0, ])
# Print the counts for patients on Empower Health intervention
cat("Number of individuals with MI events (Empower Health Intervention):", mi_count_emp, "\n")## Number of individuals with MI events (Empower Health Intervention): 58774cat("Number of individuals with Stroke events (Empower Health Intervention):", stroke_count_emp, "\n")## Number of individuals with Stroke events (Empower Health Intervention): 35638cat("Number of individuals with both MI and Stroke events (Empower Health Intervention):", both_count_emp, "\n")## Number of individuals with both MI and Stroke events (Empower Health Intervention): 6132cat("Number of individuals with first instance of MI events (Empower Health Intervention):", mi_first_count_emp, "\n")## Number of individuals with first instance of MI events (Empower Health Intervention): 29268cat("Number of individuals with first instance of Stroke events (Empower Health Intervention):", stroke_first_count_emp, "\n")## Number of individuals with first instance of Stroke events (Empower Health Intervention): 35638## Total number of deaths (Empower Health Intervention): 6569816 Plot the costs over cycles for patients on Empower Health intervention
# Clear the workspace from previous data frames for patients on Empower Health intervention
rm(m_States_Emp)
# Load the simulation results from rda file for patients on Empower Health intervention
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_Costs_Emp.rda') # On Empower Health intervention
costs_emp_df <- as.data.frame(m_Costs_Emp)
# Add individual id variable taken from rownames without the prefix "ind_"
costs_emp_df$id <- gsub("ind_", "", rownames(costs_emp_df))
costs_emp_long <- melt(costs_emp_df, id.vars = "id", variable.name = "cycle", value.name = "cost")
# Convert cycle to numeric
costs_emp_long$cycle <- as.numeric(gsub("cycle_", "", costs_emp_long$cycle))
# Plot the costs over cycles for patients on Empower Health intervention
ggplot(costs_emp_long, aes(x = cycle, y = id, fill = cost)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "black") +
labs(title = "Costs Over Cycles (Empower Health Intervention)",
x = "Cycle",
y = "Individual ID",
fill = "Cost") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Plot mean costs per cycle for patients on Empower Health
intervention
# Calculate mean costs per cycle for patients on Empower Health intervention
mean_costs_emp <- costs_emp_long %>%
group_by(cycle) %>%
summarise(mean_cost = mean(cost, na.rm = TRUE)) %>%
ungroup()
# Plot mean costs per cycle for patients on Empower Health intervention
ggplot(mean_costs_emp, aes(x = cycle, y = mean_cost)) +
geom_line() +
labs(title = "Mean Costs per Cycle (Empower Health Intervention)",
x = "Cycle",
y = "Mean Cost") +
theme_minimal()
# Plot the DALYs over cycles for patients on Empower Health
intervention
# Clear the workspace from previous data frames for patients on Empower Health intervention
rm(m_Costs_Emp)
# Load the simulation results from rda file for patients on Empower Health intervention
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_Effs_Emp.rda') # On Empower Health intervention
effs_emp_df <- as.data.frame(m_Effs_Emp)
# Add individual id variable taken from rownames without the prefix "ind_"
effs_emp_df$id <- gsub("ind_", "", rownames(effs_emp_df))
effs_emp_long <- melt(effs_emp_df, id.vars = "id", variable.name = "cycle", value.name = "eff")
# Convert cycle to numeric
effs_emp_long$cycle <- as.numeric(gsub("cycle_", "", effs_emp_long$cycle))
# Plot the DALYs over cycles for patients on Empower Health intervention
ggplot(effs_emp_long, aes(x = cycle, y = id, fill = eff)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "black") +
labs(title = "DALYs Over Cycles (Empower Health Intervention)",
x = "Cycle",
y = "Individual ID",
fill = "DALY") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Plot mean DALYs per cycle for patients on Empower Health
intervention
# Calculate mean DALYs per cycle for patients on Empower Health intervention
mean_effs_emp <- effs_emp_long %>%
group_by(cycle) %>%
summarise(mean_eff = mean(eff, na.rm = TRUE)) %>%
ungroup()
# Plot mean DALYs per cycle for patients on Empower Health intervention
ggplot(mean_effs_emp, aes(x = cycle, y = mean_eff)) +
geom_line() +
labs(title = "Mean DALYs per Cycle (Empower Health Intervention)",
x = "Cycle",
y = "Mean DALY") +
theme_minimal()
## Repeat process for patients on Usual Care intervention
# Clear the workspace from previous data frames for patients on Empower Health intervention
rm(m_Effs_Emp)
# Load the simulation results from rda file for patients on Usual Care intervention
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_States_UC.rda') # On Usual Care intervention
# load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_Costs_UC.rda') # On Usual Care intervention
# load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_Effs_UC.rda') # On Usual Care intervention
# # load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_MI_Risks_UC.rda') # On Usual Care intervention17 Renaming the dataframes for patients on Usual Care intervention
v_states_names <- c("NE", "ST", "PS", "MI", "PM", "D")
state_occupancy_uc <- as.data.frame(m_States_UC)
# Add individual id variable taken from rownames without the prefix "ind_"
state_occupancy_uc$id <- gsub("ind_", "", rownames(state_occupancy_uc))
# Reshape the data to long format
state_occupancy_uc_long <- melt(state_occupancy_uc, id.vars = "id", variable.name = "cycle", value.name = "state")
# Convert cycle to numeric
state_occupancy_uc_long$cycle <- as.numeric(gsub("cycle_", "", state_occupancy_uc_long$cycle))
# Convert state to factor
state_occupancy_uc_long$state <- factor(state_occupancy_uc_long$state, levels = v_states_names)
head(state_occupancy_uc_long)## id cycle state
## 1 1 0 NE
## 2 2 0 NE
## 3 3 0 NE
## 4 4 0 NE
## 5 5 0 NE
## 6 6 0 NE18 Plot the state occupancy over cycles for patients on Usual Care intervention
ggplot(state_occupancy_uc_long, aes(x = cycle, y = id, fill = state)) +
geom_tile() +
scale_fill_manual(values = c("NE" = "green", "ST" = "yellow", "PS" = "orange", "MI" = "red", "PM" = "purple", "D" = "black")) +
labs(title = "Health State Occupancy Over Cycles (Usual Care Intervention)",
x = "Cycle",
y = "Individual ID",
fill = "Health State") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Plot distribution of state occupancy across the six states for each
cycle for patients on Usual Care intervention
# Compute the proportions of individuals in each state for each cycle for patients on Usual Care intervention
state_proportions_uc <- state_occupancy_uc_long %>%
group_by(cycle, state) %>%
summarise(proportion = n() / nrow(state_occupancy_uc_long)) %>%
ungroup()## `summarise()` has grouped output by 'cycle'. You can override using the
## `.groups` argument.state_proportions_uc <- state_occupancy_uc_long %>%
group_by(cycle, state) %>%
summarise(count = n(), .groups = "drop") %>%
group_by(cycle) %>%
mutate(proportion = count / sum(count))
head(state_proportions_uc)## # A tibble: 6 × 4
## # Groups: cycle [3]
## cycle state count proportion
## <dbl> <fct> <int> <dbl>
## 1 0 NE 100000 1
## 2 1 NE 98899 0.989
## 3 1 ST 222 0.00222
## 4 1 MI 179 0.00179
## 5 1 D 700 0.007
## 6 2 NE 97778 0.978# Plot the distribution of state occupancy across the six states for each cycle for patients on Usual Care intervention
ggplot(state_proportions_uc, aes(x = cycle, y = proportion, fill = state)) +
geom_bar(stat = "identity", position = "fill") +
scale_fill_manual(values = c("NE" = "green", "ST" = "yellow", "PS" = "orange", "MI" = "red", "PM" = "purple", "D" = "black")) +
labs(title = "Distribution of Health State Occupancy per Cycle (Usual Care Intervention)",
x = "Cycle",
y = "Proportion",
fill = "Health State") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Plot stacked area chart for patients on Usual Care intervention
# Plot stacked area chart for patients on Usual Care intervention
ggplot(state_proportions_uc, aes(x = cycle, y = proportion, fill = state)) +
geom_area() +
scale_fill_manual(values = c("NE" = "green", "ST" = "yellow", "PS" = "orange", "MI" = "red", "PM" = "purple", "D" = "black")) +
labs(title = "Health State Distribution Over Time (Usual Care Intervention)",
x = "Cycle",
y = "Proportion",
fill = "Health State") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Counting number of CVD events for patients on Usual Care
intervention
# Count the number of CVD events(MI and Stroke) for each individual for patients on Usual Care intervention
cvd_events_uc <- state_occupancy_uc_long %>%
filter(state %in% c("MI", "ST")) %>%
group_by(id) %>%
summarise(cvd_events = n()) %>%
ungroup()
# Print the first few rows of the cvd_events_uc data frame
head(cvd_events_uc)## # A tibble: 6 × 2
## id cvd_events
## <chr> <int>
## 1 10 1
## 2 100 1
## 3 1000 1
## 4 10000 1
## 5 100000 2
## 6 10002 119 Plot the distribution of CVD events for patients on Usual Care intervention
ggplot(cvd_events_uc, aes(x = cvd_events)) +
geom_histogram(binwidth = 1, fill = "blue", color = "black", alpha = 0.7) +
labs(title = "Distribution of CVD Events (Usual Care Intervention)",
x = "Number of CVD Events",
y = "Frequency") +
theme_minimal()
20 Count the number of CVD events(MI and Stroke) for each individual-separate stroke and MI events for patients on Usual Care intervention
# Count the number of MI events for each individual for patients on Usual Care intervention
mi_events_uc <- state_occupancy_uc_long %>%
filter(state == "MI") %>%
group_by(id) %>%
summarise(mi_events = n()) %>%
ungroup()
# Count the number of Stroke events for each individual for patients on Usual Care intervention
stroke_events_uc <- state_occupancy_uc_long %>%
filter(state == "ST") %>%
group_by(id) %>%
summarise(stroke_events = n()) %>%
ungroup()
# Print the first few rows of the mi_events_uc and stroke_events_uc data frames
head(mi_events_uc)## # A tibble: 6 × 2
## id mi_events
## <chr> <int>
## 1 10 1
## 2 100000 2
## 3 10002 1
## 4 10003 3
## 5 10004 1
## 6 1001 1## # A tibble: 6 × 2
## id stroke_events
## <chr> <int>
## 1 100 1
## 2 1000 1
## 3 10000 1
## 4 10006 1
## 5 10011 1
## 6 10014 121 Plot the distribution of MI events for patients on Usual Care intervention
ggplot(mi_events_uc, aes(x = mi_events)) +
geom_histogram(binwidth = 1, fill = "red", color = "black", alpha = 0.7) +
labs(title = "Distribution of MI Events (Usual Care Intervention)",
x = "Number of MI Events",
y = "Frequency") +
theme_minimal()
# Plot the distribution of Stroke events for patients on Usual Care
intervention
ggplot(stroke_events_uc, aes(x = stroke_events)) +
geom_histogram(binwidth = 1, fill = "green", color = "black", alpha = 0.7) +
labs(title = "Distribution of Stroke Events (Usual Care Intervention)",
x = "Number of Stroke Events",
y = "Frequency") +
theme_minimal()
# Plot the distribution of MI and Stroke events for patients on Usual
Care intervention
# Combine the MI and Stroke events data frames for patients on Usual Care intervention
cvd_events_combined_uc <- merge(mi_events_uc, stroke_events_uc, by = "id", all = TRUE)
# Replace NA values with 0
cvd_events_combined_uc[is.na(cvd_events_combined_uc)] <- 0
# Reshape the data to long format
cvd_events_long_uc <- melt(cvd_events_combined_uc, id.vars = "id", variable.name = "event_type", value.name = "event_count")
# Convert event_type to factor
cvd_events_long_uc$event_type <- factor(cvd_events_long_uc$event_type, levels = c("mi_events", "stroke_events"), labels = c("MI Events", "Stroke Events"))
# Plot the distribution of MI and Stroke events for patients on Usual Care intervention
ggplot(cvd_events_long_uc, aes(x = event_count, fill = event_type)) +
geom_histogram(binwidth = 1, position = "identity", alpha = 0.7) +
labs(title = "Distribution of MI and Stroke Events (Usual Care Intervention)",
x = "Number of Events",
y = "Frequency") +
scale_fill_manual(values = c("MI Events" = "red", "Stroke Events" = "green")) +
theme_minimal()
# Plot the number of CVD events per person over their lifetime for
patients on Usual Care intervention
# Combine the MI and Stroke events data frames for patients on Usual Care intervention
cvd_events_combined_uc <- merge(mi_events_uc, stroke_events_uc, by = "id", all = TRUE)
# Replace NA values with 0
cvd_events_combined_uc[is.na(cvd_events_combined_uc)] <- 0
# Reshape the data to long format
cvd_events_long_uc <- melt(cvd_events_combined_uc, id.vars = "id", variable.name = "event_type", value.name = "event_count")
# Convert event_type to factor
cvd_events_long_uc$event_type <- factor(cvd_events_long_uc$event_type, levels = c("mi_events", "stroke_events"), labels = c("MI Events", "Stroke Events"))
# Summarise the data to get the total number of events per person for patients on Usual Care intervention
cvd_events_summary_uc <- cvd_events_long_uc %>%
group_by(id) %>%
summarise(total_mi_events = sum(event_count[event_type == "MI Events"]),
total_stroke_events = sum(event_count[event_type == "Stroke Events"])) %>%
ungroup()
# Plot the number of CVD events per person over their lifetime for patients on Usual Care intervention
ggplot(cvd_events_summary_uc, aes(x = id)) +
geom_bar(aes(y = total_mi_events), stat = "identity", fill = "red", alpha = 0.7) +
geom_bar(aes(y = total_stroke_events), stat = "identity", fill = "green", alpha = 0.7) +
labs(title = "Number of CVD Events per Person Over Their Lifetime (Usual Care Intervention)",
x = "Individual ID",
y = "Number of Events") +
scale_y_continuous(sec.axis = sec_axis(~., name = "Number of Events")) +
theme_minimal()
hist(cvd_events_summary_uc$total_mi_events, breaks = 20, main = "Distribution of MI Events (Usual Care Intervention)", xlab = "Number of MI Events", col = "red")
# Count the number of deaths for patients on Usual Care intervention
# Count the total number of deaths for all individuals for patients on Usual Care intervention- Only take the first instance of death per individual
deaths_uc <- state_occupancy_uc_long %>%
filter(state == "D") %>%
distinct(id) %>%
mutate(deaths = 1)
# Compute the total number of deaths across all individuals for patients on Usual Care intervention
total_deaths_uc <- nrow(deaths_uc)
# Print the total number of deaths for patients on Usual Care intervention
cat("Total number of deaths (Usual Care Intervention):", total_deaths_uc, "\n")## Total number of deaths (Usual Care Intervention): 7301022 Count only the first instance of MI and stroke events per individual for patients on Usual Care intervention
# Count only the first instance of MI events per individual for patients on Usual Care intervention
mi_first_instance_uc <- state_occupancy_uc_long %>%
filter(state == "MI") %>%
group_by(id) %>%
summarise(mi_first_instance = n() > 0) %>%
ungroup()
# Count only the first instance of Stroke events per individual for patients on Usual Care intervention
stroke_first_instance_uc <- state_occupancy_uc_long %>%
filter(state == "ST") %>%
group_by(id) %>%
summarise(stroke_first_instance = n() > 0) %>%
ungroup()
# View the first few rows of the mi_first_instance_emp and stroke_first_instance_emp data frames
head(mi_first_instance_uc)## # A tibble: 6 × 2
## id mi_first_instance
## <chr> <lgl>
## 1 10 TRUE
## 2 100000 TRUE
## 3 10002 TRUE
## 4 10003 TRUE
## 5 10004 TRUE
## 6 1001 TRUE## # A tibble: 6 × 2
## id stroke_first_instance
## <chr> <lgl>
## 1 100 TRUE
## 2 1000 TRUE
## 3 10000 TRUE
## 4 10006 TRUE
## 5 10011 TRUE
## 6 10014 TRUE# Count the number of individuals with MI events for patients on Empower Health intervention
mi_first_count_uc <- nrow(mi_first_instance_uc[mi_first_instance_uc$mi_first_instance == TRUE, ])
# Count the number of individuals with Stroke events for patients on Empower Health intervention
stroke_first_count_uc <- nrow(stroke_first_instance_uc[stroke_first_instance_uc$stroke_first_instance == TRUE, ])
# Print the counts for patients on Empower Health intervention
cat("Number of individuals with first instance of MI events (Usual Care Intervention):", mi_first_count_uc, "\n")## Number of individuals with first instance of MI events (Usual Care Intervention): 34794cat("Number of individuals with first instance of Stroke events (Usual Care Intervention):", stroke_first_count_uc, "\n")## Number of individuals with first instance of Stroke events (Usual Care Intervention): 4078223 Summarize people with each event and both events for patients on Usual Care intervention
# Count the number of individuals with MI events for patients on Usual Care intervention
mi_count_uc <- nrow(cvd_events_uc[cvd_events_uc$cvd_events > 0, ])
# Count the number of individuals with Stroke events for patients on Usual Care intervention
stroke_count_uc <- nrow(stroke_events_uc[stroke_events_uc$stroke_events > 0, ])
# Count the number of individuals with both MI and Stroke events for patients on Usual Care intervention
both_count_uc <- nrow(cvd_events_combined_uc[cvd_events_combined_uc$mi_events > 0 & cvd_events_combined_uc$stroke_events > 0, ])
# Print the counts for patients on Usual Care intervention
cat("Number of individuals with MI events (Usual Care Intervention):", mi_count_uc, "\n")## Number of individuals with MI events (Usual Care Intervention): 67878## Number of individuals with Stroke events (Usual Care Intervention): 40782cat("Number of individuals with both MI and Stroke events (Usual Care Intervention):", both_count_uc, "\n")## Number of individuals with both MI and Stroke events (Usual Care Intervention): 7698cat("Number of individuals with first instance of MI events (Usual Care Intervention):", mi_first_count_uc, "\n")## Number of individuals with first instance of MI events (Usual Care Intervention): 34794cat("Number of individuals with first instance of Stroke events (Usual Care Intervention):", stroke_first_count_uc, "\n")## Number of individuals with first instance of Stroke events (Usual Care Intervention): 40782## Total number of deaths (Usual Care Intervention): 7301024 Plot the costs over cycles for patients on Usual Care intervention
# Clear the workspace from previous data frames for patients on Usual Care intervention
rm(m_States_UC)
# Load the simulation results from rda file for patients on Usual Care intervention
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_Costs_UC.rda') # On Usual Care intervention
costs_uc_df <- as.data.frame(m_Costs_UC)
# Add individual id variable taken from rownames without the prefix "ind_"
costs_uc_df$id <- gsub("ind_", "", rownames(costs_uc_df))
costs_uc_long <- melt(costs_uc_df, id.vars = "id", variable.name = "cycle", value.name = "cost")
# Convert cycle to numeric
costs_uc_long$cycle <- as.numeric(gsub("cycle_", "", costs_uc_long$cycle))
# Plot the costs over cycles for patients on Usual Care intervention
ggplot(costs_uc_long, aes(x = cycle, y = id, fill = cost)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "black") +
labs(title = "Costs Over Cycles (Usual Care Intervention)",
x = "Cycle",
y = "Individual ID",
fill = "Cost") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Plot mean costs per cycle for patients on Usual Care intervention
# Calculate mean costs per cycle for patients on Usual Care intervention
mean_costs_uc <- costs_uc_long %>%
group_by(cycle) %>%
summarise(mean_cost = mean(cost, na.rm = TRUE)) %>%
ungroup()
# Plot mean costs per cycle for patients on Usual Care intervention
ggplot(mean_costs_uc, aes(x = cycle, y = mean_cost)) +
geom_line() +
labs(title = "Mean Costs per Cycle (Usual Care Intervention)",
x = "Cycle",
y = "Mean Cost") +
theme_minimal()
# Plot the DALYs over cycles for patients on Usual Care intervention
# Clear the workspace from previous data frames for patients on Usual Care intervention
rm(m_Costs_UC)
# Load the simulation results from rda file for patients on Usual Care intervention
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/m_Effs_UC.rda') # On Usual Care intervention
effs_uc_df <- as.data.frame(m_Effs_UC)
# Add individual id variable taken from rownames without the prefix "ind_"
effs_uc_df$id <- gsub("ind_", "", rownames(effs_uc_df))
effs_uc_long <- melt(effs_uc_df, id.vars = "id", variable.name = "cycle", value.name = "eff")
# Convert cycle to numeric
effs_uc_long$cycle <- as.numeric(gsub("cycle_", "", effs_uc_long$cycle))
# Plot the DALYs over cycles for patients on Usual Care intervention
ggplot(effs_uc_long, aes(x = cycle, y = id, fill = eff)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "black") +
labs(title = "DALYs Over Cycles (Usual Care Intervention)",
x = "Cycle",
y = "Individual ID",
fill = "DALY") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Plot mean DALYs per cycle for patients on Usual Care intervention
# Calculate mean DALYs per cycle for patients on Usual Care intervention
mean_effs_uc <- effs_uc_long %>%
group_by(cycle) %>%
summarise(mean_eff = mean(eff, na.rm = TRUE)) %>%
ungroup()
# Plot mean DALYs per cycle for patients on Usual Care intervention
ggplot(mean_effs_uc, aes(x = cycle, y = mean_eff)) +
geom_line() +
labs(title = "Mean DALYs per Cycle (Usual Care Intervention)",
x = "Cycle",
y = "Mean DALY") +
theme_minimal()
24.1 Save the Plots
# Save the plots to files
# ggsave("state_occupancy_empower_health.png", width = 10, height = 6)
# ggsave("state_occupancy_usual_care.png", width = 10, height = 6)
# ggsave("cvd_events_empower_health.png", width = 10, height = 6)
# ggsave("cvd_events_usual_care.png", width = 10, height = 6)
# ggsave("costs_empower_health.png", width = 10, height = 6)
# ggsave("costs_usual_care.png", width = 10, height = 6)
# ggsave("effs_empower_health.png", width = 10, height = 6)
# ggsave("effs_usual_care.png", width = 10, height = 6)25 Comparing the results between No intervention, Empower Health and Usual Care interventions
# Combine the state occupancy data for none and both interventions
state_occupancy_combined <- rbind(
mutate(state_occupancy_emp_long, intervention = "Empower Health"),
mutate(state_occupancy_uc_long, intervention = "Usual Care"),
mutate(state_occupancy_long, intervention = "No Intervention")
)
# Plot the state occupancy over cycles for both interventions
ggplot(state_occupancy_combined, aes(x = cycle, y = id, fill = state)) +
geom_tile() +
scale_fill_manual(values = c("NE" = "green", "ST" = "yellow", "PS" = "orange", "MI" = "red", "PM" = "purple", "D" = "black")) +
facet_wrap(~ intervention) +
labs(title = "Health State Occupancy Over Cycles (All Interventions)",
x = "Cycle",
y = "Individual ID",
fill = "Health State") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# Plot the distribution of state occupancy across the six states for
each cycle for none and both interventions
# Combine the state proportions data for none and both interventions
state_proportions_combined <- rbind(
mutate(state_proportions_emp, intervention = "Empower Health"),
mutate(state_proportions_uc, intervention = "Usual Care"),
mutate(state_proportions, intervention = "No Intervention")
)
# Plot the distribution of state occupancy across the six states for each cycle for none and both interventions
state_proportions_all <- ggplot(state_proportions_combined, aes(x = cycle, y = proportion, fill = state)) +
geom_bar(stat = "identity", position = "fill") +
scale_fill_manual(values = c("NE" = "green", "ST" = "yellow", "PS" = "orange", "MI" = "red", "PM" = "purple", "D" = "black")) +
facet_wrap(~ intervention) +
labs(title = "Distribution of Health State Occupancy per Cycle (All Interventions)",
x = "Cycle",
y = "Proportion",
fill = "Health State") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# View the plot
print(state_proportions_all)
26 Plot the distribution of state occupancy across the six states for each cycle for none and both interventions-label proportions in each state for the last cycle
# Plot the distribution of state occupancy across the six states for each cycle for none and both interventions with labels for last cycle
state_proportions_labeled <- ggplot(state_proportions_combined, aes(x = cycle, y = proportion, fill = state)) +
geom_bar(stat = "identity", position = "fill") +
scale_fill_manual(values = c("NE" = "green", "ST" = "yellow", "PS" = "orange", "MI" = "red", "PM" = "purple", "D" = "black")) +
facet_wrap(~ intervention) +
labs(title = "Distribution of Health State Occupancy per Cycle (All Interventions)",
x = "Cycle",
y = "Proportion",
fill = "Health State") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank()) +
geom_text(data = subset(state_proportions_combined, cycle == max(cycle)),
aes(label = paste0(state, ": ", round(proportion * 100, 1), "%")),
position = position_stack(vjust = 0.5),
size = 3,
color = "white")
# View the plot
print(state_proportions_labeled)
27 Summarise the proportion state occupancy at the final cycle for all interventions
# Summarise the proportion state occupancy at the final cycle for all interventions
final_cycle_proportions <- state_proportions_combined %>%
filter(cycle == 50) %>%
select(intervention, state, proportion)## Adding missing grouping variables: `cycle`## # A tibble: 18 × 4
## # Groups: cycle [1]
## cycle intervention state proportion
## <dbl> <chr> <fct> <dbl>
## 1 50 Empower Health NE 0.178
## 2 50 Empower Health ST 0.0144
## 3 50 Empower Health PS 0.0705
## 4 50 Empower Health MI 0.0153
## 5 50 Empower Health PM 0.0644
## 6 50 Empower Health D 0.657
## 7 50 Usual Care NE 0.108
## 8 50 Usual Care ST 0.0117
## 9 50 Usual Care PS 0.0695
## 10 50 Usual Care MI 0.0137
## 11 50 Usual Care PM 0.0673
## 12 50 Usual Care D 0.730
## 13 50 No Intervention NE 0.0712
## 14 50 No Intervention ST 0.0118
## 15 50 No Intervention PS 0.0751
## 16 50 No Intervention MI 0.0148
## 17 50 No Intervention PM 0.0744
## 18 50 No Intervention D 0.753# Create a table of final cycle proportions-convert to percentage
final_cycle_table <- final_cycle_proportions %>%
pivot_wider(names_from = state, values_from = proportion) %>%
mutate(across(-intervention, ~ round(.x * 100, 1))) # Convert to percentage and round to 1 decimal place
# Add total column
final_cycle_table <- final_cycle_table %>%
rowwise() %>%
mutate(Total = sum(c_across(-intervention)))
# Print the final cycle table
print(final_cycle_table)## # A tibble: 3 × 9
## # Rowwise: cycle
## cycle intervention NE ST PS MI PM D Total
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 50 Empower Health 17.8 1.4 7 1.5 6.4 65.7 99.8
## 2 50 Usual Care 10.8 1.2 7 1.4 6.7 73 100.
## 3 50 No Intervention 7.1 1.2 7.5 1.5 7.4 75.3 10028 Plot stacked area chart for all interventions
# Plot stacked area chart for both interventions
state_occupancy_final <- ggplot(state_proportions_combined, aes(x = cycle, y = proportion, fill = state)) +
geom_area() +
scale_fill_manual(values = c("NE" = "green", "ST" = "yellow", "PS" = "orange", "MI" = "red", "PM" = "purple", "D" = "black")) +
facet_wrap(~ intervention) +
labs(title = NULL,
x = "Cycle",
y = "Proportion",
fill = "Health State") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
# View the plot
print(state_occupancy_final)
29 Plot the distribution of CVD events for all interventions
# Combine the CVD events data for none and both interventions
cvd_events_combined <- rbind(
mutate(cvd_events_emp, intervention = "Empower Health"),
mutate(cvd_events_uc, intervention = "Usual Care"),
mutate(cvd_events, intervention = "No Intervention")
)
ggplot(cvd_events_combined, aes(x = cvd_events, fill = intervention)) +
geom_histogram(binwidth = 1, position = "identity", alpha = 0.7) +
labs(title = "Distribution of CVD Events (All Interventions)",
x = "Number of CVD Events",
y = "Frequency") +
scale_fill_manual(values = c("Empower Health" = "green", "Usual Care" = "orange", "No Intervention" = "red")) +
theme_minimal()
# Count the number of events for each intervention
# Count the total number of Mi events for each intervention
cat("Number of individuals with MI events (No Intervention):", mi_count, "\n")## Number of individuals with MI events (No Intervention): 73109## Number of individuals with MI events (Usual Care Intervention): 67878## Number of individuals with MI events (Empower Health Intervention): 58774# Count the total number of Stroke events for each intervention
cat("Number of individuals with Stroke events (No Intervention):", stroke_count, "\n")## Number of individuals with Stroke events (No Intervention): 43704## Number of individuals with Stroke events (Usual Care Intervention): 40782cat("Number of individuals with Stroke events (Empower Health Intervention):", stroke_count_emp, "\n")## Number of individuals with Stroke events (Empower Health Intervention): 35638# Count the total number of individuals with both MI and Stroke events for each intervention
cat("Number of individuals with both MI and Stroke events (No Intervention):", both_count, "\n")## Number of individuals with both MI and Stroke events (No Intervention): 8212cat("Number of individuals with both MI and Stroke events (Usual Care Intervention):", both_count_uc, "\n")## Number of individuals with both MI and Stroke events (Usual Care Intervention): 7698cat("Number of individuals with both MI and Stroke events (Empower Health Intervention):", both_count_emp, "\n")## Number of individuals with both MI and Stroke events (Empower Health Intervention): 6132# Count the number of individuals developing MI event (Only first event)
cat("Number of individuals with first instance of MI events (No Intervention):", mi_first_count, "\n")## Number of individuals with first instance of MI events (No Intervention): 37617cat("Number of individuals with first instance of MI events (Usual Care Intervention):", mi_first_count_uc, "\n")## Number of individuals with first instance of MI events (Usual Care Intervention): 34794cat("Number of individuals with first instance of MI events (Empower Health Intervention):", mi_first_count_emp, "\n")## Number of individuals with first instance of MI events (Empower Health Intervention): 29268# Count the number of individuals developing Stroke event (Only first event)
cat("Number of individuals with first instance of Stroke events (No Intervention):", stroke_first_count, "\n")## Number of individuals with first instance of Stroke events (No Intervention): 43704cat("Number of individuals with first instance of Stroke events (Usual Care Intervention):", stroke_first_count_uc, "\n")## Number of individuals with first instance of Stroke events (Usual Care Intervention): 40782cat("Number of individuals with first instance of Stroke events (Empower Health Intervention):", stroke_first_count_emp, "\n")## Number of individuals with first instance of Stroke events (Empower Health Intervention): 35638# Count the total number of deaths for each intervention
cat("Total number of deaths (No Intervention):", total_deaths, "\n")## Total number of deaths (No Intervention): 75262## Total number of deaths (Usual Care Intervention): 73010## Total number of deaths (Empower Health Intervention): 6569830 Draw graph of MI, stroke and deaths averted for Empower Health and Usual Care interventions compared to No intervention
# Calculate the number of events averted for each intervention
mi_averted_emp <- mi_count - mi_count_emp
mi_averted_uc <- mi_count - mi_count_uc
stroke_averted_emp <- stroke_count - stroke_count_emp
stroke_averted_uc <- stroke_count - stroke_count_uc
deaths_averted_emp <- total_deaths - total_deaths_emp
deaths_averted_uc <- total_deaths - total_deaths_uc
# Create a data frame for plotting
events_averted <- data.frame(
intervention = rep(c("Empower Health", "Usual Care"), each = 3),
event_type = rep(c("MI Events Averted", "Stroke Events Averted", "Deaths Averted"), times = 2),
count = c(mi_averted_emp, stroke_averted_emp, deaths_averted_emp,
mi_averted_uc, stroke_averted_uc, deaths_averted_uc)
)
# Plot the events averted for both interventions
ggplot(events_averted, aes(x = event_type, y = count, fill = intervention)) +
geom_bar(stat = "identity", position = "dodge") +
labs(title = "CVD Events and Deaths Averted by Intervention",
x = "Event Type",
y = "Number of Events Averted",
fill = "Intervention") +
scale_fill_manual(values = c("Empower Health" = "blue", "Usual Care" = "orange")) +
theme_minimal()
# Draw graph of MI, stroke and deaths averted for Empower Health and
Usual Care interventions compared to No intervention-use onliy first
instances of events
# Calculate the number of first instance events averted for each intervention
mi_first_averted_emp <- mi_first_count - mi_first_count_emp
mi_first_averted_uc <- mi_first_count - mi_first_count_uc
stroke_first_averted_emp <- stroke_first_count - stroke_first_count_emp
stroke_first_averted_uc <- stroke_first_count - stroke_first_count_uc
# Create a data frame for plotting
first_events_averted <- data.frame(
intervention = rep(c("Empower Health", "Usual Care"), each = 2),
event_type = rep(c("First MI Events Averted", "First Stroke Events Averted"), times = 2),
count = c(mi_first_averted_emp, stroke_first_averted_emp,
mi_first_averted_uc, stroke_first_averted_uc)
)
# Plot the first instance events averted for both interventions
ggplot(first_events_averted, aes(x = event_type, y = count, fill = intervention)) +
geom_bar(stat = "identity", position = "dodge") +
labs(title = "First Instance CVD Events Averted by Intervention",
x = "Event Type",
y = "Number of First Instance Events Averted",
fill = "Intervention") +
scale_fill_manual(values = c("Empower Health" = "blue", "Usual Care" = "orange")) +
theme_minimal()
# Draw graph of MI, stroke and deaths averted for Empower Health and
Usual Care interventions compared to No intervention-use onliy first
instances of events_add deaths averted and value labels-use green,red
and purple colors
# Create a data frame for plotting with deaths averted included
first_events_averted_with_deaths <- data.frame(
intervention = rep(c("Empower Health", "Usual Care"), each = 3),
event_type = rep(c("First MI Events Averted", "First Stroke Events Averted", "Deaths Averted"), times = 2),
count = c(mi_first_averted_emp, stroke_first_averted_emp, deaths_averted_emp,
mi_first_averted_uc, stroke_first_averted_uc, deaths_averted_uc)
)
# Plot the first instance events averted with deaths for both interventions
first_events_deaths_averted <- ggplot(first_events_averted_with_deaths, aes(x = event_type, y = count, fill = intervention)) +
geom_bar(stat = "identity", position = "dodge") +
geom_text(aes(label = count), position = position_dodge(width = 0.9), vjust = -0.25) +
labs(title = NULL,
x = NULL,
y = "Number of Events Averted",
fill = "Intervention") +
scale_fill_manual(values = c("Empower Health" = "black", "Usual Care" = "green")) +
theme_minimal()
# View the plot
print(first_events_deaths_averted)
31 Plot the costs over cycles for all interventions
# Combine the costs data for none and both interventions
costs_combined <- rbind(
mutate(costs_emp_long, intervention = "Empower Health"),
mutate(costs_uc_long, intervention = "Usual Care"),
mutate(costs_long, intervention = "No Intervention")
)
# Plot the costs over cycles for both interventions
ggplot(costs_combined, aes(x = cycle, y = id, fill = cost)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "black") +
facet_wrap(~ intervention) +
labs(title = "Costs Over Cycles (All Interventions)",
x = "Cycle",
y = "Individual ID",
fill = "Cost") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
32 Plot mean costs per cycle for All interventions
mean_costs_combined <- rbind(
mutate(mean_costs_emp, intervention = "Empower Health"),
mutate(mean_costs_uc, intervention = "Usual Care"),
mutate(mean_costs, intervention = "No Intervention")
)
# Plot mean costs per cycle for both interventions
ggplot(mean_costs_combined, aes(x = cycle, y = mean_cost, color = intervention)) +
geom_line() +
labs(title = "Mean Costs per Cycle (All Interventions)",
x = "Cycle",
y = "Mean Cost") +
scale_color_manual(values = c("Empower Health" = "blue", "Usual Care" = "orange", "No Intervention" = "grey")) +
theme_minimal()
# Plot the DALYs over cycles for both interventions
# Combine the DALYs data for none and both interventions
effs_combined <- rbind(
mutate(effs_emp_long, intervention = "Empower Health"),
mutate(effs_uc_long, intervention = "Usual Care"),
mutate(effs_long, intervention = "No Intervention")
)
# Plot the DALYs over cycles for both interventions
ggplot(effs_combined, aes(x = cycle, y = id, fill = eff)) +
geom_tile() +
scale_fill_gradient(low = "white", high = "black") +
facet_wrap(~ intervention) +
labs(title = "DALYs Over Cycles (All Interventions)",
x = "Cycle",
y = "Individual ID",
fill = "DALY") +
theme_minimal() +
theme(axis.text.y = element_blank(),
axis.ticks.y = element_blank(),
panel.grid.major.y = element_blank(),
panel.grid.minor.y = element_blank())
33 Plot mean DALYs per cycle for both interventions
mean_effs_combined <- rbind(
mutate(mean_effs_emp, intervention = "Empower Health"),
mutate(mean_effs_uc, intervention = "Usual Care"),
mutate(mean_effs, intervention = "No Intervention")
)
# Plot mean DALYs per cycle for both interventions
ggplot(mean_effs_combined, aes(x = cycle, y = mean_eff, color = intervention)) +
geom_line() +
labs(title = "Mean DALYs per Cycle (All Interventions)",
x = "Cycle",
y = "Mean DALY") +
scale_color_manual(values = c("Empower Health" = "blue", "Usual Care" = "orange", "No Intervention" = "grey")) +
theme_minimal()
# Load results lists
rm(list = ls()) # Clear the environment
# Load the lists from rda file
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/res_no_Trt.rda')
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/res_Empower_Health.rda')
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/res_Usual_Care.rda')34 Plot a cost effectiveness plane for the three interventions
# Extract the mean costs and effects from the results lists
mean_costs_no_trt <- res_no_Trt$mean_Dcosts
mean_costs_empower_health <- res_Empower_Health$mean_Dcosts
mean_costs_usual_care <- res_Usual_Care$mean_Dcosts
mean_effects_no_trt <- res_no_Trt$mean_Ddalys
mean_effects_empower_health <- res_Empower_Health$mean_Ddalys
mean_effects_usual_care <- res_Usual_Care$mean_Ddalys
dalys_averted_empower_health <- mean_effects_no_trt - mean_effects_empower_health
dalys_averted_usual_care <- mean_effects_no_trt - mean_effects_usual_care
# Print the mean costs and effects for each intervention
cat("Mean Costs and Effects for No Treatment Intervention:\n")## Mean Costs and Effects for No Treatment Intervention:## [1] 2945.716## [1] 3.861605## Mean Costs and Effects for Empower Health Intervention:## [1] 9292.248## [1] 3.740573## Mean Costs and Effects for Usual Care Intervention:## [1] 8303.236## [1] 3.946795## DALYs Averted for Empower Health Intervention:## [1] 0.121032## DALYs Averted for Usual Care Intervention:## [1] -0.085189635 Plot the cost-effectiveness plane
# Create a data frame for the cost-effectiveness plane
cost_effectiveness_plane <- data.frame(
costs = c(mean_costs_no_trt, mean_costs_empower_health, mean_costs_usual_care),
effects = c(mean_effects_no_trt, mean_effects_empower_health, mean_effects_usual_care),
intervention = c("No Treatment", "Empower Health", "Usual Care")
)
# Plot the cost-effectiveness plane
ggplot(cost_effectiveness_plane, aes(x = effects, y = costs, color = intervention)) +
geom_point(size = 3) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey") +
labs(title = NULL,
x = "Mean Effects (DALYs)",
y = "Mean Costs (Costs)") +
scale_color_manual(values = c("No Treatment" = "red", "Empower Health" = "blue", "Usual Care" = "green")) +
theme_minimal()
# Draw CE plane using dalys averted and costs comparing Empower Health and Usual Care interventions
# Create a data frame for the cost-effectiveness plane comparing Empower Health and Usual Care interventions
cost_effectiveness_plane_empower_vs_usual <- data.frame(
costs = c(mean_costs_empower_health, mean_costs_usual_care),
effects = c(dalys_averted_empower_health, dalys_averted_usual_care),
intervention = c("Empower Health", "Usual Care")
)
# Plot the cost-effectiveness plane comparing Empower Health and Usual Care interventions
ce_plane_empower_vs_uc <- ggplot(cost_effectiveness_plane_empower_vs_usual, aes(x = effects, y = costs, color = intervention)) +
geom_point(size = 3) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey") +
labs(title = NULL,
x = "DALYs Averted",
y = "Costs") +
scale_color_manual(values = c("Empower Health" = "blue", "Usual Care" = "green")) +
theme_minimal()
# View the plot
print(ce_plane_empower_vs_uc)
36 Plot the cost-effectiveness plane comparing Empower Health and Usual Care interventions using DALYs averted and costs-add willingness to pay (WTP) threshold line at $1000
# Define willingness to pay (WTP) threshold
wtp_threshold <- 1000
# Plot the cost-effectiveness plane comparing Empower Health and Usual Care interventions with WTP line
ggplot(cost_effectiveness_plane_empower_vs_usual, aes(x = effects, y = costs, color = intervention)) +
geom_point(size = 3) +
geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey") +
geom_abline(slope = wtp_threshold, intercept = 0, linetype = "dotted", color = "red") +
labs(title = NULL,
x = "DALYs Averted",
y = "Costs") +
scale_color_manual(values = c("Empower Health" = "blue", "Usual Care" = "green")) +
theme_minimal()
37 Compute the incremental cost-effectiveness ratio (ICER) for the two interventions
# Compute DALYs averted
dalys_averted_empower_health <- mean_effects_no_trt - mean_effects_empower_health
dalys_averted_usual_care <- mean_effects_no_trt - mean_effects_usual_care
# Compute the incremental costs and effects for the two interventions (Empower Health versus usual care)
# Compute the ICER for the two interventions
icer_empower_health <- (mean_costs_empower_health-mean_costs_usual_care) / (dalys_averted_empower_health-dalys_averted_usual_care)
icer_empower_health## [1] 4795.86538 Alternatively
# Compute the incremental cost-effectiveness ratio (ICER) for the two interventions
icer_empower_health_vs_none <- (mean_costs_empower_health - mean_costs_no_trt) / (mean_effects_no_trt-mean_effects_empower_health)
icer_usual_care_vs_none <- (mean_costs_usual_care - mean_costs_no_trt) / (mean_effects_no_trt-mean_effects_usual_care)
icer_empower_vs_usual_care <- (mean_costs_empower_health - mean_costs_usual_care) / (mean_effects_usual_care-mean_effects_empower_health )
# Print the ICER values
cat("ICER for Empower Health Intervention vs No Intervention:", icer_empower_health_vs_none, "\n")## ICER for Empower Health Intervention vs No Intervention: 52436.79## ICER for Usual Care Intervention vs No Intervention: -62889.37## ICER for Empower Health versus Usual Care Intervention: 4795.86538.0.0.1 PROCESSING PSA RESULTS- NB- THIS RESULT IS FOR 10000 PATIENTS (THE ABOVE RESULTS ARE FOR 100000 PATIENTSu)
# Import PSA datasets
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/df_PSA_res_Empower_Health.rda')
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/df_PSA_res_No_Intervention.rda')
load('/Users/jamesoguta/Documents/James Oguta/My PhD Folder-2023-2025/Trainings/KenyaCVDModel/data/df_PSA_res_Usual_Care.rda')39 Rename datasets
40 Add strategy to the intervention datasets
41 Combine into one dataset for plotting
42 Reshape data to wide format
psa_wide <- reshape(psa_results,
timevar = "strategy",
idvar = "sim",
direction = "wide")
head(psa_wide)## sim mean_Dcosts.No_Intervention mean_Ddalys.No_Intervention
## 1 1 2208.080 4.714274
## 2 2 2044.875 6.949993
## 3 3 2071.269 6.060695
## 4 4 2054.884 6.310825
## 5 5 2042.061 6.213368
## 6 6 1909.764 7.267212
## mean_Dcosts.Empower_Health mean_Ddalys.Empower_Health mean_Dcosts.Usual_Care
## 1 6086.152 4.542679 6091.434
## 2 6013.275 4.828032 5536.341
## 3 5558.201 6.358854 5711.100
## 4 5573.383 6.251195 5635.255
## 5 6039.080 4.566123 5653.076
## 6 5566.457 6.382662 5370.031
## mean_Ddalys.Usual_Care
## 1 4.457583
## 2 6.754081
## 3 5.850002
## 4 6.096533
## 5 5.998156
## 6 7.05257243 Step 3: Calculate Incremental Cost and Effect
# Calculate incremental cost and incremental DALYs (Usual_Care as comparator)
psa_wide$delta_cost <- psa_wide$mean_Dcosts.Empower_Health - psa_wide$mean_Dcosts.Usual_Care
psa_wide$delta_dalys <- psa_wide$mean_Ddalys.Usual_Care - psa_wide$mean_Ddalys.Empower_Health
# Avoid division by zero in ICER
psa_wide$icer <- with(psa_wide, ifelse(delta_dalys == 0, NA, delta_cost / delta_dalys))
# Summary of incremental costs and DALYs
summary(psa_wide$delta_cost)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -762.66 -205.94 13.97 16.61 251.95 757.11## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -2.7147 -0.6980 0.0957 0.1010 0.9074 2.8224## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -11306.3 255.2 281.4 287.2 305.5 15717.744 Cost-Effectiveness Plane (CE Plane) Plot
# Clean NA from icers for plotting points
valid <- !is.na(psa_wide$delta_dalys) & !is.na(psa_wide$delta_cost)
plot(psa_wide$delta_dalys[valid], psa_wide$delta_cost[valid],
xlab = "Incremental DALYs Averted (Usual Care vs Empower Health)",
ylab = "Incremental Cost (USD)",
main = "Cost-Effectiveness Plane",
pch = 19, col = rgb(0, 0, 1, 0.3))
abline(h = 0, col = "gray", lty = 2) # horizontal zero cost line
abline(v = 0, col = "gray", lty = 2) # vertical zero effect line
# Histogram of ICER
hist(psa_wide$icer,
breaks = 50,
main = "Distribution of ICERs (Empower Health vs Usual Care)",
xlab = "ICER (USD per DALY averted)",
col = "lightblue",
xlim = c(min(psa_wide$icer, na.rm=TRUE), quantile(psa_wide$icer, 0.95, na.rm=TRUE)))
# Cost-Effectiveness Acceptability Curve (CEAC)
# Define range of WTP thresholds, e.g., 0 to 5000 USD per DALY averted
wtp_values <- seq(0, 5000, by = 50)
# Initialize vector to store probability cost-effective at each WTP
ceac <- numeric(length(wtp_values))
for(i in seq_along(wtp_values)) {
wtp <- wtp_values[i]
# Cost-effective if ICER <= WTP
ceac[i] <- mean(psa_wide$icer <= wtp, na.rm = TRUE)
}
# Plot CEAC
plot(wtp_values, ceac,
type = "l",
lwd = 2,
col = "darkgreen",
xlab = "Willingness-to-pay Threshold (USD per DALY averted)",
ylab = "Probability Cost-Effective",
main = "Cost-Effectiveness Acceptability Curve")
grid()
## quartz_off_screen
## 3## quartz_off_screen
## 245 Add no intervention options
# Empower_Health vs No_Intervention
psa_wide$delta_cost_EmpNoInt <- psa_wide$mean_Dcosts.Empower_Health - psa_wide$mean_Dcosts.No_Intervention
psa_wide$delta_dalys_EmpNoInt <- psa_wide$mean_Ddalys.No_Intervention - psa_wide$mean_Ddalys.Empower_Health
psa_wide$icer_EmpNoInt <- psa_wide$delta_cost_EmpNoInt / psa_wide$delta_dalys_EmpNoInt
# Usual_Care vs No_Intervention
psa_wide$delta_cost_UCNoInt <- psa_wide$mean_Dcosts.Usual_Care - psa_wide$mean_Dcosts.No_Intervention
psa_wide$delta_dalys_UCNoInt <- psa_wide$mean_Ddalys.No_Intervention - psa_wide$mean_Ddalys.Usual_Care
psa_wide$icer_UCNoInt <- psa_wide$delta_cost_UCNoInt / psa_wide$delta_dalys_UCNoInt
par(mfrow = c(1,2))
# Empower Health vs No Intervention
plot(psa_wide$delta_dalys_EmpNoInt, psa_wide$delta_cost_EmpNoInt,
xlab = "Incremental DALYs Averted (No Intervention vs Empower Health)",
ylab = "Incremental Cost (USD)",
main = "CE Plane: Empower Health vs No Intervention",
pch = 19, col = rgb(1, 0, 0, 0.3))
abline(h = 0, v = 0, col = "gray", lty = 2)
# Usual Care vs No Intervention
plot(psa_wide$delta_dalys_UCNoInt, psa_wide$delta_cost_UCNoInt,
xlab = "Incremental DALYs Averted (No Intervention vs Usual Care)",
ylab = "Incremental Cost (USD)",
main = "CE Plane: Usual Care vs No Intervention",
pch = 19, col = rgb(0, 0, 1, 0.3))
abline(h = 0, v = 0, col = "gray", lty = 2)
46 CE Plane for all the three interventions
# Calculate incremental values including Empower vs Usual
cep_data <- data.frame(
sim = psa_wide$sim,
delta_cost_emp = psa_wide$mean_Dcosts.Empower_Health - psa_wide$mean_Dcosts.No_Intervention,
delta_dalys_emp = psa_wide$mean_Ddalys.No_Intervention - psa_wide$mean_Ddalys.Empower_Health,
delta_cost_uc = psa_wide$mean_Dcosts.Usual_Care - psa_wide$mean_Dcosts.No_Intervention,
delta_dalys_uc = psa_wide$mean_Ddalys.No_Intervention - psa_wide$mean_Ddalys.Usual_Care,
delta_cost_emp_vs_uc = psa_wide$mean_Dcosts.Empower_Health - psa_wide$mean_Dcosts.Usual_Care,
delta_dalys_emp_vs_uc = psa_wide$mean_Ddalys.Usual_Care - psa_wide$mean_Ddalys.Empower_Health
)
# Convert to long format for ggplot
cep_long <- bind_rows(
data.frame(
sim = cep_data$sim,
delta_cost = cep_data$delta_cost_emp,
delta_dalys = cep_data$delta_dalys_emp,
comparison = "Empower Health vs. No Intervention"
),
data.frame(
sim = cep_data$sim,
delta_cost = cep_data$delta_cost_uc,
delta_dalys = cep_data$delta_dalys_uc,
comparison = "Usual Care vs. No Intervention"
),
data.frame(
sim = cep_data$sim,
delta_cost = cep_data$delta_cost_emp_vs_uc,
delta_dalys = cep_data$delta_dalys_emp_vs_uc,
comparison = "Empower Health vs. Usual Care"
)
)
# Plot CEP with WTP threshold of $1000/DALY
ggplot(cep_long, aes(x = delta_dalys, y = delta_cost, color = comparison)) +
geom_point(alpha = 0.4, size = 1.2) +
geom_hline(yintercept = 0, linetype = "dashed", color = "grey50") +
geom_vline(xintercept = 0, linetype = "dashed", color = "grey50") +
# Add WTP line
geom_abline(slope = 1000, intercept = 0, linetype = "dotted", color = "black", size = 1) +
annotate("text", x = max(cep_long$delta_dalys, na.rm = TRUE) * 0.9,
y = 1000 * max(cep_long$delta_dalys, na.rm = TRUE) * 0.9,
label = "WTP = $1000/DALY", hjust = 1, vjust = -0.5, angle = atan(1000)*180/pi, size = 3.5) +
labs(
title = "Cost-Effectiveness Plane at US$1000 WTP Threshold",
x = "Incremental DALYs Averted",
y = "Incremental Costs (USD)",
color = "Comparison"
) +
theme_minimal()## 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.
47 Cost-Effectiveness Acceptability Curve (All Interventions)
# Step 1: Define WTP thresholds
wtp_values <- seq(0, 20000, by = 100)
# Step 2: Pivot long to wide
library(reshape2)
psa_wide_costs <- dcast(psa_results, sim ~ strategy, value.var = "mean_Dcosts")
psa_wide_dalys <- dcast(psa_results, sim ~ strategy, value.var = "mean_Ddalys")
# Step 3: Calculate NMB for each strategy at each WTP
strategies <- c("No_Intervention", "Usual_Care", "Empower_Health")
ceac_matrix <- matrix(NA, nrow = length(wtp_values), ncol = length(strategies))
colnames(ceac_matrix) <- strategies
for (i in seq_along(wtp_values)) {
wtp <- wtp_values[i]
# Calculate NMB: NMB = DALYs averted × WTP − Cost
nmb_df <- data.frame(
No_Intervention = (max(psa_wide_dalys) - psa_wide_dalys$No_Intervention) * wtp - psa_wide_costs$No_Intervention,
Usual_Care = (max(psa_wide_dalys) - psa_wide_dalys$Usual_Care) * wtp - psa_wide_costs$Usual_Care,
Empower_Health = (max(psa_wide_dalys) - psa_wide_dalys$Empower_Health) * wtp - psa_wide_costs$Empower_Health
)
# Find which strategy has the max NMB per simulation
best_strategy <- apply(nmb_df, 1, function(x) colnames(nmb_df)[which.max(x)])
# Compute CEAC (probability each strategy is cost-effective)
ceac_matrix[i, ] <- table(factor(best_strategy, levels = strategies)) / nrow(nmb_df)
}
# Step 4: Plot CEAC using matplot (base R)
matplot(wtp_values, ceac_matrix, type = "l", lty = 1, lwd = 2,
col = c("black", "red", "blue"),
xlab = "Willingness-to-pay Threshold (USD)",
ylab = "Probability Cost-Effective",
main = "Cost-Effectiveness Acceptability Curve")
legend("topright", legend = strategies, col = c("black", "red", "blue"), lty = 1, lwd = 2)
## Summary: Empower Health vs No Intervention## delta_cost_EmpNoInt delta_dalys_EmpNoInt icer_EmpNoInt
## Min. :3123 Min. :-2.4594 Min. :-3592635
## 1st Qu.:3454 1st Qu.:-0.4759 1st Qu.: -3269
## Median :3650 Median : 0.3100 Median : 2194
## Mean :3659 Mean : 0.3217 Mean : 2843
## 3rd Qu.:3872 3rd Qu.: 1.1194 3rd Qu.: 4881
## Max. :4206 Max. : 3.0302 Max. : 1862156##
## Summary: Usual Care vs No Intervention## delta_cost_UCNoInt delta_dalys_UCNoInt icer_UCNoInt
## Min. :3429 Min. :0.1882 Min. :14233
## 1st Qu.:3525 1st Qu.:0.2076 1st Qu.:15823
## Median :3627 Median :0.2160 Median :16591
## Mean :3643 Mean :0.2207 Mean :16572
## 3rd Qu.:3769 3rd Qu.:0.2295 3rd Qu.:17284
## Max. :3893 Max. :0.2692 Max. :19614## Summary: Empower Health vs No Intervention## delta_cost delta_dalys icer
## Min. :-762.66 Min. :-2.7147 Min. :-11306.3
## 1st Qu.:-205.94 1st Qu.:-0.6980 1st Qu.: 255.2
## Median : 13.97 Median : 0.0957 Median : 281.4
## Mean : 16.61 Mean : 0.1010 Mean : 287.2
## 3rd Qu.: 251.95 3rd Qu.: 0.9074 3rd Qu.: 305.5
## Max. : 757.11 Max. : 2.8224 Max. : 15717.7