0% found this document useful (0 votes)

3 views37 pages

Patient-Level Simulation Model Code

Q: Describe the model's approach for calculating additional costs and QALYs accrued from successive events.

The model employs functions like `AddOngoing` and `AddInstant` to calculate additional costs and QALYs. These functions utilize local discount rates and event times, applying continuous and discrete discounting methods to reflect the economic impact of health events over time, ensuring precise economic evaluations by considering temporal effects .

Q: What are the implications of setting individual level patient details to be stored and output in the simulation, and when should this feature be enabled?

Storing individual level patient details increases memory usage and may slow the simulation due to the substantial data being generated, especially with a large number of patients. This feature should primarily be enabled when specific, detailed analyses per patient is required for further insight or debugging purposes, ensuring comprehensive patient-level understanding rather than summary statistics .

Q: What is the model's method of integrating intervention costs into the patient's overall economic evaluation?

The model integrates intervention costs by including them as additional ongoing costs during event processing. For instance, when an event occurs in an intervention, the `AddOngoing` function adds predetermined intervention costs (`cost.int`) alongside ongoing quality of life adjustments, reflecting the financial impact of interventions throughout the patient's simulation timeline .

Q: In the simulation model, what mechanisms are used to enforce the stochastic nature of patient outcomes, particularly regarding hip fracture mortality?

The script utilizes random number generation to simulate variability in patient outcomes. For instance, it uses `runif(1)` to generate a uniform random number to determine if a patient will die from a hip fracture, comparing it against a threshold of `mortprobhip = 0.05`, ensuring that death occurs with a 5% probability .

Q: How does the simulation model's approach affect the interpretation of QALYs and costs regarding healthcare interventions?

The simulation model provides insights into the cost-effectiveness of healthcare interventions by simulating patient outcomes and calculating total and differential costs and QALYs between intervention and no-intervention scenarios. This approach allows stakeholders to evaluate the health benefits relative to economic costs, giving a comprehensive outlook on the intervention's value proposition .

Q: How does the model account for the difference in timing of event occurrences between intervention and non-intervention cases?

The model modifies the timing of events differently for intervention and non-intervention cases. For interventions, the time to the first vertebral fracture is doubled using `vf1 <- 2 * vf1` and similar adjustments are made for hip fracture events, indicating delayed onset compared to non-intervention . This manipulation illustrates how an intervention could potentially delay adverse health outcomes.

Q: What role do comments play in the main script for running the patient-level simulation model, and why are they considered important?

Comments within the main script serve as annotations that explain what the code is doing at each step. They are considered important because they promote transparency and help avoid and identify errors in the model code development .

Q: How does the script ensure that variables do not accidentally carry data from previous simulations?

The script uses the line `rm(list=ls())` at the beginning, which removes any existing data objects in the workspace environment. This ensures that variables do not retain data from previous simulations, thus preventing contamination of results .

Q: How do the codes guarantee correctness and reliability of simulation results during debugging?

Debugging is facilitated by printing out the next event and other variables during simulation if `dbg` is set to TRUE. This output helps detect and correct errors by verifying if the event logic aligns with expected outcomes, such as preventing illogical sequences like fractures occurring posthumously .

Q: Explain the reasoning behind using different utility multipliers for hip and vertebral fractures in the simulation model.

The use of different utility multipliers for hip fractures (0.75) and vertebral fractures (0.90) in the simulation model reflects the varying impacts these conditions have on a patient's quality of life. A hip fracture generally has a more significant negative effect on utility due to increased morbidity compared to a vertebral fracture, hence a lower utility multiplier .

The document provides R code for a cost-effectiveness modeling simulation using patient-level data, specifically focusing on interventions for hip and vertebral fractures. It includes a main script for running the simulation, initializing event lists, and calculating costs and quality-adjusted life years (QALYs) for both intervention and non-intervention groups. Additionally, it outlines how to export results to Excel for further analysis.

Uploaded by

Ricardo Pietrobon

We take content rights seriously. If you suspect this is your content, claim it here.

Available Formats

Download as PDF, TXT or read online on Scribd

