The effect of missing levels of nesting in multilevel analysis

Three-level linear mixed effect model

In the context of multilevel data analysis, a LMM is considered as an analysis model that accounts for correlation due to the hierarchical structure of the data. To establish background information on the analysis model for hierarchical structure data, especially for a three-level data structure, we develop a general notation and introduce a brief overview of a three-level model.

Let i be the index of the level-3 unit with size L (i=1,···,L), j be the index of the level-2 unit with size M_i (j=1,…,M_i), and k be the index of the level-1 unit with size N_ij (k=1,…,N_ij). We consider a three-level model, where a level-1 unit is nested within a level-2 unit that is also nested within a level-3 unit. Let y_ijk be the response variable for level-1 unit k in the level-2 unit j and level-3 unit i, x_ijk be the associated covariate variable observed at level-1, x_ij be the associated covariate variable at level-2, and x_i be the associated covariate at level-3. Then a general three-level linear mixed effect model is:

where eijk~N(0,σe2) are residuals at level-1 that are independent and identically distributed (iid) and capture within-cluster residual variations, v_ij=(v_ij,v′_ij) are iid random effects at level-2 and captures between-cluster intercept differences at level-2, u_i=(u_i,u′_i) are iid random effects at level-3 and captures between-cluster intercept differences at level-3.

In this study, we assume that only the regression intercept varies across clusters and hence restrict our attention to a three-level random effect model with a random intercept for each level. It is the simplest type of linear mixed model for multilevel study, although we note that an extended model including random slope can be considered. We revisit this in the discussion. Considering a random effect model with random intercepts only, α₁₁=v′_ij=γ₁₁=u′_i=0. Therefore, the model in Eq. (1) can be rewritten as

where β₀(=γ₀₀) is the intercept, and β₁(=γ₀₁),β₂(=α₀₁), and β₃(=γ₁₀) are the slope coefficients for the fixed effects at level-3, level-2, and level-1, respectively. Considering random intercepts only in the model, we assume that level-3 and level-2 residuals, u_i and v_ij, respectively, are multivariate normal with zero means and variances of σu2 and σv2, respectively, and the within-cluster residuals are normally distributed with mean zero and constant variance σe2:

The covariance structure of the response vector y=(y₁₁₁,y₁₁₂,…,y_L,M_L,N_LM) is a block-diagonal covariance matrix given by

