A practical guide and power analysis for GLMMs: detecting among treatment variation in random effects

Linear mixed model

We begin by introducing a general linear mixed model (LMM) to illustrate the variance components we are interested in (Fig. 1) and their applications in behavioral ecology. We provide only a brief introduction to LMMs here because they have been extensively discussed in several recent reviews and textbooks (Gelman & Hill, 2006; Zuur et al., 2009; Stroup, 2012; Zuur, Hilbe & Leno, 2013; Dingemanse & Dochtermann, 2013; Bates et al., 2014; Bolker, 2015). We use the notation of Stroup (2012) to facilitate a transition to the binomial GLMM model, which is the focus of our power analyses.

A two treatment linear mixed model can be written as: (1)${\displaystyle y_{ijk}|b_{0ik},b_{1ik}\sim \mathrm{Normal}\left({\mu }_{ijk},{\sigma }_{\epsilon k}^2\right)}$ (2)${\displaystyle {\eta }_{ijk}={\beta }_{0k}+b_{0ik}+\left({\beta }_{1k}+b_{ik}\right)X_{ijk}}$ (3)${\displaystyle \text{Identity link: }{\eta }_{ijk}={\mu }_{ijk}}$ (4)${\displaystyle \left[\begin{array}{c}\hfill b_{0ik}\hfill \\ \hfill b_{1ik}\hfill \end{array}\right]\sim \text{MVN}\left(\left[\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \end{array}\right],\left[\begin{array}{cc}\hfill {\sigma }_{0k}^2\hfill & \hfill {\sigma }_{01k}\hfill \\ \hfill {\sigma }_{0lk}\hfill & \hfill {\sigma }_{1k}^2\hfill \end{array}\right]\right).}$ Here, a single phenotypic measurement y_ijk is of individual i, at level j of the covariate X (in studies of animal behavior the covariate of interest is often an environmental gradient) in treatment k. This model is composed of three components: the treatment mean in environment j (β_0k + β_1kX_ijk), the unique average response of individual i across the environmental gradient (b_0k + b_1kX_ijk), and a residual error due to the variation around the mean of individual $i\left({\sigma }_{\epsilon k}^2\right)$, which is assumed to be homogenous across X and among all individuals in treatment k, but is allowed to vary by treatment. Individuals vary from the treatment mean reaction norm in both their intercept (b_0ik) and slope (b_1ik), which together compose the total phenotypic variance attributable to among-individual variation. This individual contribution is quantified using a random intercepts and slopes model with a multivariate Normal (MVN) distribution (4). Variation among individuals in intercept and slope are ${\sigma }_{0k}^2$ and ${\sigma }_{1k}^2$ respectively; covariance between intercept and slope is given by σ_01k. In a LMM, the linear predictor directly predicts the mean, as shown by the identity link function in Eq. (3). In a GLMM, the linear predictor predicts a function of the mean g(x), which must be linearized through the use of non-identity link functions; for example, we use the standard logit (log-odds) link for Binomial GLMM.

Binomial GLMM

We assess power of a binomial GLMM for detecting differences in variation by treatment. This model can be written as: (5)${\displaystyle y_{ijk}|b_{0ik},b_{1ik},v_{ijk}\sim \text{Binomial}\left(N_{ijk},{\pi }_{ijk}\right)}$ (6)${\displaystyle {\eta }_{ijk}={\beta }_0+b_{0ik}+\left({\beta }_1+b_{1ik}\right)X_{ijk}+v_{ijk}}$ (7)${\displaystyle \text{Inverse-logit: }{\pi }_{ijk}=1/\left(1+e^{-\eta }ijk\right)}$ (8)${\displaystyle \left[\begin{array}{c}\hfill b_{0ik}\hfill \\ \hfill b_{1ik}\hfill \end{array}\right]\sim \text{MVN }\left(\left[\begin{array}{c}\hfill 0\hfill \\ \hfill 0\hfill \end{array}\right],\left[\begin{array}{cc}\hfill {\sigma }_{0k}^2\hfill & \hfill {\sigma }_{01k}\hfill \\ \hfill {\sigma }_{01k}\hfill & \hfill {\sigma }_{1k}^2\hfill \end{array}\right]\right)}$ (9)${\displaystyle v_{ijk}\sim \text{Normal }\left(0,{\sigma }_{vk}^2\right).}$

