Introduction to the “Plug-In” Method of Power Analysis in R

Introduction

Loosely speaking, statistical power can be thought of as our capacity to avoid inconclusive test results. As researchers, we want our studies to provide precise estimates and powerful tests. All else being equal, statistical power will be greater with:

Of these factors, the sample size and the experimental design are the two under the most direct control of the experimenter. Diligence and care in the conduct of the experiment can help reduce the residual dispersion to a degree, but some level of background noise is inevitable. And the true difference between the parameters being compared is, of course, entirely outside the control of the researcher. Thus, the focus during the planning stages is usually on determining the required sample size to achieve some pre-specified power level (given some experimental design), determining the power of a test based on the number of available replicates (again, assuming a particular experimental design), or determining which of several candidate designs will maximize the statistical power given some number of replicates. In all of these cases, the true magnitude of the difference and the magnitude of the residual variance is assumed known, usually being based on the results of a pilot study, on findings from previously published studies, or from educated guesswork.

For simple linear models, where the data are assumed independent and normally distributed and where the only random effect is the residual variance, there are fairly simple equations available for determining power or sample size requirements (for an excellent introduction and discussion of these equations, along with much, much else, see: Piepho et al., 2022 ). For linear mixed models (i.e., where the experimental design implies the presence of additional random effects, heteroscedasticity, or correlated residuals), or for generalized linear models (i.e., where the data are characterized by a distribution other than the normal), these calculations can become much more complex and the relationships among the different factors less simple and intuitive. For this reason, alternative approaches are needed for power analysis/sample size determination/design planning. The origins of the approach to power analysis employed here, which Piepho et al., 2022 refers to as the “plug-in approach”, are actually quite old, but a recent description (for its use in SAS) can be found in the last chapter of Generalized Linear Mixed Models: Modern Concepts and Methods. It has the advantages of being very general, capable of accommodating a wide variety of designs/models, while remaining simpler to implement and less computationally demanding than simulation-based power analysis. The basic idea is to first encode information about sample size and expected treatments means (and differences) for the proposed experiment in a fake data set, to then specify a model encoding information about the design of the proposed experiment, and finally to “fit” the model while fixing (“plugging in”) the values for the random effects and dispersion parameters. The result is that the modelling software does all the work of calculating the required standard errors which are needed to derive the power of the tests under the hypothetical conditions defined by our fake data and proposed model. By manipulating the fake data set, model, or variance parameter values, it becomes possible to explore how differences in sample size, effect size, experimental design, or random effects will impact statistical power. In R, the glmmTMB package allows one to fix parameter estimates during model fitting, making it the perfect tool for this sort of power analysis.

Housekeeping

To begin, load the required packages:

library(dplyr)
library(tidyr)
library(glmmTMB)
library(emmeans)
library(powerutilities)

Example 1: A two-sample t-test

To begin, lets consider the simplest example, in which we will establish a single factor experiment with two levels (‘A’ and ‘B’) arranged in a completely randomized design (CRD) with 8 replicates per treatment. Under these circumstances, we can actually skip the modeling stage and simply perform a two-sample t-test on this data, and for that we have tools in base R for performing power analysis. By simply leaving one of the arguments in the set n, delta, sd, sig.level and power to the value NULL, that value will be calculated on the basis of the other arguments provided. Thus, if we assume a residual standard deviation of 2 and a difference in treatment means of 3, then our power is just shy of 0.8, or 80%, which is widely considered the threshold for an adequately powered study.

power.t.test(n = 8, delta = 3, sd = 2, 
             sig.level = 0.05, power = NULL, 
             type = 'two.sample', 
             alternative = 'two.sided')
#> 
#>      Two-sample t test power calculation 
#> 
#>               n = 8
#>           delta = 3
#>              sd = 2
#>       sig.level = 0.05
#>           power = 0.7965441
#>     alternative = two.sided
#> 
#> NOTE: n is number in *each* group

We can also use the alternative “plug-in” approach, by creating a fake data set and “fitting” the model using a fixed value for the residual standard deviation. The first step is to create a data set with the correct number of rows with treatment “A” and “B”, thereby encoding the proposed sample size. Then, the response value, Y, is set such that the difference between the treatments is 3 (the actual values are unimportant here, although that will not always be true):

