Statistical power

• Statistical power is defined as the probability, before a study is performed, that a particular comparison will achieve statistical significance at some predetermined level (typically a p-value below 0.05), given some assumed true effect size

• If a certain effect of interest exists (e.g. a difference between two groups) power is the chance that we actually find the effect in a given study

• A power analysis is performed by first hypothesizing an effect size, then making some assumptions about the variation in the data and the sample size of the study to be conducted, and finally using probability calculations to determine the chance of the p-value being below the threshold

• The conventional view is that you should avoid low-power studies because they are unlikely to succeed

• There are several problems with this view, but it’s often required by research funding bodies

The uses of simulation methods

Simulating data uses generating random data sets with known properties using code (or some other method). This can be useful in various contexts.

1. To better understand our models. Probability models mimic variation in the world, and the tools of simulation can help us better understand this variation. Patterns of randomness are contrary to normal human thinking and simulation helps in training our intuitions about averages and variation
2. To run statistical analyses (e.g., simulating a null distribution against which to compare a sample)
3. To approximate the sampling distribution of data and propagate this to the sampling distribution of statistical estimates and procedures

Distribution functions

Base R functions:

• rnorm(): sampling from a normal distribution

• runif(): sampling from a uniform distribution

• rbinom(): sampling from a binomial distribution

• rpois(): sampling from a Poisson distribution

(Other distributions are also available)

• sample(): sampling elements from an R object with or without replacement

• replicate(): often plays a role in conjunction with sampling functions; it is used to evaluate an expression N number of times repeatedly

From non-base packages:

• MASS::mvtnorm(): multivariate normal; sampling multiple variables with a known correlation structure (i.e., we can tell R how variables should be correlated with one another) and normally distributed errors

Sampling from a uniform distribution

The runif function returns some number (n) of random numbers from a uniform distribution with a range from \(a\) (min) to \(b\) (max) such that \(X\sim\mathcal U(a,b)\) (verbally, \(X\) is sampled from a uniform distribution with the parameters \(a\) and \(b\)), where \(-\infty < a < b < \infty\) (verbally, \(a\) is greater than negative infinity but less than \(b\), and \(b\) is finite). The default is to draw from a standard uniform distribution (i.e., \(a = 0\) and \(b = 1\)):

# Sample a vector of ten numbers and store the results in the object `rand_unifs`
# Note that the numbers will be different each time we re-run the `runif` function above.
# If we want to recreate the same sample, we should set a `seed` number first

rand_unifs <- runif(n = 10000, min = 0, max = 1);

The first 40 numbers from the sample are:

 [1] 0.890188523 0.591166379 0.399648855 0.972892989 0.101147806 0.262360107
 [7] 0.181185795 0.474828251 0.250801761 0.511658177 0.738156238 0.628122433
[13] 0.365605151 0.045501343 0.230871157 0.460637787 0.250218366 0.492162029
[19] 0.856357203 0.566203399 0.147660543 0.383201637 0.093662443 0.001523262
[25] 0.111942857 0.936084684 0.014159835 0.679044370 0.960550269 0.389671817
[31] 0.849250339 0.568449457 0.938294313 0.012541510 0.262154853 0.101214718
[37] 0.090681904 0.984074978 0.031501658 0.954211052 0.792093656 0.284537536
[43] 0.661775433 0.281384660 0.681993441 0.760708201 0.130558664 0.793561081
[49] 0.703461216 0.427160911

To visualise the entire sample, we can plot it on a histogram:

Sampling from a normal distribution

The rnorm function returns some number (n) of randomly generated values given a set mean (\(\mu\); mean) and standard deviation (\(\sigma\); sd), such that \(X\sim\mathcal N(\mu,\sigma^2)\). The default is to draw from a standard normal (a.k.a., Gaussian) distribution (i.e., \(\mu = 0\) and \(\sigma = 1\)):

rand_norms_10000 <- rnorm(n = 10000, mean = 0, sd = 1)

print(rand_norms_10000[1:20])

 [1] -0.74432228  0.33997728  0.58149118  0.05342097 -0.65516320  0.30875826
 [7]  0.42112635  0.97564581  1.80259443 -0.64037078  0.23910519 -0.03569876
