Correct specification of design matrices in linear mixed effects models: tests with graphical representation

Linear mixed effects models (LMMs) are a popular and powerful tool for analysing grouped or repeated observations for numeric outcomes. LMMs consist of a fixed and a random component, which are specified in the model through their respective design matrices. Verifying the correct specification of the two design matrices is important since mis-specifying them can affect the validity and efficiency of the analysis. We show how to use empirical stochastic processes constructed from appropriately ordered and standardized residuals from the model to test whether the design matrices of the fitted LMM are correctly specified. We define two different processes: one can be used to test whether both design matrices are correctly specified, and the other can be used only to test whether the fixed effects design matrix is correctly specified. The proposed empirical stochastic processes are smoothed versions of cumulative sum processes, which have a nice graphical representation in which model mis-specification can easily be observed. The amount of smoothing can be adjusted, which facilitates visual inspection and can potentially increase the power of the tests. We propose a computationally efficient procedure for estimating p-values in which refitting of the LMM is not necessary. Its validity is shown by using theoretical results and a large Monte Carlo simulation study. The proposed methodology could be used with LMMs with multilevel or crossed random effects.


Introduction
We consider a (single-level) linear mixed effects model (LMM) (Laird and Ware 1982) with n groups (clusters) each having n i entries, which can be described by the following equation: (1.1) The random vector y i is a vector of dependent variables with elements y i j , j = 1, . . . , n i , which are assumed to be independent across groups but correlated within a group. The matrix X i is a given n i × m matrix of fixed effects, and β is an m × 1 vector of corresponding fixed effects coefficients. Z i is a given n i ×k matrix of random effects, and b i is a k × 1 random vector of random coefficients. The random vectors b i and i are assumed to be independent, with mean zero and covariance matrices D and σ 2 I, respectively, where b i and i are independent of X i and Z i . We will assume that the matrices X i , Z i , as well as the numbers of their rows n i , are generated randomly. More precisely, we will assume that they are determined by the independent and identically distributed (i.i.d.) random elements The random vector ν i ∈ {0, 1} n max ×1 determines the n i as n i = ν i1 + · · · + ν in max and determines which n i out of the n max rows from the random matrices X # i ∈ R n max ×m , Z # i ∈ R n max ×k , are included in the matrices X i and Z i . Here, n max ∈ N is a constant. Defining ν i allows us to apply our theory to data with groups of varying sizes and at the same time use a more accessible theory based on i.i.d. observations.
The approach proposed in this paper relies on testing two null hypotheses. The first null hypothesis is that the conditional mean of y i is equal to and thus, under H F 0 , the matrices of fixed effects X 1 , . . . , X n are correctly specified. The second null hypothesis states that the conditional mean and the conditional variance of the random vector y i are correctly specified, that is Therefore, under the null hypothesis H O 0 , the matrices of fixed and random effects X 1 , . . . , X n , Z 1 , . . . , Z n are correctly specified.
Checking whether the assumed LMM is correctly specified is important since model mis-specification affects the validity and efficiency of regression analysis. The most commonly used techniques for assessing the goodness-of-fit of LMMs are graphical tools such as residual plots (Pinheiro and Bates 2000;Wu 2009). These procedures are highly subjective and often completely uninformative. Loy et al. (2017) derived an approach based on the concept of visual p-values (Majumder et al. 2013) to make such plots less subjective. Having to rely on human experts observing the plots is, however, impractical.
While there are numerous formal tests for checking the distributional assumptions of model (1.1) (Jiang 2001;Ritz 2004;Claeskens and Hart 2009;Efendi et al. 2017), only a few tests are available for checking its choice of design matrices. Tang et al. (2014) derived a test statistic for the validity of a fixed effects design matrix. The test involves partitioning the fixed effects design matrix. The performance of the test depends on the choice of the partition and can be poor if the partition is not selected appropriately. Lee and Braun (2012) used a permutation approach for inferences regarding the inclusion or exclusion of the random effects in LMMs. Pan and Lin (2005) and Sánchez et al. (2009) proposed evaluating the choice of the design matrix for generalized LMMs (GLMMs) by considering the cumulative sum (cusum) of the ordered residuals. For LMMs, this approach has no power against alternatives in which the fixed effects design matrix is correctly specified but the random effects design matrix is not. That is, it cannot detect a mis-specification of the random effects design matrix. However, it has some appealing features, and we based our methodology on it. For example, it provides, in addition to a formal hypothesis test, an informative visual presentation that gives hints about how to improve the fit of the model (Lin et al. 2002). González-Manteiga et al. (2016) proposed a test based on the distance between two empirical error distribution functions, extending the ideas of Van Keilegom et al. (2008) for cross-sectional independent data. This approach can be used to evaluate the choice of the fixed effects design matrix (also semiparametric and generalized linear models can be considered), and it has been shown to be more powerful than the approach of Pan and Lin (2005). However, there is currently a shortage of implementations of it. Likewise, methods utilizing the link between random effects and penalized regressions have been applied to test whether the fixed effects are linear, quadratic, etc. Wood 2013Wood , 2012Huet and Kuhn 2014). To the best of our knowledge, there is no test available for checking the correct specification of both design matrices.
Note that we are not interested in the distributional assumptions of model (1.1) but only in the correct specification of the design matrices for the fixed and random effects. Extending the approach presented in Pan and Lin (2005), we construct two empirical stochastic processes that inspect H O 0 , the entire model, and H F 0 , only the fixed part of the model. We formally prove that when the fixed effects design matrix is correctly specified, the process for checking only the fixed effects design matrix will be robust against the mis-specification of the random effects design matrix. Intuitively, this should hold since the estimator of β is consistent given only the correct specification of the fixed effects design matrix (Zeger et al. 1988). Since the parameter sets of the fixed and random effects parts of the LMM are chosen separately, i.e., independently, a mis-specification of the random effects part of the model is implied when the null hypothesis that both design matrices from the model are correctly specified (H O 0 ) is rejected and the null hypothesis that the fixed effects part of the model is correctly specified (H F 0 ) is not rejected. The above observation allows one to construct a procedure with which to address the mis-specifications of both LMM design matrices one at a time. The procedure combines two tests based on empirical stochastic processes, where the process for testing the entire model is novel. The asymptotic theory for these two tests is derived based on the strong fundamental stochastic process theory presented by van der Vaart and Wellner (1996). Finally, these empirical stochastic processes may be nicely visu-alized, and the deviations from the null hypothesis for each process can easily be judged from the figures. We also introduce a smoothing parameter, which facilitates visual inspection and can improve the power of the tests and show how the amount of smoothing can be determined from the data.
The challenging part for all applications of the (empirical) stochastic processes in this context is obtaining their null distribution. Given the complexity of the problem introduced by the dependence among the residuals, the asymptotic distribution for even the most trivial test statistic is analytically intractable. In application to linear models (LMs), the null distribution is obtained by using bootstrapping (Stute et al. 1998a) or simulations (Su and Wei 1991;Lin et al. 2002). The simulation approach has also been used for marginal models (MMs) (Lin et al. 2002) and single-level GLMMs (Pan and Lin 2005). We propose a computationally efficient bootstrap approach for correctly approximating the null distribution of the proposed processes, in which refitting of the LMM is not necessary. This approach creates new residuals, which will be called modified residuals, constructs stochastic processes based on them and evaluates the test statistic. We will show that several distributions could potentially be used when creating the modified residuals, as long as some mild assumptions about their moments are satisfied, and we will investigate the finite sample performance for some obvious candidate distributions.
In Sect. 2, we outline the approach, and in Sect. 3, we introduce some additional notation. Next, in Sect. 4, we present the definitions of the proposed empirical stochastic processes and define modified residuals. A short subsection regarding the algorithmic data-driven choice of the smoothing parameter concludes the section. The assumptions under which we establish the empirical stochastic processes' weak convergence are provided in Sect. 5. We show the asymptotic validity of the proposed methodology under the null (in Sect. 6) and alternative hypotheses (in Sect. 7). In Sect. 8, we showcase the finite sample performance for a selection of Monte Carlo simulation results. An application to a real data example is given in Sect. 9. The paper concludes with a summary of the most significant findings and possibilities for future research. For brevity, we only consider single-level LMMs in detail; a possible extension to LMMs with more complex random effects structures is discussed in the supplementary material. Proofs and additional simulation results are shown in the supplementary material.

