Graduate seminar on operations research¶
Homework 20¶
Add CABG as a third treatment action. The treatment is much more efficient and expensive compared to PTCA. How does the optimal policy change with your chosen values?
Package documentation at http://search.r-project.org/R/R/library/pomdp/html/solve_POMDP.html
install.packages("pomdp")
library("pomdp")
## This model assumes the patient visits a physician once a month.
## One action can be done each time.
ihd_cost <- -30 # cost of one month with severe IHD
# costs of the actions
PTCA_cost <- -50
medication_cost <- -5
angiogram_cost <- -5
stresstest_cost <- -1
# Transition probabilities each month without treatment
waiting_transitions <- matrix(c(0.9, 0.1, 0.01, 0.99), nrow = 2 , byrow = TRUE)
Treatment model without CABG¶
# define the problem
IHD_diagnostics <- POMDP(
states = c("no_IHD", "IHD"),
actions = c("stress-test", "angiogram", "medication", "PTCA", "wait"),
observations = c("no_IHD", "IHD"),
discount = 0.95, # assuming the intervals are always the same
transition_prob = list(
"stress-test" = waiting_transitions ,
"angiogram" = waiting_transitions,
"wait" = waiting_transitions,
"medication" = matrix(c(
0.9, 0.1,
0.3, 0.7), nrow = 2 , byrow = TRUE) ,
"PTCA" = matrix(c(
0.99, 0.01,
0.9, 0.1), nrow = 2 , byrow = TRUE)
),
observation_prob = list(
"stress-test" = matrix(c(
0.85, 0.15,
0.2, 0.8), nrow = 2 , byrow = TRUE),
"angiogram" = matrix(c(
1, 0,
0, 1), nrow = 2 , byrow = TRUE),
"medication" = "uniform",
"PTCA" = "uniform",
"wait" = "uniform"
) ,
reward = rbind(
R_("medication", "*", "IHD", "*", medication_cost + ihd_cost),
R_("medication", "*", "no_IHD", "*", medication_cost),
R_("PTCA", "*", "IHD", "*", PTCA_cost + ihd_cost),
R_("PTCA", "*", "no_IHD", "*", PTCA_cost),
R_("angiogram", "*", "IHD", "*", angiogram_cost + ihd_cost),
R_("angiogram", "*", "no_IHD", "*", angiogram_cost),
R_("stress-test", "*", "IHD", "*", stresstest_cost + ihd_cost),
R_("stress-test", "*", "no_IHD", "*", stresstest_cost),
R_("wait", "*", "IHD", "*", ihd_cost),
R_("wait", "*", "no_IHD", "*", 0)
) ,
start = "uniform",
)
# solve the problem
sol <- solve_POMDP(IHD_diagnostics)
#sol$solution
#sol$solution$alpha
# visualize the solution
plot_policy_graph(sol)
plot_value_function(sol)
Treatment model with CABG¶
First define a cost of CABG (should be higher than for PTCA)
CABG_cost <- -60
Then define transition and observation probabilities for CABG, as well as reward.
Note that CABG is not an observatory action, and thus should not have better observation probabilities than guessing.
The transition probability "IHD" -> "no_IHD" should be higher with CABG than with PTCA.
It also seems reasonable that patients without IHD are less likely to get it after a CABG than before.
# define the problem
IHD_diagnostics <- POMDP(
states = c("no_IHD", "IHD"),
actions = c("stress-test", "angiogram", "medication", "PTCA", "wait", "CABG"),
observations = c("no_IHD", "IHD"),
discount = 0.95, # assuming the intervals are always the same
transition_prob = list(
"stress-test" = waiting_transitions ,
"angiogram" = waiting_transitions,
"wait" = waiting_transitions,
"medication" = matrix(c(
0.9, 0.1,
0.3, 0.7), nrow = 2 , byrow = TRUE) ,
"PTCA" = matrix(c(
0.99, 0.01,
0.9, 0.1), nrow = 2 , byrow = TRUE),
# Transition probabilities if CABG is significantly better than PTCA if the patient has IHD, but otherwise the same
"CABG" = matrix(c(
0.99, 0.01,
0.999, 0.001), nrow = 2 , byrow = TRUE)
),
observation_prob = list(
"stress-test" = matrix(c(
0.85, 0.15,
0.2, 0.8), nrow = 2 , byrow = TRUE),
"angiogram" = matrix(c(
1, 0,
0, 1), nrow = 2 , byrow = TRUE),
"medication" = "uniform",
"PTCA" = "uniform",
# We cannot observe anything, thus the observation probabilities give a uniform matrix
"CABG" = "uniform",
"wait" = "uniform"
) ,
reward = rbind(
R_("medication", "*", "IHD", "*", medication_cost + ihd_cost),
R_("medication", "*", "no_IHD", "*", medication_cost),
R_("PTCA", "*", "IHD", "*", PTCA_cost + ihd_cost),
R_("PTCA", "*", "no_IHD", "*", PTCA_cost),
R_("angiogram", "*", "IHD", "*", angiogram_cost + ihd_cost),
R_("angiogram", "*", "no_IHD", "*", angiogram_cost),
R_("stress-test", "*", "IHD", "*", stresstest_cost + ihd_cost),
R_("stress-test", "*", "no_IHD", "*", stresstest_cost),
R_("wait", "*", "IHD", "*", ihd_cost),
R_("wait", "*", "no_IHD", "*", 0),
# The cost is formed in the same way as for the other methods
R_("CABG", "*", "IHD", "*", CABG_cost + ihd_cost),
R_("CABG", "*", "no_IHD", "*", CABG_cost)
) ,
start = "uniform",
)
# solve the problem
sol <- solve_POMDP(IHD_diagnostics)
#sol$solution
#sol$solution$alpha
# visualize the solution
plot_policy_graph(sol)
plot_value_function(sol)
For these values, the overall effect is quite small, but CABG will clearly be the most beneficial option when the certainty of IHD is high. If we have CABG as an option, waiting becomes a much better alternative, since getting IHD is not as bad as before! PTCA becomes a much less relevant option.
Let's try what happens if CABG is a perfect cure, and a bit cheaper than before (not that realistic).
CABG_cost <- -55
# define the problem
IHD_diagnostics <- POMDP(
states = c("no_IHD", "IHD"),
actions = c("stress-test", "angiogram", "medication", "PTCA", "wait", "CABG"),
observations = c("no_IHD", "IHD"),
discount = 0.95, # assuming the intervals are always the same
transition_prob = list(
"stress-test" = waiting_transitions ,
"angiogram" = waiting_transitions,
"wait" = waiting_transitions,
"medication" = matrix(c(
0.9, 0.1,
0.3, 0.7), nrow = 2 , byrow = TRUE) ,
"PTCA" = matrix(c(
0.99, 0.01,
0.9, 0.1), nrow = 2 , byrow = TRUE),
# Transition probabilities if CABG is a "perfect cure"
"CABG" = matrix(c(
1, 0,
1, 0), nrow = 2 , byrow = TRUE)
),
observation_prob = list(
"stress-test" = matrix(c(
0.85, 0.15,
0.2, 0.8), nrow = 2 , byrow = TRUE),
"angiogram" = matrix(c(
1, 0,
0, 1), nrow = 2 , byrow = TRUE),
"medication" = "uniform",
"PTCA" = "uniform",
# We cannot observe anything, thus the observation probabilities give a uniform matrix
"CABG" = "uniform",
"wait" = "uniform"
) ,
reward = rbind(
R_("medication", "*", "IHD", "*", medication_cost + ihd_cost),
R_("medication", "*", "no_IHD", "*", medication_cost),
R_("PTCA", "*", "IHD", "*", PTCA_cost + ihd_cost),
R_("PTCA", "*", "no_IHD", "*", PTCA_cost),
R_("angiogram", "*", "IHD", "*", angiogram_cost + ihd_cost),
R_("angiogram", "*", "no_IHD", "*", angiogram_cost),
R_("stress-test", "*", "IHD", "*", stresstest_cost + ihd_cost),
R_("stress-test", "*", "no_IHD", "*", stresstest_cost),
R_("wait", "*", "IHD", "*", ihd_cost),
R_("wait", "*", "no_IHD", "*", 0),
# The cost is formed in the same way as for the other methods
R_("CABG", "*", "IHD", "*", CABG_cost + ihd_cost),
R_("CABG", "*", "no_IHD", "*", CABG_cost)
) ,
start = "uniform",
)
# solve the problem
sol <- solve_POMDP(IHD_diagnostics)
#sol$solution
#sol$solution$alpha
# visualize the solution
plot_policy_graph(sol)
plot_value_function(sol)
Now PTCA is not a relevant option anymore since there is a more efficient option for just a slightly higher cost. The rest of the policy looks quite similar to before. CABG becomes the chosen plocy for quite many patients in this scenario.
Other values chosen for cost or transition probabilities change this outcome abit but the behaviour will mostly look like this. Most patients are not affected by the change, but for some this addition might be a major improvement.