[13] -0.76900149 -0.54166882 -0.02412268 -0.61409177 -1.12341512 -1.54125305
[19]  0.16835280 -1.26876048

Sampling from a normal distribution

• Histograms allow us to check how samples from the same distribution might vary.

• Exercise: Compare the above distribution with a normal distribution that had a standard deviation of 2 instead of 1.

• Sample 10,000 new values in rnorm with sd = 2 instead of sd = 1 and create a new histogram with hist.

• To see what the distribution of sampled data might look like given a low sample size (e.g., 10), repeat the process of sampling from rnorm(n = 10, mean = 0, sd = 1) multiple times and look at the shape of the resulting histogram.

Sampling from a Poisson distribution

A Poisson process describes events happening with some given probability over an area of time or space such that \(X\sim Poisson(\lambda)\), where the rate parameter \(\lambda\) is both the mean and variance of the Poisson distribution (note that by definition, \(\lambda > 0\), and although \(\lambda\) can be any positive real number, data are always integers, as with count data).

Sampling from a Poisson distribution can be done in R with rpois, which takes only two arguments specifying the number of values to be returned (n) and the rate parameter (lambda). There are no default values for rpois.

rand_poissons <- rpois(n = 10, lambda = 1.5)

print(rand_poissons)

 [1] 2 0 1 1 1 0 0 1 1 2

Sampling from a Poisson distribution

A histogram of a large number of values to see the distribution when \(\lambda = 4.5\):

rand_poissons_10000 <- rpois(n = 10000, lambda = 4.5)

Sampling from a binomial distribution

• A binomial distribution describes the number of successes for some number of independent trials (\(\Pr(success) = p\)).

• The rbinom function returns the number of successes after size trials, in which the probability of success in each trial is prob.

• Sampling from a binomial distribution in R with rbinom is a bit more complex than using runif, rnorm, or rpois.

• Like those previous functions, the rbinom function returns some number (n) of random numbers, but the arguments and output can be slightly confusing at first.

Sampling from a binomial distribution

• For example, suppose we want to simulate the flipping of a fair coin 1000 times, and we want to know how many times that coin comes up heads (success). We can do this with the following code:

coin_flips <- rbinom(n = 1, size = 1000, prob = 0.5)

coin_flips

[1] 491

• The above result shows that the coin came up heads 491 times. But note the (required) argument n. This allows us to set the number of sequences to run.

• If we instead set n = 2, then this could simulate the flipping of a fair coin 1000 times once to see how many times heads comes up, then repeating the whole process a second time to see how many times heads comes up again (or, if it is more intuitive, the flipping of two separate fair coins 1000 times at the same time).

Sampling from a binomial distribution

coin_flips_2 <- rbinom(n = 2, size = 1000, prob = 0.5)

coin_flips_2

[1] 496 506

• A coin was flipped 1000 times and returned 496 heads, and then another fair coin was flipped 1000 times and returned 506 heads.

• As with the rnorm and runif functions, we can check to see what the distribution of the binomial function looks like if we repeat this process.

Sampling from a binomial distribution

• Suppose that we want to see the distribution of the number of times heads comes up after 1000 flips. We can simulate the process of flipping 1000 times in a row with 10000 different coins:

coin_flips_10000 <- rbinom(n = 10000, size = 1000, prob = 0.5)

Random sampling using sample

• Sometimes it is useful to sample a set of values from a vector or list. The R function sample is very flexible for sampling a subset of numbers or elements from some structure (x) in R according to some set probabilities (prob).

• Elements can be sampled from x some number of times (size) with or without replacement (replace), though an error will be returned if the size of the sample is larger than x but replace = FALSE (default).

• Suppose we want to ask R to pick a random number from one to ten with equal probability:

rand_number_1 <- sample(x = 1:10, size = 1)

print(rand_number_1)

[1] 6

Random sampling using sample

• We can increase the size of the sample to 10:

rand_number_10 <- sample(x = 1:10, size = 10)
print(rand_number_10)

 [1]  9  5  1 10  6  8  7  3  2  4

