Nothing
#' Title
#' This function calculates the standard deviation of the model's outputs
#' @param trans_prob_mat_drug transition matrix for the drug arm
#' @param trans_prob_mat_no_drug transition matrix for the no drug arm
#' @param util_drug utility vector for the drug arm
#' @param util_no_drug utility vector for the no drug arm
#' @param cost_drug cost vector for the drug arm
#' @param cost_no_drug cost vector for the no drug arm
#' @param time_horizon null for lifetime horizon or number of years
#' @param inital_state_vec_drug initial state vector for the drug arm
#' @param inital_state_vec_no_drug initial state vector for the no drug arm
#' @param dead_state_num_drug number of dead states in the drug arm
#' @param dead_state_num_no_drug number of dead states in the no drug arm
#' @param discount_rate discount rate in percentages
#' @param cycle_length cycle length in years
#'
#' @return list
#' @export
#'
#' @examples
#' \donttest{
#' prob_vec<-extract_prob(initial_age=45,selected_country='Ethiopia',
#' time_horizon=NULL,gender = 'Female',selected_year=2015,cycle_length = (1/12),
#' prob_data = sample_data)
#' drug_arm_5_years<-matrix(data = c(0,0.029,0.087*0.59,0,0,0.1,0,0.261*0.59,
#' 0,0,0,0,0,0.325,0,0,0,0,1,0,0,0,0,0,1),nrow = 5,ncol = 5,byrow = TRUE)
#' both_arms_after_20_years<-matrix(data = c(0,0,0,0,0,0.1,0,0.261,0,0,0,0,
#' 0,0.325,0,0,0,0,1,0,0,0,0,0,1),nrow = 5,ncol = 5,byrow = TRUE)
#' after_5_years_drug_arm_no_drug_arm_5<-matrix(data = c(0,0.029,0.087,0,0,
#' 0.1,0,0.261,0,0,0,0,0,0.325,0,0,0,0,1,0,0,0,0,0,1),nrow = 5,ncol = 5,byrow = TRUE)
#' drug_arm_5_years<-change_current_prob_to_cycle_prob(drug_arm_5_years,1,(1/12))
#' both_arms_after_20_years<-change_current_prob_to_cycle_prob(both_arms_after_20_years,1,(1/12))
#' after_5_years_drug_arm_no_drug_arm_5<-change_current_prob_to_cycle_prob(
#' after_5_years_drug_arm_no_drug_arm_5,1,(1/12))
#' drug_arm<-array(c(replicate(n=5*12,drug_arm_5_years),replicate(n=15*12,
#' after_5_years_drug_arm_no_drug_arm_5),replicate(n=40*12,
#' both_arms_after_20_years)),dim = c(5,5,60*12))
#' no_drug_arm<-array(c(replicate(n=20*12,after_5_years_drug_arm_no_drug_arm_5),
#' replicate(n=40*12,both_arms_after_20_years)),dim = c(5,5,60*12))
#' up_mat<-prepare_trans_matrix(drug_arm,prob_vec,(1/12))
#' up_mat_no_drug<-prepare_trans_matrix(no_drug_arm,prob_vec,(1/12))
#' cost_vec_ethiopia<-c(0,1323/12,1841/12,0,0)
#' cost_vec_12_months<-c(20000,1323/12,1841/12,0,0)
#' util_vec<-c(0.94,0.82,0.58,0,0)
#' total_cost_vec_no_drug_arm<-t(replicate(n=60*12,cost_vec_ethiopia))
#' total_util_vec_drug<-t(replicate(n=60*12,util_vec))
#' total_util_vec_no_drug<-total_util_vec_drug
#' total_cost_drug_arm<-t(cbind(replicate(n=1,cost_vec_12_months),
#' replicate(n=59*12+11,cost_vec_ethiopia)))
#' ini_vec<-c(1,0,0,0,0)
#' run_both_arms_sim(up_mat,up_mat_no_drug,total_util_vec_drug,total_util_vec_no_drug,
#' total_cost_drug_arm,total_cost_vec_no_drug_arm,NULL,ini_vec,ini_vec,2,2,
#' both_arms_res<-covert_discount_rate(3,1/12),(1/12))
#' SD_main_func(up_mat, up_mat_no_drug, total_util_vec_drug, total_util_vec_no_drug,
#' total_cost_drug_arm,total_cost_vec_no_drug_arm, NULL, ini_vec,ini_vec, 2, 2,
#' covert_discount_rate(3,1/12), (1/12))}
SD_main_func <- function(trans_prob_mat_drug, trans_prob_mat_no_drug, util_drug, util_no_drug, cost_drug, cost_no_drug, time_horizon, inital_state_vec_drug,inital_state_vec_no_drug, dead_state_num_drug,dead_state_num_no_drug, discount_rate, cycle_length){
check_variables_compatibility(trans_prob_mat_drug, util_drug, cost_drug, inital_state_vec_drug, dead_state_num_drug)
check_variables_compatibility(trans_prob_mat_no_drug, util_no_drug, cost_no_drug, inital_state_vec_no_drug, dead_state_num_no_drug)
res_drug <- sim_function(trans_prob_mat_drug, util_drug, cost_drug, time_horizon, inital_state_vec_drug, dead_state_num_drug, discount_rate, cycle_length)
res_drug_SD<-main_sd_function(res_drug$res_matrix,cost_drug,util_drug,trans_prob_mat_drug,dead_state_num_drug,cycle_length,res_drug$res_discount_rate)
res_no_drug <- sim_function(trans_prob_mat_no_drug, util_no_drug, cost_no_drug, time_horizon, inital_state_vec_no_drug, dead_state_num_no_drug, discount_rate, cycle_length)
res_no_drug_SD<-main_sd_function(res_no_drug$res_matrix,cost_no_drug,util_no_drug,trans_prob_mat_no_drug,dead_state_num_no_drug,cycle_length,res_no_drug$res_discount_rate)
discounted_results <- matrix(data = c(as.numeric(res_no_drug$discounted_costs), res_no_drug$discounted_QALYs,res_no_drug$discounted_LYs,as.numeric(res_drug$discounted_costs),res_drug$discounted_QALYs, res_drug$discounted_LYs, sqrt(res_no_drug_SD$cost_sum_d-res_no_drug$discounted_costs^2),sqrt(res_no_drug_SD$QALY_sum_d-res_no_drug$discounted_QALYs^2), sqrt(res_no_drug_SD$LY_sum_d-res_no_drug$discounted_LYs^2), sqrt(res_drug_SD$cost_sum_d-res_drug$discounted_costs^2), sqrt(res_drug_SD$QALY_sum_d-res_drug$discounted_QALYs^2), sqrt(res_drug_SD$LY_sum_d-res_drug$discounted_LYs^2)), nrow = 2, ncol = 6, byrow = TRUE)
non_discounted_results <- matrix(data = c(as.numeric(res_no_drug$costs), res_no_drug$QALYs,res_no_drug$LYs, as.numeric(res_drug$costs), res_drug$QALYs, res_drug$LYs,sqrt(res_no_drug_SD$cost_sum_no_d-res_no_drug$costs^2),sqrt(res_no_drug_SD$QALY_sum_no_d-res_no_drug$QALYs^2), sqrt(res_no_drug_SD$LY_sum_no_d-res_no_drug$LYs^2),sqrt(res_drug_SD$cost_sum_no_d-res_drug$costs^2),sqrt(res_drug_SD$QALY_sum_no_d-res_drug$QALYs^2), sqrt(res_drug_SD$LY_sum_no_d-res_drug$LYs^2)), nrow = 2, ncol = 6, byrow = TRUE)
colnames(discounted_results) <-colnames(non_discounted_results) <- c("No Drug Costs","No Drug QALYs", "No Drug LYs","Drug Costs","Drug QALYs", "Drug LYs")
rownames(discounted_results) <-rownames(non_discounted_results) <- c("Mean", "SD")
discounted_ICER <- (res_drug$discounted_costs - res_no_drug$discounted_costs) / (res_drug$discounted_QALYs - res_no_drug$discounted_QALYs)
non_discounted_ICER <- (res_drug$costs - res_no_drug$costs) / (res_drug$QALYs - res_no_drug$QALYs)
return(list(discounted_results = discounted_results, discounted_ICER = discounted_ICER,non_discounted_results ,non_discounted_ICER=non_discounted_ICER))
}
#' Title
#' part one of the sd calculation
#' @param prob_matrix state probabilities for each time
#' @param discount_rate discount rate in percentage
#' @param currentT current time
#' @param variable LYs/QALYs/Costs
#'
#' @return list
#'
#' @examples
f1_sum<-function(prob_matrix,discount_rate,currentT,variable){
f1_sum_no_discount<-sum(prob_matrix[currentT,]*variable[currentT,]^2)
f1_sum_discount<-sum(prob_matrix[currentT,]*(variable[currentT,]/discount_rate[currentT])^2)
return(list(f1_sum_no_discount=f1_sum_no_discount,f1_sum_discount=f1_sum_discount))
}
#' Title
#' part two of the sd calculation
#' @param prob_matrix state probabilities for each time
#' @param cost cost vector
#' @param util utility vector
#' @param trans_prob_mat probability matrix
#' @param dead_state_num number of dead states
#' @param cycle_length cycle length in years
#' @param discount_rate discount rate in percentages
#'
#' @return list
#'
#' @examples
main_sd_function<-function(prob_matrix,cost,util,trans_prob_mat,dead_state_num,cycle_length,discount_rate){
T_final<-nrow(prob_matrix)
currentT<-2
ly_val_mat<-matrix(data = 0, ncol = ncol(prob_matrix), nrow = nrow(prob_matrix))
ly_val_mat[,1:(ncol(prob_matrix)-dead_state_num)]<-(1/cycle_length^-1)
cost_sum_f_1<-f1_sum(prob_matrix,discount_rate,currentT-1,cost)
util_sum_f_1<-f1_sum(prob_matrix,discount_rate,currentT-1,util/cycle_length^-1)
ly_sum_f_1<-f1_sum(prob_matrix,discount_rate,currentT-1,ly_val_mat)
cost_sum_f_1_no_d<-cost_sum_f_1$f1_sum_no_discount
cost_sum_f_1_d<-cost_sum_f_1$f1_sum_discount
util_sum_f_1_no_d<-util_sum_f_1$f1_sum_no_discount
util_sum_f_1_d<-util_sum_f_1$f1_sum_discount
ly_sum_f_1_no_d<-ly_sum_f_1$f1_sum_no_discount
ly_sum_f_1_d<-ly_sum_f_1$f1_sum_discount
while (currentT<=T_final) {
#a(n)
cost_sum_a<-f1_sum(prob_matrix,discount_rate,currentT,cost)
QALY_sum_a<-f1_sum(prob_matrix,discount_rate,currentT,util/cycle_length^-1)
LY_sum_a<-f1_sum(prob_matrix,discount_rate,currentT,ly_val_mat)
sum_b<-cond_prob_calc(prob_matrix,discount_rate,cost,util,currentT,trans_prob_mat,dead_state_num,cycle_length,ly_val_mat)
cost_sum_f_1_no_d<-rbind(cost_sum_f_1_no_d,sum_b$final_res_cost_no_d+cost_sum_a$f1_sum_no_discount)
cost_sum_f_1_d<-rbind(cost_sum_f_1_d,sum_b$final_res_cost_d+cost_sum_a$f1_sum_discount)
util_sum_f_1_no_d<-rbind(util_sum_f_1_no_d,sum_b$final_res_util_no_d+QALY_sum_a$f1_sum_no_discount)
util_sum_f_1_d<-rbind(util_sum_f_1_d,sum_b$final_res_util_d+QALY_sum_a$f1_sum_discount)
ly_sum_f_1_no_d<-rbind(ly_sum_f_1_no_d,sum_b$final_res_ly_no_d+LY_sum_a$f1_sum_no_discount)
ly_sum_f_1_d<-rbind(ly_sum_f_1_d,sum_b$final_res_ly_d+LY_sum_a$f1_sum_discount)
currentT<-currentT+1
}
return(list(cost_sum_no_d = (sum(cost_sum_f_1_no_d)), cost_sum_d = (sum(cost_sum_f_1_d)), QALY_sum_no_d = (sum(util_sum_f_1_no_d)),QALY_sum_d = (sum(util_sum_f_1_d)), LY_sum_no_d = (sum(ly_sum_f_1_no_d)),LY_sum_d = (sum(ly_sum_f_1_d))))
}
#' Title
#' conditional probability sd calculation
#' @param prob_matrix state probabilities for each time
#' @param discount_rate discount rate in percentages
#' @param cost cost vector
#' @param util utility vector
#' @param currentT current time
#' @param trans_prob_mat transition matrix
#' @param dead_state_num number of dead states
#' @param cycle_length cycle length in years
#' @param ly_val_mat LYs matrix
#'
#' @return list
#'
#' @examples
cond_prob_calc<-function(prob_matrix,discount_rate,cost,util,currentT,trans_prob_mat,dead_state_num,cycle_length,ly_val_mat){
#E[Ct,i/currentT-1]
cond_prob_cost_1_d<-cond_prob_util_1_d<-cond_prob_ly_1_d<-matrix(data = 0,nrow = 0,ncol = ncol(prob_matrix))
cond_prob_cost_1_no_d<-cond_prob_util_1_no_d<-cond_prob_ly_1_no_d<-matrix(data = 0,nrow = 0,ncol = ncol(prob_matrix))
for (t in 1:(currentT-1)) {
for (i in 1:ncol(prob_matrix)) {
ini_vector<-vector(mode = 'numeric', length = ncol(prob_matrix))
ini_vector[i]<-1
cond_prob_1_temp<-sd_sim_function(trans_prob_mat,t+currentT-t-1,ini_vector,dead_state_num,cycle_length,t)
cond_prob_temp<-(cond_prob_1_temp[nrow(cond_prob_1_temp)-1,])
cond_prob_calc_cost_1<-cond_prob_calc_1(cost,discount_rate,prob_matrix,currentT,t,i,cond_prob_temp)
cond_prob_cost_1_no_d<-rbind(cond_prob_cost_1_no_d,cond_prob_calc_cost_1$cond_prob_calc_1_no_d)
cond_prob_cost_1_d<-rbind(cond_prob_cost_1_d,cond_prob_calc_cost_1$cond_prob_calc_1_d)
cond_prob_calc_util_1<-cond_prob_calc_1(util/cycle_length^-1,discount_rate,prob_matrix,currentT,t,i,cond_prob_temp)
cond_prob_util_1_no_d<-rbind(cond_prob_util_1_no_d,cond_prob_calc_util_1$cond_prob_calc_1_no_d)
cond_prob_util_1_d<-rbind(cond_prob_util_1_d,cond_prob_calc_util_1$cond_prob_calc_1_d)
cond_prob_calc_ly_1<-cond_prob_calc_1(ly_val_mat, discount_rate, prob_matrix, currentT, t, i, cond_prob_temp)
cond_prob_ly_1_no_d<-rbind(cond_prob_ly_1_no_d,cond_prob_calc_ly_1$cond_prob_calc_1_no_d)
cond_prob_ly_1_d<-rbind(cond_prob_ly_1_d,cond_prob_calc_ly_1$cond_prob_calc_1_d)
}
}
cond_prob_cost_2_no_d<-cond_prob_util_2_no_d<-cond_prob_ly_2_no_d<-matrix(data = 0,nrow = 0,ncol = ncol(prob_matrix))
cond_prob_cost_2_d<-cond_prob_util_2_d<-cond_prob_ly_2_d<-matrix(data = 0,nrow = 0,ncol = ncol(prob_matrix))
for (i in 1:ncol(prob_matrix)) {
ini_vector<-vector(mode = 'numeric', length = ncol(prob_matrix))
ini_vector[i]<-1
cond_prob_2_temp<-sd_sim_function(trans_prob_mat,currentT+1,ini_vector,dead_state_num,cycle_length,currentT)
cond_prob_calc_cost_2<-cond_prob_calc_2(cost,discount_rate,prob_matrix,currentT,cond_prob_2_temp)
cond_prob_cost_2_no_d<-rbind(cond_prob_cost_2_no_d,cond_prob_calc_cost_2$cond_prob_calc_2_no_d)
cond_prob_cost_2_d<-rbind(cond_prob_cost_2_d,cond_prob_calc_cost_2$cond_prob_calc_2_d)
cond_prob_calc_util_2<-cond_prob_calc_2(util/cycle_length^-1,discount_rate,prob_matrix,currentT,cond_prob_2_temp)
cond_prob_util_2_no_d<-rbind(cond_prob_util_2_no_d,cond_prob_calc_util_2$cond_prob_calc_2_no_d)
cond_prob_util_2_d<-rbind(cond_prob_util_2_d,cond_prob_calc_util_2$cond_prob_calc_2_d)
cond_prob_calc_ly_2<-cond_prob_calc_2(ly_val_mat,discount_rate,prob_matrix,currentT,cond_prob_2_temp)
cond_prob_ly_2_no_d<-rbind(cond_prob_ly_2_no_d,cond_prob_calc_ly_2$cond_prob_calc_2_no_d)
cond_prob_ly_2_d<-rbind(cond_prob_ly_2_d,cond_prob_calc_ly_2$cond_prob_calc_2_d)
}
final_res_cost_no_d<-final_res(prob_matrix,currentT,cond_prob_cost_1_no_d,cond_prob_cost_2_no_d)
final_res_cost_d<-final_res(prob_matrix,currentT,cond_prob_cost_1_d,cond_prob_cost_2_d)
final_res_util_no_d<-final_res(prob_matrix,currentT,cond_prob_util_1_no_d,cond_prob_util_2_no_d)
final_res_util_d<-final_res(prob_matrix,currentT,cond_prob_util_1_d,cond_prob_util_2_d)
final_res_ly_no_d<-final_res(prob_matrix,currentT,cond_prob_ly_1_no_d,cond_prob_ly_2_no_d)
final_res_ly_d<-final_res(prob_matrix,currentT,cond_prob_ly_1_d,cond_prob_ly_2_d)
return(list(final_res_cost_no_d=final_res_cost_no_d,final_res_cost_d=final_res_cost_d,final_res_util_no_d=final_res_util_no_d,final_res_util_d=final_res_util_d,final_res_ly_no_d=final_res_ly_no_d,final_res_ly_d=final_res_ly_d))
}
#' Title
#' part one conditional probability sd calculation
#' @param variable LYs/QALYs/Costs
#' @param discount_rate discount rate in percentages
#' @param prob_matrix state probabilities for each time
#' @param currentT current time
#' @param t time t iterator
#' @param i state i iterator
#' @param cond_prob_temp temp conditional probability
#'
#' @return list
#'
#' @examples
cond_prob_calc_1<-function(variable,discount_rate,prob_matrix,currentT,t,i,cond_prob_temp){
cond_prob_calc_1_no_d<-variable[t,i]*(1/prob_matrix[currentT-1,])*prob_matrix[t,i]*cond_prob_temp
cond_prob_calc_1_no_d[is.infinite(cond_prob_calc_1_no_d) | is.nan(cond_prob_calc_1_no_d)]<-0
cond_prob_calc_1_d<-(variable[t,i]/discount_rate[t])*(1/prob_matrix[currentT-1,])*prob_matrix[t,i]*cond_prob_temp
cond_prob_calc_1_d[is.infinite(cond_prob_calc_1_d) | is.nan(cond_prob_calc_1_d)]<-0
return(list(cond_prob_calc_1_no_d=cond_prob_calc_1_no_d,cond_prob_calc_1_d=cond_prob_calc_1_d))
}
#' Title
#' part two conditional probability sd calculation
#' @param variable LYs/QALYs/Costs
#' @param discount_rate discount rate in percentages
#' @param prob_matrix state probabilities for each time
#' @param currentT current time
#' @param cond_prob_temp temp conditional probability
#'
#' @return list
#'
#' @examples
cond_prob_calc_2<-function(variable,discount_rate,prob_matrix,currentT,cond_prob_temp){
cond_prob_calc_2_no_d<-variable[currentT,]*cond_prob_temp[nrow(cond_prob_temp),]
cond_prob_calc_2_d<-(variable[currentT,]/discount_rate[currentT])*cond_prob_temp[nrow(cond_prob_temp),]
return(list(cond_prob_calc_2_no_d=cond_prob_calc_2_no_d,cond_prob_calc_2_d=cond_prob_calc_2_d))
}
#' Title
#' sum of conditional probability one and two
#' @param prob_matrix state probabilities for each time
#' @param currentT current time
#' @param cond_prob_1 temp conditional probability 1
#' @param cond_prob_2 temp conditional probability 2
#'
#' @return number
#'
#' @examples
final_res<-function(prob_matrix,currentT,cond_prob_1,cond_prob_2){
final_res<-2*sum(prob_matrix[currentT,]*colSums(cond_prob_1)*colSums(cond_prob_2))
return(final_res)
}
#' Title
#' simulation function used in the sd calculation
#' @param trans_prob_mat transition matrix
#' @param time_horizon null for lifetime horizon or number of years
#' @param inital_state_vec initial state vector
#' @param dead_state_num number of dead states
#' @param cycle_length cycle length in years
#' @param ini_t current cycle
#'
#' @return matrix
#'
#' @examples
sd_sim_function<-function(trans_prob_mat,time_horizon,inital_state_vec,dead_state_num,cycle_length,ini_t){
num_of_states<-dim(trans_prob_mat[,,1])[1]
#creates the results matrix according to the number of states
res_matrix<-matrix(data = 0,nrow=0,ncol = num_of_states)
#adds the inital state vector to the results matrix
res_matrix<-rbind(res_matrix,inital_state_vec)
cyc<-ini_t
flag<-TRUE
#runs while the sum of dead states isn't 1 or until cycle = time horizon
while(sum(inital_state_vec[(length(inital_state_vec)-dead_state_num+1):length(inital_state_vec)])<=0.99 & flag){
if(!is.null(time_horizon)){
# if time horizon is defined and it was reached change flag to false
if(cyc>=time_horizon){
flag<-FALSE
}
}
#multiply the state vector with the transition matrix
inital_state_vec<-(inital_state_vec %*% trans_prob_mat[,,cyc])
#add the result to the result matrix
res_matrix<-rbind(res_matrix,inital_state_vec)
cyc<-cyc+1
#calculate the cycle total discount rate
}
#calculate final results
return(res_matrix)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.