--- title: '3-state Markov model in R' subtitle: "With simulation-time dependency and sensitivity analysis" author: "The DARTH workgroup" output: pdf_document: default html_document: default editor_options: markdown: wrap: 72 --- Developed by the Decision Analysis in R for Technologies in Health (DARTH) workgroup. If you use this code as a basis for your own work, please cite our papers (see suggestions below) and in your acknowledgments for both the code and publications please add: **The authors thank the Decision Analysis in R for Technologies in Health (DARTH) workgroup for providing the coding framework in R.** \newpage ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, eval = TRUE) ``` Change `eval` to `TRUE` if you want to knit this document. ```{r} rm(list = ls()) # clear memory (removes all the variables from the workspace) ``` # 01 Load packages ```{r} if (!require('pacman')) install.packages('pacman'); library(pacman) # use this package to conveniently install other packages # load (install if required) packages from CRAN p_load("dplyr", "tidyr", "reshape2", "devtools", "scales", "ellipse", "ggplot2", "lazyeval", "igraph", "truncnorm", "ggraph", "reshape2", "knitr", "stringr", "diagram", "dampack", install = FALSE) # load (install if required) packages from GitHub #install_github("DARTH-git/darthtools", force = TRUE) #Uncomment if there is a newer version p_load_gh("DARTH-git/darthtools") ``` # 02 Load functions ```{r} # all functions are in the darthtools package ``` # 03 Model input ```{r} ## General setup cycle_length <- 1 # cycle length equal to one year (use 1/12 for monthly) n_cycles <- 60 / cycle_length # number of cycles v_names_cycles <- paste("cycle", 0:n_cycles) # cycle names v_names_states <- c("H", "S", "D") # state names, healthy, sick and dead n_states <- length(v_names_states) # number of health states ### Discounting factors # Update according to country guidelines d_c <- 0.03 # annual discount rate for costs d_e <- 0.015 # annual discount rate for QALYs ### Strategies v_names_str <- c("Standard of care", # store the strategy names "Treatment A", "Treatment B") n_str <- length(v_names_str) # number of strategies ## Within-cycle correction (WCC) using Simpson's 1/3 rule, alternative method = "half-cycle" v_wcc <- gen_wcc(n_cycles = n_cycles, method = "Simpson1/3") ### Transition intensity rates r_HS_SoC <- 0.05 # rate of becoming sick when healthy, under standard of care r_HS_trtA <- 0.04 # rate of becoming sick when healthy, under treatment A r_HS_trtB <- 0.02 # rate of becoming sick when healthy, under treatment B r_SD <- 0.1 # rate of dying when sick v_r_HD <- read.csv("data/lifetable_rates.csv")$x # select lifetable mortality rates when healthy = age specific mortality if(n_cycles > length(v_r_HD)){paste("Note: the lifetable data only supports as simulation of", length(v_r_HD), "years", sep = " ")} ### Transition Probabilities ### Converting rates to probabilities # function in darthtools, rate_to_prob based on p = 1 - exp( -r * cycle_length) p_HS_SoC <- rate_to_prob(r = r_HS_SoC, time = cycle_length) # probability of becoming sick when healthy, under SoC p_HS_trtA <- rate_to_prob(r = r_HS_trtA, time = cycle_length) # probability of becoming sick when healthy, under treatment A p_HS_trtB <- rate_to_prob(r = r_HS_trtB, time = cycle_length) # probability of becoming sick when healthy, under treatment B p_SD <- rate_to_prob(r = r_SD, time = cycle_length) # probability of dying when sick v_p_HD <- rate_to_prob(r = v_r_HD, time = cycle_length) # probability of dying when healthy (vector) v_p_HD <- rep( v_p_HD, each = 1 / cycle_length) ### State rewards #### Costs c_H <- 400 # cost of one year in healthy state c_S <- 1000 # cost of one year in sick state c_D <- 0 # cost of one year in dead state c_trtA <- 800 # cost of treatment A (per year) in healthy state c_trtB <- 1500 # cost of treatment B (per year) in healthy state #### Utilities u_H <- 1 # utility when healthy u_S <- 0.5 # utility when sick u_D <- 0 # utility when dead ### Discount weight for costs and effects v_dwc <- 1 / ((1 + (d_c * cycle_length)) ^ (0:n_cycles)) v_dwe <- 1 / ((1 + (d_e * cycle_length)) ^ (0:n_cycles)) ``` # 04 Construct state-transition models ```{r} m_P_diag <- matrix(0, nrow = n_states, ncol = n_states, dimnames = list(v_names_states, v_names_states)) m_P_diag["H", "H" ] = "" m_P_diag["H", "S" ] = "" m_P_diag["H", "D" ] = "" m_P_diag["S", "S" ] = "" m_P_diag["S", "D" ] = "" m_P_diag["D", "D" ] = "" layout.fig <- c(2, 1) plotmat(t(m_P_diag), t(layout.fig), self.cex = 0.5, curve = 0, arr.pos = 0.8, latex = T, arr.type = "curved", relsize = 0.85, box.prop = 0.8, cex = 0.8, box.cex = 0.7, lwd = 1) ``` ## 04.1 Initial state vector ```{r} # All starting healthy v_m_init <- c("H" = 1, "S" = 0, "D" = 0) v_m_init ``` ## 04.2 Initialize cohort traces ```{r} ### Initialize cohort trace for SoC. Set all values to zero # dimensions are number of cycles + 1 (for initial state) and number of states m_M_SoC <- matrix(0, nrow = (n_cycles + 1), ncol = n_states, dimnames = list(v_names_cycles, v_names_states)) # Store the initial state vector in the first row of the cohort trace m_M_SoC[1, ] <- v_m_init ## Initialize cohort traces for treatments A and B # Structure and initial states are the same as for SoC m_M_trtA <- m_M_trtB <- m_M_SoC ``` ## 04.3 Create transition probability arrays ```{r} ## Create transition probability arrays for strategy SoC ### Initialize transition probability array for strategy SoC a_P_SoC <- array(0, # Create 3-D array dim = c(n_states, n_states, n_cycles), dimnames = list(v_names_states, v_names_states, v_names_cycles[-length(v_names_cycles)])) # name the dimensions of the array ### Fill in array ## Standard of Care # from Healthy a_P_SoC["H", "H", ] <- 1 - p_HS_SoC - v_p_HD a_P_SoC["H", "S", ] <- p_HS_SoC a_P_SoC["H", "D", ] <- v_p_HD # from Sick a_P_SoC["S", "S", ] <- 1 - p_SD a_P_SoC["S", "D", ] <- p_SD # from Dead a_P_SoC["D", "D", ] <- 1 ## Treatment A a_P_trtA <- a_P_SoC # Store array, only replace values that are different in the treatment a_P_trtA["H", "H", ] <- 1 - p_HS_trtA - v_p_HD a_P_trtA["H", "S", ] <- p_HS_trtA ## Treatment B a_P_trtB <- a_P_SoC a_P_trtB["H", "H", ] <- 1 - p_HS_trtB - v_p_HD a_P_trtB["H", "S", ] <- p_HS_trtB ## Check if transition array and probabilities are valid # Check that transition probabilities are in [0, 1] check_transition_probability(a_P_SoC, verbose = TRUE) check_transition_probability(a_P_trtA, verbose = TRUE) check_transition_probability(a_P_trtB, verbose = TRUE) # Check that all rows sum to 1 check_sum_of_transition_array(a_P_SoC, n_states = n_states, n_cycles = n_cycles, verbose = TRUE) check_sum_of_transition_array(a_P_trtA, n_states = n_states, n_cycles = n_cycles, verbose = TRUE) check_sum_of_transition_array(a_P_trtB, n_states = n_states, n_cycles = n_cycles, verbose = TRUE) ``` # 05 Run Markov model ```{r} # Iterative solution of age-dependent cSTM for(t in 1:n_cycles){ ## Fill in cohort trace # For SoC m_M_SoC[t + 1, ] <- m_M_SoC[t, ] %*% a_P_SoC[, , t] # For strategy A m_M_trtA[t + 1, ] <- m_M_trtA[t, ] %*% a_P_trtA[, , t] # For strategy B m_M_trtB[t + 1, ] <- m_M_trtB[t, ] %*% a_P_trtB[, , t] } ## Store the cohort traces in a list l_m_M <- list(SoC = m_M_SoC, A = m_M_trtA, B = m_M_trtB) names(l_m_M) <- v_names_str ``` # 06 Plot Outputs ## 06.1 Plot the cohort trace for strategies SoC ```{r, message = FALSE} plot_trace(m_M = m_M_SoC) ``` ## 06.2 Overall Survival (OS) Print the overall survival for the Standard of Care ```{r} v_os_SoC <- 1 - m_M_SoC[, "D"] # calculate the overall survival (OS) probability v_os_SoC <- rowSums(m_M_SoC[, 1:2]) # alternative way of calculating the OS probability plot(v_os_SoC, type = 'l', ylim = c(0, 1), ylab = "Survival probability", xlab = "Cycle", main = "Overall Survival") # create a simple plot showing the OS # add grid grid(nx = n_cycles, ny = 10, col = "lightgray", lty = "dotted", lwd = par("lwd"), equilogs = TRUE) ``` ## 06.2.1 Life Expectancy (LE) ```{r} le_SoC <- sum(v_os_SoC) # summing probability of OS over time (i.e. life expectancy) le_SoC ``` ## 06.2.2 Disease prevalence ```{r} v_prev <- m_M_SoC[, "S"]/v_os_SoC plot(v_prev, ylim = c(0, 1), ylab = "Prevalence", xlab = "Cycle", main = "Disease prevalence") ``` # 07 State Rewards ```{r} ## Scale by the cycle length # Standard of Care # vector of state QALYs accrued each cycle v_u_SoC <- c(H = u_H, S = u_S, D = u_D) * cycle_length # vector of state costs accrued each cycle v_c_SoC <- c(H = c_H, S = c_S, D = c_D) * cycle_length # Treatment A # vector of state QALYs accrued each cycle v_u_trtA <- c(H = u_H, S = u_S, D = u_D) * cycle_length # vector of state costs accrued each cycle v_c_trtA <- c(H = c_H + c_trtA, S = c_S, D = c_D) * cycle_length # Treatment B # vector of state QALYs accrued each cycle v_u_trtB <- c(H = u_H, S = u_S, D = u_D) * cycle_length # vector of state costs accrued each cycle v_c_trtB <- c(H = c_H + c_trtB, S = c_S, D = c_D) * cycle_length ## Store state rewards # Store the vectors of state utilities for each strategy in a list l_u <- list(SoQ = v_u_SoC, A = v_u_trtA, B = v_u_trtB) # Store the vectors of state cost for each strategy in a list l_c <- list(SoQ = v_c_SoC, A = v_c_trtA, B = v_c_trtB) # assign strategy names to matching items in the lists names(l_u) <- names(l_c) <- v_names_str ``` # 08 Compute expected outcomes ```{r} # Create empty vectors to store total utilities and costs v_tot_qaly <- v_tot_cost <- vector(mode = "numeric", length = n_str) names(v_tot_qaly) <- names(v_tot_cost) <- v_names_str ## Loop through each strategy and calculate total utilities and costs for (i in 1:n_str) { v_u_str <- l_u[[i]] # select the vector of state utilities for the i-th strategy v_c_str <- l_c[[i]] # select the vector of state costs for the i-th strategy ### Expected QALYs and costs per cycle ## Vector of QALYs and Costs # Apply state rewards v_qaly_str <- l_m_M[[i]] %*% v_u_str # sum the utilities of all states for each cycle v_cost_str <- l_m_M[[i]] %*% v_c_str # sum the costs of all states for each cycle ### Discounted total expected QALYs and Costs per strategy and apply within-cycle correction if applicable # QALYs v_tot_qaly[i] <- t(v_qaly_str) %*% (v_dwe * v_wcc) # Costs v_tot_cost[i] <- t(v_cost_str) %*% (v_dwc * v_wcc) } ``` # 10 Deterministic Sensitivity Analysis (DSA) ```{r} ## Load model, CEA and PSA functions source('functions/Functions_cSTM_3state_time.R') ``` ## 10.1 Model input for SA ```{r} l_params_all <- list( # Transition probabilities # probability of dying r_HS_SoC = 0.05, # rate of becoming sick when healthy, under standard of care r_HS_trtA = 0.04, # rate of becoming sick when healthy, under treatment A r_HS_trtB = 0.02, # rate of becoming sick when healthy, under treatment B r_SD = 0.1, # rate of dying when sick v_r_HD = v_r_HD, # save vector of rated for dying when healthy # Costs c_H = 400, # cost of one cycle in healthy state c_S = 1000, # cost of one cycle in sick state c_D = 0, # cost of one cycle in dead state c_trtA = 800, # cost of treatment A (per cycle) in healthy state c_trtB = 1500, # cost of treatment B (per cycle) in healthy state # Utilities u_H = 1, # utility when healthy u_S = 0.5, # utility when sick u_D = 0, # utility when dead # Discount rates d_e = 0.03, # discount rate per cycle equal discount of costs and QALYs by 3% d_c = 0.03, # discount rate per cycle equal discount of costs and QALYs by 3% # Time horizon n_cycles = 60 # simulation time horizon (number of cycles) ) ``` **Test model functions** Two functions were defined in the `Functions_markov_3state_time.R` file. The first was the `decision_model()` function which takes in model parameters and outputs the Markov trace for all three strategies. The second is the `calculate_ce_out()` function which calls `decision_model()` and uses the resulting Markov trace to compute the total costs, QALYs, and net monetary benefit (NMB) for each strategy, returning a data frame of costs and QALYs. ```{r} # Try the calculate_ce_out() function df_ce <- calculate_ce_out(l_params_all, n_wtp = 50000 ) df_ce # Get strategies names (will be used to label plots) v_names_str <- df_ce$Strategy n_str <- length(v_names_str) ``` ```{r} ## Incremental cost-effectiveness ratios (ICERs) df_cea <- calculate_icers(cost = df_ce$Cost, effect = df_ce$Effect, strategies = v_names_str) df_cea ``` ## 10.2 One-way sensitivity analysis (OWSA) ```{r} options(scipen = 999) # disabling scientific notation in R # dataframe containing all parameters, their base case values, and the min and # max values of the parameters of interest df_params_owsa <- data.frame(pars = c("c_trtA", "r_HS_SoC", "u_S"), min = c(300 , 0.001, 0.32), # min parameter values max = c(3200, 0.09, 0.75) # max parameter values ) owsa_nmb <- run_owsa_det(params_range = df_params_owsa, # dataframe with parameters for OWSA params_basecase = l_params_all, # list with all parameters nsamp = 100, # number of parameter values FUN = calculate_ce_out, # function to compute outputs outcomes = c("NMB"), # output to do the OWSA on strategies = v_names_str, # names of the strategies n_wtp = 5000) # extra argument to pass to FUN ``` ```{r} plot(owsa_nmb, txtsize = 10, n_x_ticks = 4, facet_scales = "free") + theme(legend.position = "bottom") ``` ### 10.2.1 Optimal strategy with OWSA ```{r} owsa_opt_strat(owsa = owsa_nmb, txtsize = 10) ``` ## 10.3 Two-way sensitivity analysis (TWSA) ## 10.3.1 Cost example (TWSA) ```{r} # dataframe containing all parameters, their basecase values, and the min and # max values of the parameters of interest df_params_twsa <- data.frame(pars = c("c_trtA", "c_trtB"), min = c(300, 500), # min parameter values max = c(3200, 2000) # max parameter values ) twsa_nmb_c <- run_twsa_det(params_range = df_params_twsa, # dataframe with parameters for TWSA params_basecase = l_params_all, # list with all parameters nsamp = 40, # number of parameter values FUN = calculate_ce_out, # function to compute outputs outcomes = "NMB", # output to do the TWSA on strategies = v_names_str, # names of the strategies n_wtp = 5000) # extra argument to pass to FUN ``` ### 10.3.1 Plot TWSA ```{r} plot(twsa_nmb_c) ``` ## 10.3.1 Cost example (TWSA) ```{r} # dataframe containing all parameters, their basecase values, and the min and # max values of the parameters of interest df_params_twsa <- data.frame(pars = c("r_HS_trtA", "r_HS_trtB"), min = c(0.01, 0.01), # min parameter values max = c(0.06, 0.11) # max parameter values ) twsa_nmb_r <- run_twsa_det(params_range = df_params_twsa, # dataframe with parameters for TWSA params_basecase = l_params_all, # list with all parameters nsamp = 40, # number of parameter values FUN = calculate_ce_out, # function to compute outputs outcomes = "NMB", # output to do the TWSA on strategies = v_names_str, # names of the strategies n_wtp = 5000) # extra argument to pass to FUN plot(twsa_nmb_r) ``` Add a tornado plot ```{r} tornado <- owsa_tornado_new(owsa_nmb, params_basecase = l_params_all, select_str = c("Treatment.A","Treatment.B"), outcome_name = "NMB", n_wtp = 5000) tornado ``` # 11 Probabilistic Sensitivity Analysis (PSA) ## 11.1 Model input ```{r} # Store the parameter names into a vector v_names_params <- names(l_params_all) ## Test functions to generate CE outcomes and PSA dataset # Test function to compute CE outcomes calculate_ce_out(l_params_all) # Test function to generate PSA input dataset generate_psa_params(10) ## Generate PSA dataset # Number of simulations n_sim <- 1000 # Generate PSA input dataset df_psa_input <- generate_psa_params(n_sim = n_sim) # First six observations head(df_psa_input) ### Histogram of parameters ggplot(melt(df_psa_input, variable.name = "Parameter"), aes(x = value)) + facet_wrap(~Parameter, scales = "free") + geom_histogram(aes(y = ..density..)) + ylab("") + theme_bw(base_size = 16) + theme(axis.text = element_text(size = 6), axis.title.x = element_blank(), axis.title.y = element_blank(), axis.text.y = element_blank(), axis.ticks.y = element_blank()) ``` ## 11.2 Run PSA ```{r} # Initialize data.frames with PSA output # data.frame of costs df_c <- as.data.frame(matrix(0, nrow = n_sim, ncol = n_str)) colnames(df_c) <- v_names_str # data.frame of effectiveness df_e <- as.data.frame(matrix(0, nrow = n_sim, ncol = n_str)) colnames(df_e) <- v_names_str # Conduct probabilistic sensitivity analysis # Run Markov model on each parameter set of PSA input dataset n_time_init_psa_series <- Sys.time() for (i in 1:n_sim) { # i <- 1 l_psa_input_i <- as.list(df_psa_input[i,]) l_psa_input <- update_param_list(l_params_all, l_psa_input_i) # Outcomes l_out_ce_temp <- calculate_ce_out(l_psa_input) df_c[i, ] <- l_out_ce_temp$Cost df_e[i, ] <- l_out_ce_temp$Effect # Display simulation progress if (i/(n_sim/100) == round(i/(n_sim/100), 0)) { # display progress every 5% cat('\r', paste(i/n_sim * 100, "% done", sep = " ")) } } n_time_end_psa_series <- Sys.time() n_time_total_psa_series <- n_time_end_psa_series - n_time_init_psa_series print(paste0("PSA with ", scales::comma(n_sim), " simulations run in series in ", round(n_time_total_psa_series, 2), " ", units(n_time_total_psa_series))) ``` # 11.3 Visualize PSA results for CEA ```{r} ### Create PSA object l_psa <- make_psa_obj(cost = df_c, effectiveness = df_e, parameters = df_psa_input, strategies = v_names_str) l_psa$strategies <- v_names_str colnames(l_psa$effectiveness) <- v_names_str colnames(l_psa$cost) <- v_names_str # Vector with willingness-to-pay (WTP) thresholds. v_wtp <- seq(0, 30000, by = 1000) ``` ## 11.3.1 Cost-Effectiveness Scatter plot ```{r} ### Cost-Effectiveness Scatter plot txtsize <- 13 gg_scattter <- plot_psa(l_psa, txtsize = txtsize) + ggthemes::scale_color_colorblind() + ggthemes::scale_fill_colorblind() + scale_y_continuous("Cost (Thousand $)", breaks = number_ticks(10), labels = function(x) x/1000) + xlab("Effectiveness (QALYs)") + guides(col = guide_legend(nrow = 2)) + theme(legend.position = "bottom") gg_scattter ``` ## 11.3.2 Incremental cost-effectiveness ratios (ICERs) with probabilistic output ```{r} ### Incremental cost-effectiveness ratios (ICERs) with probabilistic output # Compute expected costs and effects for each strategy from the PSA df_out_ce_psa <- summary(l_psa) df_cea_psa <- calculate_icers(cost = df_out_ce_psa$meanCost, effect = df_out_ce_psa$meanEffect, strategies = df_out_ce_psa$Strategy) df_cea_psa ``` ## 11.3.3 Plot cost-effectiveness frontier with probabilistic output ```{r} ### Plot cost-effectiveness frontier with probabilistic output plot_icers(df_cea_psa, label = "all", txtsize = txtsize) + expand_limits(x = max(df_cea_psa$QALYs) + 0.1) + theme(legend.position = c(0.8, 0.3)) ``` ## 11.3.4 Cost-effectiveness acceptability curves (CEACs) and frontier (CEAF) This figure shows the Cost-effectiveness acceptability curves. The CEAC shows the probability that a specific strategy is cost-effective (based on all results from the PSA iterations). The cost-effectiveness frontier consist of the set of points corresponding to treatment alternatives that are considered cost-effective at different values of the WTP threshold. ```{r} ### Cost-effectiveness acceptability curves (CEACs) and frontier (CEAF) ceac_obj <- ceac(wtp = v_wtp, psa = l_psa) # Regions of highest probability of cost-effectiveness for each strategy summary(ceac_obj) # CEAC & CEAF plot gg_ceac <- plot_ceac(ceac_obj, txtsize = txtsize, xlim = c(0, NA), n_x_ticks = 14) + ggthemes::scale_color_colorblind() + ggthemes::scale_fill_colorblind() + theme(legend.position = c(0.8, 0.48)) gg_ceac ``` ## 11.3.5 Expected Loss Curves (ELCs) While the CEAC and efficiency frontier only show the probability of being cost effective, it says nothing about how good or bad a strategy is compared to the other strategy in case it is not cost-effective. The Expected Loss Curves combine the probability of being cost-effective with its consequences. It account for the probability that a strategy is not cost-effective and how drastic the consequences are in those scenarios. ```{r} ### Expected Loss Curves (ELCs) elc_obj <- calc_exp_loss(wtp = v_wtp, psa = l_psa) elc_obj # ELC plot gg_elc <- plot_exp_loss(elc_obj, log_y = FALSE, txtsize = txtsize, xlim = c(0, NA), n_x_ticks = 14, col = "full") + ggthemes::scale_color_colorblind() + ggthemes::scale_fill_colorblind() + # geom_point(aes(shape = as.name("Strategy"))) + scale_y_continuous("Expected Loss (Thousand $)", breaks = number_ticks(10), labels = function(x) x/1000) + theme(legend.position = c(0.4, 0.7),) gg_elc ``` ## STAR = 11.3.6 Expected value of perfect information (EVPI) In health economics, decisions are often made under uncertainty, such a the uncertainty around the implementation of strategy A and B. A Value of Information (VOI) analysis helps us assess whether it is worthwhile to reduce that uncertainty through further research before making the decision. Specifically, the VOI quantifies how much value we could gain if we had perfect information to make the best possible decision in every situation. In the analysis, the difference between the expected value under current information and the expected value with perfect information represents the Expected Value of Perfect Information (EVPI). This value tells us the maximum amount we should be willing to invest in additional data collection or studies to eliminate uncertainty. A high EVPI suggests that current uncertainty may lead to sub-optimal decisions, making further research potentially cost-effective. Conversely, a low EVPI implies that our current evidence is sufficient and that additional research may not be necessary. Therefore, VOI analysis supports evidence-based decision-making by helping prioritize research efforts where they are most needed. ```{r} ### Expected value of perfect information (EVPI) evpi <- calc_evpi(wtp = v_wtp, psa = l_psa) # EVPI plot gg_evpi <- plot_evpi(evpi, effect_units = "QALY", txtsize = txtsize, xlim = c(0, NA), n_x_ticks = 14) + scale_y_continuous("EVPI (Thousand $)", breaks = number_ticks(10), labels = function(x) x/1000) gg_evpi ``` We kindly request you to add the following Acknowledgement paragraph to your further work where DARTH code formed the basis. We also like to remind you that you can add other sources of reference to this paragraph to acknowledge code you got from others. # Acknowlegdement For this work we made use of the template developed by the Decision Analysis in R for Technologies in Health (DARTH) workgroup: . The notation of our code is based on the following provided framework and coding convention: Alarid-Escudero, F., Krijkamp, E., Pechlivanoglou, P. et al. A Need for Change! A Coding Framework for Improving Transparency in Decision Modeling. PharmacoEconomics 37, 1329–1339 (2019). . - Alarid-Escudero F, Krijkamp EM, Enns EA, Yang A, Hunink MGM Pechlivanoglou P, Jalal H. An Introductory Tutorial on Cohort State-Transition Models in R Using a Cost-Effectiveness Analysis Example. Medical Decision Making, 2023; 43(1). (Epub). - Alarid-Escudero F, Krijkamp EM, Enns EA, Yang A, Hunink MGM Pechlivanoglou P, Jalal H. A Tutorial on Time-Dependent Cohort State-Transition Models in R using a Cost-Effectiveness Analysis Example. Medical Decision Making, 2023; 43(1). Other work from DARTH can be found on the website: # Copyright for assignment work Copyright 2017, THE HOSPITAL FOR SICK CHILDREN AND THE COLLABORATING INSTITUTIONS. All rights reserved in Canada, the United States and worldwide. Copyright, trademarks, trade names and any and all associated intellectual property are exclusively owned by THE HOSPITAL FOR Sick CHILDREN and the collaborating institutions. These materials may be used, reproduced, modified, distributed and adapted with proper attribution.