This vignette assumes you have already familiarized yourself with the
basic operations of the “plug-in” approach to power analysis and the
powerutilities package. If you have not, please refer to
the “introduction-powerutilities” vignette.
Load all required packages:
library(dplyr)
library(tidyr)
library(glmmTMB)
library(powerutilities)
library(emmeans)This example is inspired by the work of Ataide et al (2025). The hypothetical planned experiment is a 3x2 factorial of insecticide treatments (spinosad and sesame oil, along with a water control) and application method (prophylactic or curative). Plants are arranged in a randomized complete block design, and then thrips will be introduced into each plant to observe the mortality rate. For illustration purposes, lets imagine we have a total of 240 thrips and the question is whether to arrange them such that there are 4 blocks with 10 thrips on each plant or 8 blocks with 5 thrips on each plant. We start by creating our treatment dictionary containing all combinations of our treatments:
ex1_trt_dict = expand.grid(Insecticide = c('Water', 'Spinosad', 'Sesame Oil'),
Application = c('Curative', 'Prophylactic'))
ex1_trt_dict#> Insecticide Application
#> 1 Water Curative
#> 2 Spinosad Curative
#> 3 Sesame Oil Curative
#> 4 Water Prophylactic
#> 5 Spinosad Prophylactic
#> 6 Sesame Oil ProphylacticWe then assign mean response values (as proportions) to that treatment dictionary. We then create two separate data sets corresponding to our two proposed designs by merging the treatment dictionary to incorporate blocks and adding a column containing the total number of thrips on each plant:
ex1_trt_dict$Mortality = c(.05, 0.20, 0.15,
.05, 0.35, 0.25)
ex1a = ex1_trt_dict |>
merge(data.frame(Block = factor(1:4))) |>
mutate(Ttl_Thrips = 10)
ex1b = ex1_trt_dict |>
merge(data.frame(Block = factor(1:8))) |>
mutate(Ttl_Thrips = 5)Now, note that as discrete proportions, we will model this as binomial data, and the binomial distribution has no dispersion parameter: the variance is equal to kp(1-p), which means its entirely determined by the proportion parameter, p, and the number of trials, k. This would seem to simplify the process of power estimation, since we have 1 less “plug-in” value to determine. However, we now also have the issue of a link function (here, the logit), which complicates the issue of specifying a block variance, since this is the variance on the logit scale. Where an investigator is required to rely on educated guesswork to find an appropriate plug-in value, it is highly recommended they employ some simulation to get a feel for what the implications of different block (or other random effect) variances have on the original scale. Here, we’ll simply take the block variance as 0.3 (again, on the logit scale).
Fit the model to the data with fewer blocks and more thrips per plant:
ex1a_mod = set_glmm(Mortality ~ Insecticide*Application + (1|Block),
data = ex1a, re_terms = c(0.3), REML = F,
family = binomial(link = 'logit'),
weights = ex1a$Ttl_Thrips)#> Warning in eval(family$initialize): non-integer #successes in a binomial glm!VarCorr(ex1a_mod)#>
#> Conditional model:
#> Groups Name Std.Dev.
#> Block (Intercept) 0.3Fit the same model to the data with more blocks and fewer thrips per plant:
# A warning is produced about non-integer number of successes, but we
# can safely ignore this.
ex1b_mod = set_glmm(Mortality ~ Insecticide*Application + (1|Block),
data = ex1b, family = binomial(link = 'logit'),
weights = ex1b$Ttl_Thrips, REML = F,
re_terms = 0.3)#> Warning in eval(family$initialize): non-integer #successes in a binomial glm!VarCorr(ex1b_mod)#>
#> Conditional model:
#> Groups Name Std.Dev.
#> Block (Intercept) 0.3Before proceeding any further with the power analysis, let’s first examine the marginal means:
emmeans(ex1a_mod, ~ Insecticide:Application, type = 'response')#> Insecticide Application prob SE df asymp.LCL asymp.UCL
#> Water Curative 0.0489 0.0345 Inf 0.0119 0.180
#> Spinosad Curative 0.1975 0.0673 Inf 0.0967 0.361
#> Sesame Oil Curative 0.1477 0.0590 Inf 0.0646 0.303
#> Water Prophylactic 0.0489 0.0345 Inf 0.0119 0.180
#> Spinosad Prophylactic 0.3482 0.0830 Inf 0.2069 0.522
#> Sesame Oil Prophylactic 0.2475 0.0739 Inf 0.1314 0.417
#>
#> Confidence level used: 0.95
#> Intervals are back-transformed from the logit scaleWhy do these estimates not coincide with the ones we specified in the “treatment dictionary”? The answer lies in the fact that what we specified is marginal proportion, while what we are estimating in the conditional proportion. This distinction, and its implications for modeling and estimation, are both subtle and important. Unfortunately, they are beyond the scope of this current work. Interested readers (and you should all be interested) will find Section 3.5 of Stroup, Ptukhina and Garai (2024) an excellent jumping-off point for understanding this topic which does not arise in either generalized linear models or linear mixed models, but is distinct to generalized linear mixed models.
Now, turn to the power of the F-tests. First, we’ll consider
the calculated power when using asymptotic Wald tests. You will notice
that the differences in power, which derive exclusively from the small
variation in the estimated coefficients, are negligibly small. These
Wald tests are the default when using car::Anova() or
emmeans::joint_tests, but they are known to
“anti-conservative”, which means that they have a tendency to result in
inflated type I error rates.
power_ftest(ex1a_mod, ddf = 'asymptotic')
power_ftest(ex1b_mod, ddf = 'asymptotic')Alternative approaches to calculating the degrees of freedom, in contrast, can lead to more pronounced differences in the power:
Residual DF:
power_ftest(ex1a_mod, ddf = 'df.residual')#> Warning in eval(family$initialize): non-integer #successes in a binomial glm!#> Warning in .warn_ddf_glmm(ddf): performance (and theoretical justification) of ddf='df.residual' for GLMMs is poorly understoodpower_ftest(ex1b_mod, ddf = 'df.residual')#> Warning in eval(family$initialize): non-integer #successes in a binomial glm!
#> Warning in eval(family$initialize): performance (and theoretical justification) of ddf='df.residual' for GLMMs is poorly understoodThese differences might not be huge, but its important to be aware of them. The take-home message is that Wald \(\chi^2\)-tests and, equivalently, F-test with Inf denominator degrees of freedom provide, at best, an upper limit for your statistical power and, during actual testing, your p-values.
Turning now to the contrasts, we start by calculating the marginal means for the models fit to each of our two candidate designs:
(ex1a_emm = emmeans(ex1a_mod, ~Insecticide:Application, type = 'response'))#> Insecticide Application prob SE df asymp.LCL asymp.UCL
#> Water Curative 0.0489 0.0345 Inf 0.0119 0.180
#> Spinosad Curative 0.1975 0.0673 Inf 0.0967 0.361
#> Sesame Oil Curative 0.1477 0.0590 Inf 0.0646 0.303
#> Water Prophylactic 0.0489 0.0345 Inf 0.0119 0.180
#> Spinosad Prophylactic 0.3482 0.0830 Inf 0.2069 0.522
#> Sesame Oil Prophylactic 0.2475 0.0739 Inf 0.1314 0.417
#>
#> Confidence level used: 0.95
#> Intervals are back-transformed from the logit scale(ex1b_emm = emmeans(ex1b_mod, ~Insecticide:Application, type = 'response'))#> Insecticide Application prob SE df asymp.LCL asymp.UCL
#> Water Curative 0.0486 0.0340 Inf 0.0120 0.177
#> Spinosad Curative 0.1968 0.0650 Inf 0.0986 0.354
#> Sesame Oil Curative 0.1471 0.0574 Inf 0.0657 0.297
#> Water Prophylactic 0.0486 0.0340 Inf 0.0120 0.177
#> Spinosad Prophylactic 0.3477 0.0795 Inf 0.2115 0.514
#> Sesame Oil Prophylactic 0.2469 0.0711 Inf 0.1341 0.410
#>
#> Confidence level used: 0.95
#> Intervals are back-transformed from the logit scaleNext, we will focus on two contrasts specifically, not because they are necessarily very interesting in the context of our hypothetical study but because they highlight an additional complication that arises specifically for generalized linear models and generalized linear mixed models which are not present for normally distributed data.
The key item to remember when considering these two contrasts is that the difference in proportions between the treatments being compared is the same for both contrasts: 0.1.
ex1_contr_list = list('Sesame Curative - Control' = c(-1, 0, 1, 0, 0, 0),
'Sesame Prophylactic - Curative' = c(0, 0, -1, 0, 0, 1))power_contrast(ex1a_emm, ex1_contr_list)#> Setting `type = "emmeans"`: power calculations must be performed on the link scale.We are warned that power calculations must be performed on the link scale, and that those differences are not consistent between the two comparisons. As a result, the power of these two tests differ.