nrep = 8 # Number of replicates
nfac = 2 # Number of treatments

# Expand grid creates a data set
# with all combinations of treatments
ex1 = expand.grid(Rep = factor(1:nrep),
                  Trt = factor(LETTERS[1:nfac])) 

# Set the response values to correspond 
# to the expected treatment differences
ex1$Y = ifelse(ex1$Trt == 'A', 0, 3)

head(ex1, n = 10)
#>    Rep Trt Y
#> 1    1   A 0
#> 2    2   A 0
#> 3    3   A 0
#> 4    4   A 0
#> 5    5   A 0
#> 6    6   A 0
#> 7    7   A 0
#> 8    8   A 0
#> 9    1   B 3
#> 10   2   B 3

Now, we construct a model in which the residual standard deviation is set to 2 rather than being estimated from the data. This is done using set_glmm, which is essentially just a wrapper for glmmTMB modified to somewhat simplify the syntax and model random effects parameterizations:

ex1_mod = set_glmm(Y ~ Trt, data = ex1, disp = 2)

# Confirm it worked:
sigma(ex1_mod)
#> [1] 2

On the basis of those results, we can calculate power of the F-test (i.e., ANOVA table) using the power_ftest function from the powerutilities package:

power_ftest(ex1_mod)

Note first, that the power calculated here is identical to the one we got from the power.t.test() function (F-tests are generalizations of t-tests to more than one mean, so an F-test for a set of 2 groups is identical to a t-test). Also note that the output includes a degree-of-freedom method, the default being “df.residual” for linear models (i.e., models without random effects and where the data are assumed normally distributed). We will not go into detail here, but the choice of DF method will effect the p-value and power:

power_ftest(ex1_mod, ddf = 'asymptotic')
#> For gaussian fixed-effects models, it is recommended to use ddf = "df.residual"

In addition to the F-test, we are usually also interested in assessing the power of specific contrasts. To do this, we first need to calculate the estimated marginal means for the factor levels we want to compare:

(ex1_emm = emmeans(ex1_mod, ~Trt)) # Estimated treatment means
#>  Trt emmean    SE  df asymp.LCL asymp.UCL
#>  A    2e-07 0.707 Inf     -1.39      1.39
#>  B    3e+00 0.707 Inf      1.61      4.39
#> 
#> Confidence level used: 0.95

We then define our contrasts of interest as a named list, each item of which is a vector of weights. For more details on how and why we specify contrasts, refer to vignette('contrasts-powerutilities'). For this example there is essentially only one meaningful contrast: testing for a difference between treatments A and B.

contr_list = list('Trt A - Trt B' = c(1, -1))

And then use the power_contrast function from powerutilities to calculate the power of that test, again assuming the proposed study matched the fake model and data we defined earlier.

power_contrast(ex1_emm, contr_list)

There are several things to note about this output. Firstly, we note that p-value and power shown here are identical to those calculated previously, as they should be. Also note that besides statistical power, this output also provides two other important metrics: the probability of a sign error (“TypeS”) and the expected ratio between a significant estimate and the true value (magnitude, or “TypeM” error). The former is the probability that, assuming you conducted this experiment, performed the test and found a statistically significant difference, that you’d get the actual relationship backwards (i.e., your estimates would have A larger than B). Naturally, we want this value as close as possible to zero. The type M error is the expected ratio between a statistically significant estimate and the true value: underpowered studies will tend, when they do happen to achieve statistically significant findings, to produce estimates which are markedly biased. We want this value as close to 1 as possible. Further details on the measures can be found in Gelman & Carlin.

(ex1_emm = emmeans(ex1_mod, ~Trt))
#>  Trt emmean    SE  df asymp.LCL asymp.UCL
#>  A    2e-07 0.707 Inf     -1.39      1.39
#>  B    3e+00 0.707 Inf      1.61      4.39
#> 
#> Confidence level used: 0.95
power_contrast(ex1_emm, contr_list)

Two additional options are worth noting here: first, we can also change the degrees of freedom calculation method, if needed, for power_contrast by simply using ddf argument as before. We also have the option to get confidence intervals on our contrasts by using infer = c(T, T):

ex1_emm = emmeans(ex1_mod, ~Trt)
power_contrast(ex1_emm, contr_list, ddf = 15, infer = c(T, T)) |> 
  print(digits = 3)
