Simulations are not only important for statistics but also for a variety of other areas, e.g. Markov chain Monte Carlo, bootstrap or stochastic disease transmission models. All simulations have in common that there is a need to introduce randomness. Typically, simulation require repeated random sampling, random number generation or both.
Therefore, to perform simulations we need to know the following:
How to select random samples?
How to generate random numbers?
How to execute R code repeatedly?
How to track the results from each simulation run?
R's function sample can perform sampling with and without replacement.
sample(1:6, 3) selects 3 random numbers from 1 to 6 without replacement, i.e. each element can occur only once in the sample.
sample(1:6, 3, replace = T) will also select 3 random numbers but with replacement. Therefore, in the sample a certain number might occure more than once. In contrast to sampling without replacement, the sample size can be higher than the number of items to be sampled from, e.g. you can simulate rolling a six-sided die 10 times with:
sample(1:6, 10, replace = T)
Often, we will not sample from a sequence of numbers but from observed data which is done in the same way, e.g. sample 50 values from variable Temp in dataset airquality:
sample(airquality$Temp, 50)
For many applications, you need to select many random samples. Therefore, you need to execute some R code several times. The R function replicate(n_times, any_r_code) will execute 'any_r_code' exactly 'n_times'.
Example: Poker hands
A deck of poker cards has 52 card with 13 different ranks (2 to 10, J, Q, K & ace). If you draw 5 random cards, What is the probabaility of getting 3 or 4 of a kind? To see how many cards of the same kind in our hand we can use the function table. It is a bit tricky to store the entire table from all simulation runs. Therfore, we store only the highest value with max.
Example: Block randomization
In randomized clinical trials, a simple random allocation scheme (flipping a coin for each participant) might by chance result in an unequal number of individuals per trial arm. Block randomization is a commonly used technique to ensure a balanced number of participants in each treatment arm. In a clinical trial with 2 treatment arms (A and B), a block size of 4 will result in 6 possible permutations AABB, ABAB, BABA, ABBA, BAAB, BBAA.
A colleague asks you to generate a random allocation sequnece for a clinical trial with 2 arms, a block size of 4 and a total sample size of 200. Generating 1 random permutation is done by:
sample(c("A","A","B","B"), 4)
All we have to do is to repeat the procedure 50 times.
Example: Bootstrap confidence intervals
Sometimes you would like to estimate the 95% confidence intervals of a parameter for which no suitable formula exists. For example, to estimate the 95% CI for the median is much more complicated compared to the interval for the arithmetic mean.
A good way to solve this problem is to apply bootstrap resampling. A bootstrap sample of a variable or dataset is a sample of the same size of the original data where the observations are sampled with replacement. Bootstrap samples of the numbers 1 to 8 might look as in the table below:
| Sample # | Value 1 | Value 2 | Value 3 | Value 4 | Value 5 | Value 6 | Value 7 | Value 8 |
|---|---|---|---|---|---|---|---|---|
| original | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| 1 | 2 | 3 | 4 | 4 | 4 | 5 | 7 | 7 |
| 2 | 1 | 1 | 3 | 6 | 6 | 6 | 7 | 8 |
| 3 | 3 | 3 | 3 | 5 | 5 | 6 | 6 | 8 |
| 4 | 1 | 3 | 4 | 4 | 4 | 5 | 7 | 8 |
| 5 | 1 | 1 | 1 | 1 | 4 | 6 | 7 | 8 |
Confidence interval estimation via bootstrap resampling requires the following steps:
1. Generate a bootstrap resample, i.e. random sample of your data of the same size with relacement.
2. Estimate the parameter oif interest.
3. Replcate steps 1 & 2 many times.
3. The 2.5% and 97.5% percentiles of all estimates represent the 95% CI.
The dataset airquality (daily air measurements in New York, 1973) has the variable Temp (temperature in Fahrenheit). We would like to estimate the 95% CI of the mean via bootstrap resampling. Note: The one sample t.test estimates the 95% CI at 76.4 - 79.4. Let's see if we get similar results.
Exercise: Birthday paradoxon
You have a group of 23 randomly chosen persons. What is the probability that (at least) 2 people share the same birthday?
We could calculate the probability with the hypergeometric distribution but we don't have the formula at hand and would like to simulate it. To simulate the birthday from 23 persons, we sample randomly 23 numbers from 1 to 365 with replacement.
To check if there are any duplicated numbers in each sample we can wrap sample in the function anyDuplicated which returns 0 if there are no duplicated values. Finally, we replicate this 1000 times and count the 0's.
Don't be surprised: the proportion is roughly 50%. Complete the code below:
Sampling from probability distributions
Rather than sampling from a set of numbers, we might want to get samples from a probability distributions, e.g. from the normal distribution. To get a sample of size 10 from a normal distribution with mean of 5 and standard deviation of 2 we can use the function rnorm as follows:
rnorm(10, mean = 5, sd = 2)
The colleague from the previous chapter visits you again and has another question: A trial should be designed and they assume that in the control group the persons will have a normally distributed outcome with mean 5 and standard deviation 4 whereas in the treatment group people are expected to show an outcome with mean 3 and standard deviation 1. Sample size is fixed at 100 persons (no refusals or loss to follow-up) and randomisation should be done with a ratio of 1:2. You are asked to determine the statistical power if the assumptions are correct.
You know the sample size formula for Student's t-test but you are not sure if the formula holds if standard deviation and sample size per group differs. Therefore, you decide to simulate it.
Note: if you have more than 1 line of code to replicate you can enclose them in curly brackets:
replicate(100, {
#... some R code ...
#... more R code ...
}
)Let's start the simulation to estimate power wich is the proportion of simulation runs with P-values below 0.05.
Similarly, we have the functions rbinom and runif to genrate binomial and uniform distributed random values.
For instance to simulate 10 coin flips use:
rbinom(10, 1, 0.5)
[1] 0 1 1 1 1 0 1 1 0 1The 10 indicates the number of values to generate, the 1 the number of trials (which is 1 for binary data) and the 0.5 the probability of success which is 50% for a coin flip.
Exercise: Binomial distribution
The binomial distribution is discrete but can approximate the normal distribution if B(n, p) n is large and p is close to 50%. Let's see if this is true. To get n large we will not flip a single coin 500 times but we will flip our entire purse (which has 20 coins inside) 500 times. Each time we will record the number of heads of the 20 coins. Finally, will plot a histogram to see the distribution. Complete the code below:
The function replicate executes the same R code several times. Unfortunately, in most situations you would like not to execute exactly the same code several times but run each time a slightly different code.
R, like all programming languages, has loops for this pourpose. The most important loop is the for loop. In loops you specify an iterator (often i) and a sequence. R will execute the code within the loop body in exactly the same way. The only difference is that the iterator gets a new value in each loop - always the next value from the elements in the sequence. The R code has to be enclosed in curly brackets.
for(i in 1:3){
print(i)
}
[1] 1
[1] 2
[1] 3Note: R usually prints all output in the console but not if the code is enclosed in curly brackets. In this case you have to specify print if you would like to display the output.
Examples: Loop over variables
Loops can also be used to loop over variables. The code below can be used to print tables of the first 6 variables of the dataset infert (case control study on risk factors for infertility).
Example: Storing results from each loop
In the introduction was mentioned that it is important to track the results from each simulation run? A good way is to create a placeholder which will be updated in each loop. Let us translate the bootstrap example into a loop:
Exercise: Euler's number
The number e can be calculate in several ways.
One way is
e = 1/0! + 1/1! + 1/2! ... + 1/n!
Where the ! operator represents the factorial. In R we calculate 5! with factorial(5). Complete the code below to use a loop to calculate 1/0! + ... + 1/10!:
Sometimes you wish that parts of your code within a loop will only be executed if a certain condition is true. In the example below the loop is executed 5 times but only if i is greater than 3, the value will be displayed.
for(i in 1:5){
if(i > 3) print(i)
}
[1] 4
[1] 5
Of course, if there is an if than there is an else, too. Assume you would like to loop over all variables of a dataset and you would like to perform one action if the variable has no missing values and another action if the variable has missing values:
To put all together we will replicate a simulation study done by Freedman in 1983.
https://amstat.tandfonline.com/doi/abs/10.1080/00031305.1983.10482729
Researchers often collect hundreds of variables and it is impossible to run a multiple regression with all of them. Often, in a first step variables will be identified which are associated with the outcome (e.g. p < 0.25) and only those will be included in the final multivariable model. It is well known that this approach leads often to spurious results.
Freedman published in 1983 the results from a simulation to show the extend of the problem.
The simulation can be broken down into 3 main steps:
1) A data set was created with 51 variables and 100 observations each. All values were pure noise, i.e. random numbers drawn from the standard normal distribution. The last variable was taken as outcome variable, the other 50 as predictors.
2) 50 bivariable linear regressions were run: outcome and each of the 50 predictors.
3) All predictor variables that were associated with p values < 0.25 were included in the final multiple regression model.
Take a minute to think about how to approach each step.
The first step seems to be easy. We can use rnorm to generate random numbers. If we combine it with matrix. We have almost a dataset:
df <- matrix(rnorm(5100, 0, 1), nrow=100)
We just need to transform it into a true dataset and we need to find out the variable names.
The 2nd step is also solvable. We run 50 bivariable regressions using always the last variable (X51) as outcome and one by one of the 50 others as predictor. A loop can easily perform this operation. There is only one tricky thing: we need to store the p-values. Extracting P-values from regression models is not really intutive but can be do ne as follows:
lin.1 <- lm(outcome ~ predictor)
P <- coef(summarize(lin.1))[2,4]
For the 3rd step we simply keep only those variables with p < 0.25 and to run the multivariable model with the remaining ones.
Example: sample with the probs argument
Sometimes you would like to select a random sample but the elements should have different selection probabilities.
Nate Silver's internet page fivethirtyeight.com had a nice probability puzzle some time ago:
https://fivethirtyeight.com/features/who-wants-to-be-a-riddler-millionaire/
You are in the game show “Who Wants to Be a Millionaire?” and you’ve made it to the final question. Out of the four choices, A, B, C and D, you’re 70 % sure the correct answer is B, and none of the remaining choices looks more plausible than another. You use your 50:50 lifeline which leaves you with two possible answers, one of them correct. Fortunately, B remains an option! How confident are you now that B is the correct answer?
On the first view this looks like a Monty Hall problem but it is in fact not. It can be easily solved using Bayes' theorem but you can also solve it easily via simulations. Note: the function setdiff(set1, set2) returns the elements of set1 which are not in set2, e.g.
setdiff(1:5, 2:4) returns 1 and 5
Example: Nested loops & Bayes' Billiards Problem
Bayes' invented a game in 1763 which Bayesianists love because it is very easy to solve with Bayesian statistics but not at all with a frequentists approach. The rules are as follows:
Alice and Bob are playing a game:
A casino has a pool table that A and B can't see.
The Casino rolls an initial ball onto the table, which comes to rest at a random position.
Than, the Casino rolls additional balls randomly. If the ball rests above the initial ball, Alice scores a point; below the initial ball, Bob scores a point.
The Casino reveals only who won each point.
##############################################
## ##
# Alice's area #
# #
# --- (0) <- initial ball --------------- #
# #
# Bob's area #
# #
# (1) <-1st new ball (Bob +1pt) #
## ##
##############################################
The first with 6 points wins.
Imagine Alice is leading 5 to 3. What is the probability that Bob is winning?
If you are a hero in Bayesian statistics or not, it is possible to answer the question via simulations. In the example below we use a loop within another loop. The implementation of nested loops is possible but the name of the iterator has to be different (i and j in this example).