Overview of the proposed approach
Our approach is based on the repeated application of Algorithm 1. In each application of the algorithm, we use residuals and fitted values from a fitted LMM with assumed design matrices X i and Z i , but we define the residuals and the fitted values in two different ways depending on which hypothesis, Construct a process based on the modified residuals and fitted values from 1. 7: Calculate the test statistic based on the process from the previous step. 8: end for 9: Calculate the p-value based on Steps 3 and 7. 10: return the p-value and visualize the processes from Steps 2 and 6 The theoretical setting under which we prove the validity of our approach assumes that there is some stochastic mechanism that is generating the groups of data y i . This mechanism is also assumed to be generating the matrices X i , Z i and is assumed to be generating groups of data that may have different sizes. The parameters β, D and σ 2 are assumed to be fixed. All our asymptotic results assume that the number of groups n grows to infinity.
The advantage of this theoretical setting is that under some assumptions given later, the proposed approach can be shown to be valid on the type of data to which linear mixed models are very often applied. Another advantage of this theoretical approach is that we do not have to deal with independent nonidentically distributed data such as in the theoretical approach considered in Hagemann (2017), which introduces various technical problems regarding measurability, making the theory less accessible.

Notation and general definitions
We use span(M), M ∈ R m 1 ×m 2 to denote the column span or the image of matrix M. Given a vector subspace V ⊂ R n , we denote its orthogonal complement by V ⊥ . As usual, a hat over a symbol will be used to denote an estimate of this quantity. For example, the estimatorV i of the matrix V i , defined as We define the following quantities: The random vectorsb i are known as BLUPs. The motivation for their use can be found in Demidenko (2005). Furthermore, we define Here, Z i Z i + denotes the Moore-Penrose inverse of the matrix Z i Z i . The Therefore, it also holds that var e F i = I. The conditional covariance matrix of e O i is an orthogonal projection matrix. Hence, the absolute values of the conditional covariances of the elements of e O i are bounded by 1. These results slightly simplify the theoretical results.
When dealing with empirical stochastic processes, we mostly use the same notation as van der Vaart and Wellner (1996). We denote the probability space as (X , A, P). The random elements X i ∈ X are defined as Recall that the random vector ν i and random matrices X # i , Z # i determine the random element D i . The random vector e # i ∈ R n max ×1 is defined so that its rows, determined by ν i , determine the random vector e i . The random element X i therefore completely determines every known quantity in (1.1).
When discussing the independence of random vectors and of matrices and vectors that are implicitly or explicitly determined by the random elements {X i } i∈N , we will always refer to the independence conditional on ν i . We do not require any additional conditions regarding the independence of the rows of matrices X i and Z i within a group.
Let δ X i denote the Dirac measure on the space (X , A). Let X 1 , . . . , X n be n random elements in (X , A). Define the empirical measure P n as P n = (δ X 1 + · · · + δ X n )/n. Define the empirical stochastic process based on n observations as: As in van der Vaart and Wellner (1996), we will index (empirical) stochastic processes by functions. For example, Let m : X → R m be some function that is linear in e i , i.e., where i is a matrix that is not dependent on e i and has elements with finite variance. For example, when using MLE or REML, i is equal to (Demidenko 2005 The symbol denotes convergence in distribution, and * denotes convergence in distribution conditional on the data, that is, conditional on X 1 , X 2 , . . . To establish * , we will use the theory from Chapter 2.9 of van der Vaart and Wellner (1996) denoting multipliers as ξ i .

Definition of the empirical stochastic processes
Let λ ∈ (0, ∞) be a constant. Define a function σ λ : R → [0, 1], which is a continuous approximation of the indicator function that we will use in constructing the empirical stochastic processes. More precisely, let p : [0, 1] → [0, 1] be a twice-differentiable monotonous function with p(0) = 1 and p(1) = 0. Then, we define the function σ λ as: Define the function p so that as λ increases, the function σ λ becomes a better approximation for the indicator function. Define σ ∞ : R → [0, 1] to be equal to where 1 A denotes the indicator function of the set A. When the function σ λ is applied to each component of a vector x ∈ R d , we write it as σ λ (x). For a random vector (4.1) and they are obtained by replacing the unknown quantities in (4.1) by their respective consistent estimators. We will show later that when mis-specifying the fixed effects design matrix, the mean function for the processŴ F n is nonzero for some t, while it is zero for all t when the fixed effects design matrix is correctly specified, regardless of the (in)correct specification of the random effects design matrix. Moreover, we will show that the mean function for the processŴ O n is nonzero for some t if the fixed and/or random effects design matrix is mis-specified. These deviations from a zero-mean stochastic process can then also easily be observed visually by plotting the process against t. In addition to the processesŴ O n andŴ F n , we define the procesŝ where X i j,l is the lth column and jth row, j = 1, . . . , n i , of the design matrix for the fixed effects X i , i = 1, . . . , n, and the sum extends only over some subset of the columns of the fixed effects design matrix. UsingŴ F S n (t) allows one to test a possible lack of fit due to the specified subset of the fixed effects covariates. By inspectinĝ W F S n (t), it is also possible to detect the omission of an important interaction effect of two or more variables. For example, if the p-value based on W F S n (t), defined as a subset of two variables, is significant at some level α and the two p-values based on W F S n (t), defined as a subset of only one variable, are both insignificant at level α, this implies that an important interaction effect between the two variables was omitted. Since the processŴ F S n is very similar to the processŴ F n , it will not be thoroughly examined separately.
We calculate p-values using either the Kolmogorov-Smirnov (KS) or Cramér-von Mises (CvM) type statistic, i.e., Instead of obtaining the theoretical probabilities for a value of a test statistic, we rely on the repeated evaluation of these test statistics on the processes defined using modified residuals, which are defined as: where the random elements X * i are equal to i are the components of the random element X i . We will show that the choice of the distribution of the multipliers, ξ i , does not matter asymptotically, as long as they are independent, zero-mean, unit-variance variables with some additional restrictions on the existence of higher-order moments, which will be presented in the next section. In small samples, the choice of the distribution can make a difference, as will be illustrated in the simulation study.
The processesŴ O * n (t) andŴ F * n (t), empirical versions of W O * n and W F * n , which are defined as: are then obtained by replacing the unknown quantities with their respective consistent estimators. Note that the underlying data do not change, but the data from the modified residuals do. Therefore, we will later prove the convergence of these empirical stochastic processes conditionally on the data, which is in line with the bootstrap approach presented in Chapter 3.6 of van der Vaart and Wellner (1996). For a very large λ, the proposed processes behave similarly to commonly used processes constructed as cumulative sum(s) of the model's residuals for LMs (Christensen and Lin 2015; Lin et al. 2002;Stute et al. 1998a;Diebolt and Zuber 1999;Su and Wei 1991;Fan and Huang 2001;Stute et al. 1998b), MMs (Lin et al. 2002) and single-level GLMMs (Pan and Lin 2005). A similar process to the W F n (t) process of (4.1) was considered by Pan and Lin (2005), but they used different residuals when constructing the process and only considered the situation in which λ = ∞. Using the function σ λ facilitates the visual inspection of the processes by smoothing them. At the same time, the choice of λ < ∞ can potentially increase the power.

Data-driven choice of the smoothing parameter
One potential data-driven choice of the smoothing parameter λ is to define a grid of different λs, selecting the one that gives the largest value of the test statistic. To obtain the correct size and to avoid the need for multiplicity correction, this is also done for the processes using the modified residuals. Since λ affects the scale of the process, and hence the scale of the test statistic, the test statistics need to be standardized to account for the difference in scaling. For this purpose, we suggest using the estimated standard error of the test statistic, which can be obtained from the processes using the modified residuals. The entire procedure is presented in Algorithm 2.
Algorithm 2 Data-driven choice of the smoothing parameter λ 1: Define a grid of λs: {λ 1 , . . . , λ K }. 2: for k = 1 : K do 3: Using λ k , calculate the test statistic for the original process, T k , and the processes using the modified residuals, T * k . 4: Estimate the standard error of T k , denoted by s k , from T * k .

5:
Define the standardized test statistic for the original process as T s k = T k /s k , and similarly define the processes using the modified residuals: T * s k = T * k /s k . 6: end for 7: Calculate the p-value using max It is reasonable to assume that when the choice of the design matrices of a certain model is assessed, this model has already been fitted to the data. Our assumptions, therefore, while seeming strict at first glance, mainly require that the model has been fitted using parameter estimators that satisfy certain conditions. The assumptions required to establish the convergence of the proposed empirical stochastic processes are listed below. A brief discussion of the assumptions is given at the end of this section.
to be at most n max , where n max is some constant that is at least k + 1. That is, X i ∈ R n i ×m , Z i ∈ R n i ×k and e i ∈ R n i ×1 . Matrices X i ∈ R n i ×m and Z i ∈ R n i ×k have full rank almost surely, and their elements have bounded second moments. (A2) We assume that we have consistent estimatorsβ,D andσ of the unknown parameters β ∈ R m , D ∈ R k×k and σ > 0, respectively. The matrix D is positive definite. The estimatorβ is asymptotically independent of the estimatorŝ D andσ . Furthermore, We assume that P m = 0 and that every element of P mm is smaller than ∞.  (2005), the model is identifiable. Hence, we could ease the requirements in assumption (A1), but we would then need to assume additional conditions regarding the estimation process.
Note that the assumption (A2) holds when the errors are distributed normally and we use MLE, REML or the Fisher-Scoring algorithm to estimate the parameters of the LMM. More information on this can be found in Chapters 2 and 3 in Demidenko (2005). Without the normality assumption, (A2) holds when using the estimator proposed by Peng and Lu (2012) under some additional regularity conditions (see Peng and Lu (2012) for more details).
Assumption (A3) helps us establish that the stochastic processes are equicontinuous. The most apparent use of this assumption is in the proof of Theorem 2, when λ = ∞.
The last condition in (A4) implies that ξ i 2,1 = ∞ 0 √ P(|ξ i | > t)dt < ∞, by exercise 2.9.1 in van der Vaart and Wellner (1996). This condition is a required assumption in the multiplier central limit theorem that we use in the proof of Theorem 3. This condition is satisfied, for example, in the case in which ξ i are Rademacher or standard normal random variables.
The assumption (A5) will be needed to show that under H 0 , W O n is a zero-mean process. Without this assumption, this result would still hold if we assumed that Z i = (1, . . . , 1) holds for i = 1, . . . , n. The reason why this result could otherwise not be established is presented in more detail in the supplementary information, where a counterexample is constructed.
Note that the dimensions m and k are assumed to be fixed. With some additional assumptions, such as those found in Jiang (1996), or by showing that certain estimators of β, D and σ satisfy the conditions in He and Shao (2000), we could extend our setting to settings with increasing m and k. In this case, we would use the same ideas, with some additional references to Sect. 2.11.3 of van der Vaart and Wellner (1996).

Asymptotic behaviour under H 0
In this section, we present a part of the theoretical justification for our method when the model is correctly specified. More precisely, we prove that the stochastic processes on which our approach is based, under the assumptions from the previous section, converge in distribution as the number of groups of data n grows to infinity. We also prove that the stochastic processes based on the modified residuals converge weakly conditionally on the data and that their limits are equal to the limits of the processes based on the original residuals. Here we only summarize the main ideas of the proofs of the presented theorems. Detailed proofs can be found in the supplementary material.
First, we justify the use of the function σ λ . Let V(F) denote the Vapnik-Chervonenkis index (VC index) of a class of functions F, as defined in van der Vaart and Wellner (1996, pp.141).
Define the classes of functions F O and F F as: (6.1) Note that P f F t = 0 for every t ∈ R. This follows from the fact that the random vector e F i has zero mean and is independent of y F i , since it is defined so that it does not depend on the random matrix X i .
To show that P f O t = 0 for every t ∈ R, we additionally assume that (A5) holds. Note the equality e O i = P i e i = P i i . P i is the orthogonal projection such that Therefore, the empirical stochastic processes W O n and W F n at t can be written as: The next theorem ensures that under some assumptions, the processes W O n and W F n converge.
Theorem 1 Under the assumption (A1), the families of functions F F defined as in (6.1) are P-Donsker. Under the assumptions (A1) and (A5), the families of functions F O defined as in (6.1) are P-Donsker.

Short proof of Theorem 1
The assumption (A5) ensures that P f O t = 0 for t ∈ R. Furthermore, the classes of functions F O and F F are both VC classes, with VC indices equal to 2. The reason for this is that the functions f O t , for t ∈ R, are constructed so that the subgraph of f O s is contained in f O t for every s ≤ t. The same holds for the functions f F t , t ∈ R. We can therefore use Theorem 2.6.7 from van der Vaart and Wellner (1996), which implies that the uniform entropy condition in Theorem 2.5.2 from van der Vaart and Wellner (1996) is satisfied. We proceed by constructing two square measurable envelopes to satisfy all of the requirements in Theorem 2.5.2, which proves this theorem. We can do this because of assumption (A1). Now, we define the following two families of functions: (6.2) We define the limit empirical stochastic process G similar to G n for f ∈ F O ∪ F F as n G n f .
Note the equality [m(X 1 ) + · · · + m(X n )]/ √ n = G n m. The next result establishes the convergence ofŴ O n andŴ F n .
Short proof of Theorem 2 We can use the assumption (A2) to show that the processeŝ W O n andŴ F n can be written as sums of either W O n or W F n , and in the case ofŴ O n , From this, we can use the continuity of the function t → σ λ (·, t) for λ < ∞ or the Assumption (A3) when λ = ∞, in combination with Lemma 2.12 from van der Vaart (1998), to show that the second summand in the above process is equicontinuous. Then, we apply the same reasoning as in the proof of Theorem 19.23 from van der Vaart (1998) to show that the processŴ O n converges weakly. The proof of the weak convergence ofŴ F n is very similar.
Theorem 3 shows that the processes based on the modified residuals, conditionally on the data, converge weakly and that their limits are equal to the limits of the processes based on the original residuals.

Theorem 3 Assume that the Assumptions (A1), (A2) and (A4) hold. ForŴ O *
n , additionally assume that (A5) holds. Then, if either λ < ∞ or λ = ∞ and Assumption (A3) holds, the processesŴ O * n andŴ F * n converge in distribution conditionally on the data toŴ Short proof of Theorem 3 Since the i.i.d. random variables ξ i have zero mean and variance 1, the distributions of finite collections of marginals ofŴ O n andŴ O * n are the same, and the distributions of finite collections of marginals ofŴ F n andŴ F * n are the same.
Furthermore, assumption (A4) allows us to use the central multiplier theorem 2.9.6 from van der Vaart and Wellner (1996). This, together with the same kind of argument as that given at the end of the previous theorem, completes the proof.

Asymptotic behaviour under some alternative hypotheses
In this section, we prove that our method rejects the null hypothesis under one of the three alternative hypotheses stated below.
We assume that we have access only to the random vectors y i from the original data, but instead of the true matrices u X i and u Z i , which determine the outcome vector y i , i = 1, . . . , n, by (where u b i and u i are independent sequences of i.i.d. random vectors, which are independent of u X i , u Z i and D i and have mean zero and variances var( u b i ) = u D and var(u i ) = u σ 2 I, where u D is some positive definite matrix and u σ > 0), we have access to some possibly different matrices X i and Z i . The numbers of columns in matrices X i and Z i may be different than the numbers of columns in u X i and u Z i . The i.i.d. sequences of random elements to which we have access will be denoted by {X i } i∈N and {D i } i∈N , and the i.i.d. sequence of random elements that generated the random vectors y i will be denoted by Note that the processes will be based on {X i } i∈N .
We will use the following two statements: (S1) We say that a sequence of random matrices X i , i ∈ N is correctly specified when the equation does not hold for at most finitely many i. (S2) We say that a sequence of random matrices Z i , i ∈ N is specified correctly when, conditionally on u D i and D i , the equation does not hold for at most a finite number of i.
Note that the negation of the two statements (S1) and (S2) is that the equations (7.1) and (7.2), respectively, do not hold for an infinite number of i but not necessarily all i ∈ N.
The three alternative hypotheses are then given as follows: (H 1 ) The random matrices X i , i ∈ N are not specified correctly, and the matrices Z i , i ∈ N are specified correctly.
(H 2 ) The random matrices X i , i ∈ N are specified correctly, and the matrices Z i , i ∈ N are not specified correctly.
(H 3 ) The random matrices X i , i ∈ N are not specified correctly, and the matrices Z i , i ∈ N are not specified correctly.
In all three alternative hypotheses, we also assume the following: (B1) The random elements X i are i.i.d., and their elements have bounded second moments. (B2) The estimatorsβ,D andσ are obtained in the same way as the estimators based on the original data in (A2). We assume that they converge to the limits β, D and σ . We can therefore write We also assume that the estimatorD and its limit D are positive definite.
Note that if the assumption (B2) is not satisfied-more specifically, if one of the estimatorsβ,D orσ does not converge-then the process of fitting the LMM may fail.
A probable cause of the violation of (S1) is that the column span of X (n) , does not contain a part of the column span of u X (n) , for every n. This happens, for example, when the matrix X (n) does not contain a certain column of u X (n) .
More concretely, assume that we use either the MLE or REML method of estimating the parameters and that the matrix u X (n) , which has full rank almost surely, is generated as in Sect. 8, i.e., that its entries are i.i.d. with the possible presence of a column of ones. Furthermore, assume that there is at least one column of u X (n) that is not present in the matrix X (n) . If the intercept is not present in the matrices X (n) and u X (n) , then (S1) does not hold. However, if the intercept is present in both matrices X (n) and u X (n) , then (S1) still holds, because the intercept causes the residuals to be centred around 0. But when the missing columns are a function of the columns that are present in both matrices X (n) and u X (n) , (S1) does not hold.
Observe that since the matrix P i , i ∈ N is defined so that it maps any vector from span(Z i ) to zero, in the case when (7.2) does not hold, span( P i u Z i ) = {0}. Therefore, the random vectors e O i and y O i are correlated conditionally on u D i and D i , and for every t in some open interval I . Additionally, because the random elements X i and u D i are assumed to be i.i.d., this interval I is the same for every i. An example of this arises when the matrices Z i do not contain a column of u Z i that is not a linear combination of other columns of u Z i , for i ∈ N. In the next theorem and proofs, we use the same quantities-random vectors, random matrices and stochastic processes-as in the previous section. In this case, however, they will not be based on random elements { u D i } i∈N but will be based on the available data {X i } i∈N .

If hypothesis (H 2 ) holds,
and the processesŴ F n andŴ F * n converge weakly to the same limit.

Short proof of Theorem 4
Note that (ξ 1 + · · · + ξ n ) P f O t / √ n and (ξ 1 + · · · + ξ n ) P f F t / √ n are zero mean. We can then use the same arguments as in Theorems 2 and 3 to prove that the processesŴ O * n andŴ F * n converge weakly conditionally on the data {X i } i=N . On the other hand, √ n P f O t and √ n P f F t do not converge under any of (H 1 ), (H 2 ) and (H 3 ), except in the case of (H 2 ) and √ n P f F t .
Theorem 4 (cases 1 and 3) implies that when the matrices X i are not correctly specified, the mean functions of the processesŴ O n (t) andŴ F n (t) will be nonzero for some t, while the mean functions of the respective bootstrapped processes,Ŵ O * n (t) andŴ F * n (t), will be zero for all t, regardless of the (in)correct specification of the matrices Z i . In contrast, when X i are correctly specified but Z i are not, this only holds for the processŴ O n (t), whereas the mean functions ofŴ O * n (t),Ŵ F n (t) andŴ F * n (t) will be zero for all t.

Simulation results
In the simulation study below, we report the results forŴ F n (t), the F-process, and for W O n (t), the O-process. The analysis was performed in R (R version 3.4.3) using the R package gofLMM (available on GitHub, https://github.com/rokblagus/gofLMM). The LMMs were fitted by using the function lme from the nlme package (Pinheiro and Bates 2000). The variance parameters were estimated by REML. The p-values were estimated by using M = 500 random realizations of null approximations simulating ξ i , i = 1, . . . , n, from a standard normal distribution (Pan and Lin 2005) (normal), a Rademacher distribution (i.e., the sign-flipping approach (Winkler et al. 2014), SF) that attaches a mass of 0.5 to the points −1 and 1 and a Mammen two-point distribution (Stute et al. 1998a) (Mammen) that attaches masses ( √ 5 + 1)/2 √ 5 and ( √ 5 − 1)/2 √ 5 to the points −( √ 5 − 1)/2 and ( √ 5 + 1)/2; these distributions have been shown to work well in regression settings (Hagemann 2017). Note that all distributions satisfy E(ξ ) = 0 and E(ξ 2 ) = 1, while E(ξ 3 ) = 0 in the case of the standard normal and Rademacher distributions and E(ξ 3 ) = 1 for the Mammen distribution, thus satisfying the conditions in Assumption (A4) in our theoretical investigation. In the definition of σ λ , two different options for specifying the function p were considered, and considering different values of λ = {0.5, 1, 2, 4, 6, 8, 10, ∞} (note that for λ = ∞, the function σ λ is the indicator function irrespective of the definition of p); these values of λ were also specified in the grid when using the data-driven approach for choosing λ presented in Sect. 4.1. The differences between the two functions p are not substantial, and hence, only the results for the function defined in (8.1) are shown here (see the supplementary material for the results obtained with the other definition of p for some values of λ). For computational reasons, σ λ is only evaluated at distinct fitted values. The p-values were simulated 5000 times; the simulation margin of errors is thus ±0.003, ±0.006 and ±0.008 for α = 0.01, 0.05 and 0.1, respectively. For the Fprocess, we compare the performance of our proposed tests with the approach proposed by Pan and Lin (2005) (equivalent to using the standard normal distribution when simulating ξ i and setting λ = ∞) and the restricted likelihood ratio test (RLRT) of Greven et al. (2008). The RLRT was performed using the function exactRLRT from the R package RLRsim (Scheipl et al. 2008).
We omit the KS test statistics from the results, since they are less powerful than those of CvM. Throughout, the results were similar in terms of size for different choices of λ, and hence, only the results for the data-driven choice of λ are shown when considering the performance of the tests under H 0 . The results for various values of α when considering the power of the tests were as expected, with larger power obtained with larger values of α; hence, only the results for α = 0.05 are shown in illustrating the performance of the tests when mis-specifying the model.
The outcome was simulated from a linear mixed model with random intercepts and slopes where j = 1, . . . , n i and i = 1, . . . , n. The number of observations for all n groups was the same. All quantities on the right side of (8.2) were simulated independently from each other: the covariates X i j,1 ∼ U(0, 1) and X i j,2 ∼ U(0, 1), random effects b i,0 ∼ N (0, 0.25) and b i,1 ∼ N (0, σ 2 b,1 ), error term i j ∼ N (0, 0.5), β 3 , and σ 2 b,1 were set according to the scenarios in the next subsections.

Example I: Size under normal errors and normal random effects
The outcome was simulated from (8.2), specifying the variance of the random effect b i,1 as σ 2 b,1 = 0.25 and β 3 = 0. The fitted model was the same as the simulated model. The simulations were performed for n = 50, 75 and n i = 5, 10, 20.
The empirical sizes of the tests were close to nominal levels for both processes and all distributions (Fig. 1). The exceptions were the situations with a smaller sample size, where the tests for the O-process based on the Mammen and standard normal distributions were slightly too conservative; the difference between the empirical rejection rate and the nominal level was, however, not substantial, and it diminished with increasing sample size. Importantly, this was not a consequence of our data-driven approach for Fig. 2 Power under a mis-specified random effects design matrix for different λ (columns) and processes (rows) using the CvM test statistics; α = 0.05. Each figure shows the empirical rejection rate (y axis) for a different variability of the missed random effect b i,1 (x axis) and a different number of groups (left, n = 50; right, n = 75). The colours specify different proposed distributions used in the proposed bootstrap approach. The shaded areas are the simulation margins of error under the null hypothesis. RLRT-restricted likelihood ratio test of Greven et al. (2008). For a normal distribution and λ = ∞ (σ λ is the indicator function), the results for the F-process are equivalent to using the approach of Pan and Lin (2005); opt-the λ chosen using the data-driven approach presented in Sect. 4.1 (colour figure online) choosing λ, since similar results were also obtained when considering fixed values of λ that formed the grid (data not shown). The RLRT was too conservative, rarely rejecting the null hypothesis regardless of the sample size. Very similar results were obtained when relaxing the normality assumption in the error and random effect terms, showing the robustness of our approach against non-normality (see the supplementary material).

Example II: Mis-specified random effects design matrix
The outcome was simulated from (8.2), β 3 = 0, varying the dispersion of b i,1 by σ 2 b,1 = 0.5, 1, 1.5 and considering n = 50, 75 and n i = 10. The fixed effects part of the fitted model was correctly specified, but the random effects part included only a random intercept.
As suggested by our theoretical results, the empirical sizes of the tests for the fixed effects part of the model were close to the nominal level for all values of λ, demonstrating robustness against a mis-specification of the random effects design matrix (Fig. 2). As in the previous examples, the RLRT did not perform well. The tests for the O-process rejected the null hypothesis more often than the nominal level. The rejection rates were larger with smaller λ and larger n and/or σ 2 b,1 , with the SF approach being the most powerful. (This was more evident with large values of λ.) The data-driven approach for choosing λ performed well, with a power that was only slightly smaller than that obtained with the λ that yielded the largest power amongst all the values considered in the grid. Fig. 3 Power under a mis-specified fixed effects design matrix for different λ (columns) and processes (rows) using the CvM test statistics; α = 0.05. Each figure shows the empirical rejection rate (y axis) for different fixed effects β 3 (x axis) and different numbers of groups (left, n = 50; right, n = 75). The colours specify different proposed distributions used in the proposed bootstrap approach. The shaded areas are the simulation margins of error under the null hypothesis. RLRT-restricted likelihood ratio test of Greven et al. (2008). For the normal distribution and λ = ∞ (σ λ is the indicator function), the results for the F-process are equivalent to using the approach of Pan and Lin (2005); opt-the λ chosen using the data-driven approach presented in Sect. 4.1 (colour figure online)

Example III: Mis-specified fixed effects design matrix
The outcome was simulated from (8.2) with b i,1 variance σ 2 b,1 = 0.25 and varying β 3 = 0.5, 1, 1.5, considering n = 50, 75 and n i = 10. The random effects part of the fitted model was correctly specified, but the fixed effects part included only the linear effects of the covariates.
The empirical rejection rates of all tests were larger than the nominal level, showing that the tests are powerful against this alternative (Fig. 3). The rejection rates when using the F-process were larger than those when using the O-process. The test based on the SF approach was the most powerful (this was more obvious for the O-process and large values of λ). The power increased with larger n and β 3 . The power with the O-process increased with smaller λ. In contrast, for the F-process, the power was generally smaller for smaller values of λ. The data-driven approach for choosing λ performed well for both processes, yielding power comparable to what could be obtained by using the value of λ that obtained the largest power. For the F-process, our approach was more powerful than the RLRT.

Example IV: Mis-specified fixed and random effects design matrices
The outcome was simulated from (8.2) with β 3 = 1 and a varying dispersion of b i,1 with σ 2 b,1 = 0.25, 0.5, 1, 1.5, considering n = 50, 75 and n i = 10. The fixed effects part of the fitted model included only the linear effects of the covariates, and the random effects part included only a random intercept. Fig. 4 Power under mis-specified fixed and random effects design matrices for different λ (columns) and processes (rows) using the CvM test statistics; α = 0.05. Each figure shows the empirical rejection rate (y axis) for a different variability of the missed random effect b i,1 (x axis) and a different number of groups (left, n = 50; right, n = 75). The colours specify different proposed distributions used in the proposed bootstrap approach. The shaded areas are the simulation margins of error under the null hypothesis. RLRTrestricted likelihood ratio test of Greven et al. (2008). For the normal distribution and λ = ∞ (σ λ is the indicator function), the results for the F-process are equivalent to those when using the approach of Pan and Lin (2005); opt-the λ chosen using the data-driven approach presented in Sect. 4.1 (colour figure online) The empirical rejection rates of all tests were larger than the nominal level, showing that the tests are powerful against this alternative (Fig. 4). The rejection rates when using the O-process were larger than those when using the F-process. As expected, the power of the F-process decreased with increasing σ 2 b,1 . Similar to the other examples, the SF approach was the most powerful, but the differences between the approaches were small in this particular example. The power of the O-process increased with smaller values of λ. In contrast, the power of the F-process was the smallest with λ = 0.1 and was comparable with the other values of λ. The proposed data-driven approach again performed well, exhibiting good power. For the F-process, our approach was more powerful than the RLRT.

Application
We apply the proposed methodology to cross-sectional data from the 2004 American National Election Study (ANES 2022) (see the supplementary material for an application to longitudinal data). The ANES is a series of surveys on voters' opinions before and after elections in the USA. The 2004 ANES data, with the outcome variable feeling thermometer reading for George W. Bush (a variable with values from 0 to 100, with higher values indicating a more positive feeling towards Bush) and a large set of possible predictor variables, were used as a real data example in Peng and Lu (2012).
Since variable effects tend to be mediated by social and cultural context at the state level, the natural way to handle this data was to fit a linear mixed effects model in which individuals are sampled from states. Peng and Lu (2012) used the data as an illustration for their proposed model selection procedure that uses iterative penalized regression to extract important variables to be included in the fixed and random effects design matrices. On the ANES data set, it yielded a final model with 11 variables for the fixed effects design matrix and two variables for the random effects design matrix, omitting the intercept in the random effects design matrix. We applied our methodology to their final model, which is a good starting point for the fine-tuning that can be done with our proposed methodology.
We prepared the 2004 ANES data in the same way as in Peng and Lu (2012) (the 1212 respondents were decreased to 1156 individuals from 24 states, the number of subjects in each state ranged from 19 to 136, five states were excluded, and the included variables were recoded in the same fashion). The regression coefficients of the nlme fit for the model (9.1), f eeling i j = β 0 +β 1 ·age i j +β 2 · educ i j +β 3 ·christian i j +β 4 · black i j +β 5 · other i j + β 6 · liberal i j + +β 7 · de f ence i j + β 8 · death i j + β 9 · democrat i j + β 10 · indep i j + β 11 · Iraq i j + b i,1 · gender i j + b i,2 · christian i j , (9.1) differ only slightly from the ones published in Peng and Lu (2012) due to different data set versions (see the supplementary material for more details). The only numerical variable in the model is age, and all the others are dichotomous. The codes 0 and 1 are used to represent not having and having a characteristic, respectively (see the supplementary material for the exact meaning of each variable).
In the first row of Fig. 5, we show our proposed empirical stochastic processes (black) for the final Peng-Lu fit of the model. The M = 500 generated null F-and O-processes obtained by using the SF approach are shown in grey. The fit is not good (there is a low p-value for the F-process), with a sequence of (mainly) negative ordered residuals in the middle of the plot following a sequence of (mainly) positive residuals.
The addition of six two-way interactions between dichotomous variables (defence-Iraq, liberal-Iraq, democrat-Iraq, indep-Iraq, democrat-black and christian-black) to the fixed effects design matrix yields a reasonable fit for the fixed effects (see the left plot in the second row of Fig. 5 for the F-process; see the longitudinal example presented in the supplementary material for an illustration of how subset F-processes can be used to detect important omitted interaction effects). However, the p-value for the O-process is still significant at the 5% level, implying that the random effects design matrix is mis-specified. The final improvement of fit comes with the inclusion of black, Iraq, and two-way interactions for christian-black and christian-Iraq in the random effects design matrix, which also cause the p-value of the O-process to become insignificant at the 5% level (Fig. 5, right plot in the third row).
By fine-tuning the model obtained by Peng and Lu (2012) (with the addition of interactions to the fixed and random effects design matrices), the interpretation of the final model yields additional insight. It can be observed, for example, that democrats' expected feeling thermometer towards Bush is very low, but this changes significantly if the person supports the war in Iraq. for the three different models (a Rademacher distribution was used to obtain the 500 simulated null processes, shown as grey lines). Row 1 corresponds to model (9.1), row 2 to the same model enhanced with 6 two-way interactions and row 3 to the final model for the 2014 ANES data set. The p-values above each graph correspond to the CvM test statistic with λ = 1200 chosen from the grid {5, 10, 50, 200, 500, 1000, 1200} using the proposed data-driven procedure

Discussion and conclusions
We proposed a novel procedure for testing the assumed specification of the design matrices of fitted LMMs, for which an asymptotic theory based on the strong fundamental stochastic process theory presented by van der Vaart and Wellner (1996) was derived; its validity in terms of size was demonstrated, and its power against several alternatives was showcased. The approach is based on inspecting different empirical stochastic processes that are constructed as smoothed cumulative sums of appropriately standardized and ordered model residuals: the O-process, the F-process, and the subset F-process. Investigating the O-process allows us to test the correct specification of both design matrices. In contrast, investigating the F-process (or the subset F-process) allows us to investigate the assumed fixed effects design matrix (or some subset thereof).
We showed that the O-process is expected to fluctuate around zero when the fixed and random effects design matrices are correctly specified, while it is not when either (or both) of the fixed and random effects design matrices are mis-specified. In contrast, it was shown that the F-process and subset F-process fluctuate around zero when the fixed effects design matrix is correctly specified, regardless of the (in)correct specification of the random effects design matrix. While we did not construct a process that would target only the random component, a mis-specification of the random effects design matrix can still be detected. That is, we showed that when there is no evidence that the fixed effects design matrix is mis-specified but enough evidence to deduce that the entire LMM is mis-specified, this implies a mis-specification of the random effects design matrix.
Any observed fluctuations or deviations from zero can be evaluated visually and by means of p-values by using the proposed computationally efficient approach, which does not require re-estimating the LMM. Several different multiplier distributions could be used when estimating the p-values. Asymptotically, the choice of the multiplier distribution is not important as long as Assumption (A4) holds. However, it could make a difference in finite samples. In our simulation study, we considered three options: a standard normal distribution, a Rademacher distribution (i.e., a signflipping approach) and a Mammen two point distribution. While the three distributions performed similarly in terms of size, attaining the nominal level, the approach based on the Rademacher distribution yielded better results in terms of power. Therefore, even though the differences between the three distributions in terms of power were not substantial, we would recommend using the Rademacher distribution, especially with a smaller sample size.
The indicator function is usually used in the context of goodness-of-fit testing (equivalently to using a cumulative sum of ordered residuals). We proposed replacing the indicator function with a continuous function of λ (and the model's fitted values), which for large values of λ is a close approximation of the indicator function. The constant λ can be seen as a smoothing parameter of the cusum process, facilitating the visual inspection of the plots and making it easier to identify potential improvements in the model's fit. In our simulations, smoothing improved the power of the O-process. For the F-process, the differences in power for most considered values of the smoothing parameter were small in our examples. The exceptions were very small values of the smoothing parameter (i.e., a very large amount of smoothing), where the power was reduced due to (excessive) smoothing. However, we have also identified examples, where smoothing can improve the power for the F-process (one example is shown in the supplementary material). We proposed a straightforward data-driven approach for choosing the amount of smoothing. This approach performed well in our simulations, obtaining power comparable to what could potentially be obtained by specifying a single value of λ that yielded the largest power amongst all the values forming the grid. There might still be room to further improve the power, e.g., by not relying on a pre-specified grid, which we think represents an interesting subject for future research. By using a different multiplier distribution (and smoothing) we were able to improve the power of the approach proposed by Pan and Lin (2005) for checking the correct specification of the fixed effects design matrix. González-Manteiga et al. (2016) proposed an approach which has greater power than the Pan and Lin (2005) approach and it would be interesting to compare it to our method; however, readily available code makes this difficult.
While we only considered single-level LMMs in detail, the proposed methodology could be adapted to multiple nested levels of random effects, but at a cost of notational inconvenience (see the supplementary information). In principle, the methodology presented here could be extended to GLMMs. However, further extensions to nonlinear link GLMMs could be problematic when trying to distinguish between the reasons for a (possible) lack of fit due to fixed or random effects design matrices since the parameter sets of fixed and random effects cannot be chosen separately, i.e., independently. Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.