# Chapter 8 Probability & Simulation

## 8.1 Learning Objectives

### 8.1.1 Basic

- Understand what types of data are best modeled by different distributions
- Generate and plot data randomly sampled from the above distributions
- Test sampled distributions against a null hypothesis
- exact binomial test
- t-test (1-sample, independent samples, paired samples)
- correlation (pearson, kendall and spearman)

- Define the following statistical terms:
- p-value
- alpha
- power
- smallest effect size of interest (SESOI)
- false positive (type I error)
- false negative (type II error)
- confidence interval (CI)

- Calculate power using iteration and a sampling function

### 8.1.2 Intermediate

- Generate 3+ variables from a multivariate normal distribution and plot them

### 8.1.3 Advanced

- Calculate the minimum sample size for a specific power level and design

## 8.2 Resources

- Stub for this lesson
- Distribution Shiny App
- Simulation tutorials
- Chapter 21: Iteration of
*R for Data Science* - Improving your statistical inferences on Coursera (week 1)
- Faux package for data simulation
- Simulation-Based Power-Analysis for Factorial ANOVA Designs (Lakens and Caldwell 2019)
- Understanding mixed effects models through data simulation (DeBruine and Barr 2019)

## 8.3 Distributions

Simulating data is a very powerful way to test your understanding of statistical concepts. We are going to use simulations to learn the basics of probability.

```
# libraries needed for these examples
library(tidyverse)
library(MASS)
set.seed(8675309) # makes sure random numbers are reproducible
```

### 8.3.1 Uniform Distribution

The uniform distribution is the simplest distribution. All numbers in the range have an equal probability of being sampled.

Take a minute to think of things in your own research that are uniformly distributed.

#### 8.3.1.1 Sample continuous distribution

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

Use `runif()`

to sample from a continuous uniform distribution.

```
u <- runif(100000, min = 0, max = 1)
# plot to visualise
ggplot() +
geom_histogram(aes(u), binwidth = 0.05, boundary = 0,
fill = "white", colour = "black")
```

#### 8.3.1.2 Sample discrete distribution

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

Use `sample()`

to sample from a discrete distribution.

You can use `sample()`

to simulate events like rolling dice or choosing from a deck of cards. The code below simulates rolling a 6-sided die 10000 times. We set `replace`

to `TRUE`

so that each event is independent. See what happens if you set `replace`

to `FALSE`

.

```
rolls <- sample(1:6, 10000, replace = TRUE)
# plot the results
ggplot() +
geom_histogram(aes(rolls), binwidth = 1,
fill = "white", color = "black")
```

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] "ferret" "bird" "cat" "dog" "ferret" "cat" "ferret" "ferret"
## [9] "ferret" "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] "ferret" "dog" "ferret" "bird" "dog" "bird" "cat" "ferret"
## [9] "dog" "dog"
```

### 8.3.2 Binomial Distribution

The binomial distribution is useful for modeling binary data, where each observation can have one of two outcomes, like success/failure, yes/no or head/tails.

#### 8.3.2.1 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 heads to 1 and tails to 0.

```
# 20 individual coin flips of a fair coin
rbinom(20, 1, 0.5)
```

`## [1] 0 1 1 1 0 0 1 0 0 0 0 1 1 1 0 0 1 1 0 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 1 1 1 1 1 1 1 0 1 0 1 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] 18`

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

.

`rbinom(10, 20, 0.5)`

`## [1] 9 7 11 9 9 9 13 11 10 10`

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`

.

#### 8.3.2.2 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)
binom.test(biased_coin, n, p = 0.5)
```

```
##
## Exact binomial test
##
## data: fair_coin and n
## number of successes = 5, number of trials = 10, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
## 0.187086 0.812914
## sample estimates:
## probability of success
## 0.5
##
##
## Exact binomial test
##
## data: biased_coin and n
## number of successes = 4, 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.1215523 0.7376219
## sample estimates:
## probability of success
## 0.4
```

Run the code above several times, noting the p-values for the fair and biased coins. Alternatively, you can simulate coin flips online and build up a graph of results and p-values.

- How does the p-value vary for the fair and biased coins?
- What happens to the confidence intervals if you increase n from 10 to 100?
- What criterion would you use to tell if the observed data indicate the coin is fair or biased?
- How often do you conclude the fair coin is biased (false positives)?
- How often do you conclude the biased coin is fair (false negatives)?

#### 8.3.2.3 Statistical terms

