Power Analysis

This is for conducting a power analysis for experimental design.

So! You're probably here because you want to do a power analysis!

No? You just found this page? Well, you should want to do a power analysis!

A power analysis is important for good experimental design, as it allows you to estimate the sample size you will need to detect an effect of a given size. To read more about the justification, check out this page: https://www.statmethods.net/stats/power.html

To do this power analysis, we will be using an R package called simr. This page demonstrates how to use previous data, either from a pilot or from a similar study to assess your needs for your next experiment. Given the longitudinal design of much of our work, we will be analyzing our data with linear mixed effects models. To check out the original vignettes this work is based on, go here: https://besjournals.onlinelibrary.wiley.com/doi/full/10.1111/2041-210X.12504 OR https://humburg.github.io/Power-Analysis/simr_power_analysis.html

Use R version 3.5.1, simR version 1.0.4

First, load these packages that we will be using:

library(lme4)
library(simr)
library(lmerTest)

lme4 and lmerTest are for making your lmers, and simr is for simulating the data.

In this example, I will be looking at two things:

1. the test-retest reliability of our study. If we have variance between two timepoints that we would not expect based on our experimental design, it may be the effect of the scanner, or our acquisition methodologies. That means any effect smaller than this cannot be reliably ascribed to the experiment (it may just be an artifact). We will use the effect size from the test-retest as a lower bound for effects we can reliably distinguish.
2. an a priori effect size. For the sake of our lmers, I will be discussing effect size in standardized betas (the linear regression coefficients given for effect size). To perform this power analysis, you must have an idea of what size of an effect you are looking for.

So, first I will calculate the test-retest reliability of our study. I will be using MRIs at 2 timepoints during the mouse adulthood (~post-natal day [pnd] 60, and ~pnd 90) for saline-treated mice. Please note, this study was not designed for a test-retest study, so it is possible that there is still some development between pnd 60 and 90 that may increase the effect size of time that we see. However, for now, it will serve as a stand in for a proper test-retest study.

Having subset my dataset to only include adult mice and controls, I first run a lmer with sex and age as fixed effects and subject id as the random effect.

retest_mod <- lmer(scale(left_CA1Rad) ~ sex + scale(age) + (1|mouse_id), data = sal_adults, REML = FALSE)

In this example, I am just looking at results for a single subfield of the hippocampus. In due time I would like to be able to perform the same analyses on all brain regions, but for the sake of this example we will focus on this area.

You can extract the coefficient representing the change between the timepoints like this:

`fixef(retest_mod)["scale(age)"]

Now, we will look at all of our data for the power analysis itself. Create a lmer including all of your interactions of interest--this should be the most complex model that you will want to fit when you have your actual data.

fm1 <- lmer(scale(left_CA1Rad) ~ sex + condition * scale(age) + (1|mouse_id), data = my_data, REML = FALSE)
summary(fm1)

Great, so the coefficients from this model and the variance will serve to help create the simulation data.

Right now, you can assess the effect size of condition, for example like this:

fixef(fm1)["conditionSAL"]

For the simulation, choose the dimension along which you would like to expand your current data set. For me, I want to test how having a higher sample size will increase our power, so I'm going to extend along the variable "mouse_id".

model_ext_mice <- extend(fm1, along = "mouse_id", n = 200) This increased my sample size to 200. If you would like to make sure that your new data preserve the distribution and variance across conditions, sexes, etc. as your initial data, take the time to make some quick histograms.

This extracts your new data:

newData <- getData(model_ext_mice)
hist_orig <- hist(my_data$left_CA1Rad, main = "Volume of left CA1 Rad", xlab = "Volumes (mm)")
hist_new <- hist(newData$left_CA1Rad, main = "Volume of left CA1 Rad (simulated)", xlab = "Volumes (mm)")

Now on to the fun stuff!

At the start, I like to assess my current power to detect an effect of a given size with the real sample size. Let's choose a smallish effect of 0.15. First, we set this to be our effect size.

fixef(fm1)["conditionSAL"] <- 0.15 IMPORTANT. Here, we actually changed the effect in our model!

We can assess the power with this function:

powerSim(fm1, test = fixed("conditionSAL", "t"))

This may take a minute, but then it outputs your given power. Mine is 73.20%. Not our 80% threshold. So, let's turn to the simulated data.

Here we are calculating a power curve with # of subjects on the x-axis and power on the y. This will tell us at what N we achieve 80% power.

fixef(model_ext_mice)["conditionSAL"] <- 0.15
set.seed(0)
p_curve_cond <- powerCurve(model_ext_mice, test = fixed("conditionSAL", "t"), along = "mouse_id", breaks = seq(1, 100, 20), nsim = 100, alpha = 0.5, progress = TRUE)
plot(p_curve_cond)

When the curved line crosses the horizontal 80% line, that is the size your total sample should be! For me, it's about 70 mice! divided by 4, that's 18 mice (rounding up).

If you're interested, you can also replace the effect size with the one derived from your test-retest lmer.

fixef(model_ext_mice)["conditionSAL"] <- fixef(retest_mod)["age"]

Then, you can calculate your curve again. If your effect size is miniscule, like mine (0.00087), then you'll have to boost your number of mice a lot.

p_curve_retest <- powerCurve(model_ext_mice, test = fixed("conditionSAL", "t"), along = "mouse_id", breaks = seq(1, 200, 50), nsim = 500, alpha = 0.05, progress = TRUE)
plot(p_curve_retest)

In this function, you can change the tested numbers with breaks (so you can plot a point every 10 mice, say, instead of 50). You can also change the number of simulations with nsim to tighten up your error bars (more is better, but more computationally intensive). alpha represents the significance levels.

Probably, like me, you will find that your model never reaches 80% power for such a tiny effect size.

Stay tuned for the same analysis on every brain structure!