1.1 Notation and Terms

In the manuscript, we focus on using R code (2023) and specifically lme4 (Bates, Machler, et al., 2015) and nlme (Pinheiro et al., 2022) syntax to illustrate how models are specified. We provide R code in the text (and supplement) to facilitate implementation as R is free to use, and all code and data are also provided in our online repository. We focus on code to make the conceptual importance of random effects clearer to readers who may be less familiar with mathematical notation and matrix algebra.

In the example below, we have a single within-subject variable (W1), so we need to account for within-subject clustering of the data. Including a random intercept for each subject (1|subject) removes between-subjects variance from the outcome (DV), so that our independent variable is explaining within-subject variance after between-subjects variance has been accounted.

# Note that R uses the ‘#' symbol to indicate programmer comments

example_model <- lmer(DV~ 
                        # Fixed-effects 
                        1+ W1 +  
                        # Random-effects 
                        (1|subject), 
                        data=DATA, REML=FALSE)

# Note that DV is the name of a continuous outcome, W1 is a factorial within-
# subject variable, and subject is a factorial subject identifier.

1.2 Defining Fixed and Random Effects

There are different definitions of fixed and random effects in the literature (B. M. Bolker et al., 2009; Gelman, 2005; Snijders & Bosker, 2012), but our definitions are consistent with the frequentist mixed-effect estimation employed in lme4 (Bates, Machler, et al., 2015; Kuznetsova et al., 2017) and many domains of research (B. M. Bolker et al., 2009; Long, 2012). Fixed effects are those effects whose levels are fully represented in the data and beyond which we do not want to generalize. Random effects are those sampled from a larger population and from which we do want to generalize. For instance, in the code below, if we measured subjects’ heart rate at rest and again following a bout of exercise, then condition (at rest or following exercise) would be a fixed effect. We are only interested in the difference between that specific type/duration of exercise and rest. In contrast, subject would be a random effect, because we have a random sample of subjects, and we want to generalize from that sample back to the larger population. Programmatically, that could be written:

example_mod <- lmer(heart_rate~ 
                        # Fixed-effects 
                        1+ condition +  
                        # Random-effects 
                        (1|subject), 
                        data=DATA, REML=FALSE)

# Based on the data below:
> DATA
  subject condition heart_rate
1      s1       ctl         60
2      s1        ex         72
3      s2       ctl         42
4      s2        ex         50

1.3 Nested and Crossed Random Effects

Readers entering the world of mixed-effect models are likely to have heard of random-effects being either nested or crossed. Nested refers to situations where levels of one factor are represented at only one level of another factor. Such nesting occurs in fully hierarchical models (Raudenbush & Bryk, 2002; Singer & Willett, 2003). For instance, children are nested within families, students might be nested within classrooms, and classrooms may be nested within schools. Note that nesting is a product of the research design/environment in which a study is conducted and may not be immediately apparent in how variables are coded. Figure 3A shows an example of unambiguously coded nested factors: students (S1-6) are nested within different classrooms (C1-3) with no student belonging to multiple classrooms. As seen in the cross tabulation, each student has only one observation in only one of the classrooms. Figure 3B shows an example of an ambiguously coded nested factor. The same two levels of student (S1-2) show up at each level of each level of classroom (C1-3). However, student S1 in class C1 is not the same student as S1 in C2 or S1 in C3. Although it is best practice to code unambiguously, to make the nesting clear, we can concatenate the classroom and student labels (classroom:student) to get a unique identifier (e.g., C1:S2). The cross tabulation then shows the correct nesting: one observation from each student in only one class. If we did not concatenate these labels, we would end up with the wrong cross tabulation with 3-observations each for S1 and S2, one in each classroom.

Crossed refers to situations when different levels of one factor are represented at different levels of another factor. Two variables are said to be fully crossed if all levels of one factor are represented at all levels of the other factor. For instance, let’s say that two subjects (S1-2) had to provide responses to video stimuli (V1-3) in our two different conditions: at rest and 1-hr post exercise. Thus, in a balanced design each subject has two observations for every video, one for each condition. Our focus is on the fixed effect of condition, but our randomly selected subjects and video stimuli need to be accounted for as random effects. If there are no missing data, we would get the fully crossed design shown in Figure 2C; each subject has three observations, one for each video for which a response was measured. If there are some missing data, these factors might become partially crossed as in Figure 2D (S1 is missing an observation for V3, S2 is missing an observation for V2).

