A Tutorial on Tailored Simulation-Based Sample-Size Planning for Experimental Designs With Generalized Linear Mixed Models

Abstract

When planning experimental research, determining an appropriate sample size and using suitable statistical models are crucial for robust and informative results. The recent replication crisis underlines the need for more rigorous statistical methodology and adequately powered designs. Generalized linear mixed models (GLMMs) offer a flexible statistical framework to analyze experimental data with complex (e.g., dependent and hierarchical) data structures. However, available methods and software for a priori sample-size planning for GLMMs are often limited to specific designs. Tailored data-simulation approaches offer a more flexible alternative. Based on a practical case study in which we focus on a binomial GLMM with two random intercepts and discrete predictor variables, the current tutorial equips researchers with a step-by-step guide and corresponding code for conducting tailored a priori sample-size planning with GLMMs. We not only focus on power analysis but also explain how to use the precision of parameter estimates to determine appropriate sample sizes. We conclude with an outlook on the increasing importance of simulation-based sample-size planning.

When planning experimental research, it is essential to determine an appropriate sample size and use appropriate statistical models to analyze the data to ensure that the results are robust and informative (Lakens, 2022b). The recent replication crisis in psychology and other disciplines has illustrated many challenges surrounding the reproducibility and reliability of study findings (Yarkoni, 2022). As a result, there is a growing need for more rigorous statistical methodology and the adoption of adequately powered experimental designs. Multiple easy-to-use software solutions exist for simple statistical models and experimental designs (Champely, 2020; Lakens & Caldwell, 2021). A priori sample-size planning for more complex research designs, such as flexible generalized linear mixed models (GLMMs), is not covered by standard software solutions. Researchers willing to use this framework will need to use data simulation. In the present work, we provide a tutorial on how to determine adequate sample sizes by performing tailored simulation-based sample-size planning for GLMMs. After introducing some theoretical background on sample-size planning, we review existing software solutions in R and discuss under which circumstances tailored data simulations are necessary. We proceed by describing the relevant steps and decisions involved in tailored data simulation, illustrated in a case study in which we focus on a binomial GLMM with two random intercepts and discrete predictor variables.

To benefit most from this tutorial article, we recommend readers familiarize themselves with basic statistical concepts such as hypothesis tests (HTs) and their statistical power, and confidence intervals (CIs) and their precision (Kumle et al., 2021; Lakens, 2022a; Riesthuis, 2024). Some knowledge of causal inference is beneficial but not necessary (Deffner et al., 2022; Lundberg et al., 2021). In addition, readers should have an understanding of how to conduct statistical analyses with R (Wickham et al., 2023) and how to simulate data (DeBruine & Barr, 2021; Hallgren, 2013; Lee et al., 2020). For data simulation, we use functions from tidyverse (Wickham et al., 2019) and the faux package (DeBruine, 2023). Finally, readers should have a basic understanding of regression modeling and GLMMs (Brown, 2021). In this tutorial, we simulate data by manually specifying the model equation of a GLMM that represents our assumed data-generating process (DeBruine & Barr, 2021). It is not necessary to understand the technical details of how GLMMs are estimated. However, it is crucial to understand the structure of a basic GLMM (e.g., logistic regression with random intercepts) and how the model assumes that the dependent variable’s values are determined by the predictor variables and the random effects.

Theoretical Background

Planning for statistical power or precision

Conducting research with insufficiently large sample sizes can have many negative consequences (Button et al., 2013). First, experiments may yield inconclusive or misleading results, hindering the accumulation of knowledge. Second, studies that are doomed never to find a postulated effect waste resources by consuming time, effort, and funding. For these reasons, many journals and funding bodies now require a sample-size justification in study protocols and grant proposals, recognizing its relevance in ensuring robust and meaningful findings. Although sample sizes can be justified with resource constraints or general heuristics, statistical arguments based on power or precision are the “gold standard” (Lakens, 2022b).

Most empirical studies in psychology and other social sciences apply hypothesis testing. Consequently, the dominant approach for determining an adequate sample size is based on power analysis (i.e., planning for power; Lakens, 2022b; Maxwell et al., 2008). Statistical power is defined as the probability that an HT has a significant p value when analyzing repeated samples from a population with a true effect of some prespecified size (Cohen, 1992). Less formally, power is described as the probability that an HT correctly rejects the null hypothesis when the alternative hypothesis is true. If the sample size (i.e., the number of participants and/or stimuli) is insufficient to detect the effects or relationships being investigated with high probability, the study is considered “underpowered.” When planning for power, a target is set for the statistical power of an HT of interest. Assuming an effect size of interest and a desired significance level, a minimum sample size can be determined that, on average, would guarantee reaching this target. The most prominent heuristic is to target a power of $1-\beta =0.8$ in combination with a Type I error rate of $\alpha =.05$ (Lakens, 2022b). However, depending on the research goals or resource constraints, there are often good reasons to move away from this standard (Benjamin et al., 2017; Lakens, Adolfi, et al., 2018). ¹

In contrast to power analysis, sample-size planning can also be based on the precision of parameter estimates (i.e., planning for precision or planning for accuracy; Lakens, 2022b; Maxwell et al., 2008). Not all research questions are best answered by hypothesis testing. It has been argued that basic research rarely requires discrete decisions on whether some effect has been “discovered” and should thus shift from hypothesis testing toward an estimation framework (Cumming, 2014; Kruschke & Liddell, 2018; McElreath, 2020), although this view is not without critique (Uygun Tunç et al., 2023). When no HTs are conducted, power analysis is not relevant for sample-size planning. In the precision framework, the target quantity commonly used for sample-size planning is the expected width of a CI (Kelley & Rausch, 2006; Lakens, 2022b; Maxwell et al., 2008). A CI with a confidence level of 0.95 provides the smallest interval with the property that 95% of individual CIs would include the true quantity of interest upon repeated sampling. Thus, a narrow CI with fewer plausible values for the quantity of interest is more informative about the size of the true effect than a wide CI. Apart from the confidence level, the width of a CI depends on the sample size. Because bigger samples carry more information, they lead to smaller CIs. When planning for precision, a target can be set for the expected width of a CI of interest. Assuming some effect size of interest and a certain confidence level, a minimum sample size can be determined that would guarantee reaching the targeted expected width. Because planning for precision is still rare, there are no common heuristics for choosing the desired width of the CI (Lakens, 2022b).

GLMMs

As study designs become more complex, psychological researchers require more sophisticated statistical models to capture their nuanced relationships and grouping structures (Yarkoni, 2022). GLMMs (also called “multilevel models”) are gaining popularity because they offer great flexibility when applied carefully (Brown, 2021; Matuschek et al., 2017; Meteyard & Davies, 2020). GLMMs are an extension of linear mixed models (LMMs), which are, in turn, extensions of linear regression models that account for correlated data, including hierarchical structures (Bolker, 2015; Fahrmeir et al., 2021). In this context, correlated data means that the value in the outcome variable for one observation may be related to the value of another observation in a systematic way that is not already accounted for by the usual (fixed) predictor variables (e.g., the age of participants). This correlation can arise for various reasons. For instance, responses to some stimuli from some participants might be more similar because the same person was measured multiple times (repeated measurements), participants belong to the same group (clustering), or participants responded to the same stimulus (stimulus effects). Thus, modeling such correlations is important whenever the data have a clear structure, whereas the grouping variables can be hierarchically nested (e.g., grouping variables students and schools: each student belongs to exactly one school) or cross-classified (e.g., grouping variables students and math exercises: each student is presented with several math exercises). LMMs are used when the outcome variable is continuous and follows a normal distribution (after conditioning on all predictor variables). They allow for the modeling of fixed effects, which capture the relationships between the usual predictors and the outcome, and random effects, which account for the different types of correlation structure and grouping effects. Random effects are typically assumed to follow a normal distribution with a mean of zero and a variance that quantifies the heterogeneity across groups. Correlated random effects can be assumed in models that contain both random intercepts and random slopes. GLMMs extend the LMM framework to accommodate nonnormally distributed continuous and categorical outcome variables. GLMMs involve a link function that connects the linear combination of predictor variables to the expected value of the outcome variable. The link function allows for modeling the relationship between predictors and the outcome in a nonlinear way that is appropriate for the specific distribution family of the outcome variable.

Simulation-Based Sample-Size Planning With GLMMs

