Probability and simulating data

Sample continuous distribution

runif(n, min=0, max=1)

Use runif() to sample from a continuous uniform distribution.

runif(10, min = 0, max = 1)

##  [1] 0.4629114 0.7160526 0.3703031 0.3297495 0.6740526 0.8923179 0.4092995
##  [8] 0.1899743 0.9876491 0.1418342

Sample discrete distribution

sample(x, size, replace = FALSE, prob = NULL)

Use sample() to sample from a discrete distribution.

Simulate a single roll of a 6-sided die.

sample(6, 1)

## [1] 5

Simulate 10 rolls of a 6-sided die. Set replace to TRUE so each roll is independent. See what happens if you set replace to FALSE.

sample(6, 10, replace = TRUE)

##  [1] 3 5 6 3 2 6 1 3 3 5

You can also use sample to sample from a list of named outcomes.

pet_types <- c("cat", "dog", "ferret", "bird", "fish")
sample(pet_types, 10, replace = TRUE)

##  [1] "cat"    "fish"   "fish"   "ferret" "dog"    "fish"   "cat"   
##  [8] "fish"   "dog"    "fish"

Ferrets are a much less common pet than cats and dogs, so our sample isn’t very realistic. You can set the probabilities of each item in the list with the prob argument.

pet_types <- c("cat", "dog", "ferret", "bird", "fish")
pet_prob <- c(0.3, 0.4, 0.1, 0.1, 0.1)
sample(pet_types, 10, replace = TRUE, prob = pet_prob)

##  [1] "cat"    "dog"    "dog"    "dog"    "ferret" "cat"    "dog"   
##  [8] "dog"    "dog"    "dog"

Sample distribution

rbinom(n, size, prob)

The rbinom function will generate a random binomial distribution.

• n = number of observations
• size = number of trials
• prob = probability of success on each trial

Coin flips are a typical example of a binomial distribution, where we can assign head to 1 and tails to 0.

20 individual coin flips of a fair coin

rbinom(20, 1, 0.5)

##  [1] 1 1 1 1 0 1 0 0 1 1 0 0 1 1 0 0 1 0 1 1

20 individual coin flips of a baised (0.75) coin

rbinom(20, 1, 0.75)

##  [1] 1 1 1 1 1 0 1 1 0 1 1 1 1 0 1 0 1 0 0 1

You can generate the total number of heads in 1 set of 20 coin flips by setting size to 20 and n to 1.

rbinom(1, 20, 0.75)

## [1] 14

You can generate more sets of 20 coin flips by increasing the n.

rbinom(10, 20, 0.5)

##  [1]  7  8 11 11  7 11 10 15  8  8

You should always check your randomly generated data to check that it makes sense. For large samples, it’s easiest to do that graphically. A histogram is usually the best choice for plotting binomial data.

flips <- rbinom(1000, 20, 0.5)

ggplot() +
  geom_histogram(
    aes(flips), 
    binwidth = 1, 
    fill = "white", 
    color = "black"
  )

Run the simulation above several times, noting how the histogram changes. Try changing the values of n, size, and prob.

Exact binomial test

binom.test(x, n, p)

You can test a binomial distribution against a specific probability using the exact binomial test.

• x = the number of successes
• n = the number of trials
• p = hypothesised probability of success

Here we can test a series of 10 coin flips from a fair coin and a biased coin against the hypothesised probability of 0.5 (even odds).

n <- 10
fair_coin <- rbinom(1, n, 0.5)
biased_coin <- rbinom(1, n, 0.6)

binom.test(fair_coin, n, p = 0.5)

## 
##  Exact binomial test
## 
## data:  fair_coin and n
## number of successes = 3, number of trials = 10, p-value = 0.3438
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.06673951 0.65245285
## sample estimates:
## probability of success 
##                    0.3

binom.test(biased_coin, n, p = 0.5)

## 
##  Exact binomial test
## 
## data:  biased_coin and n
## number of successes = 6, number of trials = 10, p-value = 0.7539
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.2623781 0.8784477
## sample estimates:
## probability of success 
##                    0.6

False positives & negatives

The probability that a test concludes the fair coin is biased is called the false positive rate (or Type I Error Rate). The alpha is the false positive rate we accept for a test. This is traditionally set at 0.05, but there are good arguments for setting a different criterion in some circumstances.

