9.2 Solution methods for implementation of tailored policies and elimination of infeasible interventions
We demonstrate how the tailored and universal blood safety portfolio problems can be solved, using the solution methods described in Sections 3.1 and 3.2 of the paper. We include an implementation in R, for which we use the following packages and source files:
The file portfolio_functions.R contains the expected cost function for the optimal portfolio model.
cost_portfolio_main(z_i, M_li, A_ji,
c_k, P_ik, d_i, w_i, n_i, c_test_j, R_jk, Q_jk,
g, c_mod_l, H_lk,
full_output = 0)
The first three function parameters for this R function are the decision variables:
z_i(\(\textbf{z}\)), a vector whose entries indicate whether donor group \(i\) is defered (0) or accepted (1).M_li(\(\textbf{M}\)), a matrix whose entries indicate whether modification \(l\) is applied to donor group \(i\) (1 if applied).A_ji(\(\textbf{A}\)), a matrix whose entries indicate whether test \(j\) is applied to donor group \(i\) (1 if applied).
The next 11 function parameters are model parameters:
c_k(\(\textbf{c}\)), a matrix whose entries indicate the expected net present monetary cost of releasing a donation infectious with disease \(k\).P_ik(\(\textbf{P}\)), a matrix whose entries indicate the prevalence of disease \(k\) among donors from group \(i\).d_i(\(\textbf{d}\)), a vector whose entries indicate the cost of replacing a deferred donor from group \(i\).w_i(\(\textbf{w}\)), a vector whose entries indicate the cost of processing a donation from donor group \(i\).n_i(\(\textbf{n}\)), a vector whose entries indicate the number of donors in donor group \(i\).c_test_j(\(\boldsymbol\phi\)), a vector whose entries indicate the per-donation cost of applying test \(j\).R_jk(\(\textbf{R}\)), a matrix whose entries indicate the sensitivity of test \(j\) for disease \(k\) (set to 0 if test \(j\) cannot detect disease \(k\)).Q_jk(\(\textbf{Q}\)), a matrix whose entries indicate the specificity of test \(j\) for disease \(k\) (set to 1 if test \(j\) cannot detect disease \(k\)).g(\(g\)) the cost of removing a donation that tests positive, including confirmatory testing and donor notification costs.c_mod_l(\(\boldsymbol\psi\)), a vector whose entries indicate the per-donation cost of applying modification \(l\).H_lk(\(\textbf{H}\)), a matrix whose entries indicate the risk-reducing multiplier for modification \(l\) on disease \(k\) (set to 1 if modification \(l\) does not reduce risk for disease \(k\)).
Finally, the function parameter full_output can be set to 1 if the user wishes to return 5 additional performance metrics in addition to the value of the cost function. Those are the testing cost, modification cost, net monetary cost of released infectious donations, and the average residual risk of releasing infectious donations. The last two metrics are reported separately for each disease and are returned as \(1 \times K\) vectors.
9.2.1 Example problem
To demonstrate the solution methods, we construct an example problem with hypothetical donor groups, diseases, tests, and modifications:
I = 6 #donor groups
J = 4 #4 tests
K = 5 #5 diseases
L = 2 #2 modifications
c_k = matrix(c(120000, 20000, 60000, 39000, 84000),
ncol = 1)
P_ik = matrix(c(1e-4, 1e-5, 1e-3, 1e-6, 1e-7,
5e-4, 7e-5, 1e-3, 1e-6, 1e-7,
1e-4, 1e-5, 1e-3, 1e-5, 1e-7,
2e-3, 1e-5, 9e-3, 1e-6, 1e-7,
1e-4, 1e-5, 4e-3, 1e-6, 1e-7,
1e-4, 1e-5, 1e-3, 1e-6, 1e-3),
nrow = I,
ncol = K,
byrow = TRUE
)
d_i = matrix(rep(90, I),
ncol = 1)
w_i = matrix(rep(20, I),
ncol = 1)
n_i = matrix(c(7736, 600, 450, 580, 534, 100),
ncol = 1)
c_test_j = matrix(c(10, 9, 4, 7),
ncol = 1)
R_jk = matrix(c(0, 0, 0, .97, 0, # Disease 4
.95, .95, 0, .95, 0, # Diseases 1, 2, and 4 multiplex
0, 0, .91, 0, 0, #Disease 3 less sensitive
0, 0, .98, 0, 0), # Disease 3 more sensitive
nrow = J,
ncol = K,
byrow = TRUE
)
Q_jk = matrix(c(1, 1, 1, .97, 1, # Disease 4
.95, .95, 1, .95, 1, # Diseases 1, 2, and 4 multiplex
1, 1, .99, 1, 1, #Disease 3 less sensitive
1, 1, .96, 1, 1), # Disease 3 more sensitive
nrow = J,
ncol = K,
byrow = TRUE
)
g = 60
c_mod_l =matrix(c(220, 75),
ncol = 1)
H_lk = matrix(c(.6, .6, .8, 1, .75,
1, 1, 1, 1, .05),
nrow = L,
ncol = K,
byrow = TRUE
)As an example, we can calculate the expected cost function (and related performance metrics) with no intervention:
z_i = matrix(rep(1, I), ncol=1)
M_li = matrix(rep(0, I*L), nrow = L, ncol = I, byrow = TRUE)
A_ji = matrix(rep(0, I*J), nrow = J, ncol = I, byrow = TRUE)
noint_metrics <- cost_portfolio_main(z_i, M_li, A_ji,
c_k, P_ik, d_i, w_i, n_i, c_test_j, R_jk, Q_jk, g, c_mod_l, H_lk,
full_output = 1)
kable(noint_metrics)| x | |
|---|---|
| obj_cost | 1.467311e+06 |
| yield | 1.000000e+04 |
| test_cost | 0.000000e+00 |
| mod_cost | 0.000000e+00 |
| downstream_net_mon_cost | 1.267311e+06 |
| avg_resid_risk1 | 2.342000e-04 |
| avg_resid_risk2 | 1.360000e-05 |
| avg_resid_risk3 | 1.624200e-03 |
| avg_resid_risk4 | 1.400000e-06 |
| avg_resid_risk5 | 1.010000e-05 |
We can calculate the same performance metrics for a policy of using all tests and modifications:
z_i = matrix(rep(1, I), ncol=1)
M_li = matrix(rep(1, I*L), nrow = L, ncol = I, byrow = TRUE)
A_ji = matrix(rep(1, I*J), nrow = J, ncol = I, byrow = TRUE)
all_tests_mods_metrics <- cost_portfolio_main(z_i, M_li, A_ji,
c_k, P_ik, d_i, w_i, n_i, c_test_j, R_jk, Q_jk, g, c_mod_l, H_lk,
full_output = 1)
kable(all_tests_mods_metrics)| x | |
|---|---|
| obj_cost | 3.776313e+06 |
| yield | 7.877982e+03 |
| test_cost | 3.000000e+05 |
| mod_cost | 2.950000e+06 |
| downstream_net_mon_cost | 1.023505e+04 |
| avg_resid_risk1 | 7.000000e-06 |
| avg_resid_risk2 | 4.000000e-07 |
| avg_resid_risk3 | 2.300000e-06 |
| avg_resid_risk4 | 0.000000e+00 |
| avg_resid_risk5 | 4.000000e-07 |
Comparing these two strategies, we see that ‘All tests and modifications’ decreases the yield, leads to high test and modification costs, and reduces the downstream net monetary cost by 99.2%, by reducing the residual risk for each of the diseases. However, the combined test and modification costs exceed $3 million, leading to a higher overall objective cost. Most likely, the optimal portfolio uses a more limited combination of deferrals, testing, and modifications.
We consider two classes of the optimal portfolio problem. In the tailored portfolios problem, policymakers are willing to apply a different set of tests and/or modifications to each of the donor groups that are not deferred. The number of feasible policies can be calculated as \(\sum_{i=0}^I \binom{I}{i} 2^{i(L+J)}\):
n_pol = 0
for (i in 1:I){
n_pol = n_pol + choose(I, i)*2^(i*(L+J))
}
print(paste0(n_pol, " policies to evaluate"))## [1] "75418890624 policies to evaluate"Policymakers may prefer not to deal with the complexity of differential testing by donor group. In that case, we can solve the problem as a uniform test/modification policy problem, and the number of feasible policies can be calculated as \(1+\sum_{i=1}^I \binom{I}{i} 2^{(L+J)}\)
n_pol = 1
for (i in 1:I){
n_pol = n_pol + choose(I, i)*2^(L+J)
}
print(paste0(n_pol, " policies to evaluate"))## [1] "4033 policies to evaluate"9.2.2 Reducing solution space for tailored policies
When tailored policies are allowed, we can take advantage of the fact that each donor group can be linearly separated in the objective function. To do so, we solve the problem for each donor group separately and then take the best policy from each group. Doing this allows us to evaluate just \(I (1+ 2^{L+J})\) policies:
## [1] "390 policies to evaluate"9.2.3 Eliminating infeasible interventions
We can reduce the search space further by eliminating some tests or modifications from consideration before running the binary integer program. To do so, we calculate \(\gamma_{ji}^{test}\) and \(\gamma_{li}^{mod}\) as described in Section A above. We define the functions inc_obj_cost_single_test and inc_obj_cost_single_mod to calculate \(\gamma^\text{test}_{j,i}\) and \(\gamma^\text{mod}_{l,i}\), respectively:
#For specific test, assess the incremental objective cost compared to 'no intervention
# separately for each donor group
inc_obj_cost_single_test <- function(c_k, P_ik, d_i, w_i, c_test, r_k, q_k, g){
vals = matrix(NA, nrow = nrow(P_ik))
for (i in 1:nrow(vals)){
#Calc percent donations removed by testing
s <- 1 - prod( q_k*(1-matrix(P_ik[i , ])) + (1- r_k) * matrix(P_ik[i , ]))
#calc per-donor incremental cost of single test vs no intervention
vals[i] <- c_test + g*s - (1-s)*t(r_k*matrix(P_ik[i , ])) %*% c_k
}
return(vals)
}
#For specific mod, assess the incremental objective cost compared to 'no intervention
# separately for each donor group
inc_obj_cost_single_mod <- function(c_k, p_k, d_i, w_i, c_mod, h_k, g){
vals = matrix(NA, nrow = nrow(P_ik))
for (i in 1:nrow(vals)){
#Calc per-donor incremental cost of single modification vs. no intervention
vals[i] = c_mod - t((1-h_k)*matrix(P_ik[i , ])) %*% c_k
}
return(vals)
}Next, we use these functions to calculate \(\gamma^\text{test}_{j,i}\) and \(\gamma^\text{mod}_{l,i}\) for our example problem:
gamma_test_ji <- matrix(NA, nrow = J, ncol = I)
row.names(gamma_test_ji) <- paste0("Test ", 1:J)
colnames(gamma_test_ji) <- paste0("Group ", 1:I)
#Loop through tests
for (j in 1:J){
vals <- inc_obj_cost_single_test(c_k = c_k,
P_ik,
d_i,
w_i,
c_test = c_test_j[j],
r_k = matrix(R_jk[j, ]),
q_k = matrix(Q_jk[j, ]),
g)
print(paste0("Test ", j,
" lowers costs in ", length(vals[vals < 0]),
" of ", length(vals),
" donor groups.",
ifelse(length(vals[vals < 0]) == 0,
paste0(" Can eliminate test ", j, "."),"")
)
)
gamma_test_ji[j, ] <- vals
}## [1] "Test 1 lowers costs in 0 of 6 donor groups. Can eliminate test 1."
## [1] "Test 2 lowers costs in 2 of 6 donor groups."
## [1] "Test 3 lowers costs in 6 of 6 donor groups."
## [1] "Test 4 lowers costs in 6 of 6 donor groups."gamma_mod_li <- matrix(NA, nrow = L, ncol = I)
row.names(gamma_mod_li) <- paste0("Mod ", 1:L)
colnames(gamma_mod_li) <- paste0("Group ", 1:I)
#Loop through mods
for (l in 1:L){
vals <- inc_obj_cost_single_mod(c_k = c_k,
P_ik,
d_i,
w_i,
c_mod = c_mod_l[l],
h_k = matrix(H_lk[l, ]),
g)
print(paste0("Mod ", l, " lowers costs in ", length(vals[vals < 0]),
" of ", length(vals),
" donor groups.",
ifelse(length(vals[vals < 0]) == 0,
paste0(" Can eliminate mod ", l, "."),"" )
)
)
gamma_mod_li[l, ] <- vals
}## [1] "Mod 1 lowers costs in 0 of 6 donor groups. Can eliminate mod 1."
## [1] "Mod 2 lowers costs in 1 of 6 donor groups."We obtain the following values for \(\gamma^\text{test}_{ji}\):
| Group 1 | Group 2 | Group 3 | Group 4 | Group 5 | Group 6 | |
|---|---|---|---|---|---|---|
| Test 1 | 11.8 | 11.8 | 11.4 | 11.8 | 11.8 | 11.8 |
| Test 2 | 7.6 | -32.4 | 7.3 | -177.6 | 7.6 | 7.6 |
| Test 3 | -49.4 | -49.4 | -49.4 | -477.4 | -210.6 | -49.4 |
| Test 4 | -46.9 | -46.9 | -46.9 | -493.6 | -215.3 | -46.9 |
and for \(\gamma^{\text{mod}}_{li}\):
| Group 1 | Group 2 | Group 3 | Group 4 | Group 5 | Group 6 | |
|---|---|---|---|---|---|---|
| Mod 1 | 203.1 | 183.4 | 203.1 | 15.9 | 167.1 | 182.1 |
| Mod 2 | 75.0 | 75.0 | 75.0 | 75.0 | 75.0 | -4.8 |
We see that test 1 and modification 1 can be eliminated for all donor groups. Additionally, for the tailored portfolio problem, we can eliminate all policies that use test 2 in donor groups 1, 3, 5, and 6, and we can eliminate policies that use modification 2 in all donor groups but group 6.
If the number of donor groups is large relative to the number of tests and modifications, one can construct a hypothetical ‘maximum prevalence’ donor group by taking the maximum prevalence for each disease across each donor group (i.e., \(\textbf{p}^{\tilde{i}}_k = \max_i \textbf{P}_{i, k}\)) and assessing each intervention on this hypothetical donor group. The advantage is that this computes in \(L+J\) iterations intead of \(I(L+J)\), but the disadvantage is that it may fail to eliminate some interventions that impact risk for multiple diseases (multiplex tests, pathogen reduction).
In our example problem, we were able to eliminate two infeasible interventions (test 1 and modification 1), decreasing our search space by approximately a factor of 4 (from 390 to 102 tailored policies if assessed for each donor group separately; from 4033 to 1009 uniform policies).
9.2.4 Solving the example portfolio problem
We start by solving the tailored portfolio problem where we allow each non-deferred donor group to receive a unique combination of tests and modifications. We construct a matrix whose rows correspond to each possible policy for a single donor group:
#Enumerate all tailored policies for each donor group
pols_tailored_per_donor_group <- as.matrix(rbind(expand.grid(z_i = 1, a1 = 0, a2 = 0:1, a3 = 0:1, a4 = 0:1, m1 = 0, m2 = 0:1),
c(rep(0, 7))))
kable(pols_tailored_per_donor_group)| z_i | a1 | a2 | a3 | a4 | m1 | m2 |
|---|---|---|---|---|---|---|
| 1 | 0 | 0 | 0 | 0 | 0 | 0 |
| 1 | 0 | 1 | 0 | 0 | 0 | 0 |
| 1 | 0 | 0 | 1 | 0 | 0 | 0 |
| 1 | 0 | 1 | 1 | 0 | 0 | 0 |
| 1 | 0 | 0 | 0 | 1 | 0 | 0 |
| 1 | 0 | 1 | 0 | 1 | 0 | 0 |
| 1 | 0 | 0 | 1 | 1 | 0 | 0 |
| 1 | 0 | 1 | 1 | 1 | 0 | 0 |
| 1 | 0 | 0 | 0 | 0 | 0 | 1 |
| 1 | 0 | 1 | 0 | 0 | 0 | 1 |
| 1 | 0 | 0 | 1 | 0 | 0 | 1 |
| 1 | 0 | 1 | 1 | 0 | 0 | 1 |
| 1 | 0 | 0 | 0 | 1 | 0 | 1 |
| 1 | 0 | 1 | 0 | 1 | 0 | 1 |
| 1 | 0 | 0 | 1 | 1 | 0 | 1 |
| 1 | 0 | 1 | 1 | 1 | 0 | 1 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 |
We iterate over each donor group, evaluating our objective function for each policy and identifying the universal portfolio policy that minimizes our objective function.
opt_tailored_pol <- list(); #opt_tailored_metrics <- list()
for (i in 1:I){
min_obj_cost <- 1e20
for (pol in 1:nrow(pols_tailored_per_donor_group)){
obj_cost <- cost_portfolio_main(z_i = matrix(pols_tailored_per_donor_group[pol, 1]),
A_ji = matrix(pols_tailored_per_donor_group[pol, 2:5]),
M_li = matrix(pols_tailored_per_donor_group[pol, 6:7]),
c_k = c_k,
P_ik = matrix(P_ik[i, ], nrow = 1),
d_i = matrix(d_i[i]),
w_i = matrix(w_i[i]),
n_i = matrix(n_i[i]),
c_test_j = c_test_j,
R_jk = R_jk,
Q_jk = Q_jk,
g = g,
c_mod_l = c_mod_l,
H_lk = H_lk,
full_output = 0)
if (obj_cost < min_obj_cost){
min_obj_cost = obj_cost
opt_tailored_pol[[paste0("z", i)]] <- pols_tailored_per_donor_group[pol, ]
}
}
}
opt_tailored_pol <- data.table(matrix(unlist(opt_tailored_pol), byrow = TRUE, nrow=I))
colnames(opt_tailored_pol) <- colnames(pols_tailored_per_donor_group)
opt_tailored_pol[,"i" := row.names(opt_tailored_pol)]| i | z_i | a1 | a2 | a3 | a4 | m1 | m2 |
|---|---|---|---|---|---|---|---|
| 1 | 1 | 0 | 0 | 1 | 0 | 0 | 0 |
| 2 | 1 | 0 | 1 | 1 | 0 | 0 | 0 |
| 3 | 1 | 0 | 0 | 1 | 0 | 0 | 0 |
| 4 | 1 | 0 | 1 | 1 | 1 | 0 | 0 |
| 5 | 1 | 0 | 0 | 0 | 1 | 0 | 0 |
| 6 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
From this we conclude that the optimal tailored policy is to use test 2 in donor groups 2 and 4; test 3 in donor groups 1, 2, 3, and 4; test 4 in donor groups 4 and 5, and to defer donor group 6.
Next we assess the optimal universal policy, where any donor group not deferred receives the same set of tests and modifications. Again we start by creating a list of all feasible donor groups:
#Enumerate all universal policies
# Specify all combinations of deferral policies (0=defer)
defer_options <- expand.grid(replicate(I, 0:1, simplify = FALSE))
names(defer_options) <- paste0("z",1:I)
tests_and_mods <- pols_tailored_per_donor_group[-1,-1]
pols_universal <- matrix(NA,
nrow = 0, ncol = I + J + L)
colnames(pols_universal) <- c(names(defer_options), colnames(tests_and_mods))
for (row in 1:nrow(defer_options)){
if(sum(defer_options[row,]) == 0 ){
pols_universal <- rbind(pols_universal,
rep(0, ncol(pols_universal)))
} else {
pols_universal <- rbind(pols_universal,
cbind(
defer_options[row, ],
tests_and_mods
)
)
}
}
#print first 3 rows of policy matrix to show structure
kable(head(pols_universal, 3)) | z1 | z2 | z3 | z4 | z5 | z6 | a1 | a2 | a3 | a4 | m1 | m2 |
|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 | 0 |
| 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
## [1] 1009Next, we iterate over each policy to identify the policy that minimizes the objective function value:
min_obj_cost <- 1e20
for (pol in 1:nrow(pols_universal)){
output <- cost_portfolio_main(z_i = matrix(unlist(pols_universal[pol, 1:6]), ncol = 1),
A_ji = outer(unlist(pols_universal[pol, 7:10]),
unlist(pols_universal[pol, 1:6])),
M_li = outer(unlist(pols_universal[pol, 11:12]),
unlist(pols_universal[pol, 1:6])),
c_k = c_k,
P_ik = P_ik,
d_i =d_i,
w_i = w_i,
n_i = n_i,
c_test_j = c_test_j,
R_jk = R_jk,
Q_jk = Q_jk,
g = g,
c_mod_l = c_mod_l,
H_lk = H_lk,
full_output = 1)
if (output["obj_cost"] < min_obj_cost){
min_obj_cost = output["obj_cost"]
opt_universal_metrics <- output
opt_universal_policy <- pols_universal[pol, ]
}
}
#Optimal policy
kable(opt_universal_policy)| z1 | z2 | z3 | z4 | z5 | z6 | a1 | a2 | a3 | a4 | m1 | m2 | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 323 | 1 | 0 | 1 | 0 | 1 | 0 | 0 | 0 | 1 | 0 | 0 | 0 |
We conclude from this that the optimal universal policy is to defer donor groups 2, 4, and 6 and to apply test 3 to all collected donations.
Finally, we can compare performance metrics for our optimal tailored and universal portfolio to the two benchmark policies we created above: the ‘all tests and modifications’ scenario, and ‘no intervention scenario’:
opt_tailored_metrics <- cost_portfolio_main(
z_i = as.matrix(opt_tailored_pol[, 2]),
A_ji = t(as.matrix(opt_tailored_pol[, 3:6])),
M_li = t(as.matrix(opt_tailored_pol[, 7:8])),
c_k = c_k,
P_ik = P_ik,
d_i =d_i,
w_i = w_i,
n_i = n_i,
c_test_j = c_test_j,
R_jk = R_jk,
Q_jk = Q_jk,
g = g,
c_mod_l = c_mod_l,
H_lk = H_lk,
full_output = 1)
metrics_to_row <- function(metrics, pol_name){
cbind(pol = pol_name, data.frame(lapply(metrics, function(x) t(data.frame(x)))))
}
dt_metrics <- rbind(metrics_to_row(noint_metrics, "No intervention"),
metrics_to_row(all_tests_mods_metrics, "All test/mods"),
metrics_to_row(opt_universal_metrics, "Opt. universal"),
metrics_to_row(opt_tailored_metrics, "Opt. tailored"))
kable(dt_metrics[, 1:6], row.names = FALSE)| pol | obj_cost | yield | test_cost | mod_cost | downstream_net_mon_cost |
|---|---|---|---|---|---|
| No intervention | 1467311.1 | 10000.000 | 0 | 0 | 1267311.11 |
| All test/mods | 3776312.8 | 7877.982 | 300000 | 2950000 | 10235.05 |
| Opt. universal | 499825.4 | 8623.510 | 34880 | 0 | 162694.08 |
| Opt. tailored | 473016.2 | 9579.677 | 55882 | 0 | 166346.65 |
Compared to the optimal universal policy, the optimal tailored policy has a lower objective function value, a higher donation yield (because no donor groups are deferred), a slightly higher test cost, and a higher net monetary cost due to released infectious donations. For both portfolios, the objective function value is a fraction of the objective function value from either of the benchmarks.