Using lme4 in R, researchers have several different options for inputting nested and crossed random effects as illustrated in Code Snippet 3.

# For students nested within classrooms we could use
(1|classroom) + (1|classroom:student)
 
# Or equivalent short cut using a slash "/":
(1|classroom/student)
# this syntax will work for ambiguous and unambiguous coding, see Figure 2

# If you have **unambiguously** coded nested factors, you can also use:
(1|classroom) + (1|student)


#For fully or partially crossed random effects, you would use:
(1|subject) + (1|video)

2.1 Random Slopes and Intercepts

To illustrate some of the modeling concerns with longitudinal data, let us consider the hypothetical study depicted in Figure 4A. This plot shows data from three different groups of participants who had different grades of spinal cord injury (high tetraplegia, cervical vertebral levels 1-4; low tetraplegia, cervical vertebral levels 5-8; or paraplegia, vertebral levels thoracic 1- sacral 5) and reported their functioning in the activities of daily living at approximately monthly intervals for 1.5 years as they underwent physical and occupational therapy. Functioning is an umbrella term for how body structures/functions (e.g., level of spinal cord injury) interact with personal and environmental factors (e.g., accessibility), to shape activity limitations and participation restrictions (Imrie, 2004; Üstun et al., 2003). Functioning in this case is the opposite of disability. For ease of interpretation, we constructed a fictitious measure of functioning on which scores range from 0 (very high disability) to 100 (very high function) and has interval/ratio properties. (Health science researchers might be more familiar with this construct as independence, but we use functioning to avoid confusion with the statistical meaning of independence.)

As shown in Figure 4A, these data will likely violate the assumptions of the ordinary least squares regression model in which errors are assumed to be independent and identically distributed. This violation can be seen most clearly when we consider that these residuals belong not only to each time point but are also unique to each participant. Consider the lowest scoring participant in the C1-4 group in the left panel of Figure 4A. In a very simple model predicting the mean value for everyone (i.e., intercept only model), this participant is going to have negative residuals (i.e., the difference between the model’s predictions and the observed value). In contrast, a high scoring participant will have mostly positive residuals. The exact value of the residuals will depend on the model, but in general we have two sources of dependence (scores from the same person tend to be more similar and scores from neighboring time points tend to be more similar) that will lead to correlations in the residuals (negatives with negatives, positives with positives), violating the assumption of independence. Violating this assumption is not trivial and can lead to considerable bias in the results of hypothesis tests (Dunlop, 1994; Sainani, 2010).

To account for this statistical dependence, the first step is to recognize that time points (time) are nested within subjects (subject) in this design. To do this, we can add a random intercept of subject to our model, as shown mathematically,

raneff_int <- lmer(functioning~ 
                        # Fixed effects 
                        1+   
                        # Random effects 
                        (1| subject), 
                        data=DATA, REML=FALSE)

raneff_lin_slope <- lmer(functioning~ 
                        # Fixed effects 
                        1+time+   
                        # Random effects 
                        (1+time| subject), 
                        data=DATA, REML=FALSE)

raneff_quad_slope <- lmer(functioning~ 
                        # Fixed effects 
                        1+time+I(time^2)+   
                        # Random effects 
                        (1+time+I(time^2)| subject), 
                        data=DATA, REML=FALSE)

2.2 Random-Effects in Non-Linear Models

The quadratic random slopes model above is what we refer to as a curvilinear model. That is, the model ultimately predicts a nonlinear trajectory for each person, but the model itself is linear in its components: \(\beta_0\) is simply added to \(\beta_1 \cdot time_i\) and added to \(\beta_2 \cdot time^2_i\) as in Eq. 4A. Forcing this nonlinearity via polynomials has strengths and weaknesses. Strengths include the ease of interpreting polynomials (e.g., the average reader can more easily conceive of time^2 compared to the log of time or time exponentiated) and they often correspond to hypotheses that researchers might have (e.g., the typical researcher might be interested in any evidence of curvature, rather than a specific power law or exponential law). However, a major weakness is that few biological or social systems follow a quadratic function (i.e., being u- or n-shaped and symmetric about a single inflection point). For instance, in Figure 4D, it is likely that participants were approaching an asymptote that could be better captured by a truly nonlinear model such as an exponential or logistic model.