To our knowledge, existing approaches for sample-size planning for GLMMs have exclusively focused on planning for power. In Table 1, we review available software packages that can be used to perform power analysis for multilevel models in R (Version 4.3.3; R Core Team, 2024). Power-analysis methods can be categorized into formula-based, summary-statistics-based, and simulation-based methods (Murayama et al., 2022). Formula-based methods rely on exact formulas to calculate power directly. Summary-statistics-based methods use statistical theory to approximate power based on formula-based methods developed for simple t tests. Simulation-based methods rely on repeatedly simulating data with a known true effect size and estimating power empirically, that is, what percentage of simulated data sets produces a significant p value. Available formula-based and summary-statistics-based software packages for multilevel models often do not include GLMMs or are limited to simple designs (Murayama et al., 2022; Westfall et al., 2014), making it necessary to build data simulations tailored specifically to the study design. A number of tutorials have been published describing how to perform such simulation-based power analysis for multilevel models (Arend & Schäfer, 2019; Brysbaert & Stevens, 2018; DeBruine & Barr, 2021; Green & MacLeod, 2016; Johnson et al., 2015; Kain et al., 2015; Kumle et al., 2021; Lafit et al., 2021; Zimmer et al., 2023). However, many of these tutorials focus on LMMs and the most common study designs (for a tutorial that also covers more advanced settings, see Kumle et al., 2021). This narrow focus provides limited guidance for researchers using more complex study designs, especially when little prior knowledge about plausible effect sizes is available (see the discussion in Kumle et al., 2021). Simulation-based power analysis with GLMMs requires making a range of assumptions about the model structure that should align with the characteristics of the data being analyzed. Existing tutorials often rely on heuristics for specifying variance components (e.g., the standard deviation of random intercepts) or assume that results from meta-analyses or data from pilot studies are available to determine plausible values for all model parameters. However, in practice, knowledge about those parameters from prior studies is often limited, which makes specifying assumptions a practical challenge (see the discussions in Maxwell et al., 2008, and Kumle et al., 2021). We discuss a number of strategies on how to specify model parameters for application-specific, tailored data simulations in a later section.

Table 1.

Overview of Software Packages for Power Analysis for (Generalized) Linear Mixed Models in R

Package name	Summary	GUI available	GLMM support	Design limitations	Reference
glmmrBase	Formula- and simulation-based power analysis; nonlinear fixed effects and flexible covariance functions	No	Yes	None	Watson (2023)
longpower	Formula-based power analysis; focus on longitudinal data	Yes	No	Limited set of study designs with two levels	Iddi & Donohue (2022)
mixedpower	Simulation-based Power analysis; based on lme4 package	No	Yes	None	Kumle et al. (2021)
mlmpower	Simulation-based Power analysis; missing data mechanisms	No	No	Limited to two levels	Enders et al. (2023)
pamm	Simulation-based power analysis; based on lme4 package	No	No	None	Martin et al. (2011)
pass.lme	Formula-based power analysis; limited documentation	No	No	None	Yu (2019)
PowerAnalysisIL	Simulation-based power analysis; focus on longitudinal data	Yes	No	Limited set of study designs with two levels	Lafit et al. (2021)
powerlmm	Simulation-based power analysis; missing data mechanisms; focus on longitudinal data	Yes	No	Limited set of study designs with two or three levels	Magnusson et al. (2018)
simglm	Simulation-based power analysis; missing data mechanisms	Yes	Yes	Limited to two or three levels	LeBeau (2019)
simr	Simulation-based power analysis; based on lme4 package	No	Yes	None	Green & MacLeod (2016)
summary_statistics_based_power	Summary-statistics-based power analysis; based on pwr package	Yes	Yes	Limited set of study designs with two levels	Murayama et al. (2022)
WebPower	Formula-based power analysis; based on longpower package; focus on longitudinal data	Yes	No	Limited set of study designs with two levels	Zhang & Yuan (2018)

Note: GUI available = whether a graphical user interface is provided (usually a Shiny app that is available online or can be downloaded to run locally); GLMM support = whether power analysis can be performed for generalized linear mixed model (otherwise only linear mixed models are supported); design limitations = whether only a limited selection of study designs can be specified (otherwise arbitrary designs within the model class are supported).

Reasons to Use Tailored Data Simulation

Performing tailored simulation-based sample-size planning is more complicated and time-consuming than using the existing software tools outlined in Table 1. Circumstances under which tailored simulation-based sample-size planning is necessary include (a) complex study designs, (b) complex statistical hypotheses, (c) planning for precision, and (d) lack of prior studies or pilot data.

First, requirements of real-world studies are often more complex than the simplified designs assumed by many user-friendly software packages for sample-size planning. One frequent issue in applied data analysis is missing data (Little & Rubin, 2014). Missings can be completely random (e.g., an electronic measurement device randomly failed for technical reasons), depend on some attributes also measured in the data set (e.g., older subjects refuse answering a question on income), or be caused by the measured attribute itself (e.g., wealthy people refuse to report their income). Moreover, many experimental designs contain conditions in which values of the predictor variables are missing by design. This can make data analysis more complicated because predictors have to be coded in specific ways that prevent the estimated GLMM from becoming unidentified. Whether missing data have an effect on the sample-size planning depends on the theoretical assumptions on how the missingness is caused. However, it is often challenging to decide whether missing data can be safely ignored in the data analysis and sample-size planning process based on a merely theoretical approach (Gomila & Clark, 2022). Tailored simulation-based approaches offer the possibility to include the assumed process of how data become missing in the data simulation, thereby determining the required sample size based on simulated data sets that contain missing values (e.g., see Lane & Hennes, 2018). As a by-product, the simulated data sets can also be used to test whether the intended data analysis provides the expected (unbiased) results despite the missing data. Although GLMMs can handle a large variety of outcome variables, researchers are becoming increasingly aware that many data sets might profit from even more sophisticated models. Common examples are zero-inflated outcomes, censoring, and nonlinear predictor effects that can be modeled with the R packages glmmTMB (Brooks et al., 2017) or brms (Bürkner, 2018). Tailored simulation-based approaches do not share the same limitations as the existing software solutions for power analysis that focus exclusively on GLMMs. As long as there is a software package available to estimate the model of interest, it is always possible to perform tailored simulation-based sample-size planning.

Second, the most common hypotheses tested in psychological research are of the type $H_0:\beta =0$ , where $\beta$ is a slope or intercept of a regression model. However, many research questions in psychology actually require testing more complex statistical hypotheses. In the new era of preregistration and Registered Reports (Chambers & Tzavella, 2022), most research questions should be tested with directed hypotheses because good theories at least postulate whether some psychological effect of interest is positive or negative. Even better theories should be able to specify the smallest effect sizes of interest (SESOIs) that must be exceeded if the effect has any practical relevance (Lakens, Scheel, & Isager, 2018). This might require a test such as $H_0:\beta \le 0.1$ . More elaborate research questions often require testing hypotheses that consist of a combination of model parameters, for example, testing simple slopes (Preacher et al., 2006) with a hypothesis such as $H_0:{\beta }_0+{\beta }_1\le 0$ . If the research question consists only of a single hypothesis of this sort, it might be possible to reduce the hypothesis to a single regression coefficient by clever coding and/or centering of predictor variables. However, research questions often consist of combined hypotheses that consist of more than one separate statistical hypothesis (for a tutorial on contrast analysis in GLMMs, see Schad et al., 2020). For example, a combined null hypothesis $H_0$ might consist of two single null hypotheses $H_{01}:{\beta }_1\le 0$ and $H_{02}:{\beta }_0+{\beta }_1\le 0$ . For some research questions, the combined null hypothesis $H_0$ would be rejected if both $H_{01}$ and $H_{02}$ are rejected. For other research questions, the combined null hypothesis $H_0$ would be rejected if $H_{01}$ or $H_{02}$ or both are rejected. If the global hypothesis $H_0$ is combined with a logical or, the p values of the single hypotheses must be corrected for multiple testing to avoid $\alpha$ inflation for the global hypothesis (Dmitrienko & D’Agostino, 2013). However, if the global hypothesis $H_0$ is combined with a logical and, a correction for multiple testing is not necessary but rather a mistake that unnecessarily reduces the power of the global HT. None of the software packages for sample-size planning in Table 1 can handle combined hypotheses, and only some can handle directed hypotheses. In contrast, our case study demonstrates how one can test directed combined hypotheses with tailored simulation-based sample-size planning.

Third, all available software packages for sample-size planning with GLMMs are based on power analysis and do not support planning for precision. Thus, researchers that want to apply an estimation strategy instead of testing statistical hypotheses (Cumming, 2014) currently cannot use the software packages outlined in Table 1. However, tailored simulation-based sample-size planning can easily handle the planning for precision approach (Maxwell et al., 2008). The only change in procedure is that instead of computing HTs for each simulated data set and estimating statistical power across repetitions, CIs are computed for each simulated data set, and their expected width is estimated.