with an ith block-diagonal matrix of Gi=σe2I(ni)+σv2J(ni), where I_(ni is an n_i×n_i identity matrix, J_(ni) is an n_i×n_i matrix of 1s, and and n_i=M_i×N_iMi is a size of a block-diagonal matrix associates with ith level-3 cluster. The off-diagonal elements are zero.

Analysis of incomplete multilevel data using single imputation

To handle the missing level-2 units in a three-level hierarchical data, we adopt a single imputation method introduced in Sanders’ study [19] that replaces the missing clusters by its nesting lower-level unit’s measurements and considers the imputed level-2 unit as a cluster containing a single observation. Specifically, it constructs imaginary level-2 clusters for level-1 measurements whose level-2 clusters were missing entirely. As a result, the single imputation method adds level-2 clusters as many as the number of level-1 units directly under the missing level-2 clusters. By applying the method, we can have a complete three-level dataset without any missing clusters at level-2 units, allowing us to apply a three-level model. The single imputation method allows us to preserve a hierarchical data structure in multilevel data analysis.

The motivating case study CATS has a three-level structure with the components of school, individual, and its repeated measures per individual, and some measurements at the individual level (level-2; academic numeracy score at baseline) were partially not observed and recorded as missing in the dataset. There were 22.4% missing academic numeracy scores at baseline at level-2, and they were filled in by single imputation that impute a corresponding academic numeracy score measured at level-1 to handle the components with missingness. In other words, we treated each level-1 unit within a missing level-2 unit is nested in a level-2 unit that only contains the level-1 unit itself [19].

Childhood to Adolescence Transition Study

Motivated by the case study of CATS, we focus on an analysis of the simulated data mimicking the CATS data [21]. The CATS is a longitudinal cohort study collected through multiple waves, and it collects mental health to investigate the effect of early depressive symptoms on academic outcomes in children from puberty through adolescence. More details about the data collection and study protocol can be found in Mundy et al. [20].

The simulated data mimicking the CATS have a three-level hierarchical structure. Individuals are nested within schools, and repeated measures within individuals were collected. The data consist of demographic, educational, and social outcomes as well as mental health outcomes for the CATS study: depressive symptoms, National Assessment Programme – Literacy and Numeracy results (NAPLAN) academic numeracy score, socio-economic status, age, and sex are collected for the study [21].

The NAPLAN numeracy score is an outcome of interest and is observed at level 2 at wave 1 (potential baseline) and at level 1 in the following waves with missing cases. As we are interested in a case where missingness happens at level-2 only, we discard cases where the NAPLAN numeracy score has missingness at level-1, and the size of data that we used is N = 2,592. The data are unbalanced as the size of level-2 differs, and the number of level-1 units within a level-2 unit varies across individual. The dataset has 163 schools; 54 of 163 schools have a single individual with a single measurement per individual, and the other 109 schools have multiple individuals with a different number of repeated measures. There were 1,142 individuals in the data, and 256 of 1,142 (22.4%) NAPLAN numeracy scores at the individual level are missing.

Simulation data

We conducted a simulation study to evaluate the performance of a three-level LMM with single imputation. In the simulation, we considered complete three-level data sets and incomplete data with various missing rates of level-2 clusters (Table 1) and compared the performance of the three-level LMM with single imputation with those of a single-level model and 2-level LMM.

Table 1

Simulation schemes to generate complete 3-level data and incomplete data with missing in level-2 units

We generated complete data sets without missing subjected to the model in Eq. (2), where we assumed the fixed coefficients of β₀=0.5, β₁=β₂=β₃=0.3, the random effects at level-2 and level-3 of ui~N(0,σu2) and vij~N(0,σv2) with σ_u=1.4 and σ_v=1, respectively, and eijk~N(0,σe2) with σ_e=1.4 for i=1,…,L, j=0,1,…,M and k=1,…,N_ij. We considered various cluster sizes at each level (Table 1). The number of level-3 clusters varied as L=30,50, and 100. For each level-3 cluster, we assumed M=10,30 clusters. The number of level-1 clusters within each level-2 cluster was randomly determined as N_ij~Unif[15,40]. For all cases, the covariates were generated as x_i~N(0,1), x_ij~exp(1)+1, and x_ijk~N(1,1.2²). Intraclass correlations were set to 0.398 and 0.602 for level-3 and level-2, respectively. For each case, we generated 500 replicates.

To evaluate the performance of the three-level LMM in the presence of missingness at the level-2 units, we generated missingness from each of the complete 3-level data sets (Table 1). We set the complete data to be missing in M level-2 clusters and its related variables such as x_ij. The missing data rate in level-2 clusters was set to 10, 25, 50, and 75% of M level-2 clusters.

Childhood to Adolescence Transition Study

We use three analysis models to determine the effect of ignorance of hidden levels in the hierarchically structured CATS data: (1) a single-level model (ordinary least square regression) that completely ignores the structure of the hierarchy and treats the data collected at level-1 only (M1), (2) a two-level LMM with a random effect for school (level-3) only and ignoring an individual (level-2) random effect (M2), and (3) a three-level LMM with single imputation considering all levels of hierarchy in the CATS study units (M3). The LMMs we fit for M2 and M3 only have random intercepts, and random slopes are not considered in the data analysis. Tables 2 and 3 report estimates of the fixed and random effects, respectively, across the three models.

Table 2

Estimates of coefficients for CATS data inferred by the three models: single-level model (M1), two-level LMM (M2), and three-level LMM (M3)

Table 3

Estimated variance components of random effects and residuals for CATS data inferred by the three models: single-level model (M1), two-level LMM (M2), and three-level LMM (M3)

Table 2 presents estimated coefficients of fixed effects in each model, and Table 3 presents the estimated variance components of random effects and residuals. Note that M2 has one random effect at school and M3 has two random effects at school and individual.

Based on a significance level of 5%, depressive symptoms appeared to be significantly meaningful in supporting the effect of mental health on academic outcomes in children from puberty through adolescence across the three models (p < 0.001). Age and NAPLAN numeracy scores observed at baseline were meaningful covariates for explaining the relationship with academic outcome in addition to the mental health outcome. That is, when children are young, have fewer depressive symptoms, and higher NAPLAN numeracy scores at baseline, academic outcomes are likely to be higher. Noticeably, the estimated coefficients of the important covariates do not differ much across the three analysis models.

Table 3 presents the estimated variance components of random effects and residuals in the multilevel data across the three analysis models with different hierarchy levels. The estimated residual variance quantifies the variation within repeated measures and decreases when random effects are included in the analysis model (M2 and M3). The estimated variance components of the random effects explain variations across schools and individuals that were not explained by M1, and the sum of the estimated variance components in M3 is larger than in M2. This implies that M3 captures more unexplained variation in the hierarchical structure data, especially the variation among individuals, which cannot be ignorable. Hence, the analysis without accounting for the variance among the unobserved level of nesting units in data analysis resulted in a larger estimation of the residual variance as shown in Table 3. It might lead to an inaccurate estimation of the total variation in the multilevel data.

Simulation data

Considering 30 scenarios (Table 1), we applied a single-level model (M1) and 2-level LMM (M2) and 3-level LMM with single imputation (M3). M1 and M2 estimate four fixed coefficients, and the standard deviation σe2 of the error term but does not estimate some of the variances of the random effects, while M3 estimates all parameters described in Eq. (2). Since M1 ignores the hierarchical structure of the data, it does not contain the random effect terms of u_i and v_ij. M2 considers level-1 and level-3 clusters by assuming that level-1 units are directly connected to level-3 clusters. Therefore, M2 does not include the random effect term of v_ij. In summary, M3 only estimates σv2, and M2 and M3 models estimates σu2. The estimates of all fixed parameters were compared to the true values of the parameters used to simulate the data. In particular, we calculated the mean squared error (MSE) across 500 replications, and the coverage probability was estimated by the proportion of replications with 95% confidence intervals containing the true value.

We expected that M3 performed well in the presence of unobserved level-2 clusters, whereas M2 would poorly estimate the coefficient β₂ of the level-2 associated covariate with missing clusters since M2 ignored the level-2 clusters and partially incorporated the date structure of level-1 and level-3 clusters. Since the M1 did not consider the hierarchical data structure, we expected that it poorly estimated β₁ and β₂.

In general, the M3 performed well and better estimated the fixed coefficients in Eq. (2) than the M1 and M2 models. As expected, in particular, we observed that the M1 and M2 models poorly estimated the coefficient β₂ of the level-2 associated covariate with missing clusters. Figs. 1 and 2 show the MSEs and the coverage probabilities for β₂, respectively, and compared those values by the three models. The MSEs by the M3 ranged from 0.00039 to 0.00557, while those by the M1 and M2 were 0.0011–0.0639 and 0.0004–0.037, respectively (Fig. 1). The coverage probabilities for β₂ by the M3 were approximately 0.814–0.964 (Fig. 2). The M1 and M2 models did poorly estimate β₂ and their coverage probabilities low as 0.27–0.406 and 0.392–0.498, respectively (Fig. 2). The MSEs for β₂ seemed to increase with the missing rate (Fig. 1); however, the coverage probabilities did not change as the missing rate increased.

Fig. 1

The MSEs for β₂, the coefficient of the level-2 associated covariate, when estimated by the single-level model (···), 2-level LMM (---), or 3-level LMM (─). For each plot of L (=30, 50, 100) level-3 clusters and M (=10, 30) level-2 clusters, each point represents the MSE across 500 simulated data sets over a range of missing rates in level-2 clusters. MSE, mean squared error; LMM, linear mixed effect model.

Fig. 2

The coverage probability for β₂, the coefficient of the 2-level associated covariate, when estimated by the single-level model (···), 2-level LMM (---), or 3-level LMM (─). For each plot of L (=30, 50, 100) level-3 clusters and M (=10, 30) level-2 clusters, each point represents the coverage probability across 500 simulated data sets over a range of missing rates in level-2 clusters. LMM, linear mixed effect model.

We observed MSEs for β₀, β_1, and β₃ similar to Fig. 1. The differences in MSEs among the three models were small, but M3 performed the best overall. Moreover, MSEs tended to increase with missing rates. Similar to Fig. 2, the coverage probabilities for β₀ and β₁ were around 0.95 from the M3, while we observed small coverage probabilities from M1 (e.g., 0.06–0.302 for β₁) and comparable coverage probabilities from M2 (e.g., 0.926–0.96 for β₁). The three models yielded similar coverage probabilities for β₃.

The three models contain different variance components. M3 can estimate the three standard deviations, σ_u, σ_v and σ_e. However, M2 does not estimate the variance σv2 among level-2 clusters, and M1 can estimate the variance σe2 among individuals but not σ_u and σ_v. We expected that the M1 and M2 models could overestimate variance components by failing to explain the 3-level data structure.

Figs. 3 and 4 show that the MSEs and coverage probabilities for σ_v from M3 worsened as the missing rate increased for each case of L and M values. At the same missing rate, the larger value of M, the better the coverage. However, more level-3 clusters did not improve MSE and coverage probability. This is because each level-3 cluster contains unobserved level-2 clusters, and hence the number of unobserved level-2 clusters increases as L increases.

Fig. 3

The MSEs for σ_v from the 3-level LMM (M3) when the level-3 cluster size changes as L = 30 (···), L = 50 (---), or L = 100 (─). For each case of M = 10 or 30 level-2 cluster sizes, 500 datasets were simulated over a range of missing rates in level-2 clusters. MSE, mean squared error; LMM, linear mixed effect model.

Fig. 4

The coverage probability for σ_v from the 3-level LMM (M3) when the level-3 cluster size changes as L = 30 (···), L = 50 (---), or L = 100 (─). For each case of M = 10 or 30 level-2 cluster sizes, 500 datasets were simulated over a range of missing rates in level-2 clusters. LMM, linear mixed effect model.

We observed that the estimates for σ_u and σ_e from M1 were unbiased with small MSEs and their coverage probabilities were close to the nominal level of 95% (0.918–0.96 for σ_u; 0.936–0.974 for σ_e). However, the coverage probabilities of M1 and M2 models were very low (e.g., 0–0.344 for σ_u and 0–0.006 for σ_e from the 2-level model), and σ_u and σe were overestimated by M1 and M2 models (see Fig. 5 for the average estimates for σe).

Fig. 5

The average estimates for σ_e when estimated by the single-level model (···), 2-level LMM (---), or 3-level LMM (─). For each plot of L (=30, 50, 100) level-3 clusters and M (=10, 30) level-2 clusters, each point represents the mean of estimates across 500 replicates, and vertical lines indicate the standard errors. The true value σ_e is 1.4. LMM, linear mixed effect model.