The probability that a test concludes the biased coin is fair is called the false negative rate (of Type II Error Rate). The power of a test is 1 minus its false negative rate (i.e., the true positive rate). Power depends on how biased the coin is and how many samples we take.

Sampling function

To estimate these rates, we need to repeat the sampling above many times. A function is ideal for repeating the exact same procedure over and over. Set the arguments of the function to variables that you might want to change. Here, we will want to estimate power for:

• different sample sizes (n)
• different coin biases (bias)
• different hypothesised probabilities (p, defaults to 0.5)

sim_binom_test <- function(n, bias, p = 0.5) {
  coin <- rbinom(1, n, bias)
  btest <- binom.test(coin, n, p)
  
  btest$p.value
}

Once you’ve created your function, test it a few times, changing the values.

sim_binom_test(100, 0.6)

## [1] 0.01203298

Calculate power

Then you can use the replicate() function to run it many times and save all the output values. You can calculate the power of your analysis by checking the proportion of your simulated analyses that have a p-value less than your alpha (the probability of rejecting the null hypothesis when the null hypothesis is true).

my_reps <- replicate(1e4, sim_binom_test(100, 0.6))

mean(my_reps < 0.05)

## [1] 0.4595

1e4 is just scientific notation for a 1 followed by 4 zeros (10000). When youre running simulations, you usually want to run a lot of them and it’s a pain to keep track of whether you’ve typed 5 or 6 zeros (100000 vs 1000000) and this will change your running time by an order of magnitude.

Sample distribution

rnorm(n, mean, sd)

We can simulate a normal distribution of size n if we know the mean and standard deviation (sd). A density plot is usually the best way to visualise this type of data if your n is large.

dv <- rnorm(1e5, 10, 2)

ggplot() +
  geom_density(aes(dv), fill = "white") +
  geom_vline(xintercept = mean(dv), color = "red") +
  geom_vline(xintercept = quantile(dv, .5 - (.6827/2)), color = "darkgreen") +
  geom_vline(xintercept = quantile(dv, .5 + (.6827/2)), color = "darkgreen") +
  geom_vline(xintercept = quantile(dv, .5 - (.9545/2)), color = "blue") +
  geom_vline(xintercept = quantile(dv, .5 + (.9545/2)), color = "blue") +
  geom_vline(xintercept = quantile(dv, .5 - (.9973/2)), color = "purple") +
  geom_vline(xintercept = quantile(dv, .5 + (.9973/2)), color = "purple") +
  scale_x_continuous(
    limits = c(0,20), 
    breaks = seq(0,20)
  )

Run the simulation above several times, noting how the density plot changes. What do the vertical lines represent? Try changing the values of n, mean, and sd.

T-test

t.test(x, y, alternative, mu, paired)

Use a t-test to compare the mean of one distribution to a null hypothesis (one-sample t-test), compare the means of two samples (independent-samples t-test), or compare pairs of values (paired-samples t-test).

You can run a one-sample t-test comparing the mean of your data to mu. Here is a simulated distribution with a mean of 0.5 and an SD of 1, creating an effect size of 0.5 SD when tested against a mu of 0. Run the simulation a few times to see how often the t-test returns a significant p-value (or run it in the shiny app).

sim_norm <- rnorm(100, 0.5, 1)
t.test(sim_norm, mu = 0)

## 
##  One Sample t-test
## 
## data:  sim_norm
## t = 3.9531, df = 99, p-value = 0.0001449
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  0.2048914 0.6178703
## sample estimates:
## mean of x 
## 0.4113809

Run an independent-samples t-test by comparing two lists of values.

a <- rnorm(100, 0.5, 1)
b <- rnorm(100, 0.7, 1)
t_ind <- t.test(a, b, paired = FALSE)
t_ind

## 
##  Welch Two Sample t-test
## 
## data:  a and b
## t = -3.4277, df = 197.94, p-value = 0.0007404
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.6423706 -0.1731738
## sample estimates:
## mean of x mean of y 
## 0.3417206 0.7494927

The paired argument defaults to FALSE, but it’s good practice to always explicitly set it so you are never confused about what type of test you are performing.

Sampling function

We can use the names() function to find out the names of all the t.test parameters and use this to just get one type of data, like the test statistic (e.g., t-value).

names(t_ind)

## [1] "statistic"   "parameter"   "p.value"     "conf.int"    "estimate"   
## [6] "null.value"  "alternative" "method"      "data.name"