Fourth, all frameworks for sample-size planning require the user to make assumptions about the expected effect size. Assuming the true effect is of this size (or greater), one can compute the (minimum) power of an HT or the (maximum) expected width of a CI. Existing software packages for sample-size planning for GLMMs usually require users to provide the assumed effect in the unit of some standardized measure of effect size. When the researcher has access to similar studies or pilot data, providing such standardized effect sizes is feasible. However, choosing effect sizes based on small pilot studies is generally not recommended because those estimates can be heavily biased (Albers & Lakens, 2018; Lakens, 2022b). Providing an informed standardized effect size can be an almost impossible challenge when no prior studies or pilot data are available. This problem is further exacerbated by the fact that GLMMs are so flexible that general heuristics of what should be considered a small effect do not exist or are difficult to defend. In the absence of prior evidence, using domain knowledge to construct a tailored data simulation is considered the best solution to determine plausible effect sizes (see Scenario 3 in Kumle et al., 2021). It would be possible to use these tailored simulations to compute standardized effect sizes that could then be inserted in existing software packages for sample-size planning. However, we would argue that when tailored data simulations are necessary to determine effect sizes anyway, performing the whole sample-size planning in a customized way is preferred over using the existing software packages.

General Steps in Tailored Simulation-Based Sample-Size Planning