0% found this document useful (0 votes)

3 views37 pages

Patient-Level Simulation Model Code

Uploaded by

Ricardo Pietrobon

We take content rights seriously. If you suspect this is your content, claim it here.

Available Formats

Download as PDF, TXT or read online on Scribd

TSD15: Cost‐effectiveness modelling using patient‐level simulation

Appendix B: R code for DES example model

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Appendix A1: Main script for running the model

#NOTES

# Main script for running model and outputting results.

# Text to the right of the "#" sign are comments and are not processed.

# Comments are an essential part of model code and should be used to describe what the code is doing.

# This promotes transparency and helps avoid and identify errors.

###############################################################################

#Start of script###############################################################

###############################################################################

# remove previous contents

rm(list=ls())
# Set the working directory

# NOTE: Use the / rather than \\ symbol

# NOTE: Put the working directory in "" symbols (to make it a character vector)

setwd("G:/ DSU PLS TSD R model")

# load in DES functions in the file DSU_DesFunctions.R

# NOTE: Must be in the same directory

source("DSU_DesFunctions.R")

################################################################################

# GLOBAL VARIABLES #############################################################

################################################################################
# Global variables ‐ placed in base environment

[Link] <‐ 0.7 # initial utility

drc <‐ 0.035 # discount rate for costs

drq <‐ 0.035 # discount rate for QALYs

#npats <‐ 50000 # number of patients

npats <‐ 500 # use this if running through line by line

# utility multipliers for hip fractures and vert fractures

[Link] <‐ 0.75

[Link] <‐ 0.90

[Link] <‐ 7000

[Link] <‐ 3000

[Link] <‐ 500 # annual cost of intervention

mortprobhip <‐ 0.05 # probability of dying from a hip fracture

###################################################################################

#set options########################################################################

###################################################################################

# Debug mode flag: this will print out the event list and accrued costs and qalys

# as the simulation is being run

dbg <‐ T

# Flag for whether individual level patient details should be stored and output

ind <‐ T

if (ind==T){

PatData <‐ vector("list", length=npats) # empty list with 50000 elements

}
##################################################################################

# Run the simulation##############################################################

###################################################################################

# to fix simulation results ‐ use [Link]. Useful for debugging as removes random variation

[Link](1) # or any integer

RunSim() # run simulation

##################################################################################

#examine results ##################################################################

###################################################################################

# To get summary stats

[Link] # total (discounted) costs for intervention arm

[Link] # total (discounted) qalys for intervetnion arm

[Link] # total (discounted) costs for no intervention arm

[Link] # total (discounted) qalys for no intervention arm

[Link] # Total difference in costs between intervention and no intervention pairs

[Link] # Total difference in qalys between intervention and no intervention pairs

# means

[Link]/npats

# With ind mode: lapply function to count the number of

# events
# Number of deaths not relating to hip fractures in intervention arm

sum(sapply(PatData, function(x) x$int$deathother))

# Number of deaths from hip fractures in intervention arm

sum(sapply(PatData, function(x) x$int$deathhip))

# Total number of hip fractures in intervention arm

sum(sapply(PatData, function(x) x$int$nhip))

# Total number of vert fractures in intervention arm

sum(sapply(PatData, function(x) x$int$nvert))

# For non intervention arm equivalents of above, change 'int' to 'noint'

sum(sapply(PatData, function(x) x$noint$deathother))

sum(sapply(PatData, function(x) x$noint$deathhip))

sum(sapply(PatData, function(x) x$noint$nhip))

sum(sapply(PatData, function(x) x$noint$nvert))

# To access individual patient information for a single patient

# Say the 37th patient...

PatData[[37]]$noint # No intervention copy of the patient

PatData[[37]]$int # Intervention copy of the same patient

# To access the event list for a particular patient

PatData[[37]]$int$evtlist

PatData[[37]]$noint$evtlist

##################################################################################

#output results to excel###########################################################

###################################################################################
# Create object for exporting to Excel