t_ind$statistic

##        t 
## -3.42771

Alternatively, use broom::tidy() to convert the output into a tidy table.

broom::tidy(t_ind)

If you want to run the simulation many times and record information each time, first you need to turn your simulation into a function.

sim_t_ind <- function(n, m1, sd1, m2, sd2) {
  v1 <- rnorm(n, m1, sd1)
  v2 <- rnorm(n, m2, sd2)
  t_ind <- t.test(v1, v2, paired = FALSE)
  
  return(t_ind$p.value)
}

Run it a few times to check that it gives you sensible values.

sim_t_ind(100, 0.7, 1, 0.5, 1)

## [1] 0.8879239

Now replicate the simulation 1000 times.

my_reps <- replicate(1e4, sim_t_ind(100, 0.7, 1, 0.5, 1))

alpha <- 0.05
power <- mean(my_reps < alpha)
power

## [1] 0.2882

Run the code above several times. How much does the power value fluctuate? How many replications do you need to run to get a reliable estimate of power?

Compare your power estimate from simluation to a power calculation using power.t.test(). Here, delta is the difference between m1 and m2 above.

power.t.test(n = 100, delta = 0.2, sd = 1, sig.level = alpha, type = "two.sample")

## 
##      Two-sample t test power calculation 
## 
##               n = 100
##           delta = 0.2
##              sd = 1
##       sig.level = 0.05
##           power = 0.2902664
##     alternative = two.sided
## 
## NOTE: n is number in *each* group

You can plot the distribution of p-values.

ggplot() + 
  geom_histogram(
    aes(my_reps), 
    binwidth = 0.05, 
    boundary = 0,
    fill = "white", 
    color = "black"
  )

What do you think the distribution of p-values is when there is no effect (i.e., the means are identical)? Check this yourself.

Make sure the boundary argument is set to 0 for p-value histograms. See what happens with a null effect if boundary is not set.

Bivariate Normal

n <- 100 # number of random samples

x <- rnorm(n, 0, 1)
y <- x + rnorm(n, 0, 1)

cor(x, y)

## [1] 0.6836834

cor(x, y, method = "spearman")

## [1] 0.569769

n <- 1000 # number of random samples
rho <- 0.5 # population correlation between the two variables

mu <- c(10, 20) # the means of the samples
stdevs <- c(5, 6) # the SDs of the samples

# correlation matrix
cor_mat <- matrix(c(  1, rho, 
                    rho,   1), 2) 

# create the covariance matrix
sigma <- (stdevs %*% t(stdevs)) * cor_mat

# sample from bivariate normal distribution
bvn <- mvrnorm(n, mu, sigma) 

cor(bvn) # check correlation matrix

##           [,1]      [,2]
## [1,] 1.0000000 0.4864211
## [2,] 0.4864211 1.0000000

bvn %>%
  as_tibble() %>%
  ggplot(aes(V1, V2)) +
    geom_point(alpha = 0.5) + 
    geom_smooth(method = "lm") +
    geom_density2d()

Sample distribution

n <- 200 # number of random samples
rho1_2 <- 0.5 # correlation betwen v1 and v2
rho1_3 <- 0 # correlation betwen v1 and v3
rho2_3 <- 0.7 # correlation betwen v2 and v3

mu <- c(10, 20, 30) # the means of the samples
stdevs <- c(8, 9, 10) # the SDs of the samples

# correlation matrix
cor_mat <- matrix(c(     1, rho1_2, rho1_3, 
                    rho1_2,      1, rho2_3,
                    rho1_3, rho2_3,      1), 3) 

sigma <- (stdevs %*% t(stdevs)) * cor_mat
bvn3 <- mvrnorm(n, mu, sigma)

cor(bvn3) # check correlation matrix

##             [,1]      [,2]        [,3]
## [1,]  1.00000000 0.4757010 -0.05956941
## [2,]  0.47570101 1.0000000  0.67172523
## [3,] -0.05956941 0.6717252  1.00000000

3D Plots

You can use the plotly library to make a 3D graph.

library(plotly)

marker_style = list(
    color = "#ff0000", 
    line = list(
      color = "#444", 
      width = 1
    ), 
    opacity = 0.5,
    size = 5
  )

bvn3 %>%
  as_tibble() %>%
  plot_ly(x = ~V1, y = ~V2, z = ~V3, marker = marker_style) %>%
  add_markers()