• Note that all numbers from 1 to 10 have been sampled, but in a random order. This is because the default is to sample without replacement, meaning that once a number has been sampled for the first element in rand_number_10, it is no longer available to be sampled again.

• We can change this and allow for sampling with replacement:

rand_number_10_r <- sample(x = 1:10, size = 10, replace = TRUE)

print(rand_number_10_r)

 [1] 10  9  4  7  9  4  6 10  2  1

• Note that the numbers {4, 9, 10} are now repeated in the set of randomly sampled values above.

Random sampling using sample

• So far, because we have not specified a probability vector prob, the function assumes that every element in 1:10 is sampled with equal probability

• Here’s an example in which the numbers 1-5 are sampled with a probability of 0.05, while the numbers 6-10 are sampled with a probability of 0.15, thereby biasing sampling toward larger numbers; we always need to ensure that these probabilities need to sum to 1.

prob_vec      <- c( rep(x = 0.05, times = 5), rep(x = 0.15, times = 5))

rand_num_bias <- sample(x = 1:10, size = 10, replace = TRUE, prob = prob_vec)

print(rand_num_bias)

 [1] 10  8  1 10  7  8  4 10 10  5

Sampling random characters from a list

• We can also sample characters from a list of elements; it is no different than sampling numbers

• For example, if we want to create a simulated data set that includes three different species of some plant or animal, we could create a vector of species identities from which to sample:

species <- c("species_A", "species_B", "species_C");

• We can then sample from these three possible categories. For example:

sp_sample <- sample(x = species, size = 24, replace = TRUE, 
                    prob = c(0.5, 0.25, 0.25))

• What did the code above do?

 [1] "species_C" "species_B" "species_A" "species_B" "species_A" "species_B"
 [7] "species_B" "species_B" "species_A" "species_B" "species_A" "species_B"
[13] "species_A" "species_C" "species_A" "species_A" "species_B" "species_A"
[19] "species_A" "species_A" "species_B" "species_B" "species_A" "species_B"

Simulating data with known correlations

• We can generate variables \(X_{1}\) and \(X_{2}\) that have known correlations \(\rho\) with with one another.

• For example: two standard normal random variables with a sample size of 10000, and with correlation between them of 0.3:

N   <- 10000
rho <- 0.3
x1  <- rnorm(n = N, mean = 0, sd = 1)
x2  <- (rho * x1) + sqrt(1 - rho*rho) * rnorm(n = N, mean = 0, sd = 1)

• These variables are generated by first simulating the sample \(x_{1}\) (x1 above) from a standard normal distribution. Then, \(x_{2}\) (x2 above) is calculated as

\(x_{2} = \rho x_{1} + \sqrt{1 - \rho^{2}}x_{rand}\),

where \(x_{rand}\) is a sample from a normal distribution with the same variance as \(x_{1}\).

Simulating data with known correlations

• We can generate variables \(X_{1}\) and \(X_{2}\) that have known correlations \(\rho\) with with one another.

• For example: two standard normal random variables with a sample size of 10000, and with correlation between them of 0.3:

N   <- 10000
rho <- 0.3
x1  <- rnorm(n = N, mean = 0, sd = 1)
x2  <- (rho * x1) + sqrt(1 - rho*rho) * rnorm(n = N, mean = 0, sd = 1)

• Does the correlation equal rho (with some sampling error)?

cor(x1, x2)

[1] 0.2967485

Simulating data with known correlations

• There is a more efficient way to generate any number of variables with different variances and correlations to one another.

• We need to use the MASS library, which can be installed and loaded as below:

• In the MASS library, the function mvrnorm can be used to generate any number of variables for a pre-specified covariance structure.

Simulating a regression design

• We can use simulation to test rather complex study designs

• Imagine you are interested in students attitude towards smoking and how it depends on the medium of the message and the focus of the message

• We want to know whether people’s attitude is different after seeing a visual anti-smoking message (these pictures on the package) vs. a text-message (the text belonging to that picture)

• We are interested in whether the attitude that people report is different after seeing a message that regards the consequences on other people (e.g. smoking can harm your loved ones) as compared to yourself (smoking can cause cancer)