#>        Contrast Estimate  SE  LCL  UCL NumDF DenDF  Fval Fcrit   Pval  Power  TypeS TypeM
#> 1 Trt A - Trt B     -3.0 1.0 -5.1 -0.9     1    15 9.000 4.543 0.0090 0.8006 <.0001 1.115
#> 
#> Degrees-of-freedom method: user-specified
#> α = 0.05

Now let’s introduce a small additional complication: suppose we have reason to expect that samples given the “A” treatment in our study will have somewhat lower residual variance than samples given the “B” treatment. Will this effect the power of our test? There is no option in the power.t.test function to accomodate heteroscedasticity, but with set_glmm it’s quite simple:

ex1_mod2 = set_glmm(Y ~ Trt, data = ex1, disp = c(1.75, 2.25), dispformula = ~ Trt - 1)

In the code above, we’ve supplied two values for the dispersion parameter, and added the argument dispformula = ~ Trt - 1. This formula tells glmmTMB to model two dispersion parameters, one corresponding to level of ‘Trt’ (we added “-1” so that we could provide the values as the residual standard deviations, rather than as “baseline residual std. dev.” and “difference from baseline residual std. dev.”). We can now see that the estimated marginal means have different standard errors, reflecting their different residual standard deviations:

(ex1_emm2 = emmeans(ex1_mod2, ~Trt)) 
#>  Trt emmean    SE  df asymp.LCL asymp.UCL
#>  A        0 0.619 Inf     -1.21      1.21
#>  B        3 0.795 Inf      1.44      4.56
#> 
#> Confidence level used: 0.95

Here, we can see that this results in a very minor loss of power, from .7965 to .7904.

power_contrast(ex1_emm2, contr_list) |> 
  print(digits = 3)
#>        Contrast Estimate  SE NumDF DenDF  Fval Fcrit   Pval  Power  TypeS TypeM
#> 1 Trt A - Trt B     -3.0 1.0     1    14 8.862 4.600 0.0100 0.7904 <.0001 1.119
#> 
#> Degrees-of-freedom method: df.residual
#> α = 0.05

In fact, the math here is not so difficult, and we could have worked this out from first principles. For two means calculated from samples of size \(n_1\) and \(n_2\), respectively, drawn from populations with residual standard deviations of \(s_1\) and \(s_2\), respectively, the formula for the standard error of the difference is:

\[SED = \sqrt{\frac{s_1^2}{n_1} + \frac{s_2^2}{n_2}}\] In situations where the n’s and s’s are equal in the two groups, this simplifies to:

\[SED = \sqrt{\frac{2s^2}{n}} \]

Which when \(s = 2\) and \(n = 8\) for both groups becomes:

\[\sqrt{\frac{2*2^2}{8}} = 1\] Whereas keeping \(n = 8\) for both groups but setting \(s_1 = 1.75\) and \(s_2 = 2.25\), we get:

\[\sqrt{\frac{1.75^2}{8} + \frac{2.25^2}{8}} \approx 1.0078\] The effect is small, but real: heteroscedasticity has increased the standard error of the difference, which will reduce power. We’ve now seen that for this simple example the “plug-in” method produces the same results as the basic power equations do. However, this “plug-in” approach will also work for a much wider variety of designs/models than we have tools for in base R or that are convenient to work though by hand or from first principles.

Example 2: A Split-Plot Randomized Complete Block Design

