Home page

This code is developed by Joke Bilcke and Lander Willem for illustration and teaching purposes only. For further information, questions or comments please contact us.

Goal

Define the uncertainty distributions the following parameters of our flu model:

Set working directory

The first step when creating a new R file is to specify a working directory. You need to tell the R session where you want it to work. To do this you use the command to “set the working directory” or setwd().

You can label this working directory as ‘home’, and have other locations stored for plotting or for getting data. This location should be an existing folder on your computer. For example

home <- "~/Documents/Modelling_intro/" ## on a OS X
home <- "C:\\Documents\\Modelling_intro\\" ## on windows
setwd(home)

Set number of samples

nsamples = 10000

FYI: random number seed

set.seed(5698) # use this to get exactly the same results each time you run the code (e.g. for code checking)
sample(10)
##  [1]  5  1  4 10  3  6  8  9  2  7
sample(10)
##  [1]  9  2 10  3  8  5  6  4  7  1
# start over...
set.seed(5698) # use this to get exactly the same results each time you run the code (e.g. for code checking)
sample(10)
##  [1]  5  1  4 10  3  6  8  9  2  7

note: use set.seed() only once per script!

Normal distributions

Let’s first explore the parameter distribution if we assume the sample size large enough for the central theorem and there are no restrictions: sample from a normal distribution to include parameter uncertainty for the mean.

set.seed(5698) # use this to get exactly the same results each time you run the code (e.g. for code checking)

# if probability (p) and sample size (n) given: se = sqrt(p*(1-p)/n
force_of_infection_norm      = rnorm(nsamples,0.04, sqrt(0.04 * (1-0.04) / 100))    
p_hospital_norm   = rnorm(nsamples,0.1, sqrt(0.1 * (1-0.1) / 100))        
case_fatality_ratio_norm      = rnorm(nsamples,0.03, sqrt(0.03 * (1-0.03) / 1000)) 

# if mean and se given (and n unknown)
life_expectancy_norm = rnorm(nsamples,12,2) 

# if CI given: se = (CI_high - CI_low)/(2*1.96)
vaccine_efficacy_norm       = rnorm(nsamples,0.6, (0.7 - 0.5) / (2 * 1.96))         

# if mean and se given
adm_cost_norm       = rnorm(nsamples, 50, 4)        
cost_hosp_norm      = rnorm(nsamples, 2000, 100)
cost_non_hosp_norm  = rnorm(nsamples, 100, 5)   
QALYloss_hosp_case_norm      = rnorm(nsamples, 0.018, 0.0018) 
QALYloss_non_hosp_case_norm  = rnorm(nsamples, 0.0082, 0.0018)

Other distributions

For some parameters, the sample size is rather small or they should be strictly positive, hence other distribution might be preferred. Let’s first define some help functions to facilitate sampling using a given mean, se, proportion, samples size or confidence interval:

sample_gamma <- function(nsamples, mean, se){
  rgamma(nsamples, shape = mean^2 / se^2, rate = mean / se^2)
}

sample_beta_event <- function(nsamples, events, sample_size){
  return(rbeta(nsamples, shape1 = events , shape2 = sample_size - events))
}

sample_beta_mean <- function(nsamples, mean, se){
  alpha = mean^2 * (1 - mean) / se^2 
  beta  = alpha * (1 - mean) / mean
  return(rbeta(nsamples, shape1 = alpha, shape2 = beta))
}

Let’s sample from an adjusted distribution:

# probability
force_of_infection_beta    = sample_beta_event(nsamples,4,100) # = rbeta(nsamples,4,96)

...

Explore density plots

To check if sample size is sufficiently large to use normal distribution:

plot(density(force_of_infection_norm), main = 'force of infection', ylim = c(0,25))
lines(density(force_of_infection_beta),col = 2)
abline(v = mean(force_of_infection_norm), lty = 2)
grid()
legend('topright',
      c('Norm()','Beta()', 'mean'),
      col = c(1,2,1),
      lty=c(1,1,2),
      lwd=1)

Summary

# over to you...