Although the details differ for each study, every tailored simulation-based sample-size planning requires a series of steps and decisions. We introduce each step in a theoretical section, followed by the practical application based on a case study. To keep the tutorial manageable, we focus on a binomial GLMM with two random intercepts and discrete predictor variables. All code in this article and simulation results are available on OSF (https://osf.io/dhwf4/) and on our tutorial website (https://timo-ko.github.io/glmm_simulation_tutorial/).

Step 1: define the estimand

Theory

The first step in every research process is a clear definition of the theoretical estimand (Lundberg et al., 2021), that is, the theoretical quantity that is necessary to answer a specific research question. The estimand consists of a quantity that can be described for each unit under investigation and a clear definition of the target population for which the quantity is of interest. For example, an estimand might be the probability that a clinical psychologist makes the correct diagnosis for a psychiatric patient with major depression averaged across all clinical psychologists and depressed patients in psychiatric institutions in a given country.

The estimand should always be defined outside of any statistical model because there are usually a range of statistical methods that can be used to estimate the same estimand depending on the study design (e.g., a randomized experiment) that will produce the observed data in the planned study. For many common research questions in psychology, it is possible to express the estimand as a statistical quantity that can be estimated with a regression model, for example, a single $\beta$ coefficient. However, this is not possible for all estimands, which is why the literature discusses many estimation strategies beyond regression (Deffner et al., 2022; Lundberg et al., 2021).

Practice

In the present case study, we consider the effectiveness of feedback provided by an artificial intelligence (AI) embedded in a clinical-decision support system. The context is a clinical setting in which task experts (i.e., radiologists) and nonexpert (i.e., medical students) must detect bleedings based on head scans from computer tomography (CT). The AI model can provide initial diagnostic advice, which can be used as guidance by the humans, who are required to make the final diagnostic decision. The research goal is to validate the effectiveness of the AI-enabled advice. We consider the AI-enabled advice effective if the following pattern holds: We expect that for both radiologists and medical students, correct AI advice leads to a higher probability of accurately diagnosing a CT scan compared with no AI advice presented, and we expect that for both task experts and nonexperts, incorrect advice leads to a lower probability of accurately diagnosing a CT scan compared with no advice presented.

It becomes clear that our estimand consists of four comparisons between experimental conditions (Lundberg et al., 2021). However, the verbal description is still somewhat vague, which is why we try to give a more precise expression for each comparison:

$$\begin{array}{l}P\left(correctdiagnosis|correctadvice,averageexpert,averagescan\right)-\\ P\left(correctdiagnosis|noadvice,averageexpert,averagescan\right)\end{array}$$

$$\begin{array}{l}P\left(correctdiagnosis|noadvice,averageexpert,averagescan\right)-\\ P\left(correctdiagnosis|incorrectadvice,averageexpert,averagescan\right)\end{array}$$

$$\begin{array}{l}P\left(correctdiagnosis|correctadvice,averagestudent,averagescan\right)-\\ P\left(correctdiagnosis|noadvice,averagestudent,averagescan\right)\end{array}$$

$$\begin{array}{l}P\left(correctdiagnosis|noadvice,averagestudent,averagescan\right)-\\ P\left(correctdiagnosis|incorrectadvice,averagestudent,averagescan\right)\end{array}$$

For example, the first expression is concerned with the difference between the probability that a correct diagnosis is made if correct AI advice is presented and the probability that a correct diagnosis is made if no AI advice is presented. This contrast is quantified for a hypothetical typical expert and a typical scan, where typical is usually defined as an average score on all attributes of the expert or scan. ² To complete our definition of the estimand, we have to define our target population that consists of persons, stimuli, and AI advice: With respect to persons, we are interested only in experts (i.e., radiologists) and nonexperts (i.e., medical students). With respect to stimuli, we are interested only in the head CT scans made from subjects that do or do not suffer from intracerebral hemorrhage. Finally, we are interested only in AI advice given by a specific AI-enabled system.

Although the estimand is initially defined outside of a statistical model, it is useful only if we find a way to estimate it based on observed data. For our exemplary research question, it is possible to construct an experimental study in which all participants are confronted with the same set of head CT scans but the kind of AI advice given for each scan is randomly assigned within participants. This random intervention allows us to produce an empirical estimate of our estimand, although in reality, each person receives only one kind of AI advice (correct advice, incorrect advice, no advice at all) for each scan. We show later how each of the probability expressions in our estimand can be modeled with the same GLMM. Estimating this GLMM based on the data observed in our planned study will produce an estimate for each probability, and these estimates can then be combined to compute an estimate for each of the four probability contrasts. For pedagogical reasons, we skip the concrete definition of our estimate until we have discussed how to simulate data based on a concrete GLMM in the next section.

Step 2: simulate the data-generating process

Theory

When the estimand has been defined, the next step in the research process is to write code that simulates the data-generating process of the planned study. This requires specifying a generative process for all predictor variables used in the final data analysis. Realistic assumptions can be quite challenging for observational studies or continuous predictor variables, which is beyond the scope of this tutorial. However, this step is much easier for experimental studies with only categorical predictor variables because the distribution of predictors is fixed by the study design. When all predictor variables have been simulated, one can use the structure of a suitable GLMM to simulate the dependent variable. To simulate the GLMM, one requires plausible values for all model parameters. We discuss strategies on how these values can be obtained later. Because we have full control over the data-generating process in a tailored simulation, it is possible to model specific aspects of the planned study, such as data missing by design or assuming that some subjects might drop out. The quality of the results of the sample-size planning crucially depends on the plausibility of the simulated data-generating process. If the simulated data-generating process is less complex than the true process, one can expect that sample-size planning underestimates the required sample size and that the planned analysis has inflated Type I error rates (Matuschek et al., 2017). However, Matuschek et al. (2017) also showed in simulations with LMMs that although fitting maximal models (i.e., that include all possible random effects) safeguards against inflated Type I error rates, this can lead to a great loss in power if the variance in the random effects is actually small. Thus, we argue that even a simplified data-generating process (e.g., only a small number of interaction effects, only random intercepts and no random slopes, assuming that data are missing completely at random) that is plausible only under idealized circumstances can yield informative results and is preferred over performing no systematic sample-size planning.

Practice

In our case study, radiologists and students review a series of head CT scans to assess the presence of a bleeding. An AI model provides initial diagnostic advice to assist their decision-making. In the control condition, no AI advice is presented. When AI advice is given, this advice can be either correct or incorrect. The type of advice (no advice, incorrect advice, correct advice) is randomized within subjects across CT scans. After reviewing a CT scan, participants deliver a diagnosis (bleeding or no bleeding), which may be either accurate or inaccurate. This experimental design introduces some missing values by design because the advice is neither correct nor incorrect when no advice is present, which must be taken into account when simulating and analyzing the data. In this example, recruiting task experts (i.e., radiologists) is more challenging because of their limited availability, whereas nonexperts (i.e., students) are more readily accessible. The goal of simulation-based sample-size planning is to determine how many task experts and nonexperts must be recruited to achieve sufficient statistical power or precision in the planned experiment.

Our specific GLMM

In a GLMM, the expected value of the dependent variable $Y$ conditioned on the vector of predictor variables $\mathbf{X}$ and random effects $\mathbf{U}$ , transformed by a link function $g()$ is modeled as a linear combination $\eta$ of the predictor variables $\mathbf{X}$ , the random effects $\mathbf{U}$ , and the model parameters $\beta$ (Fahrmeir et al., 2021):

$$g\left(E\left(Y|\mathbf{X}=\mathbf{x},\mathbf{U}=\mathbf{u}\right)\right)=\eta \mathrm{.}$$

.

Equivalently, the conditional expected value is modeled as the linear combination η, transformed by the inverse link function $g^{-1}()$ :

$$E\left(Y|\mathbf{X}=\mathbf{x},\mathbf{U}=\mathbf{u}\right)=g^{-1}\left(\mathrm{\eta }\right).$$

.

If the dependent variable (i.e., diagnostic decision) $Y$ is a binary variable with values 0 (i.e., inaccurate) or 1 (i.e., accurate), the conditional expected value is equivalent to the probability:

$$P_{si}:=P\left(Y=1|\mathbf{X}=\mathbf{x},\mathbf{U}=\mathbf{u}\right).$$

.

In our case study, $P_{si}$ is the conditional probability that subject $s$ gives the correct response to item (i.e., CT scan) $i$ .

In such a setting, we model this probability as

$$P_{si}=inverse\_logit\left({\eta }_{si}\right),$$

with the inverse-logit link $g^{-1}\left({\eta }_{si}\right)=inverse\_logit\left({\eta }_{si}\right)=\frac{exp\left({\eta }_{si}\right)}{1+exp\left({\eta }_{si}\right)}$ or equivalently,

$$logit\left(P_{si}\right)={\eta }_{si},$$

with the logit link $g\left(P_{si}\right)=logit\left(P_{si}\right)=ln\left(\frac{P_{si}}{1-P_{si}}\right)$ .

In our case study, the probability of making an accurate diagnostic decision is assumed to depend on the following predictors:

• $advise\_presen{t}_{si}:\mathrm{whether}\mathrm{subject}\mathrm{was}\mathrm{presented}\mathrm{with}\mathrm{AI}\mathrm{advice}\left(\mathrm{1}\right)\mathrm{or}\mathrm{not}\left(0\right)\mathrm{when}\mathrm{asked}\mathrm{to}\mathrm{assess}\mathrm{item}\mathrm{i};$

• $advise\_correc{t}_{si}:\mathrm{whether}\mathrm{this}\mathrm{advice}\mathrm{was}\mathrm{correct}\left(\mathrm{1}\right)\mathrm{or}\mathrm{not}\left(0\right);$

• $exper{t}_s:\mathrm{whether}\mathrm{subject}\mathrm{was}\mathrm{a}\mathrm{task}\mathrm{expert}\left(\mathrm{1}\right)\mathrm{or}\mathrm{not}\left(0\right);$

and the following random effects:

• •  $u_{0s}:\mathrm{the}\mathrm{deviation}\mathrm{of}\mathrm{subject}s\mathrm{from}\mathrm{the}\mathrm{average}\mathrm{ability}\mathrm{to}\mathrm{solve}\mathrm{an}\mathrm{item}\left(\mathrm{i}.\mathrm{e}.,\mathrm{CT}\mathrm{scan}\right)$ with average difficulty; assumed to be distributed as $u_{0s}\sim N(0,{\sigma }_S^2)$ ;

• •  $u_{oi}:thedeviationofitem(i.e.,CTscan)ifromtheaveragedifficultytobesolvedbyapersonwithaverageability;assumedtobedistributedas{u}_{0}i\sim N(0,{\sigma }_I^2).$

In total, we assume the model

$$\begin{array}{l}logit\left[P_{si}\right]=({\beta }_0+u_{0s}+u_{0i})+{\beta }_a\cdot advice\_presen{t}_{si}\\ +{\beta }_c\cdot advice\_correc{t}_{si}+{\beta }_e\cdot exper{t}_s+{\beta }_{ea}\cdot exper{t}_s\cdot \\ advice\_presen{t}_{si}+{\beta }_{ec}\cdot exper{t}_s\cdot advice\_correc{t}_{si}\end{array}$$

or equivalently,

$$\begin{array}{l}{P}_{si}=inverse\_logit[\left({\beta }_0+u_{0s}+u_{0i}\right)+{\beta }_a\cdot advice\_presen{t}_{si}\\ +{\beta }_c\cdot advice\_correc{t}_{si}+{\beta }_e\cdot exper{t}_s+{\beta }_{ea}\cdot exper{t}_s\cdot \\ advice\_presen{t}_{si}+{\beta }_{ec}\cdot exper{t}_s\cdot advice\_correc{t}_{si}]\end{array}$$

with model parameters, ${\beta }_e$ , ${\beta }_a$ , ${\beta }_c$ , ${\beta }_{ea}$ , ${\beta }_{ec}$ , ${\sigma }_S$ , and ${\sigma }_I$ .

In the GLMM literature, this would be called a binomial GLMM with two random intercepts (for subjects and items), two Level 1 predictors ( $advice\_present$ , $advice\_correct$ ), one Level 2 predictor ( $expert$ ), and two cross-level interactions ( $expert\cdot advice\_present$ , $expert\cdot advice\_correct$ ). To limit complexity, we do not consider random slopes, additional predictors, or higher-level interactions.

Simulation function in R

The following R function simulates a full data set structured according to the design of our case study:

In the first line of the function definition, we set some default parameter values (which we explain in the next section) and load the packages we use to manipulate and simulate data. ³ In our case study, each subject (n_subjects in total) is assumed to respond to each item (i.e., CT scan; n_items in total). Thus, the add_random command creates a fully crossed data.frame with n_subjects × n_items rows. We add a between-subjects effect with the add_between command, simulating that about $25\%$ of subjects are experts. The next two lines simulate that in two-thirds of trials, subjects will be presented with AI advice, and if advice is presented, the advice will be correct in about $80\%$ of cases (the variable advice_correct is always 0 when no advice is presented). Next, we simulate one random effect for each subject (u0s) and for each item (u0i). As assumed by standard GLMMs, the add_ranef function draws the random effects from a normal distribution with a mean of 0 and a standard deviation specified by the user. With all design variables done, we are ready to simulate our model equation outlined in the last section. The linear predictor variable linpred (η in the GLMM model equations) combines the predictor variables, random effects, and model parameters as assumed by our model. We then transform the linear predictor with the inverse-link function to compute y_prob, the probability that the subject correctly solved the item (in R, the inverse-logit link is computed with plogis, and the logit link is computed with qlogis). In the final step, we simulate the binary dependent variable y_bin (i.e., whether the subject makes an accurate diagnostic decision for the CT scan) by—for each trial—drawing from a Bernoulli distribution with success probability y_prob.

Step 3: specify the population parameters

Theory

Population parameters are all model parameters estimated in a GLMM, in particular, the regression coefficients of the fixed effects and the standard deviation of the random effects (and the correlation between random effects in more complicated models). In the absence of previous studies with the same design or pilot data, strategies to specify population parameters will always require access to some domain knowledge. In our experience, research teams for projects in which tailored simulation-based sample-size planning is necessary typically consist of “analysts” (who are primarily responsible for study design and statistical analysis) and “domain experts” (who often provide the original research question and are responsible for data collection). In smaller projects, analysts and domain experts might be the same persons. In these research settings, a broad review of the available literature typically plays a more important role. Because most domain knowledge can be expressed only in unstandardized measurement units of a specific application, we argue that unstandardized effect sizes are usually preferable over standardized effect sizes for tailored simulation-based sample-size planning. The basic idea of all strategies is to quantify or visualize the data-generating process implied by certain values of population parameters in an intuitive way that enables calibration of population parameters based on the available domain knowledge. Although we use frequentist model estimation in our tutorial, many strategies demonstrated in this article are inspired by research on monitoring the plausibility of model assumptions (Gelman et al., 2020; Schad et al., 2021) and eliciting prior information (Bockting et al., 2024; Hartmann et al., 2020; Mikkola et al., 2023; Stefan et al., 2022) in applied Bayesian statistics.

Practice

When introducing the simulation function for our case study, we have used theoretically plausible values as defaults for all model parameters ( ${\beta }_0$ , ${\beta }_e$ , ${\beta }_a$ , ${\beta }_c$ , ${\beta }_{ea}$ , ${\beta }_{ec}$ , ${\sigma }_S$ , and ${\sigma }_I$ ) but have discussed where these numbers came from. The starting point for all parameter values was based on results from distantly related study designs in the literature. In addition, we had repeated discussions with our affiliated domain experts in radiology to check whether our assumptions regarding participants’ diagnostic performance seemed plausible.

We now outline our main strategy to determine plausible parameter values for the fixed effects ( $\beta$ parameters): Unfortunately, the model parameters in a binomial GLMM are hard to interpret in isolation because first, the parameters are connected to the modeled probability via the nonlinear inverse-logit link and second, we also have to consider the random effects. The most simple interpretation, which allows us to ignore the random effects for now, works by imagining a subject with average ability ( $u_{0s}=0$ ) responding to an item (i.e., CT scan) with average difficulty ( $u_{0i}=0$ ). Then, the model-implied probability that such a person solves such an item accurately is given by the following:

$$\begin{array}{l}P\left(Y=1|X=x,U=0\right)\\ =inverse\_logit[{\beta }_0+{\beta }_a\cdot advice\_presen{t}_{si}+{\beta }_c\cdot \\ advice\_correc{t}_{si}+{\beta }_e\cdot exper{t}_s+{\beta }_{ea}\cdot exper{t}_s\cdot \\ advice\_presen{t}_{si}+{\beta }_{ec}\cdot exper{t}_s\cdot advice\_correc{t}_{si}]\end{array}$$

In fact, we would need only the full equation if the subject is an expert and correct advice is presented. In all other experimental conditions, some terms drop from the equation because they are multiplied by 0. For example, the probability that a student with average ability solves an item with average difficulty when no advice is presented requires only the intercept:

$$\begin{array}{l}P\left(\begin{array}{l}Y=1|advice\_present=0,advice\_correct\\ =0,expert=0,u_{0s}=0,u_{0i}=0\end{array}\right)\\ =inverse\_logit\left[{\beta }_0\right].\end{array}$$

We can revert this perspective by choosing plausible probability values based on domain knowledge and deriving the parameter values implied by these probabilities for each experimental condition.

Table 2 shows our set of assumptions concerning the probability that an average subject solves an average item for each experimental condition and the corresponding equations implied by the model. Table 2 can be used to compute the implied values for the $\beta$ parameters, starting with the first equation and reinserting the computed $\beta$ values in all following equations (b_0 stands for the intercept ${\beta }_0$ , b_e stands for the slope ${\beta }_e$ , etc.):

Table 2.

Assumed Probabilities That an Average Subject Solves an Average Item in Each Experimental Condition

Experimental condition	$P\left(Y=1|X=x,U=0\right)$	Implied equation
No advice, student	0.70	$logit\left(0.70\right)={\beta }_0$
No advice, expert	0.90	$logit\left(0.90\right)={\beta }_0+{\beta }_e$
Incorrect advice, student	0.40	$logit\left(0.40\right)={\beta }_0+{\beta }_a$
Incorrect advice, expert	0.85	$logit\left(0.85\right)={\beta }_0+{\beta }_e+{\beta }_a+{\beta }_{ea}$
Correct advice, student	0.90	$logit\left(0.90\right)={\beta }_0+{\beta }_a+{\beta }_c$
Correct advice, expert	0.95	$logit\left(0.95\right)={\beta }_0+{\beta }_e+{\beta }_a+{\beta }_c+{\beta }_{ea}+{\beta }_{ec}$

Note: Implied equations are derived based on the model equations and setting all random intercept terms to 0.

It is always possible to double-check these computations by transforming the parameter values back to probabilities, for example,

$$\begin{array}{r}\hfill P\left(\begin{array}{l}Y=1|expert=1,advice\_present=1,\\ advice\_correct=1,u_{0s}=0,u_{0i}=0\end{array}\right)\\ \hfill =inverse\_logit\left[{\beta }_0+{\beta }_e+{\beta }_a+{\beta }_c+{\beta }_{ea}+{\beta }_{ec}\right]\end{array}$$

which we compute in R as

This leaves the question of how to determine plausible values for the two remaining model parameters ( ${\sigma }_S$ and ${\sigma }_I$ ) that are the standard deviations for the random intercepts. For this, we introduce two more strategies in the following sections.

Examine insightful descriptive statistics

Theory

The mathematical structure of GLMMs determines which patterns in data would be produced by the model if a specific set of values for the population parameters is specified. The knowledge of how to simulate from a GLMM enables us to compute insightful descriptive statistics that can be compared with available domain knowledge much more easily than the opaque values of model parameters. For example, domain experts might not be able to directly choose plausible values for the coefficients in a logistic regression model (which are measured on the log-odds scale). However, they should be able to reason about the expected ratio of the binary dependent variable in different experimental conditions, that is, which relative frequency they expect to observe. The job of the analyst who is familiar with the mathematical structure of the GLMM is to produce the model-implied value of the insightful descriptive statistic that is expected by the domain expert. Although insightful descriptive statistics usually depend on the model parameters in a nonlinear way, it is not necessary to solve the exact relationship mathematically. Instead, one can simply adjust the population parameters by trial and error until the model-implied quantities produce the desired result.

Practice

In the last section, we showed how we can derive the model-implied probability that a subject with average ability solves an item with average difficulty for each experimental condition. Although these derivations are straightforward, it is important not to misinterpret their implications: In binomial GLMMs, the average probability to solve an item (averaged across persons of varying ability and items of varying difficulty) is not equal to the probability that a person with average ability solves an item with average difficulty (Fahrmeir et al., 2021). The first perspective implies a so-called marginal interpretation, and the second one implies a conditional interpretation.

For example, we determined the $\beta$ parameters in a way that corresponds to a desired conditional probability of 0.95 that an expert with average ability solves an item with average difficulty when presented with correct advice (the conditional perspective). However, even if the model assumptions were true, the relative frequency of experts that solve items presented with correct advice would not be 0.95 in a large sample of subjects drawn from their natural distribution of ability and items drawn from their natural distribution of difficulty (the marginal perspective). How much the two probabilities differ depends on the standard deviations of the random intercepts (the two probabilities are equal only if both standard deviations are zero). We want to use the model-implied observed proportion of correct diagnoses in each experimental condition as an insightful descriptive statistic to determine plausible values for the random-effect standard deviations. We simulate a large data set (for which the observed values of the descriptive statistic will be close to their model-implied true values) and simply compute the relative frequency of correct diagnoses for each experimental condition:

We tried using these descriptive statistics to judge together with domain experts whether our chosen values for the random-effect standard deviations would produce data that aligned with their expectations. Although the result was deemed plausible, these statistics were not informative enough to determine a final set of plausible parameter values (e.g., doubling the standard deviations from 0.5 to 1 produces only minor changes in relative frequency). For this reason, we additionally look at insightful model-based quantities.

Examine insightful model-based quantities

Theory

Because GLMMs are complicated models, descriptive statistics alone are usually not enough to specify plausible values for all model parameters. This is especially true for the standard deviation of random effects that have nonlinear (and often unexpected) effects on the model-implied results. An important advantage of data simulation (in which one has full control over parameter values and sample sizes) is that one can produce insightful model-based quantities that can never be directly observed in an actual empirical data set. For example, in a logistic model with random intercepts for participants, one can produce a visualization of the implied distribution of the probability that a participant, on average, solves a cognitive task. Although domain knowledge will probably not suffice to specify this distribution completely, it should be possible to rule out implausible boundary conditions. For example, the domain expert might deem it implausible that the 5% most able participants have a probability of more than 0.99 to solve the difficult cognitive task.

Practice

The discussed inequality of conditional and marginal effects in GLMMs (Fahrmeir et al., 2021) makes their interpretation more difficult. One must be careful when specifying parameter values based on previous studies or pilot data that use the marginal interpretation (e.g., a pilot study providing an estimate of how often radiologists make an accurate diagnosis based on brain scans). However, this does not mean that we cannot use the marginal interpretation (average probability across persons and items) to inform plausible parameter values: When parameter values have been selected, we can compute the implied marginal distributions and compare this information with our domain knowledge. Then, we can iteratively adjust the parameter values until we are satisfied with the implied distributions. In the last section, we simulated a large data set and computed descriptive statistics, the relative frequencies of correct diagnoses, for each experimental condition. We now use the model-implied probability of each simulated data point (stored in the variable y_prob) to visualize the whole model-implied marginal distribution of correct diagnoses for each experimental condition:

Figure 1 shows the model-implied marginal distributions, including means and 66% and 95% confidence intervals. We can see that indeed, the average probabilities (black dots) slightly differ from the probabilities of average subjects and items considered previously. This difference increases with the variability of the random effects. We can use plots like Figure 1 as a useful tool to determine whether the specified standard deviations of the subject and item random intercepts ( ${\sigma }_S$ and ${\sigma }_I$ ) are reasonable by comparing the ranges and overlap between conditions to domain knowledge.

Fig. 1.

Marginal distributions including means and 66% and 95% confidence intervals for all experimental conditions.

For the next plot, we have set the item standard deviation to almost zero ( ${\sigma }_I=0.01$ ). This gives us a better way to see the variability between persons. As an example, Figure 2 reveals a number of implicit assumptions about the comparison between experts and students: With incorrect advice, virtually all experts have a higher probability of making a correct diagnosis compared with students when considering only items with average difficulty. In contrast, there is considerable overlap in probability between experts and students with no advice and even higher overlap with correct advice. Patterns such as these should be considered carefully and discussed with the domain experts. Parameter values ( $\beta$ parameters and ${\sigma }_S$ ) should be adjusted if the implications do not seem reasonable.

Fig. 2.

Marginal distributions including means and 66% and 95% confidence intervals for all experimental conditions with the standard deviation of item random intercepts set to 0.01.

The final plot demonstrates that these plots are also useful for spotting standard deviations that were specified too high. For Figure 3, we have set ${\sigma }_S=3$ and ${\sigma }_I=3$ . This implies that in each experimental condition, the probabilities that a subject solves an item are overwhelmingly close to either 0 or 1 and nothing in between, which is not a plausible assumption. These high standard deviations do not account for the inherent variability and complexity of human performance. For example, the expectation that an expert with low ability and incorrect advice would solve a difficult item with a probability close to zero is not convincing.

Fig. 3.

Marginal distributions including means and 66% and 95% confidence intervals for all experimental conditions with the standard deviation of subject and item random intercepts set to 3.

Iterate with domain expertise

Theory

Tailored simulation-based sample-size planning always requires access to domain knowledge. However, gathering domain knowledge and the specification of population parameter values is not always straightforward and can be better described as an iterative process in which the plausibility of simulated data is repeatedly compared with the available domain expertise. As a first step, dedicated domain experts (if available) can be interviewed to elicit their domain knowledge about how the data of the planned study are expected to look. Because most domain experts are no experts in statistical modeling and GLMMs, they often struggle without further guidance to communicate their knowledge in a way that is useful when specifying the parameters for data simulation. For this reason, we suggest that after an initial interview, the analyst who is familiar with the structure of the GLMM selects an initial set of insightful descriptive statistics and model-based quantities. Then, they reenter into an iterative discussion with the domain experts in which some set of population values are selected and the plausibility of resulting implied quantities is evaluated. The population parameters are updated until the domain experts are satisfied with the result. During this process, the monitored descriptive statistics and model-based quantities can be updated or extended to capture as much available domain knowledge as possible. Even if no dedicated domain experts are available, the basic principles stay the same except that initial parameter values are often derived from a more extensive literature review. Descriptive statistics and model-based quantities are still used to evaluate whether the implied assumptions about the data-generating process are plausible. Often, this judgment evolves over several iterations as the analyst becomes more familiar with the model and the study design. ⁴

Practice

All parameter values in our present case study have been determined based on repeated discussions with domain experts in radiology to validate our assumptions. Initially, we reviewed the literature to establish a reasonable baseline performance rate for examining head CT scans for intracranial hemorrhage. Existing studies indicate that radiologists typically demonstrate high accuracies around 90%, whereas interns have been shown to perform below 80%, and medical students fall even shorter. For simplicity, we assumed plausible probability values of .90 for experts and .70 for students, respectively. Our experts confirmed that these values are realistic baselines for reviewing diverse head CT images without AI assistance. Subsequently, we consulted several published articles that investigated the effect of correct and incorrect advice on decision-making performance in other settings. From their findings, we inferred that both experts and students should benefit from correct advice and suffer losses from incorrect advice. However, the magnitude of these effects should be substantially greater for students given their stronger reliance on advice compared with experts. We further validated the plausibility of our estimated gains and losses with the collaborating radiologists until we settled on the probabilities shown in Table 2.

Step 4: estimate the statistical model

Theory

At this point, the researcher is capable of producing a simulated data set that is comparable with the actual data set to be collected in the planned study. The next step is to specify how the statistical model shall be estimated in the actual study collected later. This requires the selection of a statistical framework, a software package that is capable of estimating the model class of interest, an estimation algorithm, and the specific model structure including all fixed effects, random effects, and the model family of the dependent variable.

Note that this does not always mean that one will specify the same GLMM that was used when specifying the data-generating process. On the one hand, using a simpler model for data simulation than for model estimation can be a useful strategy in scenarios in which making plausible assumptions for complicated random-effect structures and interactions is not feasible or the omitted components are expected to have small effects (Matuschek et al., 2017). On the other hand, using a more complex model for data simulation than for model estimation can be a useful strategy in scenarios in which one has specific domain knowledge about aspects of the data-generating process that are still difficult to estimate with the current state of the art in multilevel modeling. In fact, applying models that are more parsimonious than the true data-generating process can provide a better trade-off between Type I error rates (or bias) and power (or precision) in many scenarios (Matuschek et al., 2017).

Practice

In our case study, we use the lme4 R package (Bates et al., 2015), which is a state-of-the-art tool for fitting frequentist GLMMs. ⁵ For the current example, we simulate data according to our model, in which 100 subjects respond to 50 items (we use set.seed to make the simulation reproducible). However, for the sake of the exercise, we can imagine that this would be real data resulting from our future experiment and think about how we would analyze these data:

The lme4 package uses a special syntax for model specification. Our specific GLMM is represented by the following formula:

The first three lines look similar to any linear model in R (general intercept indicated by 1, main effects indicated by variable names in the data set, interactions indicated by variable1:variable2). The fourth line specifies a random intercept for each subject (1|subject) and for each item (1|item). The complete set of rules for the syntax is outlined in Bates et al. (2015) and in the documentation of the lme4 package.

In lme4, a GLMM is fitted with the glmer function. By setting family = "binomial", we request a binomial GLMM appropriate for our binary dependent variable y_bin (the binomial GLMM uses the canonical logit link by default), which is defined as an accurate (1) versus inaccurate (0) diagnosis. We use the default estimation algorithm (for a list of alternative options, see ?glmerControl):

We can inspect the estimates for all model parameters with the summary command:

In the model output, the Estimate column in the Fixed effects table contains the estimates for the $\beta$ parameters, and the Std.Dev. column in the Random effects table contains the estimates for ${\sigma }_S$ and ${\sigma }_I$ .

Step 5: compute the estimate

Theory

In previous steps, we have defined the theoretical estimand, written a data-simulation function, and specified how to estimate a GLMM using simulated data. The next step is to specify how to compute a concrete point estimate of the theoretical estimand within the framework of the fitted GLMM. For some research questions, the estimate corresponds with a single regression coefficient. In more complicated scenarios, the estimate is computed from a combination of coefficients. Beyond computing the point estimate, we have already discussed that both hypothesis testing and interval estimation could be used to answer the research question. The decision on testing or estimating is then followed by selecting the specific statistical method that shall be applied to compute the HTs or CIs (e.g., compute HTs and CIs with the marginaleffects R package using the delta method).

Practice

In the estimand section, we translated a verbal description of our research question into four probability statements that are specified outside of any specific statistical model. For a concrete estimate within the context of our specified GLMM, we must compute the following probability contrasts:

$$\begin{array}{l}P\left(Y=1|advice\_present=1,advice\_correct=1,expert=1,u_{0s}=0,u_{0i}=0\right)-\\ P\left(Y=1|advice\_present=0,advice\_correct=0,expert=1,u_{0s}=0,u_{0i}=0\right)\end{array}$$

$$\begin{array}{l}P\left(Y=1|advice\_present=0,advice\_correct=0,expert=1,u_{0s}=0,u_{0i}=0\right)-\\ P\left(Y=1|advice\_present=1,advice\_correct=0,expert=1,u_{0s}=0,u_{0i}=0\right)\end{array}$$

$$\begin{array}{l}P\left(Y=1|advice\_present=1,advice\_correct=1,expert=0,u_{0s}=0,u_{0i}=0\right)-\\ P\left(Y=1|advice\_present=0,advice\_correct=0,expert=0,u_{0s}=0,u_{0i}=0\right)\end{array}$$

$$\begin{array}{l}P\left(Y=1|advice\_present=0,advice\_correct=0,expert=0,u_{0s}=0,u_{0i}=0\right)-\\ P\left(Y=1|advice\_present=1,advice\_correct=0,expert=0,u_{0s}=0,u_{0i}=0\right).\end{array}$$

We have already discussed how to compute the involved probabilities in the section on specifying population parameters. Plugging in the GLMM model equation produces an equation for computing each contrast if all model parameters were known. When we want to estimate the above contrasts based on observed data, the only difference is that model parameters are not known, and we instead use the corresponding parameter estimates.

We could use our knowledge of the structure of the GLMM to determine the exact formula needed to compute the contrasts of interest and then plug in the parameter estimates manually from the summary(fit) output. However, this would be tedious, and we can use R to compute this contrast without doing the math. Using the first contrast, (correct advice, expert) – (no advice, expert), as our example, we could apply the predict function of the lme4 package to compute the predicted probability for a correct diagnosis based on our fitted model, plug in the two sets of predictor values, and compute the difference between the two estimated probabilities.

The argument type = "response" specifies that predictions are made on the probability scale (instead of the log-odds scale of the $\beta$ parameters), whereas re.form = NA sets all random effects to 0. We could use this method to compute point estimates for all four contrasts that are part of our estimand. However, depending on whether we are interested in hypothesis testing or parameter estimation, we also need a method to compute HTs or CIs. The marginaleffects package (Arel-Bundock et al., in press) is a very flexible, increasingly popular package to compute HTs and CIs for contrasts with a variety of statistical models, including GLMMs estimated with lme4. First, we specify a grid of all combinations of predictor variables and then compute estimated probabilities for all experimental conditions in our experiment with the predictions function:

Estimate	Std. Error	z	Pr(>|z|)	S	2.5 %	97.5 %	advice_present	advice_correct	expert
0.738	0.02134	34.6	<0.001	867.2	0.696	0.779	0	0	0
0.902	0.01739	51.9	<0.001	Inf	0.868	0.936	0	0	1
0.975	0.00421	231.6	<0.001	Inf	0.966	0.983	0	1	0
0.954	0.01454	65.6	<0.001	Inf	0.925	0.982	0	1	1
0.423	0.03221	13.1	<0.001	128.6	0.360	0.486	1	0	0
0.874	0.02793	31.3	<0.001	711.4	0.819	0.928	1	0	1
0.909	0.00967	94.1	<0.001	Inf	0.890	0.928	1	1	0
0.939	0.01091	86.1	<0.001	Inf	0.918	0.961	1	1	1
Type: response
Columns: rowid, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, advice_present, advice_correct, expert, y_bin

The point estimates for all experimental conditions are reported in the Estimate column. Note that the output also contains the two missing-by-design conditions that will never be observed in the actual study ( $advice\_present=0,advice\_correct=1,expert=1$ and $advice\_present=0,advice\_correct=1,expert=0$ ). This is no problem as long as we never interpret those estimates. Next, we use the estimated probabilities to compute the four specific contrasts that are part of our estimand. For this, we must specify which rows in preds have to be subtracted from each other. We use the hypotheses function to compute our four contrasts of interest together with HTs and CIs. We use the default inference options of the marginaleffects package that compute HTs and CIs based on the approximate delta method:

Term	Estimate	Std. Error	z	Pr(>|z|)	S	2.5 %	97.5 %	p (NonSup)	p (NonInf)	p (Equiv)
b8=b2	0.0374	0.0162	2.31	0.021	5.6	0.00563	0.0691	0.989	0.0105	0.989
b2=b6	0.0282	0.0279	1.01	0.312	1.7	–0.02653	0.0830	0.844	0.1562	0.844
b7=b1	0.1717	0.0173	9.93	<0.001	74.8	0.13780	0.2056	1.000	<0.001	1.000
b1=b5	0.3145	0.0280	11.24	<0.001	95.0	0.25965	0.3693	1.000	<0.001	1.000
Type: response
Columns: term, estimate, std.error, statistic, p.value, s.value, conf.low, conf.high, statistic.noninf, statistic.nonsup, p.value.noninf, p.value.nonsup, p.value.equiv

The expression "b8 = b2" is special syntax to subtract the estimate in Row 8 from the estimate in Row 2 in the preds-output. The argument equivalence = c(0, 0) can be used to compute one-sided p values, testing whether the contrast in the population is smaller than 0 (p (NonSub) column) or greater than 0 (p (NonInf) column). The point estimates for four contrasts are reported in the Estimate column. Note that to facilitate interpretation, we arranged the contrasts in a way that we theoretically expect positive values for all four of them.

Hypothesis testing

If we chose hypothesis testing for our case study, we would test a combined null hypothesis $H_0$ that consists of four separate null hypotheses:

$$\begin{array}{l}{H}_{01}:P\left(Y=1|advice\_present=1,advice\_correct=1,expert=1,u_{0s}=0,u_{0i}=0\right)\le \\ P\left(Y=1|advice\_present=0,advice\_correct=0,expert=1,u_{0s}=0,u_{0i}=0\right)\\ H_{02}:P\left(Y=1|advice\_present=0,advice\_correct=0,expert=1,u_{0s}=0,u_{0i}=0\right)\le \\ P\left(Y=1|advice\_present=1,advice\_correct=0,expert=1,u_{0s}=0,u_{0i}=0\right)\\ H_{03}:P\left(Y=1|advice\_present=1,advice\_correct=1,expert=0,u_{0s}=0,u_{0i}=0\right)\le \\ P\left(Y=1|advice\_present=0,advice\_correct=0,expert=0,u_{0s}=0,u_{0i}=0\right)\\ H_{04}:P\left(Y=1|advice\_present=0,advice\_correct=0,expert=0,u_{0s}=0,u_{0i}=0\right)\le \\ P\left(Y=1|advice\_present=1,advice\_correct=0,expert=0,u_{0s}=0,u_{0i}=0\right).\end{array}$$

The combined null hypothesis $H_0$ should be rejected only if all individual null hypotheses are rejected (i.e., intersection-union setting; Dmitrienko & D’Agostino, 2013). In such cases, the error probabilities do not accumulate, and we would waste power when correcting for multiple tests.

With a standard significance level of $\alpha =0.05$ , we would not reject all four null hypotheses (the p value in the p (NonInf) column for the second hypothesis is not significant) and therefore also not reject the combined null hypothesis for this particular (simulated) data set. Note that this decision would be wrong because we have simulated the data such that the combined alternative hypothesis $H_1$ is actually true in the population.

Interval estimation

If we chose parameter estimation for our case study, we would focus on the two-sided CIs of the four contrasts of interest. With a standard confidence level of $1-\alpha =0.95$ plausible values are clearly in the positive range for the first, third, and fourth contrasts, whereas both negative and positive values seem plausible for the second contrast. Note that because of the constrained range of the probability scale, the width of the CI differs between the four contrasts (which is the expected behavior for binomial GLMMs). The smallest width is observed for the first contrast (expert with correct advice vs. expert without advice), in which both underlying probabilities are close to 1. The largest width is observed for the fourth contrast (student with incorrect advice vs. student without advice), in which both underlying probabilities are closer to 0.5.

Step 6: perform repeated simulations

Theory

Conducting all previous steps enables the analyst to (a) simulate a data set, (b) estimate a GLMM, and (c) compute HTs or CIs for estimands of interest, mirroring the analysis that will later be conducted on the actual data set of the planned study. To produce an estimate of power or precision, the last step is to perform these previous steps repeatedly. On a conceptual level, we first require a function that takes a sample size and a full set of population parameter values as input. When planning for power, the function should return the p value(s) of the HT(s) of interest when conducted on the simulated data set. When planning for precision, the function should return the width of the CI(s) of interest. Second, we run this function repeatedly with the same sample-size and population parameters. Because fitting GLMMs can quickly become time-consuming, it is recommended to use parallel computing, that is, running simulations on multiple cores of the computer at the same time to reduce total run time. Third, the results of the repeated simulation are collected and aggregated. When planning for power, we compute the relative frequency of a significant p value(s) across repeated simulations. When planning for precision, we compute the average width of the CI(s). Finally, we repeat the complete simulation for different sample sizes to determine how big the sample must be to achieve the targeted power or precision.

Practice

Wrapping the simulate function already constructed earlier, the helper function sim_and_analyse performs all previous steps (simulate a data set, fit a GLMM, compute p values and CIs) in a single command:

We use the future (Bengtsson, 2021) and furrr (Vaughan & Dancho, 2022) packages to perform computations in parallel. First, we enable parallelization with the plan function and specify how many parallel cores (“workers”) of our computer to use (users can find out the maximum number of cores on their computer with the command parallel::detectCores()) and set a seed to make the simulation reproducible:

The next code chunk specifies a simulation grid with different settings for both the number of subjects (n_subjects) and the number of items (n_items), each combination being repeated rep times. The more repetitions, the more accurately power and precision can be estimated. We chose 300 repetitions for the data simulation at hand to strike a balance between achieving a robust estimate and remaining computationally feasible. With the current settings, this simulation takes several hours on a MacBook Pro from 2020 with an M1 chip and 16 GB working memory. If you want to quickly experiment with the code yourself, a setting with workers = 4 and rep = 5 should finish in less than 5 min, even on smaller machines:

The result of this computation is a data frame that contains the p values and CIs of all specified contrasts for each simulated data set. In some iterations (predominantly in conditions with small sample sizes), model estimation did not converge with the lme4 package. When the model fails to converge, it means that the statistical model being fitted to the data failed to reach a stable or valid solution during the estimation process. We do not remove these results because nonconvergence can also happen when analyzing the real data we plan to collect; thus, we want to factor in this possibility to keep our simulation more realistic.

Power results

For our exemplary combined hypothesis, power is defined as the (long-run) percentage of simulations in which all four p values of our individual hypotheses are significant at the $\alpha =.05$ level. Based on our simulation outcomes, we compute a power estimate for each combination of n_subjects × n_items (including 95% CIs) and visualize the results with the following code: ⁶

As should be the case, power estimates in Figure 4 increase with both the number of subjects and the number of items. The CIs reported here indicate how precisely power was estimated by our simulation. Higher precision (which would be reflected in narrower CIs) could be obtained by increasing the number of repetitions (rep) in the simulation. In practice, data simulations are often run multiple times with adjusted combinations of sample sizes. When running for the first time, it might be revealed that power is way too low (or much higher than required) for some combinations of n_subjects and n_items. When narrowing down the best combination that achieves sufficient power while at the same time striking a good balance of how many subjects and items are practically feasible, later rounds of data simulation will typically include a smaller grid of sample sizes combined with a higher number of repetitions. This will assure high precision for the final power estimates, which are then used for the sample-size justification of the planned study.

Fig. 4.

Simulation-based power estimates including 95% confidence intervals of the case study for different numbers of subjects and items, based on a significance level of .05.

When target power has been specified, the number of subjects and the number of items in our study design can be traded against each other based on practical considerations. Although it is recommended to justify the power (Lakens, 2022b), we adopt the common heuristic to target a power of 0.8 to detect an effect of the expected size implied by our data simulation. This could be achieved by collecting data from 250 subjects (about 25% of which will be experts), each completing the same 50 items (with advice present in about 67% of cases, which is correct in about 80% of cases with present advice). If collecting data from 250 subjects is not feasible, an alternative would be to recruit 200 subjects but increase the length of the experiment to 70 items. However, 70 items might take too long to complete for the radiologists participating in the study, who have a busy schedule. The simulation suggests that it might also be possible to plan a shorter experiment with only 30 items if it is feasible to recruit an even higher number of subjects (to be determined by additional rounds of power analysis). Design parameters that also affect power and could be investigated in the simulation to find a more optimal trade-off are the ratio of experts to students, the frequency of whether advice is presented at all, and whether it is correct or incorrect. ⁷

Precision results

When planning for precision, one could monitor the width of all four CIs at the same time. However, because the CIs of the four contrasts strongly differ in width, it is not trivial to decide which width one should target when deciding on the appropriate sample size. In contrast to planning for power, there are no common standards on how to specify the targeted precision (Lakens, 2022b). For our example, we use a simple heuristic, but we strongly encourage readers to think about better alternatives that are appropriate in their own applications. Our simulations show that the smallest CI can be expected for the first contrast (expert with correct advice vs. expert without advice). The true contrast in probability for an average expert and an average item in this condition is plogis(b_0 + b_e + b_a + b_c + b_ea + b_ec) – plogis(b_0 + b_e) = 0.05. We want the width of this CI to be smaller than 0.1. This would mean that if the point estimate happens to be close to the true value, the plausible values inside of a 95% CI would all be positive.

Thus, in our example, precision is defined as the (long-run) average width of a 95% CI for the probability contrast between experts with correct advice and experts without advice. Of course, a lower width implies better precision. Based on our simulation outcomes, we compute the precision estimate for each combination of n_subjects × n_items (including 95% CIs) and visualize the results with the following code:

As should be the case, precision estimates in Figure 5 increase (i.e., average width of CI decreases) with the number of included subjects and items. The CIs reported here indicate how precisely the expected width of the CI for our focal contrast was estimated by our simulation. Applying our simple heuristic of targeting an expected CI width that is smaller than 0.1, we see the same trade-off between the number of subjects and the number of items as with planning for power. We could either choose 100 subjects and 30 items or 250 subjects and 10 items. Note that our simple heuristic for determining sample size in the planning for precision scenario was quite liberal. This is reflected in the result that we would need a smaller sample size than in the planning for power scenario. With a more conservative precision target, the result is generally the opposite: As a rule, precise parameter estimates usually require bigger samples than null hypothesis testing.

Fig. 5.

Simulation-based precision estimates (expected width of confidence intervals) including 95% confidence intervals of the case study for different numbers of subjects and items, based on a confidence level of .95.

Step 7 (optional): conduct sensitivity analysis

In our case study, we have performed simulation-based sample-size planning from a single set of parameter values that reflect our assumptions of an expected effect size. Instead of extracting this expected effect size from meta-analyses or pilot data, which has been the main focus of previous tutorials (e.g., Kumle et al., 2021), we have demonstrated some strategies to determine plausible parameter values in GLMMs based on domain knowledge. When sample sizes are chosen based on the results of our simulation-based power analysis, a future study will be informative to reject the null hypothesis if an effect of our expected size is present (or estimate the effect with satisfying precision). However, if the true effect is indeed smaller, power (or precision) will be lower, and the study might not be sufficiently informative. A common, more conservative strategy for sample-size justification is to perform sample-size planning for the SESOI. An effect smaller than the SESOI would be considered too small to be interesting or practically meaningful even if the effect is not actually zero (King, 2011; Lakens, Scheel, & Isager, 2018). Unfortunately, specifying a plausible SESOI is a challenging task (for strategies on how to specify a SESOI, see Lakens, 2022a or 2022b; Riesthuis, 2024). When domain knowledge or formal theories about the research topic of interest are too vague to specify a meaningful SESOI, it is still recommended to demonstrate power or precision for different effect sizes in a sensitivity power analysis (Lakens, 2022b). By simulating power (or precision) for different effect sizes (in addition to the different number of subjects and items), one can make sure that power (or precision) would still be sufficient to detect smaller effect sizes than the expected effect or at least get an impression of how strongly power (or precision) depends on the size of the true effect. For our case study that investigates combined hypotheses in a GLMM modeling framework, sensitivity analysis would require manually specifying additional sets of plausible parameter values that reflect scenarios with smaller or larger differences between groups with respect to our specific research question. Power (or precision) could then be simulated for several of these scenarios (across different numbers of subjects and items, as considered earlier). In Steps 2 and 4 of the tutorial, we briefly discussed scenarios in which the applied statistical model is more (or less) parsimonious than the data-generating process. Sensitivity analysis can be used to assess the consequences of a mismatch by investigating different combinations of data-generating processes and statistical models. Such extended simulations can be a great resource to make an informed decision about the concrete model to estimate for the planned study in which the true data-generating process is unknown. Although recommended, we do not present a sensitivity analysis for our case study to keep the length of the tutorial manageable.

Conclusion and Outlook

In this tutorial, we provided a step-by-step guide on how to perform tailored simulation-based sample-size planning for GLMMs based on a concrete case study. To conclude, we want to give an outlook on five developments regarding the future role of simulation-based sample-size planning in experimental research.

First, as experimental designs become more complex and the appropriate flexible statistical frameworks, such as GLMMs, become more popular, there is an ever-growing need for simulation-based sample-size planning in experimental research.

Second, the ability to conduct simulation-based sample-size planning becomes an increasingly valuable skill that should be taught to experimental researchers. By incorporating such training into research methods courses and workshops, researchers can improve the quality of their experimental designs and enhance the rigor of their studies. The need to reason about how to simulate plausible data that is in line with the research hypothesis while not violating domain expertise on how plausible data should look might also contribute to planning more insightful studies that can answer more precise research questions (Yarkoni, 2022).

Third, given the significant disconnect between the amount of effort required to perform tailored simulation-based sample-size planning and the perceived effort estimated by researchers and collaborators in experimental research, there is a need to address the mismatch in effort perception. Many researchers request simulation-based power analyses from statisticians or methodological experts without fully comprehending the complexity and time-consuming nature of these tailored simulations. Therefore, it is crucial to raise awareness about the effort involved to ensure realistic expectations and effective collaboration with methodological experts.

Fourth, tailored data simulations and power analyses are not mere technicalities; they are valuable research contributions that deserve adequate recognition in experimental research. Their importance can be reflected by highlighting the simulation work in a publication or even allocating them a separate publication or incorporating them as a significant component of Stage 1 preregistered reports (Chambers & Tzavella, 2022).

Finally, simulation-based sample-size planning aligns well with the principles of open science and preregistration and could be further integrated. When researchers have access to simulated data based on their prespecified model, analyzing the collected data set becomes straightforward and unambiguous. By preregistering their simulation-based sample-size plan, researchers enhance the transparency and accountability of their experimental procedures, contributing to the credibility and reproducibility of research.