Here, y_ijk is the number of “successes” in N_ijk observations of the ith individual in treatment k at the jth sampling occasion. When an environmental covariate (X) is present, we assume one sampling occasion occurs at each level of the covariate j. Here, in the absence of an environmental covariate, the linear predictor reduces to η_ijk = β₀ + b_0ik + v_ijk and the jth occasion is simply a repeated sampling occasion in the same conditions. Note, when N_ijk = 1 there is only 1 observation per sampling occasion j, making y_ijk a Bernoulli response variable (see Supplemental Information 1). When y_ijk is Bernoulli, overdispersion (v_ijk) and thus within-individual variation is not identifiable.

In this model π_ijk describes the underlying probability of individual i in treatment k at occasion j exhibiting a behavior. Variation in π isdetermined by the linear combination of predictors on the logit (log-odds) scale: group intercept (β₀), group slope (β₁), individual unique intercept (b_0ik), slope (b_1ik), and observation level overdispersion that decrease predictive power at each observation (v_ijk). This linear predictor is transformed with the inverse-logit link to produce π_ijk, which follows a logit-Normal-binomial mixed distribution.

We use an observation-level random effect to model additive overdispersion (Browne et al., 2005), which models increased variance (following a Normal distribution with variance ${\sigma }_{vk}^2\left.\right)$ in the linear predictor on the link scale (Nakagawa & Schielzeth, 2010). Overdispersion is used to quantify within-individual variation because it models variation in π between each sampling occasion j for each individual. Here the magnitude of overdispersion is allowed to vary by treatment (for an example of multiple data sets where this occurs see Hinde & Demétrio, 1998), which is a focus of our power analysis.

The transformation through the inverse-logit function makes each of the three target variance components difficult to visualize with a concise figure. However, because the binomial GLMM model follows similar patterns as the LMM, we present power analyses for the binomial GLMM using the visual aid presented for the LMM (Fig. 1). Finally, we simulate data for a fully balanced design without losing generality. See Martin et al. (2011) and Van de Pol (2012) for a discussion on experimental designs where individuals are assayed in partially overlapping environments and when only single measurements are obtained for some individuals.

Simulations

All data were simulated in the R statistical programming environment using newly developed simulation capabilities of the lme4 package (>version 1.1–6, Bates et al., 2014). Guidelines for parameterizing the GLMMs and running data simulations and power analyses are provided in Supplemental Information 1. For a given total sample size, we present simulations for determining the optimal ratio of total number of individuals versus the number of repeated measures within individuals needed to provide power to detect a difference among treatments 80% of the time. We conducted simulations for multiple ratios of individuals to total observations within individuals, varying both sampling occasions (j) and Bernoulli observations within sampling occasions (n). Next, we describe simulations that evaluate how increasing “noise” (variation in non-target random effects) affects power to detect differences in targeted variance comparisons.

For both scenarios we simulate data with biologically relevant parameter values that illustrate common trends in power. At extreme parameter values the trends presented here may not hold due to interactions between the variance components that arise at the boundaries of binomial space. We do not dwell on these exceptions since they are unrealistic for most empirical data sets, but suggest exploration of these exceptions with code provided in Supplemental Information 1.

We ran 2,800 simulations for each combination of parameter values. The significance of a given random effect was assessed using likelihood ratio tests (LRTs) between models with and without the focal random effect. To correct for the known conservatism of the LRT when testing for σ² = 0 (due to a null value on the boundary of parameter space), we adopted the standard correction of dividing all p-values by 2 (Pinheiro & Bates, 2001; Verbeke & Molenberghs, 2000; Fitzmaurice, Laird & Ware, 2004; Zuur et al., 2009). This correction was appropriate for all p-values because each LRT compared models that differed in only a single degree of freedom. Power is estimated as the percentage of simulations that provide a corrected p-value smaller than 0.05. We insured the validity of a nominal p-value of 0.05 by confirming that 2,800 simulations of a scenario with equivalent standard deviations in both treatments did not result in rejecting the null hypothesis more than 5% of the time. Under extremely low numbers of individuals (∼2–4) power to detect differences in the null case exceeded 5% (∼10–15%), possibly inflating power in these cases. Regardless, random effects cannot be reliably estimated with such low sample sizes and therefore in most cases such experimental designs should be avoided.

Scenario 1: Determining the optimal sampling scheme

Most researchers face limitations imposed by time, money and access to samples, and are therefore confronted with the question of how resources should be divided between individuals and measures within individuals. To investigate the optimal allocation of sampling effort between the number of individuals and number of observations per individual, we simulated two data sets for each variance comparison (see Table 1 for a summary of all simulations).