A detailed discussion of nonlinear models lies outside of our focus on random effects, but see Chapter 6 of Singer & Willett (2003), and Pinheiro & Bates (2006). To summarize however, (1) regardless of the form the nonlinear model takes, there is the possibility to include random effects for each parameter, and that (2) adding random effects can often lead to issues of over-fitting, so the random effects included in the final model should be justified through a combination of theoretical rationale and empirical model fit.

Here, we present an illustrative example of a negative exponential model. First, let us consider a three-parameter exponential function:

set.seed(100)
neg_exp_rand_mod <- nlme(functioning ~ b_0i + 
                           (b_1i)*(exp(b_2i * time)),
                         data = DATA,
                         fixed = b_0i + b_1i + b_2i ~ 1,
                         random = list(b_0i ~ 1, b_1i ~ 1, b_2i ~1),
                         groups = ~ subject,
                         start = c(80, -70, -1),
                         na.action = na.omit)

3.0.1 Code Snippet 8

# Head of the data file showing the first 12 observations (subject S1 and S2)

   subject altitude condition heart_rate
1       S1      low      rest         47
2       S1      low       imm         63
3       S1      low     delay         48
4       S1     high      rest         54
5       S1     high       imm         64
6       S1     high     delay         62
7       S2      low      rest         62
8       S2      low       imm         74
9       S2      low     delay         60
10      S2     high      rest         65
11      S2     high       imm         80
12      S2     high     delay         70
# … 

# Option A: Random Intercept of Subject Only
heart_rate ~ 1 + condition*altitude + 
             (1|subject)


# Option B: Random Intercepts for All Factors
heart_rate ~ 1 + condition*altitude + 
             (1|subject)+
             (1|condition)+
             (1|altitude)


# Option C: Random Intercepts of Subject:Repeated Measures
heart_rate ~ 1 + condition*altitude + 
             (1|subject) + 
             (1|subject:condition) + 
             (1|subject:altitude)


# Option D: Random slopes of condition and altitude within each subject
heart_rate ~ 1 + condition*altitude + 
             (1+condition+altitude|subject) 


# Option E: Random slopes including the interaction term
heart_rate ~ 1 + condition*altitude + 
             (condition*altitude|subject) 

When we analyze these data using a repeated measures ANOVA, we obtain main effects of Altitude, F(1,9) = 58.36, Condition, F(2,18)=142.15, and the Altitude × Condition interaction, F(2,18)=8.63. We can use these values as a benchmark to evaluate the performance of our mixed-effect models. Using these values as our gold standard assumes that the factorial ANOVA is the model that researchers want given their design. We acknowledge that this may not always be the case, but we think it is generally best to choose random effects based on the study design and given that most researchers would analyze data from this study with a two-way repeated-measures ANOVA, it is reasonable to use that benchmark.

Before we can evaluate these different random effects structures, it is important to briefly explain two different methods of estimation: Maximum Likelihood (ML) and Restricted Maximum Likelihood (REML). Without getting into the details of how these estimators actually work (Bates, Machler, et al., 2015; Pinheiro & Bates, 2006; Singer & Willett, 2003) a key difference is that ML will base estimates on all of the available information (i.e., fixed and random effects) but does this as an iterative procedure. As a conceptual approximation, for each iteration of ML the fixed components are estimated first, and then random components estimated second given those fixed effects. Mathematically, this happens simultaneously, but conceptually we can think about the fixed effects being needed first, as the sample mean needs to be calculated to then calculate the sample variance. As with the calculation of the sample variance, this means that variances in ML tend to be under-estimated. However, this bias in ML versus REML should only substantially affect the estimation of the variance when the number of observations, \(n\), is small relative to the parameters being estimated, \(p\). Readers are likely familiar with this issue when calculating the variance of a sample, where we divide the sum of squared errors by \(n-1\) (analogous to REML) rather than simply n (analogous to ML). Thus, the de-biasing effect of dividing by \(n-p\) relative to just dividing by n depends on the values of both \(n\) and \(p\). When the number of observations is much larger than the number of parameters (\(n>>p\)), then the choice between REML and ML will usually have a trivial effect on the variance components. Importantly, however, REML does not correct this bias simply by changing the denominator of the variance components, but by using a fundamentally different algorithm for estimation.