The **effect** is some measure of your data. This will depend on the type of data you have and the type of statistical test you are using. For example, if you flipped a coin 100 times and it landed heads 66 times, the effect would be 66/100. You can then use the exact binomial test to compare this effect to the **null effect** you would expect from a fair coin (50/100) or to any other effect you choose. The **effect size** refers to the difference between the effect in your data and the null effect (usually a chance value).

{#p-value} The **p-value** of a test is the probability of seeing an effect at least as extreme as what you have, if the real effect was the value you are testing against (e.g., a null effect). So if you used a binomial test to test against a chance probability of 1/6 (e.g., the probability of rolling 1 with a 6-sided die), then a p-value of 0.17 means that you could expect to see effects at least as extreme as your data 17% of the time just by chance alone.

{#alpha} If you are using null hypothesis significance testing (**NHST**), then you need to decide on a cutoff value (**alpha**) for making a decision to reject the null hypothesis. We call p-values below the alpha cutoff **significant**. In psychology, alpha is traditionally set at 0.05, but there are good arguments for setting a different criterion in some circumstances.

{#false-pos}{#false-neg} The probability that a test concludes there is an effect when there is really no effect (e.g., concludes a 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. The probability that a test concludes there is no effect when there really is one (e.g., concludes a biased coin is fair) is called the **false negative rate** (or *Type II Error Rate*). The **beta** is the false negative rate we accept for a test.

The false positive rate is not the overall probability of getting a false positive, but the probability of a false positive *under the null hypothesis*. Similarly, the false negative rate is the probability of a false negative *under the alternative hypothesis*. Unless we know the probability that we are testing a null effect, we can't say anything about the overall probability of false positives or negatives. If 100% of the hypotheses we test are false, then all significant effects are false positives, but if all of the hypotheses we test are true, then all of the positives are true positives and the overall false positive rate is 0.

{#power}{#sesoi} **Power** is equal to 1 minus beta (i.e., the **true positive rate**), and depends on the effect size, how many samples we take (n), and what we set alpha to. For any test, if you specify all but one of these values, you can calculate the last. The effect size you use in power calculations should be the smallest effect size of interest (**SESOI**). See (Lakens, Scheel, and Isager 2018)(https://doi.org/10.1177/2515245918770963) for a tutorial on methods for choosing an SESOI.

Let's say you want to be able to detect at least a 15% difference from chance (50%) in a coin's fairness, and you want your test to have a 5% chance of false positives and a 10% chance of false negatives. What are the following values?

- alpha =
- beta =
- false positive rate =
- false negative rate =
- power =
- SESOI =

{#conf-int} The **confidence interval** is a range around some value (such as a mean) that has some probability (usually 95%, but you can calculate CIs for any percentage) of containing the parameter, if you repeated the process many times.

A 95% CI does *not* mean that there is a 95% probability that the true mean lies within this range, but that, if you repeated the study many times and calculated the CI this same way every time, you'd expect the true mean to be inside the CI in 95% of the studies. This seems like a subtle distinction, but can lead to some misunderstandings. See (Morey et al. 2016)(https://link.springer.com/article/10.3758/s13423-015-0947-8) for more detailed discussion.

#### 8.3.2.4 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 effects (
`bias`

) - different hypothesised probabilities (
`p`

, defaults to 0.5)

```
sim_binom_test <- function(n, bias, p = 0.5) {
# simulate 1 coin flip n times with the specified bias
coin <- rbinom(1, n, bias)
# run a binomial test on the simulated data for the specified p
btest <- binom.test(coin, n, p)
# returun the p-value of this test
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.9204108`

#### 8.3.2.5 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))
alpha <- 0.05 # this does not always have to be 0.05
mean(my_reps < alpha)
```

`## [1] 0.4663`

`1e4`

is just scientific notation for a 1 followed by 4 zeros (`10000`

). When you're running simulations, you usually want to run a lot of them. 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.

### 8.3.3 Normal Distribution

#### 8.3.3.1 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)
# proportions of normally-distributed data
# within 1, 2, or 3 SD of the mean
sd1 <- .6827
sd2 <- .9545
sd3 <- .9973
ggplot() +
geom_density(aes(dv), fill = "white") +
geom_vline(xintercept = mean(dv), color = "red") +
geom_vline(xintercept = quantile(dv, .5 - sd1/2), color = "darkgreen") +
geom_vline(xintercept = quantile(dv, .5 + sd1/2), color = "darkgreen") +
geom_vline(xintercept = quantile(dv, .5 - sd2/2), color = "blue") +
geom_vline(xintercept = quantile(dv, .5 + sd2/2), color = "blue") +
geom_vline(xintercept = quantile(dv, .5 - sd3/2), color = "purple") +
geom_vline(xintercept = quantile(dv, .5 + sd3/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`

.

#### 8.3.3.2 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 = 4.988, df = 99, p-value = 2.608e-06
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## 0.3087384 0.7166283
## sample estimates:
## mean of x
## 0.5126833
```

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 = -2.9732, df = 193.23, p-value = 0.003323
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
## -0.6300477 -0.1275070
## sample estimates:
## mean of x mean of y
## 0.3850962 0.7638735
```

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.

#### 8.3.3.3 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)
t_ind$statistic
```

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

Alternatively, use `broom::tidy()`

to convert the output into a tidy table.

`broom::tidy(t_ind)`

```
## # A tibble: 1 x 10
## estimate estimate1 estimate2 statistic p.value parameter conf.low conf.high
## <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 -0.379 0.385 0.764 -2.97 0.00332 193. -0.630 -0.128
## # … with 2 more variables: method <chr>, alternative <chr>
```

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) {
# simulate v1
v1 <- rnorm(n, m1, sd1)
#simulate v2
v2 <- rnorm(n, m2, sd2)
# compare using an independent samples t-test
t_ind <- t.test(v1, v2, paired = FALSE)
# return the p-value
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.08335778`

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.2959`

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.

### 8.3.4 Bivariate Normal

#### 8.3.4.1 Correlation

You can test if two continuous variables are related to each other using the `cor()`

function.

Below is one way to generate two correlated variables: `a`

is drawn from a normal distribution, while `x`

and `y`

the sum of and another value drawn from a random normal distribution. We'll learn later how to generate specific correlations in simulated data.

```
n <- 100 # number of random samples
a <- rnorm(n, 0, 1)
x <- a + rnorm(n, 0, 1)
y <- a + rnorm(n, 0, 1)
cor(x, y)
```

`## [1] 0.3525645`

Set `n`

to a large number like 1e6 so that the correlations are less affected by chance. Change the value of the **mean** for `a`

, `x`

, or `y`

. Does it change the correlation between `x`

and `y`

? What happens when you increase or decrease the **sd** for `a`

? Can you work out any rules here?

`cor()`

defaults to Pearson's correlations. Set the `method`

argument to use Kendall or Spearman correlations.

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

`## [1] 0.3062106`

#### 8.3.4.2 Sample distribution

What if we want to sample from a population with specific relationships between variables? We can sample from a **bivariate normal distribution** using `mvrnorm()`

from the `MASS`

package.

```
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 <- MASS::mvrnorm(n, mu, sigma)
cor(bvn) # check correlation matrix
```

```
## [,1] [,2]
## [1,] 1.000000 0.539311
## [2,] 0.539311 1.000000
```

Plot your sampled variables to check everything worked like you expect. It's easiest to convert the output of `mvnorm`

into a tibble in order to use it in ggplot.

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

```
## Warning: `as_tibble.matrix()` requires a matrix with column names or a `.name_repair` argument. Using compatibility `.name_repair`.
## This warning is displayed once per session.
```

### 8.3.5 Multivariate Normal

You can generate more than 2 correlated variables, but it gets a little trickier to create the correlation matrix.

#### 8.3.5.1 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 <- MASS::mvrnorm(n, mu, sigma)
cor(bvn3) # check correlation matrix
```

```
## [,1] [,2] [,3]
## [1,] 1.00000000 0.526778 0.04238033
## [2,] 0.52677804 1.000000 0.73857497
## [3,] 0.04238033 0.738575 1.00000000
```

Alternatively, you can use the (in-development) package faux to generate any number of correlated variables. It also allows to to easily name the variables and has a function for checking the parameters of your new simulated data (`check_sim_stats()`

).

```
#devtools::install_github("debruine/faux")
library(faux)
```

```
##
## ************
## Welcome to faux. For support and examples visit:
## http://debruine.github.io/faux/
## - Get and set global package options with: faux_options()
## ************
```

```
bvn3 <- faux::rnorm_multi(
n = n,
vars = 3,
mu = mu,
sd = stdevs,
r = c(rho1_2, rho1_3, rho2_3),
varnames = c("A", "B", "C")
)
faux::check_sim_stats(bvn3)
```

```
## Warning: All elements of `...` must be named.
## Did you want `multisim_data = c(A, B, C)`?
```

```
## # A tibble: 3 x 7
## n var A B C mean sd
## <dbl> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 200 A 1 0.48 -0.14 10.4 7.88
## 2 200 B 0.48 1 0.63 21.0 8.37
## 3 200 C -0.14 0.63 1 30.6 9.28
```

#### 8.3.5.2 3D Plots

You can use the `plotly`

library to make a 3D graph.

`library(plotly)`

```
##
## Attaching package: 'plotly'
```

```
## The following object is masked from 'package:MASS':
##
## select
```

```
## The following object is masked from 'package:ggplot2':
##
## last_plot
```

```
## The following object is masked from 'package:stats':
##
## filter
```

```
## The following object is masked from 'package:graphics':
##
## layout
```

```
marker_style = list(
color = "#ff0000",
line = list(
color = "#444",
width = 1
),
opacity = 0.5,
size = 5
)
bvn3 %>%
as_tibble() %>%
plot_ly(x = ~A, y = ~B, z = ~C, marker = marker_style) %>%
add_markers()
```

## 8.4 Example

This example uses the Growth Chart Data Tables from the US CDC. The data consist of height in centimeters for the z-scores of –2, -1.5, -1, -0.5, 0, 0.5, 1, 1.5, and 2 by sex (1=male; 2=female) and half-month of age (from 24.0 to 240.5 months).

### 8.4.1 Load & wrangle

We have to do a little data wrangling first. Have a look at the data after you import it and relabel `Sex`

to `male`

and `female`

instead of `1`

and `2`

. Also convert `Agemos`

(age in months) to years. Relabel the column `0`

as `mean`

and calculate a new column named `sd`

as the difference between columns `1`

and `0`

.

`orig_height_age <- read_csv("https://www.cdc.gov/growthcharts/data/zscore/zstatage.csv") `

```
## Parsed with column specification:
## cols(
## Sex = col_character(),
## Agemos = col_character(),
## `-2` = col_double(),
## `-1.5` = col_double(),
## `-1` = col_double(),
## `-0.5` = col_double(),
## `0` = col_double(),
## `0.5` = col_double(),
## `1` = col_double(),
## `1.5` = col_double(),
## `2` = col_double()
## )
```

```
height_age <- orig_height_age %>%
filter(Sex %in% c(1,2)) %>%
mutate(
sex = recode(Sex, "1" = "male", "2" = "female"),
age = as.numeric(Agemos)/12,
sd = `1` - `0`
) %>%
dplyr::select(sex, age, mean = `0`, sd)
```

If you run the code above without putting `dplyr::`

before the `select()`

function, you might get an error message. This is because the `MASS`

package also has a function called `select()`

and, since we loaded `MASS`

after `tidyverse`

, the `MASS`

function becomes the default. When you loaded `MASS`

, you should have seen a warning like "The following object is masked from ‘package:dplyr’: select". You can use functions with the same name from different packages by specifying the package before the function name, separated by two colons.

### 8.4.2 Plot

Plot your new data frame to see how mean height changes with age for boys and girls.

```
ggplot(height_age, aes(age, mean, color = sex)) +
geom_smooth(aes(ymin = mean - sd, ymax = mean + sd), stat="identity")
```

### 8.4.3 Get means and SDs

Create new variables for the means and SDs for 20-year-old men and women.

```
height_sub <- height_age %>% filter(age == 20)
m_mean <- height_sub %>% filter(sex == "male") %>% pull(mean)
m_sd <- height_sub %>% filter(sex == "male") %>% pull(sd)
f_mean <- height_sub %>% filter(sex == "female") %>% pull(mean)
f_sd <- height_sub %>% filter(sex == "female") %>% pull(sd)
height_sub
```

```
## # A tibble: 2 x 4
## sex age mean sd
## <chr> <dbl> <dbl> <dbl>
## 1 male 20 177. 7.12
## 2 female 20 163. 6.46
```

### 8.4.4 Simulate a population

Simulate 50 random male heights and 50 random female heights using the `rnorm()`

function and the means and SDs above. Plot the data.

```
sim_height <- tibble(
male = rnorm(50, m_mean, m_sd),
female = rnorm(50, f_mean, f_sd)
) %>%
gather("sex", "height", male:female)
ggplot(sim_height) +
geom_density(aes(height, fill = sex), alpha = 0.5) +
xlim(125, 225)
```

Run the simulation above several times, noting how the density plot changes. Try changing the age you're simulating.

### 8.4.5 Analyse simulated data

Use the `sim_t_ind(n, m1, sd1, m2, sd2)`

function we created above to generate one simulation with a sample size of 50 in each group using the means and SDs of male and female 14-year-olds.

```
height_sub <- height_age %>% filter(age == 14)
m_mean <- height_sub %>% filter(sex == "male") %>% pull(mean)
m_sd <- height_sub %>% filter(sex == "male") %>% pull(sd)
f_mean <- height_sub %>% filter(sex == "female") %>% pull(mean)
f_sd <- height_sub %>% filter(sex == "female") %>% pull(sd)
sim_t_ind(50, m_mean, m_sd, f_mean, f_sd)
```

`## [1] 0.009012868`

### 8.4.6 Replicate simulation

Now replicate this 1e4 times using the `replicate()`

function. This function will save the returned p-values in a list (`my_reps`

). We can then check what proportion of those p-values are less than our alpha value. This is the power of our test.

```
my_reps <- replicate(1e4, sim_t_ind(50, m_mean, m_sd, f_mean, f_sd))
alpha <- 0.05
power <- mean(my_reps < alpha)
power
```

`## [1] 0.657`

### 8.4.7 One-tailed prediction

This design has about 65% power to detect the sex difference in height (with a 2-tailed test). Modify the `sim_t_ind`

function for a 1-tailed prediction.

You could just set `alternative`

equal to "greater" in the function, but it might be better to add the `alternative`

argument to your function (giving it the same default value as `t.test`

) and change the value of `alternative`

in the function to `alternative`

.

```
sim_t_ind <- function(n, m1, sd1, m2, sd2, alternative = "two.sided") {
v1 <- rnorm(n, m1, sd1)
v2 <- rnorm(n, m2, sd2)
t_ind <- t.test(v1, v2, paired = FALSE, alternative = alternative)
return(t_ind$p.value)
}
alpha <- 0.05
my_reps <- replicate(1e4, sim_t_ind(50, m_mean, m_sd, f_mean, f_sd, "greater"))
mean(my_reps < alpha)
```

`## [1] 0.758`

### 8.4.8 Range of sample sizes

What if we want to find out what sample size will give us 80% power? We can try trial and error. We know the number should be slightly larger than 50. But you can search more systematically by repeating your power calculation for a range of sample sizes.

This might seem like overkill for a t-test, where you can easily look up sample size calculators online, but it is a valuable skill to learn for when your analyses become more complicated.

Start with a relatively low number of replications and/or more spread-out samples to estimate where you should be looking more specifically. Then you can repeat with a narrower/denser range of sample sizes and more iterations.

```
alpha <- 0.05
power_table <- tibble(
n = seq(20, 100, by = 5)
) %>%
mutate(power = map_dbl(n, function(n) {
ps <- replicate(1e3, sim_t_ind(n, m_mean, m_sd, f_mean, f_sd, "greater"))
mean(ps < alpha)
}))
ggplot(power_table, aes(n, power)) +
geom_smooth() +
geom_point() +
geom_hline(yintercept = 0.8)
```

`## `geom_smooth()` using method = 'loess' and formula 'y ~ x'`

Now we can narrow down our search to values around 55 (plus or minus 5) and increase the number of replications from 1e3 to 1e4.

```
power_table <- tibble(
n = seq(50, 60)
) %>%
mutate(power = map_dbl(n, function(n) {
ps <- replicate(1e3, sim_t_ind(n, m_mean, m_sd, f_mean, f_sd, "greater"))
mean(ps < alpha)
}))
##ggplot(power_table, aes(n, power)) +
## geom_smooth() +
## geom_point() +
## geom_hline(yintercept = 0.8) +
## scale_x_continuous(breaks = sample_size)
```

### E References

Lakens, Daniël, and Aaron R Caldwell. 2019. “Simulation-Based Power-Analysis for Factorial Anova Designs,” May. PsyArXiv. doi:10.31234/osf.io/baxsf.

DeBruine, Lisa M, and Dale J Barr. 2019. “Understanding Mixed Effects Models Through Data Simulation,” June. PsyArXiv. doi:10.31234/osf.io/xp5cy.

Lakens, Daniël, Anne M. Scheel, and Peder M. Isager. 2018. “Equivalence Testing for Psychological Research: A Tutorial.” *Advances in Methods and Practices in Psychological Science* 1 (2): 259–69. doi:10.1177/2515245918770963.

Morey, Richard D., Rink Hoekstra, Jeffrey N. Rouder, Michael D. Lee, and Eric-Jan Wagenmakers. 2016. “The Fallacy of Placing Confidence in Confidence Intervals.” *Psychonomic Bulletin & Review* 23 (1): 103–23. doi:10.3758/s13423-015-0947-8.