ExcelData <‐ [Link](

ptnum=1:npats,

itthip=rep(NA, npats),

ittvert1=rep(NA, npats),

ittvert2=rep(NA, npats),

ittdeath=rep(NA, npats),

ihipcount=rep(NA, npats),

ivertcount=rep(NA, npats),

idthhip=rep(NA,npats),

idthall=rep(NA, npats),

icosts=rep(NA, npats),

iqalys=rep(NA, npats),

ntthip=rep(NA, npats),

nttvert1=rep(NA, npats),

nttvert2=rep(NA, npats),

nttdeath=rep(NA, npats),
nhipcount=rep(NA, npats),

nvertcount=rep(NA, npats),

ndthhip=rep(NA,npats),

ndthall=rep(NA, npats),

ncosts=rep(NA, npats),

nqalys=rep(NA, npats)

for (i in 1:npats){

itthip <‐ PatData[[i]]$int$evtlist$evttime[PatData[[i]]$int$evtlist$evtname=="hfrac"]

ntthip <‐ PatData[[i]]$noint$evtlist$evttime[PatData[[i]]$noint$evtlist$evtname=="hfrac"]

ittvert <‐ PatData[[i]]$int$evtlist$evttime[PatData[[i]]$int$evtlist$evtname%in%c("vfrac1","vfrac2")]

nttvert <‐ PatData[[i]]$noint$evtlist$evttime[PatData[[i]]$noint$evtlist$evtname%in%c("vfrac1","vfrac2")]

ittdeath <‐ PatData[[i]]$int$evtlist$evttime[PatData[[i]]$int$evtlist$evtname=="death"]

nttdeath <‐ PatData[[i]]$noint$evtlist$evttime[PatData[[i]]$noint$evtlist$evtname=="death"]

ExcelData[i,"itthip"] <‐ itthip

ExcelData[i,"ntthip"] <‐ ntthip

ExcelData[i,c("ittvert1", "ittvert2")] <‐ ittvert

ExcelData[i,c("nttvert1", "nttvert2")] <‐ nttvert

ExcelData[i,"ittdeath"] <‐ ittdeath

ExcelData[i,"nttdeath"] <‐ nttdeath

ExcelData[i,"ihipcount"] <‐ PatData[[i]]$int$nhip

ExcelData[i, "ivertcount"] <‐ PatData[[i]]$int$nvert

ExcelData[i,"nhipcount"] <‐ PatData[[i]]$noint$nhip

ExcelData[i, "nvertcount"] <‐ PatData[[i]]$noint$nvert

ExcelData[i, "idthhip"] <‐ PatData[[i]]$int$deathhip

ExcelData[i, "idthall"] <‐ PatData[[i]]$int$deathother

ExcelData[i, "ndthhip"] <‐ PatData[[i]]$noint$deathhip

ExcelData[i, "ndthall"] <‐ PatData[[i]]$noint$deathother

ExcelData[i, "icosts"] <‐ PatData[[i]]$int$thscosts

ExcelData[i, "iqalys"] <‐ PatData[[i]]$int$thsqalys

ExcelData[i, "ncosts"] <‐ PatData[[i]]$noint$thscosts

ExcelData[i, "nqalys"] <‐ PatData[[i]]$noint$thsqalys

[Link](ExcelData, file="DSU_PatData.csv")

+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

Appendix A2: Functions called by the main script

# Functions called by DSU_DesScript.R

# Text to the right of the "#" sign are comments and are not processed.

# Comments are an essential part of model code and should be used to describe what the code is doing.
# This promotes transparency and helps avoid and identify errors.

# Initialise the event list for the no intervention group

# up to two vfracs could occur

[Link] <‐ function(){

vfracs <‐ rweibull(n=2, shape=2, scale=8) # This code simulates a vector of 2 independent Weibull variates, indicating sojourn times to vertebral fracture

hfrac <‐ rweibull(n=1, shape=4, scale=10) # Similarly, this code simulates a time to hip fracture

death <‐ rnorm(n=1, mean=12, sd=3) # As above, but now time to death using normal distribution

output <‐ [Link](

evtname=c( # This code assigns a vector of names stored as strings

"hfrac", # This allows the user to refer to the corresponding event times by name

"vfrac1", # It is good practice to assign names to elements of the model

"vfrac2",

"death"

evttime=c( # This code assigns the event times to a vector, which can be referred to later in the model
hfrac,

vfracs[1], # Note that using square brackets [] selects values from a vector, in this case the "vfracs" vector

vfracs[1] + vfracs[2],

max(0, rnorm(n=1, mean=12, sd=3)) # the normal distribution used above has non‐zero probability of returning a negative value

) # Using the max function prevents this from occuring

# Sort event list by evttime

output <‐ output[order(output$evttime),]

if (ind==T){

[Link]$noint$evtlist <<‐ output

return(output)

}
# [Link] : input is output from [Link]

[Link] <‐ function(input){

output <‐ input

# find the hfrac event and double the time to it

thisrow <‐ which(output$evtname=="hfrac") # note reference to vector of names

output$evttime[thisrow] <‐ 2 * output$evttime[thisrow] # Reassign to output vector

# find the sojourn times for vfrac events and doubles the time to the first one.

thisrow <‐ which(output$evtname=="vfrac1")

vf1 <‐ output$evttime[thisrow]

thisrow <‐ which(output$evtname=="vfrac2")

vf2 <‐ output$evttime[thisrow] ‐ vf1 # This returns the sojourn time rather than the raw event time

vf1 <‐ 2 * vf1 # Time to first vfrac is doubled

vf2 <‐ vf1 + vf2 # Time to second vfrac FROM time to first vfrac is NOT doubled, it is simply added to the updated time to first vert frac
output$evttime[which(output$evtname=="vfrac1")] <‐ vf1 # Assign these fracture times to output vector

output$evttime[which(output$evtname=="vfrac2")] <‐ vf2

#sort event list by evttime

output <‐ output[order(output$evttime),]

if (ind==T){

[Link]$int$evtlist <<‐ output # This code stores individual patient output if this option has been chosen (see DesScript code)

} # Useful for debugging and validation, but increases memory requirements and may slow computation

return(output)

# AddOngoing: Calculate additional qalys and costs accrued from previous to current event

# Inputs:

# lcldrq : local version of discount rate for qalys

# lcldrc : local version of discount rate for costs

# lclprvtime: time of previous event

# lclcurtime: time of current event

# lclvalq : fixed value of qalys

# lclvalc : fixed value of costs

InstantDiscount<‐function(rate){ # For compound continuous discounting ‐ use INSTANTANEOUS rate.

log(1+rate)

AddOngoing <‐ function(lcldrq=0.035, lcldrc=0.035, lclprvtime, lclcurtime, lclvalq, lclvalc){

Instantdrq <‐ InstantDiscount(lcldrq)

Instantdrc <‐ InstantDiscount(lcldrc)

# calculate additional qalys

addqalys <‐ ((lclvalq)/(0 ‐ Instantdrq)) * (exp(lclcurtime * ( 0 ‐ Instantdrq)) ‐ exp(lclprvtime * (0 ‐ Instantdrq)))

# calculate additional costs

addcosts <‐ ((lclvalc)/(0 ‐ Instantdrc)) * (exp(lclcurtime * ( 0 ‐ Instantdrc)) ‐ exp(lclprvtime * (0 ‐ Instantdrc)))

# combine additional costs and additional qalys in a list

output <‐ list(addqalys=addqalys, addcosts=addcosts)

return(output)

#######################

# Add Instantaneous costs and qalys

# Inputs:

# lcldrq : local version of discount rate for qalys

# lcldrc : local version of discount rate for costs

# lclcurtime: time of current event

# lclvalq : fixed value of qalys

# lclvalc : fixed value of costs

AddInstant <‐ function(lcldrq=0.035, lcldrc=0.035, lclcurtime, lclvalq, lclvalc){

addinstqalys <‐ lclvalq * ((1+lcldrq)^(‐lclcurtime)) # Note use of DISCRETE TIME discounting for instantaneous costs and benefits

addinstcosts <‐ lclvalc * ((1+lcldrc)^(‐lclcurtime))

# combine additional costs and additional qalys in a list

output <‐ list(addinstqalys=addinstqalys, addinstcosts=addinstcosts)

return(output)

GetNxtEvt <‐ function(intervention=F){ # This function identifies which event is to be processed next for each patient, depending on intervention

if (intervention==F){ # It loops through events until death has occured, either due to hip fracture or other causes

# Do the following if no intervention

# if there is at least one event still to process

if (dim([Link])[1]> 0){

nextevt <‐ [Link]$evtname[1]

nextevttime <‐ [Link]$evttime[1]

# remove next event as it has now been processed

# Debugging line: If debugging is enabled (see DesScript.R) then this will print out the next event.

# This helps detect errors in model logic (e.g. fracture occurs after death).

if (dbg==T){

print([Link])

[Link] <<‐ [Link][‐1,]

output <‐ list(evt=nextevt, evttime=nextevttime)

} else {

output <‐ NULL

} else {

# do the following if intervention

if (dim([Link])[1]> 0){

nextevt <‐ [Link]$evtname[1]

nextevttime <‐ [Link]$evttime[1]

# remove next event as it has now been processed

# Debugging line:

if (dbg==T){

print([Link])

[Link] <<‐ [Link][‐1,]

output <‐ list(evt=nextevt, evttime=nextevttime)

} else {

output <‐ NULL

return(output)

# ReactEvt : react to the next event

# thisevt : a two element list containing the output from GetNextEvt

# $evt : event

# $evttime : event time

# intervention: boolean flag : TRUE or FALSE

ReactEvt <‐ function(thisevt, intervention){ # This function processes the next event (as identified in the GetNextEvt function)

evt <‐ thisevt$evt # Identify event type

prevtime <‐ curtime # Identify time of previous event

curtime <<‐ thisevt$evttime # Identify time of next event

if (intervention==F){

# No intervention logic

if (evt=="death"){

if(ind==T){

[Link]$noint$deathother <<‐ 1 # This indicates that the patient is now dead ‐ no further events

# create variable to additional ongoing costs and qalys

additionals <‐ AddOngoing(lclprvtime=prevtime, lclcurtime=curtime, lclvalq=utilmlt, lclvalc=0)

instadditionals <‐ AddInstant(lclcurtime=curtime, lclvalq=0, lclvalc=0)

thsqalys <<‐ thsqalys + additionals$addqalys + instadditionals$addinstqalys

thscosts <<‐ thscosts + additionals$addcosts + instadditionals$addinstcosts

curtime <<‐ Inf # Set current time to Infinity so patient level loop stops

} else if (evt %in% c("vfrac1","vfrac2")){

if (ind==T){ # Code chunks like this are processed only if individual patient info is being saved (ind=TRUE)

[Link]$noint$nvert <<‐ [Link]$noint$nvert + 1 # This records the number of vertebral fractures experienced AFTER event occurs

if (prvvert==F){

# ===================================

# [Logic if no previous vert fracture]

# ===================================

additionals <‐ AddOngoing(lclprvtime=prevtime, lclcurtime=curtime, lclvalq=utilmlt, lclvalc=0)

instadditionals <‐ AddInstant( lclcurtime=curtime, lclvalq=0, lclvalc=[Link])

thsqalys <<‐ thsqalys + additionals$addqalys + instadditionals$addinstqalys

thscosts <<‐ thscosts + additionals$addcosts + instadditionals$addinstcosts

# utility multiplier at previous

utilmlt <<‐ utilmlt * [Link] # Utility is set to previous utility (utilmlt) multiplied by the multiplier for a vertebral fracture

prvvert <<‐ T

} else {

# ====================================

# [Logic if a previous vert fracture]

# ====================================

additionals <‐ AddOngoing(lclprvtime=prevtime, lclcurtime=curtime, lclvalq=utilmlt, lclvalc=0)

instadditionals <‐ AddInstant( lclcurtime=curtime, lclvalq=0, lclvalc=[Link])

thsqalys <<‐ thsqalys + additionals$addqalys + instadditionals$addinstqalys

thscosts <<‐ thscosts + additionals$addcosts + instadditionals$addinstcosts

} else if (evt=="hfrac"){
if (ind==T){

[Link]$noint$nhip <<‐ [Link]$noint$nhip + 1 # Records that patient had a hip fracture

# ====================================

# Hip fracture logic

# ====================================

additionals <‐ AddOngoing(lclprvtime=prevtime, lclcurtime=curtime, lclvalq=utilmlt, lclvalc=0)

instadditionals <‐ AddInstant( lclcurtime=curtime, lclvalq=0, lclvalc=[Link])

thsqalys <<‐ thsqalys + additionals$addqalys + instadditionals$addinstqalys

thscosts <<‐ thscosts + additionals$addcosts + instadditionals$addinstcosts

utilmlt <<‐ utilmlt * [Link] # as with vertebral ‐ multiply previous utility by hip fracture modifier

# Death occurs with 5% probability

patdies <‐ runif(1) < mortprobhip # This code generates a uniform variate and compares it with the probability of hip fracture related mortality

# If the variate is less than the probability (0.05 in this example), "patdies" is set to TRUE
if (patdies) {

if (ind==T){

[Link]$noint$deathhip <<‐ [Link]$noint$deathhip + 1 # This indicates that the patient died of a hip fracture

curtime <<‐ Inf

} else {

# if this is reached, then something has gone wrong

stop("Event type not recognised") # Use debugging to identify why this occurred and fix error as appropriate

} else { #

#==============================================================

# Intervention logic ‐ follows same format as non‐intervention logic

#==============================================================

if (evt=="death"){

if (ind==T){

[Link]$int$deathother <<‐ [Link]$int$deathother + 1

}

# create variable to additional ongoing costs and qalys

additionals <‐ AddOngoing(lclprvtime=prevtime, lclcurtime=curtime, lclvalq=utilmlt, lclvalc=[Link])

instadditionals <‐ AddInstant(lclcurtime=curtime, lclvalq=0, lclvalc=0)

thsqalys <<‐ thsqalys + additionals$addqalys + instadditionals$addinstqalys

thscosts <<‐ thscosts + additionals$addcosts + instadditionals$addinstcosts

curtime <<‐ Inf # Set current time to Infinity so patient level loop stops

# create variable to additional ongoing costs and qalys

} else if (evt %in% c("vfrac1","vfrac2")){

if (ind==T){

[Link]$int$nvert <<‐ [Link]$int$nvert + 1

if (prvvert==F){
# ===================================

# [Logic if no previous vert fracture]

# ===================================

additionals <‐ AddOngoing(lclprvtime=prevtime, lclcurtime=curtime, lclvalq=utilmlt, lclvalc=[Link])

instadditionals <‐ AddInstant(lclcurtime=curtime, lclvalq=0, lclvalc=[Link])

thsqalys <<‐ thsqalys + additionals$addqalys + instadditionals$addinstqalys

thscosts <<‐ thscosts + additionals$addcosts + instadditionals$addinstcosts

# utility multiplier at previous

utilmlt <<‐ utilmlt * [Link]

prvvert <<‐ T

} else {

# ====================================

# [Logic if a previous vert fracture]

# ====================================

additionals <‐ AddOngoing(lclprvtime=prevtime, lclcurtime=curtime, lclvalq=utilmlt, lclvalc=[Link])

instadditionals <‐ AddInstant(lclcurtime=curtime, lclvalq=0, lclvalc=[Link])

thsqalys <<‐ thsqalys + additionals$addqalys + instadditionals$addinstqalys

thscosts <<‐ thscosts + additionals$addcosts + instadditionals$addinstcosts

} else if (evt=="hfrac"){

if (ind==T){

[Link]$int$nhip <<‐ [Link]$int$nhip + 1

additionals <‐ AddOngoing(lclprvtime=prevtime, lclcurtime=curtime, lclvalq=utilmlt, lclvalc=[Link])

instadditionals <‐ AddInstant(lclcurtime=curtime, lclvalq=0, lclvalc=[Link])

thsqalys <<‐ thsqalys + additionals$addqalys + instadditionals$addinstqalys

thscosts <<‐ thscosts + additionals$addcosts + instadditionals$addinstcosts

utilmlt <<‐ utilmlt * [Link]

# Death occurs with 5% probability

patdies <‐ runif(1) < mortprobhip

if (patdies) {

if (ind==T){

[Link]$int$deathhip <<‐ [Link]$int$deathhip + 1

curtime <<‐ Inf

} else {

# if this is reached, then something has gone wrong

stop("Event type not recognised")

# set to return nothing (NULL object) so earlier operations are not returned by accident

return(NULL)

}
# Enclose simulation within a function

RunSim <‐ function(){

# initialise variable of total costs and QALYs (note use of superassignment operator (<<‐))

[Link] <<‐ 0 # total QALYs accrued by all patients ‐ no intervention

[Link] <<‐ 0 # total costs accrued by all patients ‐ no intervention

[Link] <<‐ 0 # total QALYs accrued by all patients ‐ intervention

[Link] <<‐ 0 # total costs accrued by all patients ‐ intervention

[Link] <<‐ 0 # difference in QALYs between intervention and comparator

[Link] <<‐ 0 # difference in costs between intervention and comparator

# Outer loop, repeat for each patient

for (i in 1:npats){

if (ind==T){

[Link] <<‐ list(

int=list(

nvert=0,

nhip=0,

deathhip=0,

deathother=0

noint=list(

nvert=0,

nhip=0,

deathhip=0,

deathother=0

# Debugging line

if (dbg==T){

cat(paste("\n#####\n[", i, "]\n#####\n")) # The cat function prints text to console

print("No Intervention Patient")

# Generate event list ‐ no intervention

[Link] <<‐ [Link]()

# Generate event list ‐ intervention

[Link] <<‐ [Link]([Link])

# For the no intervention patient:

# current time, set as global variable to 0

curtime <<‐ 0

utilmlt <<‐ 0.7 # initially, the patient has a utility of 0.7

prvvert <<‐ F # no previous vfract

# QALYs and costs for this patient

thsqalys <<‐ 0

thscosts <<‐ 0
while(curtime < Inf){

# Get next event, process, repeat

Evt <‐ GetNxtEvt(intervention=F)

if ([Link](Evt)==F){

ReactEvt(Evt, intervention=F)

} else {curtime <<‐ Inf}

if(dbg==T){

print(paste("No Intervention, Qalys:", round(thsqalys, 2), "; cost:", round(thscosts,0) ))

if (ind==T){

[Link]$noint$thsqalys <<‐ thsqalys

[Link]$noint$thscosts <<‐ thscosts

[Link] <<‐ [Link] + thsqalys

[Link] <<‐ [Link] + thscosts

# subtracting costs and qalys as from no intervention arm, and want to know int ‐ noint

[Link] <<‐ [Link] ‐ thsqalys

[Link] <<‐ [Link] ‐ thscosts

# for the intervention patient

if (dbg==T){ print("Intervention Patient")}

# reset curtime

curtime <<‐ 0

utilmlt <<‐ 0.7 # initially, the patient has a utility of 0.7

prvvert <<‐ F # no previous vfrac

# QALYs and costs for this patient

thsqalys <<‐ 0

thscosts <<‐ 0

# if the event list has been emptied

emptylist <<‐ F

while(curtime < Inf){

Evt <‐ GetNxtEvt(intervention=T)

if ([Link](Evt)==F){

ReactEvt(Evt, intervention=T)

} else {curtime <<‐ Inf}

if(dbg==T){

print(paste("Intervention, Qalys:", round(thsqalys, 2), "; cost:", round(thscosts,0) ))

if (ind==T){

[Link]$int$thsqalys <<‐ thsqalys

[Link]$int$thscosts <<‐ thscosts

if (ind==T){

PatData[[i]] <<‐ [Link]

}

[Link] <<‐ [Link] + thsqalys

[Link] <<‐ [Link] + thscosts

# adding costs and qalys as from intervention arm, and no int has already been subtracted

[Link] <<‐ [Link] + thsqalys

[Link] <<‐ [Link] + thscosts

Common questions