In REML, the outcome is effectively transformed into residuals, variance components are then iteratively estimated, and then the fixed effects are estimated using generalized least squares (for an accessible discussion see McNeish (2017)). Thus, unlike ML where we simultaneously estimate the fixed and random effects in every iteration, with REML the fixed effects can essentially be treated as 0 at first because we are estimating residuals, and partitioning that variance based on the best fitting values of the random effects. In the final step of REML, we thus already have a variance/covariance matrix for the random effects and can estimate the best fitting fixed effects given that matrix.

These computational differences mean that REML and ML have complementary strengths and weakness. For instance, if you wanted to compare models with different fixed effects to ascertain the best fitting model, then you need to use ML because it is agnostic to the degrees of freedom lost to the different fixed effects in the different models. However, ML is also prone to under-estimate the variance components, especially in small samples. Thus, if there is no model comparison (e.g., you are a fitting a single mixed model with random effects determined by the design) or if you are comparing models with identical fixed effects (e.g., trying to determine only the preferred random effect structure), then REML will lead to less biased estimation of the variance components than ML.

3.1 Option A: Random Intercept of Subject Only – A Poor Choice.

In our experience, many researchers select Option A, including only a random intercept of subject, (1|subject). Including a random intercept for each subject is necessary but not sufficient to account for all the statistical dependencies in this within-subject design. This can be seen in Table 1, where we contrast the denominator degrees of freedom and F-values from a mixed-factorial ANOVA (Lawrence, 2016) and the five different mixed-effect models. Option A leads to greatly inflated denominator degrees of freedom for each of the effects; the denominator degrees of freedom are all 45, much larger than the corresponding 9 or 18 degrees of freedom in the repeated measures ANOVA. This spuriously increases statistical power (note the larger F-values for each test), in a loose sense making it seem as though we have more independent pieces of data than we do. Inspecting the degrees of freedom helps to quantify the conceptual problem with this random-effects structure: there a several strata in the data (Venables & Ripley, 1997) at which we could look at the effect of altitude and condition. With N=60 observations and only a random intercept of subject, we have essentially divided the data into two strata: a between subjects’ stratum (n=10, df=9 after the grand mean is estimated) and within subjects’ stratum (df=50 after between subject variance is removed). However, not all observations within a subject are independent of each other, as each subject has multiple observations at each altitude and condition.

3.2 Option B: Random Intercepts for All Factors – A Very Bad Choice.

To account for the fact that observations are clustered within subjects, altitudes, and conditions, we have often seen authors use a random-effects structure like Option B, treating altitude and condition as random effects. At first, these three random intercepts might seem reasonable, but there is a major conceptual flaw with this model: Our observations are not crossed by altitude and condition overall, they are only crossed by altitude and condition within each subject. The correct way to handle this is by nesting altitude and condition within each subject and accounting for those crossed factors at the within-subject level, which is precisely what we do in Option C. Before discussing that structure, however, let us dwell on Option B as an example of what not to do. As we will see, it is bad practice to enter a variable as both a fixed and random effect, because it is statistically non-sensical.

When estimated with REML, the main-effects for altitude and condition all but disappear with F(1, 45)=0.50 and F(2,45)=1.52, respectively. These F-values are greatly reduced because we have pulled out the variance due to altitude and condition in our random intercepts prior to the estimation of the fixed effects, leaving a trivial amount of variation to be explained. Recall that for REML, variance components are estimated from a residualized outcome, and fixed effects are then estimated using generalized least squares given the matrix of random effects. The interaction term remains similar to what we saw with Option A, F(2,45)=6.29, because that variance is not accounted for by the random intercepts.

When estimated with ML, the results of Option B are identical to Option A because now the fixed effects are iteratively estimated with the redundant random effects, which leaves no variance to be explained by the random intercepts of altitude and condition. The variance for both random intercepts is now estimated to be 0 (see footnotes for Table 1), and we are effectively left with a random intercept of subject.

These wildly varying estimates for altitude and condition help to illustrate how our different estimators work (i.e., REML vs. ML) and hopefully illustrate how incorrect Option B is. It is non-sensical to treat altitude as both a fixed effect and random effect, it needs to be one or the other. Note, how this is also different from longitudinal models where we estimated a fixed effect of time, but also had a random slope for time within each subject (i.e., (1+time|subject)). It is perfectly acceptable to estimate a fixed effect at the group-level, while accounting for individual differences in an effect through the random effects, discussed in Options C-E below.

