R script for a Markov model

Excel spreadsheet is probably the most popular medium for implementing Markov models. I learned the approach from a step-by-step guide by Damon Douglas (published through Knol, a service from Google and now discontinued). I tried to translate the guide into a script for R. Any comments or suggestions for improvement are welcome.

horizon<-10 # model horizon (years, days, whatever)
n=1000 # population size for a population matrix
sim_n<-1000 # number of simulations
n_strat<-2 # number of strategies
n_states<-3 # number of states

# create a table for simulation results;
sim_table<-rep(c(0,0,0,0), times=sim_n)
dim(sim_table)<-c(sim_n,2,n_strat)

#start simulation;
z<-1
for (z in 1:sim_n) {

#STRATEGY in cycles are defined via "t"
#The originating states in cycles are defined via "j" - columns
#The destination states in cycles are defined via "i" - rows

# creating a table for transition probabilities
# creating a supplementary table for standard deviations

trmx<-rep(0, times=n_states*n_states*n_strat)
sdmx<-rep(0, times=n_states*n_states*n_strat)
dim(trmx)<-c(n_states,n_states,n_strat)
dim(sdmx)<-c(n_states,n_states,n_strat)

# state in a column is an originating state; state in a row is a destination state;
# e.g. [3,2,1] means "transitional probability of moving from the state 2 to the state 3
# STRATEGY A
prob_A<-c(0.7, 0.2, 0.1, 0, 0.5, 0.5, 0, 0,1)
sd_A<- c(0.1,0.05,0.003,0, 0.1, 0.1, 0, 0,0)
# STRATEGY B
prob_B<-c(0.5, 0.35, 0.15, 0, 0.4, 0.6, 0, 0,1)
sd_B<-c(0.1,0.05,0.003,0, 0.1, 0.1, 0, 0,0)

dim(prob_A)<-c(n_states, n_states)
dim(sd_A)<-c(n_states, n_states)
dim(prob_B)<-c(n_states, n_states)
dim(sd_B)<-c(n_states, n_states)

trmx[,,1]<-prob_A
trmx[,,2]<-prob_B
sdmx[,,1]<-sd_A
sdmx[,,2]<-sd_B

##########################################
#RANDOMIZING PROBABILITIES IN THE MATRIX
##########################################

i<-1
j<-1
t<-1

for (t in 1:length(trmx[i,j,])) {
i<-1
j<-1
for (j in 1:length(trmx[i,,t])) {
i<-1
for (i in 1:length(trmx[,j,t])) {
a<-trmx[i,j,t]
b<-sdmx[i,j,t]
if (b>0) {
alf_rand<-(a^2*(1-a)/b^2 - a)
beta_rand<-(alf_rand*(1-a)/a)
trmx[i,j,t]<-rbeta(1,alf_rand,beta_rand)
}
else {
trmx[i,j,t]<-trmx[i,j,t]
}}}}

##################################
# ALL TRANSITIONAL PROBABILITIES
# IN EVERY COLUMN MUST SUM TO 1
##################################

i<-1
j<-1
t<-1
for (t in 1:length(trmx[i,j,])){
j<-1
for (j in 1:length(trmx[i,,t])) {
sum_j<-sum(trmx[,j,t])
i<-1
for (i in 1:length(trmx[,j,t])) {
trmx[i,j,t]<-trmx[i,j,t]/sum_j
}}}

##############################################################
# POPULATION TABLE FOR a 3 STATE MODEL
##############################################################

# the first year is for a population at start;
# this state could also be called "year 0"

poptab<-rep(c(n,0,0,rep(0, times=(horizon-1)*n_states)), times=n_strat)
dim(poptab)<-c(horizon, n_states, n_strat)

# a short explanation is worth being here;
# number of people in state 1 is calculated as follows:
# State 1, year 2 = state1_year1*prob[1->1] + state2_year1*prob[2->1] + state3_year1*prob[3->1]
# State 2, year 2 = state1_year1*prob[1->2] + state2_year1*prob[2->2] + state3_year1*prob[3->2]
# State 3, year 2 = state1_year1*prob[1->3] + state2_year1*prob[2->3] + state3_year1*prob[3->3]

t<-1
y<-1
j<-1

for (t in 1:length(poptab[y,j,])){
y<-1
for (y in 2:horizon) {
j<-1
for (j in 1:length(poptab[y,,t])){
poptab[y,j,t]<-crossprod(poptab[y-1,,t],trmx[j,,t])
}}}

#####################
# COST CALCULATIONS
#####################

# Cost of being in a state per patient across the horizon
# No discounting is implemented yet
state_cost<-c(rep(200, times=horizon), rep(1000, times=horizon), rep(0, times=horizon), rep(95, times=horizon), rep(1000, times=horizon), rep(0, times=horizon))
dim(state_cost)<-c(horizon,n_states,n_strat)
cost_table<-rep(0, times=(horizon*n_states*n_strat))
dim(cost_table)<-c(horizon,n_states,n_strat)

i<-1
j<-1
t<-1
for (t in 1:length(poptab[i,j,])){
j<-1
for (j in 1:length(poptab[i,,t])) {
i<-1
for (i in 1:length(poptab[,j,t])) {
cost_table[i,j,t]<-poptab[i,j,t]*state_cost[i,j,t]
}}}

#####################
# QALY CALCULATIONS
#####################
# no discounting on QALY
# Quality of life is 0.8, 0.5, and 0 for State 1, 2, and 3
state_qual<-rep(c(rep(0.8, times=horizon),rep(0.5, times=horizon), rep(0, times=horizon)), times=n_strat)
dim(state_qual)<-c(horizon,n_states,n_strat)
qual_table<-rep(0, times=(horizon*n_states*n_strat))
dim(qual_table)<-c(horizon,n_states,n_strat)

i<-1
j<-1
t<-1
for (t in 1:length(poptab[i,j,])){
j<-1
for (j in 1:length(poptab[i,,t])) {
i<-1
for (i in 1:length(poptab[,j,t])) {
qual_table[i,j,t]<-poptab[i,j,t]*state_qual[i,j,t]
}}}

###################
i<-1
for (i in 1:n_strat){
sim_table[z,,i]<-c(sum(cost_table[,,i]), sum(qual_table[,,i]))
}
# finish simulations
}

# ICER
ICER<-rep(0, times=sim_n)
i<-1
for (i in 1:sim_n) {
ICER[i]<- (sim_table[i,1,1] - sim_table[i,1,2])/(sim_table[i,2,1] - sim_table[i,2,2])
}
plot(ICER)

# plot(ecdf(ICER))

# C/E Acceptability Curve
CEAC<-rep(0,times=sim_n*2)
dim(CEAC)<-c(sim_n,2)
colnames(CEAC)<-c("Incremental effectiveness", "Incremental cost")
i<-1
for (i in 1:sim_n) {
CEAC[i,1]<- (sim_table[i,2,1] - sim_table[i,2,2])
CEAC[i,2]<- (sim_table[i,1,1] - sim_table[i,1,2])
}
plot(CEAC)