|
| 1 | +#------------------------------------------------------------------------------# |
| 2 | +#### Calculate cost-effectiveness outcomes #### |
| 3 | +#------------------------------------------------------------------------------# |
| 4 | +#' Calculate cost-effectiveness outcomes |
| 5 | +#' |
| 6 | +#' \code{calculate_ce_out} calculates costs and effects for a given vector of |
| 7 | +#' parameters using a simulation model. |
| 8 | +#' @param l_params_all List with all parameters of decision model |
| 9 | +#' @param n_wtp Willingness-to-pay threshold to compute net benefits |
| 10 | +#' @return A data frame with discounted costs, effectiveness and NMB. |
| 11 | +#' @export |
| 12 | +calculate_ce_out <- function(l_params_all, n_wtp = 100000 ,verbose =verbose){ # User defined |
| 13 | + with(as.list(l_params_all), { |
| 14 | + |
| 15 | + ### Run decision model to get transition dynamics array |
| 16 | + ########################### Process model inputs ########################### |
| 17 | + ## Number of cycles |
| 18 | + n_cycles <- (n_age_max - n_age_init)/cycle_length # time horizon, number of cycles |
| 19 | + ## Age-specific transition probabilities to the Dead state |
| 20 | + # compute mortality rates |
| 21 | + v_r_S1Dage <- v_r_HDage * hr_S1 # Age-specific mortality rate in the Sick state |
| 22 | + v_r_S2Dage <- v_r_HDage * hr_S2 # Age-specific mortality rate in the Sicker state |
| 23 | + # transform rates to probabilities adjusting by cycle length |
| 24 | + p_HS1 <- rate_to_prob(r = r_HS1, t = cycle_length) # constant annual probability of becoming Sick when Healthy conditional on surviving |
| 25 | + p_S1H <- rate_to_prob(r = r_S1H, t = cycle_length) # constant annual probability of becoming Healthy when Sick conditional on surviving |
| 26 | + p_S1S2 <- rate_to_prob(r = r_S1S2, t = cycle_length) # constant annual probability of becoming Sicker when Sick conditional on surviving |
| 27 | + v_p_HDage <- rate_to_prob(v_r_HDage, t = cycle_length) # Age-specific mortality risk in the Healthy state |
| 28 | + v_p_S1Dage <- rate_to_prob(v_r_S1Dage, t = cycle_length) # Age-specific mortality risk in the Sick state |
| 29 | + v_p_S2Dage <- rate_to_prob(v_r_S2Dage, t = cycle_length) # Age-specific mortality risk in the Sicker state |
| 30 | + |
| 31 | + ## Annual transition probability of becoming Sicker when Sick for treatment AB |
| 32 | + # Apply hazard ratio to rate to obtain transition rate of becoming Sicker when |
| 33 | + # Sick for treatment B |
| 34 | + r_S1S2_trtAB <- r_S1S2 * hr_S1S2_trtAB |
| 35 | + # Transform rate to probability to become Sicker when Sick under treatment AB |
| 36 | + # adjusting by cycle length conditional on surviving |
| 37 | + p_S1S2_trtAB <- rate_to_prob(r = r_S1S2_trtAB, t = cycle_length) |
| 38 | + |
| 39 | + ###################### Construct state-transition models ################### |
| 40 | + ## Initial state vector |
| 41 | + # All starting healthy |
| 42 | + v_m_init <- c(H = 1, S1 = 0, S2 = 0, D = 0) # initial state vector |
| 43 | + # Number of health states |
| 44 | + n_states <- length(v_m_init) |
| 45 | + # Health state names |
| 46 | + v_names_states <- names(v_m_init) |
| 47 | + |
| 48 | + ## Initialize cohort trace for age-dependent (ad) cSTM for strategies SoC and AB |
| 49 | + m_M_SoC <- matrix(0, |
| 50 | + nrow = (n_cycles + 1), ncol = n_states, |
| 51 | + dimnames = list(0:n_cycles, v_names_states)) |
| 52 | + # Store the initial state vector in the first row of the cohort trace |
| 53 | + m_M_SoC[1, ] <- v_m_init |
| 54 | + ## Initialize cohort trace for strategy AB |
| 55 | + # Structure and initial states are the same as for SoC |
| 56 | + m_M_strAB <- m_M_SoC |
| 57 | + |
| 58 | + ## Initialize transition dynamics arrays which will capture transitions |
| 59 | + ## from each state to another over time for strategy SoC |
| 60 | + a_A_SoC <- array(0, |
| 61 | + dim = c(n_states, n_states, n_cycles + 1), |
| 62 | + dimnames = list(v_names_states, v_names_states, 0:n_cycles)) |
| 63 | + # Set first slice of a_A_SoC with the initial state vector in its diagonal |
| 64 | + diag(a_A_SoC[, , 1]) <- v_m_init |
| 65 | + # Initialize transition-dynamics array for strategy AB |
| 66 | + # Structure and initial states are the same as for SoC |
| 67 | + a_A_strAB <- a_A_SoC |
| 68 | + |
| 69 | + ## Create transition arrays |
| 70 | + # Initialize 3-D array |
| 71 | + a_P_SoC <- array(0, dim = c(n_states, n_states, n_cycles), |
| 72 | + dimnames = list(v_names_states, v_names_states, 0:(n_cycles - 1))) |
| 73 | + ### Fill in array |
| 74 | + ## From H |
| 75 | + a_P_SoC["H", "H", ] <- (1 - v_p_HDage) * (1 - p_HS1) |
| 76 | + a_P_SoC["H", "S1", ] <- (1 - v_p_HDage) * p_HS1 |
| 77 | + a_P_SoC["H", "D", ] <- v_p_HDage |
| 78 | + ## From S1 |
| 79 | + a_P_SoC["S1", "H", ] <- (1 - v_p_S1Dage) * p_S1H |
| 80 | + a_P_SoC["S1", "S1", ] <- (1 - v_p_S1Dage) * (1 - (p_S1H + p_S1S2)) |
| 81 | + a_P_SoC["S1", "S2", ] <- (1 - v_p_S1Dage) * p_S1S2 |
| 82 | + a_P_SoC["S1", "D", ] <- v_p_S1Dage |
| 83 | + ## From S2 |
| 84 | + a_P_SoC["S2", "S2", ] <- 1 - v_p_S2Dage |
| 85 | + a_P_SoC["S2", "D", ] <- v_p_S2Dage |
| 86 | + ## From D |
| 87 | + a_P_SoC["D", "D", ] <- 1 |
| 88 | + |
| 89 | + ## Initialize transition probability matrix for strategy AB as a copy of SoC's |
| 90 | + a_P_strAB <- a_P_SoC |
| 91 | + |
| 92 | + ## Only need to update the probabilities involving the transition from Sick to Sicker, p_S1S2 |
| 93 | + # From S1 |
| 94 | + a_P_strAB["S1", "S1", ] <- (1 - v_p_S1Dage) * (1 - (p_S1H + p_S1S2_trtAB)) |
| 95 | + a_P_strAB["S1", "S2", ] <- (1 - v_p_S1Dage) * p_S1S2_trtAB |
| 96 | + |
| 97 | + ### Check if transition probability matrices are valid |
| 98 | + ## Check that transition probabilities are [0, 1] |
| 99 | + check_transition_probability(a_P_SoC, verbose = verbose) |
| 100 | + check_transition_probability(a_P_strAB, verbose = verbose) |
| 101 | + ### Check that all rows for each slice of the array sum to 1 |
| 102 | + check_sum_of_transition_array(a_P_SoC, n_states = n_states, n_cycles = n_cycles, verbose = verbose) |
| 103 | + check_sum_of_transition_array(a_P_strAB, n_states = n_states, n_cycles = n_cycles, verbose = verbose) |
| 104 | + |
| 105 | + #### Run Markov model #### |
| 106 | + ## Iterative solution of age-dependent cSTM |
| 107 | + for (t in 1:n_cycles) { |
| 108 | + ## Fill in cohort trace |
| 109 | + # For SoC |
| 110 | + m_M_SoC[t + 1, ] <- m_M_SoC[t, ] %*% a_P_SoC[, , t] |
| 111 | + # For strategy AB |
| 112 | + m_M_strAB[t + 1, ] <- m_M_strAB[t, ] %*% a_P_strAB[, , t] |
| 113 | + |
| 114 | + ## Fill in transition-dynamics array |
| 115 | + # For SoC |
| 116 | + a_A_SoC[, , t + 1] <- diag(m_M_SoC[t, ]) %*% a_P_SoC[, , t] |
| 117 | + # For strategy AB |
| 118 | + a_A_strAB[, , t + 1] <- diag(m_M_strAB[t, ]) %*% a_P_strAB[, , t] |
| 119 | + } |
| 120 | + |
| 121 | + ## Store the cohort traces in a list |
| 122 | + l_m_M <- list(m_M_SoC, |
| 123 | + m_M_strAB) |
| 124 | + names(l_m_M) <- v_names_str |
| 125 | + |
| 126 | + ## Store the transition array for each strategy in a list |
| 127 | + l_a_A <- list(a_A_SoC, |
| 128 | + a_A_strAB) |
| 129 | + names(l_m_M) <- v_names_str |
| 130 | + |
| 131 | + |
| 132 | + #### State Rewards #### |
| 133 | + ## Vector of state utilities under strategy SoC |
| 134 | + v_u_SoC <- c(H = u_H, |
| 135 | + S1 = u_S1, |
| 136 | + S2 = u_S2, |
| 137 | + D = u_D) * cycle_length |
| 138 | + ## Vector of state costs under strategy SoC |
| 139 | + v_c_SoC <- c(H = c_H, |
| 140 | + S1 = c_S1, |
| 141 | + S2 = c_S2, |
| 142 | + D = c_D) * cycle_length |
| 143 | + ## Vector of state utilities under strategy AB |
| 144 | + v_u_strAB <- c(H = u_H, |
| 145 | + S1 = u_trtAB, |
| 146 | + S2 = u_S2, |
| 147 | + D = u_D) * cycle_length |
| 148 | + ## Vector of state costs under strategy AB |
| 149 | + v_c_strAB <- c(H = c_H, |
| 150 | + S1 = c_S1 + c_trtAB, |
| 151 | + S2 = c_S2 + c_trtAB, |
| 152 | + D = c_D) * cycle_length |
| 153 | + |
| 154 | + ## Store the vectors of state utilities for each strategy in a list |
| 155 | + l_u <- list(SoC = v_u_SoC, |
| 156 | + AB = v_u_strAB) |
| 157 | + ## Store the vectors of state cost for each strategy in a list |
| 158 | + l_c <- list(SoC = v_c_SoC, |
| 159 | + AB = v_c_strAB) |
| 160 | + |
| 161 | + # assign strategy names to matching items in the lists |
| 162 | + names(l_u) <- names(l_c) <- names(l_a_A) <- v_names_str |
| 163 | + |
| 164 | + ## create empty vectors to store total utilities and costs |
| 165 | + v_tot_qaly <- v_tot_cost <- vector(mode = "numeric", length = n_str) |
| 166 | + names(v_tot_qaly) <- names(v_tot_cost) <- v_names_str |
| 167 | + |
| 168 | + ## Number of cycles |
| 169 | + n_cycles <- (n_age_max - n_age_init)/cycle_length # time horizon, number of cycles |
| 170 | + |
| 171 | + ## Discount weight for costs and effects |
| 172 | + v_dwc <- 1 / ((1 + d_e * cycle_length) ^ (0:n_cycles)) |
| 173 | + v_dwe <- 1 / ((1 + d_c * cycle_length) ^ (0:n_cycles)) |
| 174 | + |
| 175 | + ## Within-cycle correction (WCC) using Simpson's 1/3 rule |
| 176 | + v_wcc <- darthtools::gen_wcc(n_cycles = n_cycles, |
| 177 | + method = "Simpson1/3") # vector of wcc |
| 178 | + |
| 179 | + #### Loop through each strategy and calculate total utilities and costs #### |
| 180 | + for (i in 1:n_str) { # i <- 1 |
| 181 | + v_u_str <- l_u[[i]] # select the vector of state utilities for the i-th strategy |
| 182 | + v_c_str <- l_c[[i]] # select the vector of state costs for the i-th strategy |
| 183 | + a_A_str <- l_a_A[[i]] # select the transition array for the i-th strategy |
| 184 | + |
| 185 | + #### Array of state rewards #### |
| 186 | + # Create transition matrices of state utilities and state costs for the i-th strategy |
| 187 | + m_u_str <- matrix(v_u_str, nrow = n_states, ncol = n_states, byrow = T) |
| 188 | + m_c_str <- matrix(v_c_str, nrow = n_states, ncol = n_states, byrow = T) |
| 189 | + # Expand the transition matrix of state utilities across cycles to form a transition array of state utilities |
| 190 | + a_R_u_str <- array(m_u_str, |
| 191 | + dim = c(n_states, n_states, n_cycles + 1), |
| 192 | + dimnames = list(v_names_states, v_names_states, 0:n_cycles)) |
| 193 | + # Expand the transition matrix of state costs across cycles to form a transition array of state costs |
| 194 | + a_R_c_str <- array(m_c_str, |
| 195 | + dim = c(n_states, n_states, n_cycles + 1), |
| 196 | + dimnames = list(v_names_states, v_names_states, 0:n_cycles)) |
| 197 | + |
| 198 | + #### Apply transition rewards #### |
| 199 | + # Apply disutility due to transition from H to S1 |
| 200 | + a_R_u_str["H", "S1", ] <- a_R_u_str["H", "S1", ] - du_HS1 |
| 201 | + # Add transition cost per cycle due to transition from H to S1 |
| 202 | + a_R_c_str["H", "S1", ] <- a_R_c_str["H", "S1", ] + ic_HS1 |
| 203 | + # Add transition cost per cycle of dying from all non-dead states |
| 204 | + a_R_c_str[-n_states, "D", ] <- a_R_c_str[-n_states, "D", ] + ic_D |
| 205 | + |
| 206 | + #### Expected QALYs and Costs for all transitions per cycle #### |
| 207 | + # QALYs = life years x QoL |
| 208 | + # Note: all parameters are annual in our example. In case your own case example is different make sure you correctly apply . |
| 209 | + a_Y_c_str <- a_A_str * a_R_c_str |
| 210 | + a_Y_u_str <- a_A_str * a_R_u_str |
| 211 | + |
| 212 | + #### Expected QALYs and Costs per cycle #### |
| 213 | + ## Vector of QALYs and Costs |
| 214 | + v_qaly_str <- apply(a_Y_u_str, 3, sum) # sum the proportion of the cohort across transitions |
| 215 | + v_cost_str <- apply(a_Y_c_str, 3, sum) # sum the proportion of the cohort across transitions |
| 216 | + |
| 217 | + #### Discounted total expected QALYs and Costs per strategy and apply half-cycle correction if applicable #### |
| 218 | + ## QALYs |
| 219 | + v_tot_qaly[i] <- t(v_qaly_str) %*% (v_dwe * v_wcc) |
| 220 | + ## Costs |
| 221 | + v_tot_cost[i] <- t(v_cost_str) %*% (v_dwc * v_wcc) |
| 222 | + } |
| 223 | + |
| 224 | + ## Vector with discounted net monetary benefits (NMB) |
| 225 | + v_nmb <- v_tot_qaly * n_wtp - v_tot_cost |
| 226 | + |
| 227 | + ## data.frame with discounted costs, effectiveness and NMB |
| 228 | + df_ce <- data.frame(Strategy = v_names_str, |
| 229 | + Cost = v_tot_cost, |
| 230 | + Effect = v_tot_qaly, |
| 231 | + NMB = v_nmb) |
| 232 | + |
| 233 | + return(df_ce) |
| 234 | + } |
| 235 | + ) |
| 236 | +} |
| 237 | + |
| 238 | +#------------------------------------------------------------------------------# |
| 239 | +#### Generate a PSA input parameter dataset #### |
| 240 | +#------------------------------------------------------------------------------# |
| 241 | +#' Generate parameter sets for the probabilistic sensitivity analysis (PSA) |
| 242 | +#' |
| 243 | +#' \code{generate_psa_params} generates a PSA dataset of the parameters of the |
| 244 | +#' cost-effectiveness analysis. |
| 245 | +#' @param n_sim Number of parameter sets for the PSA dataset |
| 246 | +#' @param seed Seed for the random number generation |
| 247 | +#' @return A data.frame with a PSA dataset of he parameters of the |
| 248 | +#' cost-effectiveness analysis |
| 249 | +#' @export |
| 250 | +generate_psa_params <- function(n_sim = 1000, seed = 071818){ |
| 251 | + set.seed(seed) # set a seed to be able to reproduce the same results |
| 252 | + df_psa <- data.frame( |
| 253 | + # Transition probabilities (per cycle) |
| 254 | + r_HS1 = rgamma(n_sim, shape = 30, rate = 170 + 30), # constant rate of becoming Sick when Healthy conditional on surviving |
| 255 | + r_S1H = rgamma(n_sim, shape = 60, rate = 60 + 60), # constant rate of becoming Healthy when Sick conditional on surviving |
| 256 | + r_S1S2 = rgamma(n_sim, shape = 84, rate = 716 + 84), # constant rate of becoming Sicker when Sick conditional on surviving |
| 257 | + hr_S1 = rlnorm(n_sim, meanlog = log(3), sdlog = 0.01), # hazard ratio of death in Sick vs healthy |
| 258 | + hr_S2 = rlnorm(n_sim, meanlog = log(10), sdlog = 0.02), # hazard ratio of death in Sicker vs healthy |
| 259 | + |
| 260 | + # Effectiveness of treatment AaB |
| 261 | + hr_S1S2_trtAB = rlnorm(n_sim, meanlog = log(0.6), sdlog = 0.02), # hazard ratio of becoming Sicker when Sick under treatment AB |
| 262 | + |
| 263 | + ## State rewards |
| 264 | + # Costs |
| 265 | + c_H = rgamma(n_sim, shape = 100, scale = 20), # cost of remaining one cycle in state H |
| 266 | + c_S1 = rgamma(n_sim, shape = 177.8, scale = 22.5), # cost of remaining one cycle in state S1 |
| 267 | + c_S2 = rgamma(n_sim, shape = 225, scale = 66.7), # cost of remaining one cycle in state S2 |
| 268 | + c_trtAB = rgamma(n_sim, shape = 156.3, scale = 160), # cost of treatment AB (per cycle) |
| 269 | + c_D = 0, # cost of being in the death state |
| 270 | + |
| 271 | + # Utilities |
| 272 | + u_H = rbeta(n_sim, shape1 = 200, shape2 = 3), # utility when healthy |
| 273 | + u_S1 = rbeta(n_sim, shape1 = 130, shape2 = 45), # utility when sick |
| 274 | + u_S2 = rbeta(n_sim, shape1 = 230, shape2 = 230), # utility when sicker |
| 275 | + u_D = 0, # utility when dead |
| 276 | + u_trtAB = rbeta(n_sim, shape1 = 300, shape2 = 15), # utility when being treated with AB |
| 277 | + |
| 278 | + ## Transition rewards |
| 279 | + du_HS1 = rbeta(n_sim, shape1 = 11, shape2 = 1088), # disutility when transitioning from Healthy to Sick |
| 280 | + ic_HS1 = rgamma(n_sim, shape = 25, scale = 40), # increase in cost when transitioning from Healthy to Sick |
| 281 | + ic_D = rgamma(n_sim, shape = 100, scale = 20) # increase in cost when dying |
| 282 | + ) |
| 283 | + return(df_psa) |
| 284 | +} |
0 commit comments