3.3 Option C: Random Intercepts for Subject:Repeated Measures – A Good Choice.

Random effects Option C is a very good choice to appropriately account for the statistical dependencies in these data. In that model, the degrees of freedom and the F-values match the repeated measures ANOVA exactly using REML. These outputs match because not only are we accounting for the fact that data-points come from the same subject, (1|subject), but we also account for the fact that each subject can respond differently to each altitude, (1|subject:altitude), and each condition, (1|subject:condition). As shown in a truncated design matrix in Figure 6B, having a unique intercept for each subject:factor accounts for the fact that altitude and condition are crossed factors, within each subject. We have thus set up appropriate strata in the data at which we can evaluate our fixed effects (Venables & Ripley, 1997).

Note also that we cannot include an intercept of (1|subject:condition:altitude) because that would amount to a random deviate for every observation, leaving no residual error (i.e., such a three-way random intercept is confounded with the error term when we only have one observation person per condition, and lme4 would return an error).

As shown in Table 1, if we use ML instead of REML the results no longer match the ANOVA precisely. The degrees of freedom are slightly increased, as would expected given the differences in estimation, and this in turn leads to slightly increased F-values. Thus, the choice of estimation method will have an effect when it comes to statistical inference, but if we fit these models in a large data set (where \(n>>p\)), the practical difference between REML and ML would be small.

It is important to note that this random-intercepts method gets more complicated with more within-subject factors. For instance, if we were to add a third within-subject factor to a fully factorial design (call it W3), then not only would this add an additional fixed effect, but it adds considerable statistical dependency that needs to be accounted for in the random effects. In that situation, observations come not only from each subject in each condition, but from each subject in each combination of the conditions, so we would need to add intercepts to account for the two-way interactions within each subject (e.g.,(1|subject:condition:altitude) + (1|subject:altitude:W3) + (1|subject:condition:W3)). It is also important to note that our heuristic for the degrees of freedom may start to breakdown in more complex designs and although the degrees of freedom will be similar to the ANOVA, they will no longer match exactly (e.g., with sparse random-effects when there may only be two-levels of a factor (Venables & Ripley, 1997)). To save on space, we present a detailed example of a three-way within-subject design in the online supplement.

3.4 Option D: Random Slopes of Condition and Altitude within in Subject – Also a Good Choice.

Another reasonable alternative would be to include random slopes for each factor within each subject as in Option D. Conceptually, this model accounts for the within-subjects nature of the design similar to Option C, but it approaches the problem in a very different way and with more complexity that may not always be desirable. First, recall that when fitting random slopes, we do not only estimate their variance, but also their covariance. This allows us to account for potential correlations between the effect of altitude and condition across subjects, which is good, but increases the number of parameters being estimated, which is potentially problematic.

To simplify the model, one could also constrain the random slopes to be independent (e.g., using || in place of | in lme4 syntax, and see Frossard & Renaud (2019)). However, one would still run into a second issue which is that random slopes are actually a series of contrasts for categorical factors. That is, for altitude we only have two levels and therefore we need only one random slope to capture that difference. For condition we have three levels and therefore two random slopes are needed, which would also add more covariances if the model is not constrained. As the number of levels of the categorical factor/s increases, the model becomes substantially more complicated.

Third, and finally, when using this random-slopes approach, the user needs to cognizant of what kind of contrasts are being used (as with centering of continuous variables). By default, R will use treatment coding with the first level, alpha-numerically, treated as the reference group, but this may not always be desirable and will affect the estimates and interpretation of the model, particularly when interactions are present.

3.5 Option E: Random slopes including the interaction term – Only if you have replicates!

Estimating the variance around our effects is a very useful step towards understanding individual differences (e.g., how much does the effect of altitude differ between people?), so the random slopes structure in Option D is appealing. Ideally, we would like to take this a step further and get an estimate of how the interaction between altitude and condition differs between people as well, so we can add the random slopes that capture that interaction in Option E. However, we cannot actually estimate this model given the current data. With only one observation for each subject in each condition, if we try to fit a unique three-way interaction for every subject, we will have perfect prediction and no residual error (as with a random intercept of (1|subject:condition:altitude)). This can be seen in the footnotes under Table 1, we only have 60 observations and if we want random slopes for each factor and their interaction, then we end up with 60 random effects. As the error message says, Without replicates, this model is over-specified and cannot be fit to these data.

