ce_markov.R

ce_res <- ce_markov(start_pop = c(0, 1),
                    p_matrix,
                    c_matrix,
                    q_matrix)

#
ce_markov <- function(start_pop,
                      p_matrix,
                      c_matrix,
                      q_matrix,
                      n_cycles = 12,
                      s_names = NULL,
                      t_names = NULL) {
  
  n_states <- nrow(p_matrix)
  n_treat <- dim(p_matrix)[3]
  
  pop <- array(data = NA,
               dim = c(n_states, n_cycles, n_treat),
               dimnames = list(state = s_names,
                               cycle = NULL,
                               treatment = t_names))
  for (i in 1:n_states) {
    pop[i, cycle = 1, ] <- start_pop[i]
  }
  
  cycle_costs <- array(NA,
                       dim = c(n_treat, n_cycles),
                       dimnames = list(treatment = t_names,
                                       cycle = NULL))
  
  cycle_QALYs <- array(NA,
                       dim = c(n_treat, n_cycles),
                       dimnames = list(treatment = t_names,
                                       cycle = NULL))
  
  total_costs <- setNames(c(NA, NA), t_names)
  total_QALYs <- setNames(c(NA, NA), t_names)
  
  for (i in 1:n_treat) {
    for (j in 2:n_cycles) {
      
      pop[, cycle = j, treatment = i] <-
        pop[, cycle = j - 1, treatment = i] %*% p_matrix[, , treatment = i]
    }
    
    cycle_QALYs[i, ] <-
      q_matrix[treatment = i, ] %*% pop[, , treatment = i]
    
    cycle_costs[i, ] <-
      c_matrix[treatment = i, ] %*% pop[, , treatment = i]
    
    total_costs[i] <- sum(cycle_costs[treatment = i, ])
    total_QALYs[i] <- sum(cycle_QALYs[treatment = i, ])
  }
  
  list(pop = pop,
       cycle_costs = cycle_costs,
       cycle_QALYs = cycle_QALYs,
       total_costs = total_costs,
       total_QALYs = total_QALYs)
}