Study design:

DV: attitude towards smoking (0-100) IV1: medium (text vs. visual) IV2: focus (internal vs. external)

This is, there are 4 groups:

• group_TI will receive text-messages that are internal
• group_TE will receive text-messages that are external
• group_VI will receive visual messages that are internal
• group_VE will receive visual messages that are external

Example: simulating a regression design

• assume that we expect that people’s attitude will be more negative after seeing a visual rather than text message if the focus is internal (i.e. the message is about yourself) because it might be difficult to imagine that oneself would get cancer after reading a text but seeing a picture might cause fear regardless

• for the external focus on the other hand, we expect a more negative attitude after reading a text as compared to seeing a picture, as it might have more impact on attitude to imagine a loved one get hurt than seeing a stranger in a picture suffering from the consequences of second-hand smoking

• we expect that the internal focus messages will be related to lower attitudes compared to the external focus messages on average but we expect no main-effect of picture vs. text-messages

Example: simulating a regression design

• visualize some rough means that show the desired behavior that we described in words earlier and see where we are going

• we could make the overall mean of the internal focus groups (group_TI and group_VI) 20 and the mean of the external groups (group_TE and group_VE) 50 (this would already reflect the main-effect but also a belief that the smoking-attitudes are on average quite negative as we assume both means to be on the low end of the scale)

• assume that the mean of group_TI is 30 while the mean of group_VI is 10 and we could assume that the mean of group_TE is 40 and the mean of group_VE is 60

Example: simulating a regression design

• in the new example there is a difference between the two media groups on average but it is only .50 points, so arguably it is small enough to represent the assumption of no effect, as in real-life no effect in terms of a difference being actually 0 is rather rare

• come up with some reasonable standard-deviation; if we start at 50 and we want most people to be < 80, we can set the 2-SD bound at 80 to get a standard-deviation of 15 (80-50)/2.

• let’s assume that each of our groups has a standard-deviation of 15 points.

group_TI = normal(n, 43, 15)
group_VI = normal(n, 40, 15)
group_TE = normal(n, 45, 15)
group_VE = normal(n, 47, 15)

Example: simulating a regression design

n <- 1e5
group_TI <-  rnorm(n, 43, 15)
group_VI <-  rnorm(n, 40, 15)
group_TE <-  rnorm(n, 45, 15)
group_VE <-  rnorm(n, 47, 15)

participant <- c(1:(n*4))
focus <- rep(c("internal", "external"), each = n*2)
media <- rep(c("text", "visual"), each = n, times = 2)

data <- data.frame(participant = participant, focus = focus, media = media, score = c(group_TI, group_VI, group_TE, group_VE))

summary(data)

  participant       focus              media               score       
 Min.   :1e+00   Length:400000      Length:400000      Min.   :-26.82  
 1st Qu.:1e+05   Class :character   Class :character   1st Qu.: 33.41  
 Median :2e+05   Mode  :character   Mode  :character   Median : 43.69  
 Mean   :2e+05                                         Mean   : 43.72  
 3rd Qu.:3e+05                                         3rd Qu.: 54.00  
 Max.   :4e+05                                         Max.   :115.45  

Ready for power-analysis

• Some additional assumptions: suppose we have enough funding for a sizeable data collection and the aim is to ensure that we do not draw unwarranted conclusions from the research

• We should then set the alpha-level at a more conservative value (\(\alpha = .001\)); with this, we expect to draw non-realistic conclusions in the interaction effect in only about 1 in every 1,000 experiments

• We also want to be sure that we do detect an existing effect and keep our power high at 95%; with this, we expect that if there is an interaction effect, we would detect it in 19 out of 20 cases (only miss it in 1 out of 20, or 5%)

• Running the power-simulation can be very memory-demanding and the code can run a very long time to complete; it’s advised to start from various low-resolution sample-sizes (e.g. n = 10, n = 100, n = 200, etc.) to get a rough idea of where we can expect our loop to end. Then, the search can be made more specific in order to identify a more precise sample size.

Example: simulating a regression design

We can plot the results form the power-simulation:

• At roughly 900 participants we observe sufficient power