We are basing this proposed experiment on the one described in Frank Yates’ 1935 paper entitled “Complex experiments” (https://doi.org/10.2307/2983638), which was structured as a randomized complete block design (RCBD) with six blocks, a main-plot factor of genotype with three levels (‘Marvellous’, ‘Victory’ and ‘GoldenRain’), and a sub-plot factor of nitrogen application rate, at 4 levels (here, we use: 0, 20, 40, 60 kg/ha). Let’s imagine that at zero nitrogen fertilization, we expect Victory and Golden Rain to yield 1400 kg/ha and Marvellous to yield 1500 kg/ha. Meanwhile, you expect increases of 150 kg/ha at each higher level of N, except for Victory, which you expect to gain 200 kg/ha for each additional increase in N. Finally, you expect similar standard deviation among blocks, mainplots, and the residual error of around 250, 200 and 250, respectively.

The first step is to create a “treatment dictionary”, which describes the treatment structure. For a factorial treatment such like this, the expand.grid function can be used to create all combinations of the two factor levels:

ex2_trt_dict = expand.grid(gen = c('GoldenRain', 'Marvellous', 'Victory'), 
                           nitro = factor(c(0, 20, 40, 60))) |> 
  arrange(gen)
ex2_trt_dict
#>           gen nitro
#> 1  GoldenRain     0
#> 2  GoldenRain    20
#> 3  GoldenRain    40
#> 4  GoldenRain    60
#> 5  Marvellous     0
#> 6  Marvellous    20
#> 7  Marvellous    40
#> 8  Marvellous    60
#> 9     Victory     0
#> 10    Victory    20
#> 11    Victory    40
#> 12    Victory    60

The, the expected mean response values need to be added in. This is the most tedious part of the “plug-in” process, but if you keep yourself organized, doesn’t need to be too bad. Here, I’ve kept my treatments arranged in an organized manner, and am entering the expected treatment means as formulae, so that I can more easily see/compare my written expectations with the actual values. Note that if it helps, the treatment dictionary can also be created in Excel (or the spreadsheet software of your choice), and then simply imported into R.

ex2_trt_dict$yield = c(1400 + 150*0, 
                       1400 + 150*1, 
                       1400 + 150*2, 
                       1400 + 150*3, 
                       1500 + 150*0, 
                       1500 + 150*1, 
                       1500 + 150*2, 
                       1500 + 150*3, 
                       1400 + 200*0, 
                       1400 + 200*1, 
                       1400 + 200*2, 
                       1400 + 200*3)
ex2_trt_dict 
#>           gen nitro yield
#> 1  GoldenRain     0  1400
#> 2  GoldenRain    20  1550
#> 3  GoldenRain    40  1700
#> 4  GoldenRain    60  1850
#> 5  Marvellous     0  1500
#> 6  Marvellous    20  1650
#> 7  Marvellous    40  1800
#> 8  Marvellous    60  1950
#> 9     Victory     0  1400
#> 10    Victory    20  1600
#> 11    Victory    40  1800
#> 12    Victory    60  2000

We then create a data frame encoding the design. Here, it will simply be a data frame with one column for “Block”, but in other contexts we might need to use expand grid and incorporate additional columns for location and/or year, for example.

nblk = 6
ex2_des_df = data.frame(block = factor(1:nblk))
ex2_des_df
#>   block
#> 1     1
#> 2     2
#> 3     3
#> 4     4
#> 5     5
#> 6     6

Using the merge function will then replicate the whole treatment dictionary for every row (i.e., block) in our design data frame:

ex2 = merge(ex2_trt_dict, ex2_des_df)

head(ex2, 10)
#>           gen nitro yield block
#> 1  GoldenRain     0  1400     1
#> 2  GoldenRain    20  1550     1
#> 3  GoldenRain    40  1700     1
#> 4  GoldenRain    60  1850     1
#> 5  Marvellous     0  1500     1
#> 6  Marvellous    20  1650     1
#> 7  Marvellous    40  1800     1
#> 8  Marvellous    60  1950     1
#> 9     Victory     0  1400     1
#> 10    Victory    20  1600     1

We can now fit the model. Fixing the value of the residual standard deviation is done in exactly the same way as in the previous example, but here we also have two random effects standard deviations to incorporate as well: one for block and one for genotypes-within-blocks. In order to know the order that those terms need to be supplied, there is a function called theta_finder. Here, we can see that for this specification of the random effects, our values should be ordered (1) main-plot (i.e., genotype-within-block), (2) block:

theta_finder(yield ~ gen*nitro + (1|block/gen), data = ex2)
#>  Group     Name        Std.Dev. Cor. (us)
#>  gen:block (Intercept) 1                 
#> 
#>  Group     Name        Std.Dev. Cor. (us)
#>  block     (Intercept) 2

As before, the set_glmm function handles the requisite transformations, so that the random effects can simply be given as standard deviations, and any correlations present among the random effects can be given as correlations:

ex2_mod = set_glmm(yield ~ gen*nitro + (1|block/gen), data = ex2, 
                   re_terms = c(200, 250), disp = 225)

VarCorr(ex2_mod, sigma = extract_disp(ex2_mod)) # It worked
#> 
#> Conditional model:
#>  Groups    Name        Std.Dev.
#>  gen:block (Intercept) 200     
#>  block     (Intercept) 250     
#>  Residual              225

And we can now assess the power of the F-tests associated with our proposed study. Here, since we have a mixed model, the default denominator degrees of freedom is now ‘kenward-roger’, which reduces to the classical “containment” method in the balanced split-plot example, but which is broadly applicable to a range of different experimental designs containing nesting structures. One of the most important things to notice is that, because of the split-plot aspect of the experimental design, we have fewer denominator degrees of freedom for the main-plot factor (genotype) than we do for either of the sub-plot factors (nitrogen and the genotype x nitrogen interaction). All else being equal, fewer denominator degrees of freedom implies lower power. In this case, we can see that we have very high power for detecting effect of nitrogen, but very low power for detecting either the main effect of genotype and for the genotype:nitro interaction.

power_ftest(ex2_mod)

For the sake of illustration, let us imagine that Victory is something of a “standard” cultivar which we want to use as a reference or control. Let us suppose then that we are interested in comparing the performance of Golden Rain and Marvellous against the performance of Victory separately at each nitrogen level. To do this, we first calculate our marginal means for each genotype while using nitro as a grouping variable (via the “|” operator):

(ex2_emm1 = emmeans(ex2_mod, ~ gen|nitro, ddf = 'kenward-roger'))
#> nitro = 0:
#>  gen        emmean  SE   df lower.CL upper.CL
#>  GoldenRain   1400 160 16.3     1062     1738
#>  Marvellous   1500 160 16.3     1162     1838
#>  Victory      1400 160 16.3     1062     1738
#> 
#> nitro = 20:
#>  gen        emmean  SE   df lower.CL upper.CL
#>  GoldenRain   1550 160 16.3     1212     1888
#>  Marvellous   1650 160 16.3     1312     1988
#>  Victory      1600 160 16.3     1262     1938
#> 
#> nitro = 40:
#>  gen        emmean  SE   df lower.CL upper.CL
#>  GoldenRain   1700 160 16.3     1362     2038
#>  Marvellous   1800 160 16.3     1462     2138
#>  Victory      1800 160 16.3     1462     2138
#> 
#> nitro = 60:
#>  gen        emmean  SE   df lower.CL upper.CL
#>  GoldenRain   1850 160 16.3     1512     2188
#>  Marvellous   1950 160 16.3     1612     2288
#>  Victory      2000 160 16.3     1662     2338
#> 
#> Confidence level used: 0.95

We can then define our contrasts only in terms of genotype, and these contrasts are then automatically applied to each level of the grouping factor:

ex2_contr1 = list('Golden - Victory' = c(1, 0, -1), 
                  'Marvellous - Victory' = c(0, 1, -1))

First, we can see that where the true difference is zero, we have power of 0.05, which makes sense, since this is the false positive rate alpha! However, we can also see that a true difference of only 100kg/ha, given the magnitude of the variation among plots, blocks, etc., is very low at just above 0.086. Even more concerning, note that even for the pairwise comparison with the greatest difference (Golden - Victory @ 60 nitro) the power is still very low, at only about 13.2%. If such comparisons were the primary objectives of the study and we had no way to increase the sample size, we would be forced to seriously consider whether such an experiment was worth the investment.

power_contrast(ex2_emm1, ex2_contr1)

Now what would happen if employed a different design? What if we realized we didn’t need to use a split-plot design and could use a simple RCBD?

ex2_mod2 = set_glmm(yield ~ gen*nitro + (1|block), data = ex2, 
                    re_terms = 250, disp = 225)
power_ftest(ex2_mod2)

While still low, the power of our test among genotypes has jumped from about 8.7% to 26.7%! Meanwhile, there has been only very minor change in the power of the other tests. This is because the jump from 10 denominator degrees of freedom to 55 is much greater than the jump from 45 to 55. Now, what would be the implication of conducting an experiment as a split-plot RCBD but then ignoring the design during modeling and instead analysing it as an RCBD? We would expect to be artificially lowering our p-values for each of these tests by exactly the amounts shown in these tables! This is commonly known as “psuedo-replication”, and it is something to always avoid.