First, using three hypothetical total numbers of Bernoulli observations per treatment (total sample size per treatment, TSS_T), we manipulated either the ratio of individuals to sampling occasions (${\sigma }_{0k}^2$ and ${\sigma }_{1k}^2$), or the ratio of individuals to Bernoulli observations within sampling occasions $\left({\sigma }_{vk}^2\right)$. For comparisons of ${\sigma }_{0k}^2$ and ${\sigma }_{1k}^2$ we manipulated the ratio of individuals to sampling occasions, holding the number of Bernoulli observations constant at 5, because power follows a non-monotonic pattern across these ratios for ${\sigma }_{0k}^2$ and ${\sigma }_{1k}^2$(Figs. 2 and 3). Conversely, for comparisons of ${\sigma }_{vk}^2$ we manipulated the ratio of individuals to Bernoulli observations and held the number of sampling occasions constant at 5 because power follows a non-monotonic pattern across ratios of individuals to Bernoulli observations for ${\sigma }_{vk}^2$ (Fig. 4). For comparisons of ${\sigma }_{0k}^2$, and ${\sigma }_{vk}^2$ we simulated TSS_T of 600, 1,200 and 2,400, and for comparisons of ${\sigma }_{1k}^2{\mathrm{TSS}}_T$ were 300, 600, and 1,200. For example, for b_1ik with a TSS_T of 300, the most extreme ratios were 30 individuals with 2 sampling occasions and 2 individuals with 30 sampling occasions. While using only 2 samples for a grouping variable (individuals) is never suggested for a random effect, we include this combination as an illustration of the low power that results from an ill-conceived sampling scheme. For each variance comparison we simulated three different effect sizes (2, 2.5, and 3 fold difference in standard deviation by treatment).

Next, we simulated data sets with increasing numbers of Bernoulli observations for comparisons of ${\sigma }_{0k}^2$ and ${\sigma }_{1k}^2$ (Figs. 5A and 5B) and with increasing numbers of sampling occasions for comparisons of ${\sigma }_{vk}^2$ (Fig. 5C). For these simulations we used 1, 3, 5, 10 and 15 Bernoulli observations or sampling occasions. Ratios of individuals to sampling occasions (${\sigma }_{0k}^2$ and ${\sigma }_{1k}^2$) or individuals to Bernoulli observations (${\sigma }_{vk}^2\left.\right)$ followed the intermediate TSS_T from the simulations described above. For example, for comparisons of ${\sigma }_{0k}^2$ we simulated 1, 3, 5, 10 and 15 Bernoulli observations for ratios of individuals to sampling occasions ranging from 120:2 to 2:120. For all comparisons we simulated data using an effect size of a 2.5 fold difference in standard deviation by treatment.

In all Scenario 1 simulations, both β₀ and β₁ were constrained to a single value for all treatments. For comparisons of among-individual variation in intercept no environmental covariate was used causing each sampling occasion to occur in the same conditions. Additionally, ${\sigma }_{vk}^2$ was held constant among treatments. For comparisons of among-individual variation in slope we held ${\sigma }_{vk}^2$ constant. Finally, for comparisons of within-individual variation in intercept, no environmental covariate was included and ${\sigma }_{0k}^2$ was held constant among treatments. All parameter values used in simulations for both Scenarios can be found in Table S1.

Our goal in Scenario 1 was to isolate changes in a single variance parameter, but exploration of the dependence among multiple variance components and the mean may be warranted if it is relevant for a specific problem. Incorporating concurrent changes in intercept, slope and overdispersion parameters can be easily implemented with slight modifications to the code presented in the online supplement. We show initial results of relaxing some of these assumptions in Scenario 2, but full exploration of these possibilities are beyond the scope of this paper.

Scenario 2: Measuring the ratio of overdispersion to effect size

Decreasing the ratio of the variance in the target random effect to total variance influences power to detect differences in the target variance among treatments. Therefore, we simulated four levels of “noise” (magnitude of non-target random effect variance) assuming a Normal distribution with increasing standard deviations (0.1, 0.5, 1.0, 2.0) (Fig. 6). These correspond to ratios of target variance parameter effect size to non-target variance of 25:1, 5:1, 5:2, and 5:4. For comparisons of ${\sigma }_{0k}^2$ and ${\sigma }_{1k}^2$,“noise” was simulated with increasing variation in within-individual variation $\left({\sigma }_{vk}^2\right)$, while for comparisons of ${\sigma }_{vk}^2$ noise was simulated with among-individual variation in intercept $\left({\sigma }_{0k}^2\right)$. For each variance parameter ratios of individuals to repeated measures followed the largest TSS_T sampling scheme used in Scenario 1 and an ES of a 2.5× difference in standard deviation by treatment.

