| Title: | Power Analysis for Linear and Generalized Linear Mixed Models |
|---|---|
| Description: | Prospective power analysis and design planning using the "plug-in" approach described by Piepho et al. (2022) <doi:10.1017/S0021859622000466> and Stroup, Ptukhina and Garai (2024) <doi:10.1201/9780429092060-24>. Expected treatment means and the proposed experimental design are encoded in an artificial data set, and a model is "fit" to it with random effect variances, correlations and residual dispersion fixed at user-specified values rather than being estimated. The resulting standard errors are used to compute the power of F-tests and user-defined contrasts, as well as Type S (sign) and Type M (magnitude) error rates (Gelman and Carlin, 2014 <doi:10.1177/1745691614551642>). Models are fit with 'glmmTMB', so the approach accommodates non-normal responses, heteroscedasticity, nested and crossed random effects, correlated random slopes, and structured covariances such as AR(1) and spatial correlation. Helper functions convert standard deviations and correlations into the internal parameterization of 'glmmTMB' and report the order in which random effect parameters must be supplied. Denominator degrees of freedom can be computed by the Kenward-Roger, residual or asymptotic methods (or can be specified by users), and contrasts are specified through 'emmeans'. Vignettes cover estimated marginal means, contrast specification, and worked examples for common agricultural designs. |
| Authors: | Simon Riley [aut, cre] (ORCID: <https://orcid.org/0000-0002-4833-2026>) |
| Maintainer: | Simon Riley <[email protected]> |
| License: | BSD_3_clause + file LICENSE |
| Version: | 0.7.1 |
| Built: | 2026-10-22 10:18:11 UTC |
| Source: | https://github.com/simonshamusriley/powerutilities |
Data on barley yield for use in interaction contrast example.
barleybarley
barleyA data frame with 20 rows and 4 columns:
Factor. The replicate number
Factor. The nitrolime application rate.
Factor. The superphosphate application rate.
Numeric. The yield in kg/ha.
Saville, D. J., & Wood, G. R. (1991). Statistical Methods: The Geometric Approach. Springer Science & Business Media.
This is a helper function used for getting sigma() when it is not estimated
but fixed.
extract_disp(mod, ...)extract_disp(mod, ...)
mod |
A glmmTMB model. |
... |
Ignored. |
Numeric. The value of the (estimated or specified) residual dispersion in the model.
# Create synthetic data set: nrep = 8 # Number of replicates nfac = 2 # Number of treatments dat = expand.grid(Rep = factor(1:nrep), Trt = factor(LETTERS[1:nfac])) dat$Y = ifelse(dat$Trt == 'A', 0, 3) # "Fit" the model while setting residual standard deviation to 2: mod = set_glmm(Y ~ Trt, data = dat, disp = 2) # Extract residual SD: extract_disp(mod)# Create synthetic data set: nrep = 8 # Number of replicates nfac = 2 # Number of treatments dat = expand.grid(Rep = factor(1:nrep), Trt = factor(LETTERS[1:nfac])) dat$Y = ifelse(dat$Trt == 'A', 0, 3) # "Fit" the model while setting residual standard deviation to 2: mod = set_glmm(Y ~ Trt, data = dat, disp = 2) # Extract residual SD: extract_disp(mod)
Data on fodder beet yields for use in class contrasts example.
greenfeedsgreenfeeds
greenfeedsA data frame with 32 rows and 3 columns:
Factor. Beet variety.
Factor. The replicate number.
Numeric. Beet yield in tons/ha
...
Saville, D. J., & Wood, G. R. (1991). Statistical Methods: The Geometric Approach. Springer Science & Business Media.
Provides a method for outputting powertables as html
## S3 method for class 'powertable' knit_print(x, ...)## S3 method for class 'powertable' knit_print(x, ...)
x |
A |
... |
Unused. |
Invisibly returns x
This function calculates the power of contrasts.
power_contrast(emm, contr_list, ddf = NULL, alpha = 0.05, n_sims = 10000, ...)power_contrast(emm, contr_list, ddf = NULL, alpha = 0.05, n_sims = 10000, ...)
emm |
An emmGrid object associated with a |
contr_list |
A (named) list of contrast specifications. |
ddf |
Either a method for calculating denominator degrees of freedom ( currently supported options are "df.residual", "asymptotic" and "kenward-roger"), a numeric value, or NULL (the default), in which case a method is selected based on the model type. |
alpha |
Numeric. The nominal type I error rate. Defaults to 0.05. |
n_sims |
Numeric. The number of simulations to use for calculating type M error rate. If set to zero, the closed form, asymptotic calculations are used. |
... |
Other arguments passed to emmeans::contrast. |
A powertable object: a data frame with one row per contrast (or
per contrast within each level of any by variables), containing:
ContrastThe contrast name, taken from contr_list.
If emm has by variables, one column for each,
identifying the group the contrast was computed within.
EstimateThe expected value of the contrast, on the link scale.
SEIts standard error under the specified design and variance parameters.
LCL, UCL
Lower and upper confidence limits; present only if
confidence intervals are requested via infer.
NumDF, DenDF
Numerator (always 1) and denominator degrees of freedom for the test.
FvalThe expected F statistic (the squared t or z ratio).
FcritThe critical value of F at significance level alpha.
PvalThe expected p-value.
PowerThe probability of rejecting the null hypothesis of a zero contrast, given the expected estimate.
TypeSThe probability that a statistically significant estimate has the wrong sign (Gelman and Carlin, 2014).
TypeMThe expected ratio of the absolute value of a
statistically significant estimate to the true value (the
exaggeration ratio). Inf when the true contrast is zero.
The alpha and ddf attributes record the significance level and the
denominator degrees-of-freedom method used.
library(emmeans) # Create synthetic data set: nrep = 8 # Number of replicates nfac = 2 # Number of treatments dat = expand.grid(Rep = factor(1:nrep), Trt = factor(LETTERS[1:nfac])) dat$Y = ifelse(dat$Trt == 'A', 0, 3) # "Fit" the model while setting residual standard deviation to 2: mod = set_glmm(Y ~ Trt, data = dat, disp = 2) # Create emmeans object: emm = emmeans(mod, ~ Trt) # Define contrasts (here there is only one): contr = list('A-B' = c(1, -1)) # Calculate power of the contrast power_contrast(emm, contr)library(emmeans) # Create synthetic data set: nrep = 8 # Number of replicates nfac = 2 # Number of treatments dat = expand.grid(Rep = factor(1:nrep), Trt = factor(LETTERS[1:nfac])) dat$Y = ifelse(dat$Trt == 'A', 0, 3) # "Fit" the model while setting residual standard deviation to 2: mod = set_glmm(Y ~ Trt, data = dat, disp = 2) # Create emmeans object: emm = emmeans(mod, ~ Trt) # Define contrasts (here there is only one): contr = list('A-B' = c(1, -1)) # Calculate power of the contrast power_contrast(emm, contr)
This function calculates the power of F-tests for each fixed effect term in the model using one of several methods for determining the denominator degrees of freedom.
power_ftest(mod, ddf = NULL, alpha = 0.05, ...)power_ftest(mod, ddf = NULL, alpha = 0.05, ...)
mod |
A |
ddf |
Either a method for calculating denominator degrees of freedom ( currently supported options are "df.residual", "asymptotic" and "kenward-roger"), a numeric vector, or NULL (the default), in which case a method is selected based on the model type. |
alpha |
The nominal type I error rate. Defaults to 0.05. |
... |
Other values passed to emmeans. |
A powertable object: a data frame with one row per fixed-effect
term and columns Term, NumDF, DenDF, Fval (expected F statistic),
Fcrit (critical value at alpha), Pval (expected p-value) and Power.
Attributes alpha and ddf record the significance level and the
degrees-of-freedom method used.
# Create synthetic data set: nrep = 8 # Number of replicates nfac = 2 # Number of treatments dat = expand.grid(Rep = factor(1:nrep), Trt = factor(LETTERS[1:nfac])) dat$Y = ifelse(dat$Trt == 'A', 0, 3) # "Fit" the model while setting residual standard deviation to 2: mod = set_glmm(Y ~ Trt, data = dat, disp = 2) # Calculate power power_ftest(mod)# Create synthetic data set: nrep = 8 # Number of replicates nfac = 2 # Number of treatments dat = expand.grid(Rep = factor(1:nrep), Trt = factor(LETTERS[1:nfac])) dat$Y = ifelse(dat$Trt == 'A', 0, 3) # "Fit" the model while setting residual standard deviation to 2: mod = set_glmm(Y ~ Trt, data = dat, disp = 2) # Calculate power power_ftest(mod)
This function controls the formatting when printing powertable
objects.
## S3 method for class 'powertable' print(x, digits = 1, pdigits = getOption("pdigits", default = 4), ...)## S3 method for class 'powertable' print(x, digits = 1, pdigits = getOption("pdigits", default = 4), ...)
x |
A |
digits |
Integer. The number of digits to print for all columns other than the p-value and power. |
pdigits |
Integer. The number of digits to print for the p-value and power columns |
... |
Unused |
The input x (a powertable object), returned invisibly. Called
for its side effect of printing the table with p-values, power and
estimates rounded for display, followed by the degrees-of-freedom method
and significance level used.
Provides a method formatting and displaying retermslist objects
## S3 method for class 'retermslist' print(x, ...)## S3 method for class 'retermslist' print(x, ...)
x |
An |
... |
Unused. |
The input x (a retermslist object), returned invisibly. Called
for its side effect of printing one table per random-effects term, listing
the grouping factor, the name of each random effect, and the position at
which each standard deviation and correlation must be supplied in
re_terms when calling set_glmm().
This is a wrapper to help with fixing random effects terms in glmmTMB model for use in subsequent power analysis. In particular, it assists and converting random effect standard deviation and correlations into the log standard deviation and scaled log cholesky factors that glmmTMB expects.
set_glmm(formula, data, re_terms = NULL, disp = NULL, REML = TRUE, ...)set_glmm(formula, data, re_terms = NULL, disp = NULL, REML = TRUE, ...)
formula |
A two-sided model formula. |
data |
A data.frame. |
re_terms |
A vector of random effects standard deviations and, possibly, correlations. |
disp |
The value at which to fix the residual dispersion. |
REML |
Logical. Whether to use REML estimation (default) or, alternatively, ML. |
... |
Other arguments passed to glmmTMB. |
A glmmTMB object.
# Create synthetic data set: nrep = 8 # Number of replicates nfac = 2 # Number of treatments dat = expand.grid(Rep = factor(1:nrep), Trt = factor(LETTERS[1:nfac])) dat$Y = ifelse(dat$Trt == 'A', 0, 3) # "Fit" the model while setting residual standard deviation to 2: mod = set_glmm(Y ~ Trt, data = dat, disp = 2) # Confirm it worked: sigma(mod)# Create synthetic data set: nrep = 8 # Number of replicates nfac = 2 # Number of treatments dat = expand.grid(Rep = factor(1:nrep), Trt = factor(LETTERS[1:nfac])) dat$Y = ifelse(dat$Trt == 'A', 0, 3) # "Fit" the model while setting residual standard deviation to 2: mod = set_glmm(Y ~ Trt, data = dat, disp = 2) # Confirm it worked: sigma(mod)
When performing a "plug-in" power analysis, the values for the random effects must be supplied in the order in which they are saved internally in the glmmTMB model. It is not always clear (especially where there are correlated random slopes and intercepts) just what this ordering is. This helper function takes a model formula and a data set and outputs a table enumerating all the implied random effects and their ordering.
theta_finder(formula, data, ...)theta_finder(formula, data, ...)
formula |
A model formula (containing random effects). |
data |
A data frame containing all the columns included in |
... |
Other arguments passed to |
A 'retermslist' object numbering each of the random effect variances and correlations implied by the model formula.
library(dplyr) # Following the original analysis by Cornelius and Archbold, row x spacing # is treated as main-plots, and stock within main-plots as sub-plots. apple_des = agridat::archbold.apple |> mutate(across(rep:gen, factor), mp = interaction(row, spacing, sep = '-'), sp = interaction(row, spacing, stock, sep = '-')) # If we opt to treat blocks as random, then here is the # appropriate model for a split-split-plot RCBD: ssp_mod = ~ spacing*stock*gen + (1|rep/mp/sp) # The first theta value corresponds to to the standard deviation among # sub-plots, the second theta value corresponds to the of the standard deviation # main-plots and the third theta value corresponds to the standard deviation # among blocks: theta_finder(formula = ssp_mod, data = apple_des) # Purely for the sake of illustration, lets add hypothetical random slopes # for soil pH at subplot level. In the data, the values don't need to be # well simulated or even meaningful: apple_des2 = apple_des |> mutate(pH = 7) # Note that mp and sp are already explicitly defined in terms of their # nesting structure, so defining them here is terms of their interactions # isn't really necessary ssp_cov_mod = ~ spacing*stock*gen + (1|rep) + (1|rep:mp) + (pH + 1|rep:mp:sp) # Not one but two new theta values have been added to the model: we also # have the correlation between the pH slope and split-plot intercepts. theta_finder(formula = ssp_cov_mod, data = apple_des2) # Finally, be aware that redundant random effects terms, which would simply # be consolidated when fitting a real model, are problematic here. In this example, # rep is included twice: bad_mod = ~ spacing*stock*gen + (1|rep) + (1|rep/mp) try(theta_finder(formula = bad_mod, data = apple_des))library(dplyr) # Following the original analysis by Cornelius and Archbold, row x spacing # is treated as main-plots, and stock within main-plots as sub-plots. apple_des = agridat::archbold.apple |> mutate(across(rep:gen, factor), mp = interaction(row, spacing, sep = '-'), sp = interaction(row, spacing, stock, sep = '-')) # If we opt to treat blocks as random, then here is the # appropriate model for a split-split-plot RCBD: ssp_mod = ~ spacing*stock*gen + (1|rep/mp/sp) # The first theta value corresponds to to the standard deviation among # sub-plots, the second theta value corresponds to the of the standard deviation # main-plots and the third theta value corresponds to the standard deviation # among blocks: theta_finder(formula = ssp_mod, data = apple_des) # Purely for the sake of illustration, lets add hypothetical random slopes # for soil pH at subplot level. In the data, the values don't need to be # well simulated or even meaningful: apple_des2 = apple_des |> mutate(pH = 7) # Note that mp and sp are already explicitly defined in terms of their # nesting structure, so defining them here is terms of their interactions # isn't really necessary ssp_cov_mod = ~ spacing*stock*gen + (1|rep) + (1|rep:mp) + (pH + 1|rep:mp:sp) # Not one but two new theta values have been added to the model: we also # have the correlation between the pH slope and split-plot intercepts. theta_finder(formula = ssp_cov_mod, data = apple_des2) # Finally, be aware that redundant random effects terms, which would simply # be consolidated when fitting a real model, are problematic here. In this example, # rep is included twice: bad_mod = ~ spacing*stock*gen + (1|rep) + (1|rep/mp) try(theta_finder(formula = bad_mod, data = apple_des))