Replicates, as the name implies, are repeated observations at the same level of the data. For instance, assume that rather than one average heart rate for each subject in each condition, we ran our subjects through our exercise protocol on two different occasions with a sufficient washout period in between exposures to reduce order effects. We would now have two observations for each subject in each combination of altitude and condition, yielding 120 observations. Having only two observations per person may not be enough to get very stable estimates of the variance for these within-subject effects, but it would give us the observations necessary to estimate a model using Option E. Thus, having replicates is often very desirable but comes with increased complexity and often implicit assumptions about which an author might not be aware. For instance, let’s be generous and say that I have 10 observations for every combination of altitude and exercise. I could estimate a model using Option E, but that would assume that all 10 replicates are independent of each other. Depending on how those observations were obtained, independence might not be a tenable assumption. For instance, there are likely temporal order effects with neighboring observations being more similar. As a result, our model might need to become more complex, for instance estimating a covariance structure for the residuals rather than assuming that all errors are independent (Pinheiro et al., 2022; Pinheiro & Bates, 2006).

4.1 Variation in Multiple Nested Factors (e.g., Subjects within Study Sites)

Returning to our hypothetical spinal cord injury data, in Figure 7A we can see multiple sources of statistical dependency. Time points are nested within individual subjects, but now we also have subjects nested within different research sites (A-D). We can account for this additional level of nesting in our R code by using a random intercept of study site. We also include a fixed effect of spinal cord injury grade (C1-4, C5-8, T1-S5) and its interactions with time:

raneff_quad_site<-lmer(functioning~ 
                           # Fixed-effects 
                           1+time*grade+I(time^2)*grade+ 
                           # Random-effects 
                           (1+time|subject)+(1|site), data=DAT1, REML=FALSE)

4.2 Variation in Multiple Crossed Factors (e.g., Variation in Subjects and Stimuli)

In our linguistic example, we must account for the fact we have multiple observations not only for each participant, but also for each stimulus. Note that unlike the longitudinal data, there is no reliable pattern to the data within a participant, because the order of the stimuli was randomized for each person. As such, we can account for this statistical dependence with random intercepts of subject and stimulus, and a fixed effect of the modality of the stimuli (modality).

rand01 <- lmer(log(RT)~
                 # Fixed Effects 
                 1+modality+ 
                 # Random Effects
                 (1|subject)+(1|stimulus),
               data=DAT2, REML=TRUE)

DAT3 <- DAT2 %>% group_by(subject, modality) %>%
  summarize(RT = mean(RT))

agg_mod <- lmer(log(RT)~
                 # Fixed Effects 
                 1+modality+ 
                 # Random Effects
                 (1|subject),
               data=DAT3, REML=TRUE)

rand_slopes <- lmer(log(RT)~
                 # Fixed Effects 
                 1+modality+ 
                 # Random Effects
                 (1+modality|subject)+(1|stimulus),
               data=DAT2, REML=TRUE)

6.1 Data Accessibility

Data and code are available via a version controlled online repository: https://doi.org/10.5281/zenodo.10173823

These supplementary materials can also be viewed at: https://keithlohse.github.io/mixed_effects_models/

6.2 Author Contributions

• Contributed to conception and design: KRL, AJK, MJS
• Contributed to acquisition of data: NA
• Contributed to analysis and interpretation of data: KRL, AJK, MJS
• Drafted and/or revised the article: KRL, AJK, MJS
• Approved the submitted version for publication: KRL, AJK, MJS

6.3 Conflict of Interest

KRL is a member of the educational leadership team for Reproducible Rehabilitation (NIH NICHD/NCMRR R25HD105583), teaching content around statistical analysis and mixed-effect models specifically.

6.4 Funding

KRL receives salary support as a methodologist from: R01MH123723-01A1, 2 R37 HD068290-10, R01 AR081881, and NICHD R01 HD108240. The authors received no specific financial support for the research, authorship, and/or publication of this article.

6.5 Acknowledgments

The authors would like to thank many students and trainees who have asked important questions about mixed-effect regression over the years. We would also like to thank two anonymous reviewers whose insightful comments greatly improved the clarity and technical accuracy of the manuscript.

6.6 Preprint

The pre-publication version of this manuscript can be found on arxiv, https://arxiv.org/abs/2209.14349