Scenario 1: Determining the optimal sampling scheme

Power to detect differences between treatments for each variance component increases with total sample size (TSS_T) and effect size (ES) (Figs. 2–5). For a given TSS_T power depends on the ratio of the number of individuals to the number of repeated measures per individual; however, the optimal ratio of individuals to repeated measures varies depending on TSS_T and target variance parameter. For example, power to detect both ${\sigma }_{0k}^2$ and ${\sigma }_{1k}^2$ is non-monotonic across ratios of individuals to sampling occasions (Figs. 2 and 3), but is an increasing function of the number of Bernoulli observations within sampling occasions (Figs. 5A and 5B). Additionally, for each variance parameter the ratio of individuals to repeated measures within individuals that maximizes power is dependent on both TSS_T and ES. As TSS_T and ES increases, greater numbers of individuals relative to repeated measures within individuals leads to higher power for each variance parameter.

At a low sample size (TSS_T = 600) (Fig. 2A) power to detect ${\sigma }_{0k}^2$ ismaximized at a ratio of individuals to sampling occasions of 6:5 at smaller effect sizes (2×, 2.5×) and 10:3 at a large effect size (3×). Under a larger sample size and a small effect size (TSS_T = 2,400, ES = 2×) (Fig. 2C) power is maximized at a ratio of approximately 2:1, while under a large sample size and large effect size (TSS_T = 2,400, ES = 2.5×, 3×) (Fig. 2C), power is maximized at ratios ranging from approximately 5:1 to 13:1.

At a low sample size and effect size (TSS_T = 300, ES = 2×), power to detect ${\sigma }_{1k}^2$ is maximized at a ratio of 12:5 (Fig. 3A), while larger sample sizes and effect sizes (e.g., TSS_T = 600, ES = 2.5×, 3.0×; TSS_T = 1,200, ES = 2.5×, 3×) favor ratios heavily weighted towards having more individuals versus more repeated measures (ratios ranging from approximately 5:1 to 10:1; Figs. 3B and 3C). Power to detect ${\sigma }_{1k}^2$ is higher overall and less sensitive to deviations from the optimum ratio than power to detect ${\sigma }_{0k}^2$ (Fig. 3).

Power to detect ${\sigma }_{vk}^2$ follows a strikingly different pattern than ${\sigma }_{0k}^2$ and ${\sigma }_{1k}^2$. Power to detect ${\sigma }_{vk}^2$ is non-monotonic across ratios of individuals to the number of Bernoulli observations within sampling occasions (Fig. 4), and is an increasing function of the number of sampling occasions (Fig. 5C). At low sample sizes (e.g., TSS_T = 600) power to detect ${\sigma }_{vk}^2$ ismaximized by devoting nearly all of the available resources to repeated measures within individuals (ratios of approximately 1:30 to 3:40; Fig. 4A); however, at a large sample size and effect size (e.g., TSS_T = 2,400, ES = 3.0) power is maximized at a ratio of individuals to Bernoulli observations of approximately 5:6 (Fig. 4C).

Scenario 2: Power under increasing non-target random effect variance

Power to detect differences in variance components is strongly affected by the proportion of total variance that can be attributed to the target variance component (Fig. 6). Increasing variance in non-target random effects decreases power to detect differences in the target variance parameter by treatment. However, the ratio of target to non-target variance does not alter the optimal ratio of individuals to repeated measures for the target variance comparison (Fig. 6). Figure 6A demonstrates that power to detect ${\sigma }_{0k}^2$ decreases substantially as the magnitude of within-individual variation increases. Detecting differences in ${\sigma }_{1k}^2$ depends only on total random effect variation at extreme ratios of individuals to sampling occasions (e.g., 80:3) (Fig. 6B). Finally, detection of ${\sigma }_{vk}^2$ is largely independent of the magnitude of among-individual variation at large ratios of ES to non-target variance, as indicated by overlapping curves in Fig. 6C. However, when among-individual variation in intercept is very large (Fig. 6C: Red curve), power to detect ${\sigma }_{vk}^2$ decreases because individual mean responses approach 0 or 1, reducing the amount of detectable within-individual variation.

Further considerations