Maximum Likelihood Estimation of Hierarchical Linear Models from Incomplete Data: Random Coefficients, Statistical Interactions, and Measurement Error

Abstract

We consider two-level models where a continuous response R and continuous covariates C are assumed missing at random. Inferences based on maximum likelihood or Bayes are routinely made by estimating their joint normal distribution from observed data $R_{\text{obs}}$ and $C_{\text{obs}}$. However, if the model for R given C includes random coefficients, interactions, or polynomial terms, their joint distribution will be nonstandard. We propose a family of unique factorizations involving selected “provisionally known random effects” u such that $h(R_{\text{obs}},C_{\text{obs}}|u)$ is normally distributed and u is a low-dimensional normal random vector; we approximate $h(R_{\text{obs}},C_{\text{obs}})=\int h(R_{\text{obs}},C_{\text{obs}}|u)g(u)\mathit{\text{du}}$ via adaptive Gauss-Hermite quadrature. For polynomial models, the approximation is exact but, in any case, can be made as accurate as required given sufficient computation time. The model incorporates random effects as explanatory variables, reducing bias due to measurement error. By construction, our factorizations solve problems of compatibility among fully conditional distributions that have arisen in Bayesian imputation based on the Gibbs Sampler. We spell out general rules for selecting u, and show that our factorizations can support fully compatible Bayesian methods of imputation using the Gibbs Sampler. Supplementary materials for this article are available online.

KEYWORDS:

• Adaptive Gauss Hermite quadrature
• Compatibility
• Missing at random
• Provisionally known random effects
• The EM algorithm

1 Introduction

Large-scale surveys and experiments within the social and health sciences frequently meet four conditions that supply the focus of this article. First, the data typically have a hierarchical structure, with respondents nested within local organizational units such as schools and hospitals or repeated measures nested within persons. Second, missing data are pervasive. Third, partially observed covariates may be measured with error. Finally, the covariates of interest may have random coefficients, statistical interactions or polynomial terms.

These characteristics have received some attention in recent methodological research. A popular approach conceives the response variables and partially observed covariates as outcomes within a multivariate, hierarchical linear model (HLM) under the assumption that the data are missing at random (MAR; Rubin 1976), an assumption often thought reasonable given a presumably rich set of covariates (Schafer and Yucel 2002; Goldstein, Carpenter, and Browne 2014). Missing values are imputed from their posterior predictive distribution of the missing values using what has been termed a “fully conditional specification” (FCS) using the Gibbs sampler. FCS requires the analyst to impute missing values for each variable subject to missingness, conditional on all other unknowns. This approach has been shown to function well when the process generating the joint distribution of unknowns is reasonably assumed multivariate normal. Under the four conditions just described, however, multivariate normality is not possible even if the separate conditional distributions are normal. A concern involves the bias that can arise from the incompatibility between the multiple normal conditional distributions that generate the imputations and the assumed joint distribution of the observed data (Erler 2016; Enders, Mistler, and Keller 2016; Enders, Du, and Keller 2020).

In this article, we propose to address the problem of compatibility (Arnold and Press 1989; Liu et al. 2014; Bartlett et al. 2015) within the framework of maximum likelihood (ML) estimation under the ignorable missing data assumption that data are MAR and that the parameter spaces for the multivariate HLM and the missing data mechanism are distinct (Rubin 1976; Little and Rubin 2002).

We will first review currently popular methods of inference for normal-theory multilevel models given incomplete data and show how to modify such approaches in which random effects often become explanatory variables, allowing us to model fallible measurement as a source of incomplete data. We’ll then describe our approach when incompletely observed covariates have random coefficients or polynomial terms, including statistical interactions. If a carefully selected subset of these random effects is conditionally assumed known, the model of interest can plausibly follow a normal-theory specification. The “provisionally known random effects (PKREs)” thus, selected can be integrated out of the likelihood using numerical methods. Capitalizing on the invariance properties of ML estimates (MLE), we show that our reparameterized model can be translated back to the original model’s parameter space. Auxiliary variables are introduced to increase the robustness of the MAR assumption. We will illustrate application of the model by studying income inequality and mathematics achievement in U.S. elementary schools. Although we base our case study on MLE, the likelihood factorization we propose can readily be implemented within a Bayesian approach with assurance of compatibility between conditional distributions and the joint distribution of observed data elements. We elucidate general rules for selecting low-dimensional “provisionally constant” random effects for a general class of models involving random coefficients or polynomial terms, including interactions.

Section 2 reviews estimation of normal-theory hierarchical models from data MAR. Section 3 explains how to estimate an analytic hierarchical model with random coefficients and cross-level interaction effects given data MAR via a PKRE. Section 4 describes estimation of the analytic model using auxiliary covariates to strengthen the MAR assumption. Section 5 extends the model with polynomial terms including within-level interactions and spells out rules for selecting provisionally constant random effects. Section 6 illustrates analysis of income inequality in math achievement. Section 7 evaluates our estimators by simulation. Finally, Section 8 discusses the limitations and extensions.

2 Inference for Multivariate Normal HLMs from Incomplete Data Using Random Effects as Predictors

We begin with a review of the normal theory case. Following Schafer and Yucel (2002), our scientific interest focuses on the regression of a response variable $R^{*}$ on covariates $C^{*}$. The elements of $R^{*}$ and $C^{*}$, which are continuously measured, are partially observed. Ours is a two-level HLM (Lindley and Smith 1972; Dempster, Rubin, and Tsutakawa 1981) in which the response variable $R^{*}$ is a characteristic of “level-1 units” (e.g., students) who are clustered within level-2 units (e.g., schools). In contrast, the covariates may be characteristics of either level-1 units or level-2 units. In longitudinal studies, the level-1 units might be repeated occasions of measurements clustered within persons at level 2. Our scientist is primarily interested in the conditional distribution $f_1(R^{*}|C^{*})$ assumed normally distributed. However, to account for the missing values in $C^{*}$, we propose a normal linear model $f_2(C^{*})$. Because of missing data, we cannot separately estimate the parameters of f₁ and f₂ without discarding the cases with missing values on any element of $R^{*}$ or $C^{*}$. It is well known that a procedure that analyzes only the cases with complete data is prone to bias and/or loss of efficiency (Little and Rubin 2002). The cost can be particularly high when a missing item of $C^{*}$ varies at level 2, in which case all of the level-1 units in the level-2 unit having missing values are discarded along with that level-2 unit itself.

To make efficient inference possible, and following Schafer and Yucel (2002), we compose each outcome vector $Y^{*}={[R^{*}{C}^{*T}]}^T$ and write a multivariate HLM $h(Y^{*})$ (1) $$Y^{*}=X^{*}\alpha +Z^{*}b+r^{*}\sim N(X^{*}\alpha ,V^{*}=Z^{*}\Omega Z^{*T}+{\Sigma }^{*})$$(1) where $b\sim N(0,\Omega )$ and $r^{*}\sim N(0,{\Sigma }^{*})$. Here, $X^{*}$ and $Z^{*}$ are composed of completely observed covariates; α is a vector of fixed regression coefficients while b and $r^{*}$ are independent random effects that vary at levels 2 and 1, respectively. We partition the complete data (CD) $Y^{*}$ into components $Y^{*}=(Y_{\text{obs}},Y_{\text{mis}})$. In particular, if $Y^{*}$ is N by 1 but we observe only $M(\le N)$ elements $Y_{\text{obs}}=Y$ of $Y^{*}$, we construct an M-by-N matrix O in which every row contains a single entry equal to unity indicating which value of $Y^{*}$ that is observed. All other entries in the same row are 0. Our model for the observed Y is therefore (2) $$Y=X\alpha +\mathit{\text{Zb}}+r\sim N(X\alpha ,V=Z\Omega Z^T+\Sigma )$$(2) where $Y=O{Y}^{*},X=O{X}^{*},Z=O{Z}^{*},r=O{r}^{*}$ and $\Sigma =O{\Sigma }^{*}{O}^T$. Assuming data MAR, we can make efficient estimates of the parameters $\theta =(\alpha ,\Omega ,{\Sigma }^{*})$ using ML or Bayes inference from the observed data according to (1). The most common method in recent literature is a Bayesian approach based on multiple imputation (MI) that we will consider in the final section of this article (see Section 8). However, Schafer and Yucel (2002) showed how one can obtain MLE using the EM algorithm and use the estimates for MI. Shin and Raudenbush (2007), showed how to recover MLE for the analytic model $f_1(R^{*}|C^{*})$ by constructing model (1) carefully and estimating the model without a need for imputations. This approach allows some components of $Y^{*}$ to vary at level 2.

2.1 Estimation via the EM Algorithm

The EM algorithm (Dempster, Laird, and Rubin 1977) requires evaluation at each iteration m + 1 of the conditional expected CD score given the observed data and parameter estimates at iteration m. To find the CD score, our model for the CD $Y^{*}$ is a multivariate HLM (1). Write model (1) at the level of cluster j and let $\phi$ be an arbitrary scalar element of $(\Omega ,{\Sigma }^{*})$. The CD score equations are well known (see, Raudenbush and Bryk 2002, chap. 14): (3) $$\begin{array}{c}{S}_{\alpha ,\text{CD}}=\sum _j{X}_j^{*T}{V}_j^{*-1}(Y_j^{*}-X_j^{*}\alpha ),\\ S_{\phi ,\text{CD}}=\frac{1}{2}\sum _j{\left(\frac{d\text{vec}(V_j^{*})}{d\phi }\right)}^T(V_j^{*-1}\otimes V_j^{*-1})\\ \hspace{1em}\text{vec}\left[(Y_j^{*}-X_j^{*}\alpha ){(Y_j^{*}-X_j^{*}\alpha )}^T-V_j^{*}\right].\end{array}$$(3)

The conditional expected score equations given the observed Y_j thus, clearly depend on the conditional mean and variance of $Y_j^{*}$ which we readily derive from the fact that $Y_j^{*(m+1)}|Y_j,{\widehat{\theta }}^{(m)}\sim N({\widehat{Y}}_j^{*(m+1)},{\widehat{V}}_j^{*(m+1)})$ given iteration-m estimates ${\widehat{\theta }}^{(m)}=({\widehat{\alpha }}^{(m)},{\widehat{\Omega }}^{(m)},{\widehat{\Sigma }}^{*(m)})$ (Shin and Raudenbush 2007) where (4) $$\begin{array}{c}{\widehat{Y}}_j^{*(m+1)}=X_j^{*}{\widehat{\alpha }}^{(m)}+{\widehat{V}}_j^{*(m)}{O}_j^T{\widehat{V}}_j^{(m)-1}(Y_j-X_j{\widehat{\alpha }}^{(m)}),\\ {\widehat{V}}_j^{*(m+1)}={\widehat{V}}_j^{*(m)}-{\widehat{V}}_j^{*(m)}{O}_j^T{\widehat{V}}_j^{(m)-1}{O}_j{\widehat{V}}_j^{*(m)}.\end{array}$$(4)

2.2 Example 1: Contextual Effects Model

We apply this framework to the “contextual effects model” (Willms 1986) for the study of income inequality in the mathematics achievement of U.S. elementary school children. This model decomposes the association between family income and educational achievement into a within-school component and a between-school component. In its simplest form, the model is typically written as (5) $$R_{\mathit{\text{ij}}}={\gamma }_{00}+{\gamma }_{10}(C_{\mathit{\text{ij}}}-{\overline{C}}_j)+{\gamma }_{01}({\overline{C}}_j-\overline{C})+u_{0j}+e_{\mathit{\text{ij}}}$$(5) where R_ij is a measure of math achievement for student i in school j, C_ij is a measure of the income of that child’s parents, ${\overline{C}}_j$ is the sample mean of C_ij within school j, and $\overline{C}$ is the overall sample mean income for $i=1,\dots ,n_j$ and $j=1,\dots ,J$. Here γ₁₀ is known as the “within-school coefficient” while γ₀₁ is the “between school” coefficient; and $u_{0j}$ and e_ij are independent, normally distributed random effects that vary at levels 2 and 1, respectively. Of interest is the “contextual component” ${\gamma }_c={\gamma }_{01}-{\gamma }_{10}$ which, if positive, suggests that attending a school with high-income peers predicts elevated achievement net of the contribution of a student’s family income.

Two problems arise in the conventional analysis (Shin and Raudenbush 2010). First, past analyses have treated parental income as completely observed when, in fact, most surveys report substantial fractions of missing data on income. Second, the sample mean ${\overline{C}}_j$ will be a noisy proxy for the actual mean income of parents in a school if the sample size per school is modest, as is the case in most U.S. national surveys. Using random effects as explanatory variables within the MAR framework addresses both issues. Reflecting that income is partially observed, we propose the CD model (6) $$R_{\mathit{\text{ij}}}^{*}={\gamma }_{00}+{\gamma }_{10}{ϵ}_{\mathit{\text{ij}}}^{*}+{\gamma }_{01}{\nu }_j+u_{0j}+e_{\mathit{\text{ij}}}^{*},C_{\mathit{\text{ij}}}^{*}=\delta +{\nu }_j+{ϵ}_{\mathit{\text{ij}}}^{*}$$(6) where ν_j and ${ϵ}_{\mathit{\text{ij}}}^{*}$ are independent, normally distributed random effects at levels 2 and 1, respectively. The joint distribution of $R_{\mathit{\text{ij}}}^{*}$ and $C_{\mathit{\text{ij}}}^{*}$ may be written (7) $$\begin{array}{c}\left[\begin{array}{c}{R}_{\mathit{\text{ij}}}^{*}\\ C_{\mathit{\text{ij}}}^{*}\end{array}\right]=\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}{\gamma }_{00}\\ \delta \end{array}\right]+\left[\begin{array}{cc}1& 0\\ 0& 1\end{array}\right]\left[\begin{array}{c}{\gamma }_{01}{\nu }_j+u_{0j}\\ {\nu }_j\end{array}\right]\\ +\left[\begin{array}{c}{\gamma }_{10}{ϵ}_{\mathit{\text{ij}}}^{*}+e_{\mathit{\text{ij}}}^{*}\\ {ϵ}_{\mathit{\text{ij}}}^{*}\end{array}\right].\end{array}$$(7)

Stacking the equations within level-2 unit j, we have the general form of the model (8) $$Y_j^{*}=X_j^{*}\alpha +Z_j^{*}{b}_j+r_j^{*}\sim N(X_j^{*}\alpha ,V_j^{*}=Z_j^{*}\Omega Z_j^{*T}+I_{n_j}\otimes {\Sigma }^{*})$$(8) for $X_j^{*}=Z_j^{*}=1_{n_j}\otimes I_2,b_j\sim N(0,\Omega )$ and $r_j^{*}\sim N(0,I_{n_j}\otimes {\Sigma }^{*})$ where we denote $1_m$ as a vector of m unities and I_m an m-by-m identity matrix for a positive integer m. The HLM score equations for θ are familiar (Raudenbush and Bryk 2002, chap. 14; Shin and Raudenbush 2010) and will therefore not be elaborated here.

2.3 Compatibility

A set of conditional distributions $f_1(R^{*}|C^{*},{\theta }_1)$ and $f_2(C^{*}|R^{*},{\theta }_2)$ is said to be compatible if there exist a joint distribution $\{h(Y^{*}|\theta ):\theta \in \Theta \}$ and surjective maps $\{t_j:\Theta \to {\Theta }_j:j=1,2\}$ such that for each j, ${\theta }_j\in {\Theta }_j$ and $\theta \in t_j^{-1}({\theta }_j)=\{\theta :t_j(\theta )={\theta }_j\}$, we have $f_1(R^{*}|C^{*},{\theta }_1)=h(R^{*}|C^{*},\theta )$ and $f_2(C^{*}|R^{*},{\theta }_2)=h(C^{*}|R^{*},\theta )$ (Liu et al. 2014). Assuming a prior $p(\theta )$, Schafer and Yucel (2002) developed the Gibbs sampler based on multivariate normal (MN) $h(Y^{*}|\theta )$ that is compatible with an analytic HLM $f_1(R^{*}|C^{*},{\theta }_1)$ when $C^{*}$ is linearly associated with $R^{*}$. The MN $h(Y^{*}|\theta )$, however, cannot be compatible with $f_1(R^{*}|C^{*},{\theta }_1)$ when $C^{*}$ has nonlinear effects (Kim, Sugar, and Belin 2015; Enders, Du, and Keller 2020). Goldstein, Carpenter, and Browne (2014) and Enders, Du, and Keller (2020) factored $h(Y^{*}|\theta )=f_1(R^{*}|C^{*},{\theta }_1)f_2(C^{*}|\theta {\prime }_2)$ and estimated a compatible HLM $f_1(R^{*}|C^{*},{\theta }_1)$ with the nonlinearities by the Gibbs sampler via a Metropolis algorithm. We extend the HLM further with the nonlinear effects of $C^{*}$ that includes random effects as latent covariates. We estimate $h(Y^{*}|\theta )$ efficiently by ML via a PKRE and translate the estimates to the ML estimates of the compatible HLM as explained in the next section.

3 Coping with Random Coefficients and Interactions

Our scientific focus is again on $f_1(R^{*}|C^{*},\nu )$ assumed normal in distribution as in (6). However, the model now includes elements of $C^{*}$ that have nonlinearities including random coefficients, polynomial terms, or interactions. In this case, even if we can reasonably assume that $f_2(C^{*})$ is normal, the joint distribution $h(Y^{*})$ cannot be normal.

Recent research on the problem of nonlinearities has focused primarily on two widely used methods of analysis of multilevel incomplete data. Perhaps the most popular method of imputation for HLMs under the MAR assumption imputed missing values by a series of sequential univariate regression models (Raghunathan et al. 2001), which is also known as MI by fully conditional specification or FCS (van Buuren et al. 2006). However, these conditionals will not be compatible with the joint distribution of interest in the presence of the nonlinearities of interest in this article, as shown by Enders, Mistler, and Keller (2016) and Enders, Du, and Keller (2020). The main alternative approach estimates the joint distribution of the outcome and covariates subject to missingness. Missing values may be imputed based on their fully conditional distributions (Liu, Taylor, and Belin 2000; Schafer and Yucel 2002; Goldstein et al. 2009) and by maximum likelihood (Shin and Raudenbush 2007, 2010; Ren and Shin 2016). These normal-theory models were not designed to handle the nonlinearities of interest here. By means of a Gibbs sampler via the Metropolis Hastings algorithm, Goldstein, Carpenter, and Browne (2014) imputed missing values of a response and covariates including interaction and polynomial terms having fixed effects in a multilevel model where covariates and response may be continuous or categorical. Similarly, Erler (2016) took the sequential fully Bayesian approach (Ibrahim, Chen, and Lipsitz 2002) that expresses the joint distribution of variables MAR, including the outcome, into a series of univariate conditional models to handle missing values of cluster-level continuous and discrete covariates having fixed effects. Erler et al. (2019) extended the approach to imputing missing values of level-1 covariates having fixed effects. Enders, Du, and Keller (2020) showed, however, that these approaches do not guarantee compatibility when the partially observed covariates have random coefficients. Lacking a formal model for the joint distribution of interest, these approaches appear to fall short of ensuring compatibility.

3.1 Factorization of the Likelihood Based on Provisionally Known Random Effects

To cope with nonlinearities induced by partially observed covariates having random coefficients, polynomials, or interaction terms, we model the joint distribution $h(Y^{*}|\nu )p(\nu )$ induced by the scientific model of interest $f_1(R^{*}|C^{*},\nu )$ and the model for the covariates, $f_2(C^{*}|\nu )p(\nu )$. The problem is similar to the problem of estimation of generalized linear mixed models (Hedeker and Gibbons 1994; Raudenbush, Yang, and Yosef 2000). Using the notation of (2), we must evaluate (9) $$h(Y)=\int h(Y|b)p(b)\mathit{\text{db}}.$$(9)

The integral just defined does not have closed form in the presence of nonlinearities and must be approximated numerically. To facilitate the approximation, we choose a PKRE, call it u such that (10) $$h(Y|u)=\int h(Y|b,u)p(b|u)\mathit{\text{db}}$$(10)

is a normal HLM discussed in Section 2. The problem of approximation is then to evaluate (11) $$h(Y)=\int h(Y|u)g(u)\mathit{\text{du}}.$$(11)

The computational challenge is to select u such that the dimension of the analytic integral (10) is maximized while the dimension of numerical approximation (11) is minimized.

3.2 Example 2: A Partially Observed Covariate Having a Random Coefficient

We return to our example of income inequality in U.S. elementary schools. A common finding in multilevel studies of educational achievement is that relationship between student socioeconomic background and achievement varies from school to school (Raudenbush and Bryk 1986). This variation could reflect variation in school organization, composition, and resources. To assess the variation in income inequality within schools, we expand the contextual model (6) to allow for a random coefficient (12) $$\begin{array}{c}{R}_{\mathit{\text{ij}}}^{*}=({\gamma }_{00}+u_{0j}^{*})+({\gamma }_{10}+u_{1j}^{*}){ϵ}_{\mathit{\text{ij}}}^{*}+e_{\mathit{\text{ij}}}^{*},\\ C_{\mathit{\text{ij}}}^{*}=\delta +{\nu }_j+{ϵ}_{\mathit{\text{ij}}}^{*}\end{array}$$(12) where $u_{0j}^{*}\sim N(0,{\tau }_{00}),u_{1j}^{*}\sim N(0,{\tau }_{11}),{\nu }_j\sim N(0,{\tau }_{\nu \nu }),\text{cov}(u_{0j}^{*},u_{1j}^{*})={\tau }_{01},\text{cov}(u_{0j},{\nu }_j)={\tau }_{0\nu }$ and $\text{cov}(u_{1j}^{*},{\nu }_j)={\tau }_{1\nu }$. This is model (6) for $u_{0j}^{*}={\gamma }_{01}{\nu }_j+u_{0j}$ if ${\tau }_{11}=0$. Let $\text{var}(u_j^{*})=\tau$ for $u_j^{*}={[u_{0j}^{*}{u}_{1j}^{*}{\nu }_j]}^T$. Level-1 random effects $e_{\mathit{\text{ij}}}^{*}\sim N(0,{\sigma }^2)$ and ${ϵ}_{\mathit{\text{ij}}}^{*}\sim N(0,{\sigma }_{\mathit{\text{cc}}})$ are assumed independent of each other and of the level-2 random effects. We denote the parameters of the CD model (12) as ${\theta }_{(12)}^{*}=({\gamma }_{00},{\gamma }_{10},\tau ,{\sigma }^2,\delta ,{\sigma }_{\mathit{\text{cc}}})$.

Clearly $R_{\mathit{\text{ij}}}^{*}$ cannot be marginally normal in distribution because of the multiplication of the two normal random effects $u_{1j}^{*}$ and ${ϵ}_{\mathit{\text{ij}}}^{*}$. Our strategy is to select one of these two random effects to be considered “provisionally known”; we choose $u_{1j}^{*}$ for this purpose because it has lower dimension varying across schools than does ${ϵ}_{\mathit{\text{ij}}}^{*}$ which varies across students within each school. Therefore, we write (13) $$\left[\begin{array}{c}{u}_{0j}^{*}\\ {\nu }_j\end{array}\right]|u_{1j}^{*}\sim N\left(\left[\begin{array}{c}{\alpha }_{0|1}\\ {\alpha }_{\nu |1}\end{array}\right]u_{1j}^{*},\Omega =\left[\begin{array}{cc}{\tau }_{00|1}& {\tau }_{0\nu |1}\\ {\tau }_{\nu 0|1}& {\tau }_{\nu \nu |1}\end{array}\right]\right)$$(13) where ${\alpha }_{k|1}={\tau }_{k1}{\tau }_{11}^{-1}$ and ${\tau }_{\mathit{\text{kk}}\prime |1}={\tau }_{\mathit{\text{kk}}\prime }-{\alpha }_{k|1}{\tau }_{11}{\alpha }_{k\prime |1}$ for $k,k\prime =0,\nu$. We can therefore write the “provisional” joint model $h(Y_j^{*}|u_{1j}^{*})$ for $Y_j^{*}={[Y_{1j}^{*T}\cdots Y_{n_{j}j}^{*T}]}^T$ as (14) $$Y_j^{*}=X_j^{*}\alpha +Z_j^{*}{b}_j+r_j^{*}\sim N(X_j^{*}\alpha ,V_j^{*}=Z_j^{*}\Omega Z_j^{*T}+{\Psi }_j^{*})$$(14)

where $X_j^{*}=1_{n_j}\otimes [I_2{I}_2{u}_{1j}^{*}],\alpha ={[{\gamma }_{00}\delta {\alpha }_{0|1}{\alpha }_{\nu |1}]}^T,Z_j^{*}=1_{n_j}\otimes I_2,b_j={[u_{0j}^{*}-{\alpha }_{0|1}{u}_{1j}^{*}{\nu }_j-{\alpha }_{\nu |1}{u}_{1j}^{*}]}^T\sim N(0,\Omega )$ and $r_j^{*}={[r_{1j}^{*T}\cdots r_{n_{j}j}^{*T}]}^T\sim N(0,{\Psi }_j^{*}=I_{n_j}\otimes {\Sigma }_j^{*})$ for (15) $$\begin{array}{c}{r}_{\mathit{\text{ij}}}^{*}=\left[\begin{array}{cc}1& {\gamma }_{10}+u_{1j}^{*}\\ 0& 1\end{array}\right]\left[\begin{array}{c}{e}_{\mathit{\text{ij}}}^{*}\\ {ϵ}_{\mathit{\text{ij}}}^{*}\end{array}\right],\\ {\Sigma }_j^{*}=\left[\begin{array}{cc}{({\gamma }_{10}+u_{1j}^{*})}^2{\sigma }_{\mathit{\text{cc}}}+{\sigma }^2& ({\gamma }_{10}+u_{1j}^{*}){\sigma }_{\mathit{\text{cc}}}\\ ({\gamma }_{10}+u_{1j}^{*}){\sigma }_{\mathit{\text{cc}}}& {\sigma }_{\mathit{\text{cc}}}\end{array}\right].\end{array}$$(15)

The CD score is therefore familiar; see (3). To complete the EM algorithm to estimate $h(Y_j^{*}|u_{1j}^{*})\varphi (u_{1j}^{*};0,{\tau }_{11})$, we will maximize the likelihood $L(\theta )={\prod }_{j=1}^{J}h(Y_j)$ for (16) $$\begin{array}{c}h(Y_j)=\int h(Y_j|u_{1j}^{*})\varphi (u_{1j}^{*};0,{\tau }_{11})d{u}_{1j}^{*},\\ h(Y_j|u_{1j}^{*})\sim N(X_j\alpha ,V_j=Z_j\Omega Z_j^T+{\Psi }_j)\end{array}$$(16) where $h(Y_j|u_{1j}^{*})$ is from (14), $X_j=O_j{X}_j^{*},Z_j=O_j{Z}_j^{*}$ and ${\Psi }_j=O_j{\Psi }_j^{*}{O}_j^T={\oplus }_{i=1}^{n_j}{\Sigma }_{\mathit{\text{ij}}}$ for $O_j={\oplus }_{i=1}^{n_j}{O}_{\mathit{\text{ij}}}$ and ${\Sigma }_{\mathit{\text{ij}}}=O_{\mathit{\text{ij}}}{\Sigma }_j^{*}{O}_{\mathit{\text{ij}}}^T$. We use adaptive Gauss-Hermite Quadrature (AGHQ) to numerically approximate integral (16) (Naylor and Smith 1982; Pinheiro and Bates 1995; Rabe-Hesketh, Skrondal, and Pickles 2002).

An additional AGHQ step is needed to complete the E step by evaluating the expectation of a CD score component $S_{\text{CD}j}$ of cluster j from (3) (17) $$E(S_{\text{CD}j}|Y_j)=\int \int S_{\text{CD}j}f(Y_j^{*}|Y_j,u_{1j}^{*})g(u_{1j}^{*}|Y_j)d{Y}_{\text{mis}j}d{u}_{1j}^{*}$$(17) where $f(Y_j^{*}|Y_j,u_{1j}^{*})$ has the familiar form of the empirical Bayes posterior normal density with the means and covariance matrix in (4). We use AGHQ to approximate the outer integral with respect to the univariate random effect, $u_{1j}^{*}$; see Appendix C for detail. Because $g(u_{1j}^{*}|Y_j)$ is nonstandard, we use the Bayes theorem to find (18) $$g(u_{1j}^{*}|Y_j)=\frac{h(Y_j|u_{1j}^{*})\varphi (u_{1j}^{*};0,{\tau }_{11})}{h(Y_j)}.$$(18)

By the invariance property of MLE, we translate the MLE of the “provisional” parameters $\theta =(\alpha ,\Omega ,{\gamma }_{10},{\sigma }_{\mathit{\text{cc}}},{\sigma }^2,{\tau }_{11})$ of model (14) to a one-to-one transformation ${\widehat{\theta }}_{(12)}^{*}$ via model (13).

3.3 Example 3: Cross-Level Interaction Effects Involving Partially Observed Covariates

Following Lee and Bryk (1989), we wish to extend the contextual effects model in two ways to allow the level-1 covariate: (i) to have random coefficients as in the previous section; and (ii) to interact with the level-2 covariate. We therefore write the model (19) $$\begin{array}{c}{R}_{\mathit{\text{ij}}}^{*}=({\gamma }_{00}+{\gamma }_{01}{\nu }_j+u_{0j})+({\gamma }_{10}+{\gamma }_{11}{\nu }_j+u_{1j}){ϵ}_{\mathit{\text{ij}}}^{*}+e_{\mathit{\text{ij}}}^{*},\\ C_{\mathit{\text{ij}}}^{*}=\delta +{\nu }_j+{ϵ}_{\mathit{\text{ij}}}^{*}\end{array}$$(19) where $u_{0j}$ and $u_{1j}$ are, as before, bivariate normal, but conditional on ν_j with variances ${\tau }_{00|\nu }$ and ${\tau }_{11|\nu }$, respectively, and covariance ${\tau }_{01|\nu }$. Other random effects are as defined in model (12). The parameters are ${\theta }_{(19)}^{*}=({\gamma }_{00},{\gamma }_{01},{\gamma }_{10},{\gamma }_{11},{\tau }_{00|\nu },{\tau }_{01|\nu },{\tau }_{11|\nu },{\sigma }^2,\delta ,{\tau }_{\nu \nu },{\sigma }_{\mathit{\text{cc}}})$.

This model involves two products of normal theory random effects: ${\nu }_j{ϵ}_{\mathit{\text{ij}}}^{*}$ and $u_{1j}{ϵ}_{\mathit{\text{ij}}}^{*}$. Using the logic of the last section, we wish to choose a PKRE such that, conditional on that effect, our CD joint model will be a normal-theory HLM. The question is how to choose this random effect. We might provisionally hold both ν_j and $u_{1j}$ constant, but we would prefer to minimize the dimension of the provisionally constant random effects such that the computational burden of numerical approximation is a minimum. Alternatively, we could provisionally hold ${ϵ}_{\mathit{\text{ij}}}^{*}$ constant. But while ${ϵ}_{\mathit{\text{ij}}}^{*}$ is a scalar, it varies across all level-1 units, which could be a very large set. Instead, we choose to provisionally hold constant the scalar level-2 random effect (20) $$u_{1j}^{*}={\gamma }_{11}{\nu }_j+u_{1j}.$$(20)

In addition, we define $u_{0j}^{*}={\gamma }_{01}{\nu }_j+u_{0j}$ to represent the model parsimoniously as $R_{\mathit{\text{ij}}}=({\gamma }_{00}+u_{0j}^{*})+({\gamma }_{10}+u_{1j}^{*}){ϵ}_{\mathit{\text{ij}}}^{*}+e_{\mathit{\text{ij}}}$. This model is therefore equivalent to model (12) to imply (13)–(18) and a one-to-one correspondence between ${\theta }_{(12)}^{*}$ and ${\theta }_{(19)}^{*}$.

Consequently, the CD model (12) is a one-to-one transformation of the provisional model $h(Y_{\mathit{\text{ij}}}^{*}|u_{1j}^{*})\varphi (u_{1j}^{*};0,{\tau }_{11})$ in (14) by distribution (13) and, also, of $h(Y_{\mathit{\text{ij}}}^{*}|{\nu }_j)\varphi ({\nu }_j;0,{\tau }_{\nu \nu })=f_1(R_{\mathit{\text{ij}}}^{*}|C_{\mathit{\text{ij}}}^{*},{\nu }_j)f_2(C_{\mathit{\text{ij}}}^{*}|{\nu }_j)\varphi ({\nu }_j;0,{\tau }_{\nu \nu })$ in (19) by (20). We choose to estimate joint model (14) via the PKRE for efficient computation that is guaranteed to be compatible with the scientific model $f_1(R_{\mathit{\text{ij}}}^{*}|C_{\mathit{\text{ij}}}^{*},{\nu }_j)$ by the one-to-one correspondence. Whereas scientific interest focuses on ${\theta }_{(19)}^{*}$, we will be estimating the parameters $\theta =(\alpha ,\Omega ,{\gamma }_{10},{\sigma }_{\mathit{\text{cc}}},{\sigma }^2,{\tau }_{11})$ of the provisional joint model (14). We then exploit the invariance property of MLE again, translating the MLE of θ back to those of ${\theta }_{(19)}^{*}$.

The standard approach of replacing ν_j and ${ϵ}_{\mathit{\text{ij}}}^{*}$ with ${\overline{C}}_j^{*}-{\overline{C}}^{*}$ and $C_{\mathit{\text{ij}}}^{*}-{\overline{C}}_j^{*}$, respectively, in model (19) produces biased estimation of $({\gamma }_{01},{\gamma }_{11},{\tau }_{00|\nu },{\tau }_{01|\nu },{\tau }_{11|\nu })$ even if $R_{\mathit{\text{ij}}}^{*}$ and $C_{\mathit{\text{ij}}}^{*}$ were fully observed; see Appendix A. Model (19) can readily incorporate multiple covariates having random effects and also having multiple cross-level interactions. Moreover, it is straightforward to include covariates having fixed coefficients. This is important given the need to add auxiliary information to strengthen the robustness of the MAR assumption.

4 Auxiliary Covariates

Our focus is on estimation of income inequality in achievement of model (19). Because certain variables such as family income have severe missing rates, it robustifies the MAR assumption to involve auxiliary covariates (e.g., parent occupation and pretest score) correlated with missing values or patterns (Collins, Schafer, and Kam 2003). We consider two approaches to augment auxiliary covariates, partially observed or measured with error, to the CD model. One approach is to assume such covariates to be linearly associated with the outcome and income. Violation of the linearity, however, may produce biased estimation. The other approach is to augment them as responses to $R_{\mathit{\text{ij}}}^{*}$, thereby allowing them to be nonlinearly associated with the outcome and income. We then transform the MLE of the CD model to those of a nested model (19).

4.1 Linearly Associated Auxiliary Covariates

To augment auxiliary covariates that are linearly associated with the outcome and income, we extend the scalar $C_{\mathit{\text{ij}}}^{*}$ of model (12) to a vector ${\mathit{C}}_{\mathit{\text{ij}}}^{*}={[C_{\mathit{\text{ij}}}^{*}{A}_{1\mathit{\text{ij}}}^{*T}{A}_{2j}^{*T}]}^T$ consisting of income $C_{\mathit{\text{ij}}}^{*}$ and auxiliary covariates $A_{1\mathit{\text{ij}}}^{*}$ at level 1 and $A_{2\mathit{\text{ij}}}^{*}$ at level 2. We then write the CD model (21) $$\begin{array}{c}{R}_{\mathit{\text{ij}}}^{*}=({\gamma }_{00}+u_{0j}^{*})+({\gamma }_{10}+u_{1j}^{*}){ϵ}_{\mathit{\text{ij}}}^{*}+{\gamma }_{20}^T{ϵ}_{1\mathit{\text{ij}}}^{*}+e_{\mathit{\text{ij}}}^{*},\\ {\mathit{C}}_{\mathit{\text{ij}}}^{*}=\mathit{\delta }+{\mathit{\nu }}_j+{\mathit{\epsilon }}_{\mathit{\text{ij}}}^{*}\end{array}$$(21) where $\mathit{\delta }={\left[\delta {\delta }_1^T{\delta }_2^T\right]}^T,{\mathit{\nu }}_j={\left[{\nu }_j{\nu }_{1j}^T{\nu }_{2j}^T\right]}^T$ and ${\mathit{\epsilon }}_{\mathit{\text{ij}}}^{*}={\left[{ϵ}_{\mathit{\text{ij}}}^{*}{ϵ}_{1\mathit{\text{ij}}}^{*T}{0}^T\right]}^T={[{\mathit{ϵ}}_{\mathit{\text{ij}}}^{*}{0}^T]}^T$ for the means ${[{\delta }_1^T{\delta }_2^T]}^T$ and school-specific random effects ${[{\nu }_{1j}^T{\nu }_{2j}^T]}^T$ of ${[A_{1\mathit{\text{ij}}}^{*T}{A}_{2j}^{*T}]}^T$ and the child-specific random effects ${ϵ}_{1\mathit{\text{ij}}}^{*}$ of $A_{1\mathit{\text{ij}}}^{*}$. Other components γ₀₀, γ₁₀, $u_{0j}^{*}={\gamma }_{01}^T{\mathit{\nu }}_j+u_{0j},u_{1j}^{*}={\gamma }_{11}^T{\mathit{\nu }}_j+u_{1j}$ and $e_{\mathit{\text{ij}}}^{*}$ are defined in model (19) except the linear effects γ₂₀ of ${ϵ}_{1\mathit{\text{ij}}}^{*}$. Again $e_{\mathit{\text{ij}}}^{*}$ and ${\mathit{ϵ}}_{\mathit{\text{ij}}}^{*}\sim N(0,{\Sigma }_{ϵ})$ are independent of each other and $u_j^{*}\sim N(0,\tau )$ for ${\Sigma }_{ϵ}=\left[\begin{array}{cc}{\sigma }_{\mathit{\text{cc}}}& {\Sigma }_{c1}\\ {\Sigma }_{1c}& {\Sigma }_{11}\end{array}\right]$ and $u_j^{*}={[u_{0j}^{*}{u}_{1j}^{*}{\mathit{\nu }}_j^T]}^T$; let $\text{cov}(u_{0j}^{*},{\mathit{\nu }}_j)={\tau }_{0\mathit{\nu }},\text{cov}(u_{1j}^{*},{\mathit{\nu }}_j)={\tau }_{1\mathit{\nu }}$ and ${\mathit{\nu }}_j\sim N(0,{\tau }_{\mathit{\nu \nu }})$. Auxiliary predictors $({\nu }_{1j},{\nu }_{2j})$ and ${ϵ}_{1\mathit{\text{ij}}}^{*}$ are linearly associated with the outcome and income at levels 2 and 1, respectively. The parameters are ${\theta }_{(21)}^{*}=({\gamma }_{00},{\gamma }_{10},{\gamma }_{20},\tau ,{\sigma }^2,\mathit{\delta },{\Sigma }_{ϵ})$.

We select $u_{1j}^{*}$ provisionally known again. Let $Y_{\mathit{\text{ij}}}^{*}={[R_{\mathit{\text{ij}}}^{*}{C}_{\mathit{\text{ij}}}^{*}{A}_{1\mathit{\text{ij}}}^{*T}]}^T$ and $A_{2j}^{*}$ be of respective lengths p₁ and p₂. Using model (13) now for $f(u_{0j}^{*},{\mathit{\nu }}_j|u_{1j}^{*})$ and ${\alpha }_{\mathit{\nu }|1}={[{\alpha }_{\nu |1}{\alpha }_{\nu 1|1}^T{\alpha }_{\nu 2|1}^T]}^T$, we find the provisional joint model (14) this time for $Y_j^{*}={[Y_{1j}^{*T}\cdots Y_{n_{j}j}^{*T}{A}_{2j}^{*T}]}^T$ where $$\begin{array}{c}{X}_j^{*}=\text{diag}\{{[X_{11j}^{*T}\cdots X_{1n_{j}j}^{*T}]}^T,X_{2j}^{*}\},\\ \alpha ={\left[{\alpha }_1^T{\alpha }_2^T\right]}^T,\\ Z_j^{*}=\text{diag}\{1_{n_j}\otimes I_{p_1},I_{p_2}\},\\ b_j={[u_{0j}^{*}-{\alpha }_{0|1}{u}_{1j}^{*}{\mathit{\nu }}_j^T-{\alpha }_{\mathit{\nu }|1}^T{u}_{1j}^{*}]}^T,\\ r_j^{*}={[r_{1j}^{*T}\cdots r_{n_{j}j}^{*T}{0}^T]}^T,\\ {\Psi }_j^{*}=\text{diag}\{I_{n_j}\otimes {\Sigma }_j^{*},0\}\end{array}$$for $X_{1\mathit{\text{ij}}}^{*}=[I_{p_1}{I}_{p_1}{u}_{1j}^{*}],X_{2j}^{*}=[I_{p_2}{I}_{p_2}{u}_{1j}^{*}],{\alpha }_1^T=[{\gamma }_{00}\delta {\delta }_1^T{\alpha }_{0|1}{\alpha }_{\nu |1}{\alpha }_{\nu 1|1}^T],{\alpha }_2^T={[{\delta }_2^T{\alpha }_{\nu 2|1}^T]}^T$, (22) $$\begin{array}{c}{r}_{\mathit{\text{ij}}}^{*}=\left[\begin{array}{cc}1& B_{1j}^T\\ 0& I_{p_1-1}\end{array}\right]\left[\begin{array}{c}{e}_{\mathit{\text{ij}}}^{*}\\ {\mathit{ϵ}}_{\mathit{\text{ij}}}^{*}\end{array}\right],\\ {\Sigma }_j^{*}=\left[\begin{array}{cc}{B}_{1j}^T{\Sigma }_{ϵ}{B}_{1j}+{\sigma }^2& B_{1j}^T{\Sigma }_{ϵ}\\ {\Sigma }_{ϵ}{B}_{1j}& {\Sigma }_{ϵ}\end{array}\right]\end{array}$$(22) denoting $B_{1j}^T=[{\gamma }_{10}+u_{1j}^{*}{\gamma }_{20}^T]$.

We estimate ${\theta }_{(21)}^{*}$ using its one-to-one transformation $\theta =(\alpha ,\Omega ,{\gamma }_{10},{\gamma }_{20},{\Sigma }_{ϵ},{\sigma }^2,{\tau }_{11})$ of $h(Y_j^{*}|u_{1j}^{*})\varphi (u_{1j}^{*};0,{\tau }_{11})$ by the EM algorithm, computing $h(Y_j)$ and $E(S_{\text{CD}j}|Y_j)$ by AGHQ as before. See Appendix B for the E step. The scalar PKRE $u_{1j}^{*}$ yields efficient computation.

To translate ${\theta }_{(21)}^{*}$ to ${\theta }_{(19)}^{*}$, we use model (13) to let ${\beta }_{\mathit{\text{kj}}}={\gamma }_{k0}+u_{\mathit{\text{kj}}}^{*}$ and find $E({\beta }_{\mathit{\text{kj}}}|{\nu }_j)={\gamma }_{k0}+{\gamma }_{k1}^{*}{\nu }_j$ and $\text{cov}({\beta }_{\mathit{\text{kj}}},{\beta }_{k\prime j}|{\nu }_j)={\tau }_{\mathit{\text{kk}}\prime }^{*}={\tau }_{\mathit{\text{kk}}\prime }-{\gamma }_{k1}^{*}{\gamma }_{k\prime 1}^{*}{\tau }_{\nu \nu }$ for ${\gamma }_{k1}^{*}={\tau }_{k\nu }/{\tau }_{\nu \nu }$ and $k,k\prime =0,1$. We then marginalize auxiliary ${ϵ}_{1\mathit{\text{ij}}}^{*}$ out to obtain $$\begin{array}{c}E(R_{\mathit{\text{ij}}}^{*}|{ϵ}_{\mathit{\text{ij}}}^{*},{\nu }_j)={\gamma }_{00}+{\gamma }_{01}^{*}{\nu }_j+({\gamma }_{10}+{\gamma }_{11}^{*}{\nu }_j){ϵ}_{\mathit{\text{ij}}}^{*}+{\gamma }_{20}^{T}E({ϵ}_{1\mathit{\text{ij}}}^{*}|{ϵ}_{\mathit{\text{ij}}}^{*}),\\ \text{var}(R_{\mathit{\text{ij}}}^{*}|{ϵ}_{\mathit{\text{ij}}}^{*},{\nu }_j)={\tau }_{00}^{*}+2{\tau }_{01}^{*}{ϵ}_{\mathit{\text{ij}}}^{*}+{\tau }_{11}^{*}{ϵ}_{\mathit{\text{ij}}}^{*2}\\ \hspace{1em}+{\gamma }_{20}^T\text{var}({ϵ}_{1\mathit{\text{ij}}}^{*}|{ϵ}_{\mathit{\text{ij}}}^{*}){\gamma }_{20}+{\sigma }^2\end{array}$$that should be of form $({\gamma }_{00}+{\gamma }_{01}{\nu }_j)+({\gamma }_{10}+{\gamma }_{11}{\nu }_j){ϵ}_{\mathit{\text{ij}}}^{*}$ and ${\tau }_{00|\nu }+2{\tau }_{01|\nu }{ϵ}_{\mathit{\text{ij}}}^{*}+{\tau }_{11|\nu }{ϵ}_{\mathit{\text{ij}}}^{*2}+{\sigma }^2$, respectively. See Appendix D for detail. We illustrate this approach in Sections 6 and 7.

4.2 Nonlinearly Associated Auxiliary Covariates

The linearity assumption between $A_{1\mathit{\text{ij}}}^{*}$ and $(R_{\mathit{\text{ij}}}^{*},C_{\mathit{\text{ij}}}^{*})$ in model (21) may be violated to produce biased estimation of ${\theta }_{(19)}^{*}$. In that case, we augment $A_{1\mathit{\text{ij}}}^{*}$ to multivariate ${\mathit{R}}_{\mathit{\text{ij}}}^{*}={[R_{\mathit{\text{ij}}}^{*}{A}_{1\mathit{\text{ij}}}^{*T}]}^T$ of length r, allowing them to be nonlinearly associated in (23) $${\mathit{R}}_{\mathit{\text{ij}}}^{*}=({\gamma }_{00}+u_{0j}^{*})+({\gamma }_{10}+u_{1j}^{*}){ϵ}_{\mathit{\text{ij}}}^{*}+e_{\mathit{\text{ij}}}^{*},e_{\mathit{\text{ij}}}^{*}\sim N(0,{\Sigma }_e)$$(23) and ${\mathit{C}}_{\mathit{\text{ij}}}^{*}={[C_{\mathit{\text{ij}}}^{*}{A}_{2j}^{*T}]}^T$ in (21) for conformable vectors γ₀₀ and γ₁₀ of fixed effects and random vectors $u_{0j}^{*}$ and $u_{1j}^{*}$ independent of $e_{\mathit{\text{ij}}}^{*}$ as before. It is straightforward to find the provisional joint model (14) for $Y_{\mathit{\text{ij}}}^{*}=({\mathit{R}}_{\mathit{\text{ij}}}^{*},C_{\mathit{\text{ij}}}^{*})$ and $A_{2j}^{*}$, selecting r-by-1 $u_{1j}^{*}$ to be provisionally known. We estimate the joint model efficiently and, then, translate the bivariate distribution of $R_{\mathit{\text{ij}}}^{*}$ and $C_{\mathit{\text{ij}}}^{*}$ to ${\theta }_{(19)}^{*}$; see Appendix D again for the translation. Numerical approximation is now intensive with respect to vector $u_{1j}^{*}$. Finally, some of $A_{1\mathit{\text{ij}}}^{*}$ may be linearly associated with the outcome and income while others may not. We illustrate this case in Data Analysis.

With additional known auxiliary covariates, we marginalize them out first given their expectation and covariance matrix estimated from sample before the translation above.

5 Within-Level Interactions and Polynomial Terms

We write a CD model including the level-2 interaction effects γ₀₄ of ${\nu }_j{\nu }_{1j}$ (24) $$\begin{array}{c}{R}_{\mathit{\text{ij}}}^{*}=({\gamma }_{00}+{\gamma }_{01}{\nu }_j+{\gamma }_{02}^T{\nu }_{1j}+{\gamma }_{03}^T{\nu }_{2j}+{\gamma }_{04}^T{\nu }_j{\nu }_{1j}+u_{0j})\\ \hspace{1em}+{\gamma }_{10}^T{\mathit{ϵ}}_{\mathit{\text{ij}}}^{*}+e_{\mathit{\text{ij}}}^{*}\end{array}$$(24) and ${\mathit{C}}_{\mathit{\text{ij}}}^{*}$ as in model (21) where $u_{0j}\sim N(0,{\tau }_{00|\mathit{\nu }})$ and ${\mathit{ϵ}}_{\mathit{\text{ij}}}^{*}$ has fixed effects γ₁₀ for simplicity. We select one interactive term ν_j, $\le$ the other in dimension, provisionally known to find (25) $$\begin{array}{c}\left[\begin{array}{c}{u}_{0j}\\ {\nu }_{1j}\\ {\nu }_{2j}\end{array}\right]|{\nu }_j\sim N(\left[\begin{array}{c}0\\ {\alpha }_{\nu 1|\nu }{\nu }_j\\ {\alpha }_{\nu 2|\nu }{\nu }_j\end{array}\right],\\ \Omega =\left[\begin{array}{ccc}{\tau }_{00|\mathit{\nu }}& 0& 0\\ 0& {\tau }_{\nu 1\nu 1|c}& {\tau }_{\nu 1\nu 2|c}\\ 0& {\tau }_{\nu 2\nu 1|c}& {\tau }_{\nu 2\nu 2|c}\end{array}\right]).\end{array}$$(25)

Let $b_j={[u_{0j}{\nu }_{1j}^T-{\alpha }_{\nu 1|\nu }^T{\nu }_j{\nu }_{2j}^T-{\alpha }_{\nu 2|\nu }^T{\nu }_j]}^T={[u_{0j}{b}_{1j}^T{b}_{2j}^T]}^T$ to find the model given ν_j (26) $$\begin{array}{c}{R}_{\mathit{\text{ij}}}^{*}=X_{\mathit{\text{Rj}}}^T{\alpha }_R+Z_{\mathit{\text{Rj}}}^T{b}_j+{\gamma }_{10}^T{\mathit{ϵ}}_{1\mathit{\text{ij}}}^{*}+e_{\mathit{\text{ij}}}^{*}\\ \sim N(X_{\mathit{\text{Rj}}}^T{\alpha }_R,Z_{\mathit{\text{Rj}}}^T\Omega Z_{\mathit{\text{Rj}}}+{\gamma }_{10}^T{\Sigma }_{ϵ}{\gamma }_{10}+{\sigma }^2)\\ C_{\mathit{\text{ij}}}^{*}=\delta +{\nu }_j+{ϵ}_{\mathit{\text{ij}}}^{*},\\ A_{1\mathit{\text{ij}}}^{*}={\delta }_1+{\alpha }_{\nu 1|\nu }{\nu }_j+b_{1j}+{ϵ}_{1\mathit{\text{ij}}}^{*},\\ A_{2j}^{*}={\delta }_2+{\alpha }_{\nu 2|\nu }{\nu }_j+b_{2j}\end{array}$$(26) for $X_{\mathit{\text{Rj}}}^T=[1{\nu }_j{\nu }_j^2],{\alpha }_R={[{\gamma }_{00}{\gamma }_{01}+{\gamma }_{02}^T{\alpha }_{\nu 1|\nu }+{\gamma }_{03}^T{\alpha }_{\nu 2|\nu }{\gamma }_{04}^T{\alpha }_{\nu 1|\nu }]}^T$ and $Z_{\mathit{\text{Rj}}}^T=\left[1{\gamma }_{02}^T+{\gamma }_{04}^T{\nu }_j{\gamma }_{03}^T\right]$ where $C_{\mathit{\text{ij}}}^{*}$ varies within but not between clusters. As in Section 4.1, we stack these equations to find the implied provisional model (14) for $Y_j^{*}={[Y_{1j}^{*T}\cdots Y_{n_j}^{*T}{A}_{2j}^{*T}]}^T$, and compute $h(Y_j)$, setting $S_{\text{CD}j}=1$ below, and $E(S_{\text{CD}j}|Y_j)$ by (27) $$h(Y_j)E(S_{\text{CD}j}|Y_j)=\int E(S_{\text{CD}j}|{\nu }_j,Y_j)h(Y_j|{\nu }_j)\varphi ({\nu }_j;0,{\tau }_{\mathit{\text{cc}}})d{\nu }_j$$(27) numerically for $h(Y_j|{\nu }_j)$ from the provisional model.

We now consider another CD model including level-1 interaction effects (28) $$R_{\mathit{\text{ij}}}^{*}=({\gamma }_{00}+{\gamma }_{01}^T{\mathit{\nu }}_j+u_{0j})+{\gamma }_{10}{ϵ}_{\mathit{\text{ij}}}^{*}+({\gamma }_{20}^T+{\gamma }_{30}^T{ϵ}_{\mathit{\text{ij}}}^{*}){ϵ}_{1\mathit{\text{ij}}}^{*}+e_{\mathit{\text{ij}}}^{*}$$(28) and ${\mathit{C}}_{\mathit{\text{ij}}}^{*}$ in model (21) where ${ϵ}_{\mathit{\text{ij}}}^{*}$ and ${ϵ}_{1\mathit{\text{ij}}}^{*}$ have main (γ₁₀ and γ₂₀) and interaction effects (γ₃₀). We select an interactive term ${ϵ}_{\mathit{\text{ij}}}^{*},\le {ϵ}_{1\mathit{\text{ij}}}^{*}$ in dimension, provisionally known again, and find $${ϵ}_{1\mathit{\text{ij}}}^{*}|{ϵ}_{\mathit{\text{ij}}}^{*}\sim N({\alpha }_{ϵ1|c}{ϵ}_{\mathit{\text{ij}}}^{*},{\Sigma }_{1|c})$$

for ${\alpha }_{ϵ1|c}={\Sigma }_{1c}{\sigma }_{\mathit{\text{cc}}}^{-1}$ and ${\Sigma }_{1|c}={\Sigma }_{11}-{\alpha }_{ϵ1|c}{\sigma }_{\mathit{\text{cc}}}{\alpha }_{ϵ1|c}^T$. Given ${ϵ}_{\mathit{\text{ij}}}^{*}$ provisionally constant, let $u_{0j}^{*}={\gamma }_{01}^T{\mathit{\nu }}_j+u_{0j}$ and $a_{1\mathit{\text{ij}}}^{*}={ϵ}_{1\mathit{\text{ij}}}^{*}-{\alpha }_{ϵ1|c}{ϵ}_{\mathit{\text{ij}}}^{*}$ to express (29) $$\begin{array}{c}{R}_{\mathit{\text{ij}}}^{*}=X_{\mathit{\text{Rij}}}^T{\alpha }_R+u_{0j}^{*}+B_{1\mathit{\text{ij}}}^T{a}_{1\mathit{\text{ij}}}^{*}+e_{\mathit{\text{ij}}}^{*}\\ \sim N(X_{\mathit{\text{Rij}}}^T{\alpha }_R,{\tau }_{00}+B_{1\mathit{\text{ij}}}^T{\Sigma }_{1|c}{B}_{1\mathit{\text{ij}}}+{\sigma }^2),\\ C_{\mathit{\text{ij}}}^{*}|{ϵ}_{\mathit{\text{ij}}}^{*}\sim N(\delta +{ϵ}_{\mathit{\text{ij}}}^{*},{\tau }_{\mathit{\text{cc}}}),\\ A_{1\mathit{\text{ij}}}^{*}={\delta }_1+{\alpha }_{ϵ1|c}{ϵ}_{\mathit{\text{ij}}}^{*}+{\nu }_{1j}+a_{1\mathit{\text{ij}}}^{*},\\ A_{2j}^{*}={\delta }_2+b_{2j}\end{array}$$(29) where $X_{\mathit{\text{Rij}}}^T=[1{ϵ}_{\mathit{\text{ij}}}^{*}{ϵ}_{\mathit{\text{ij}}}^{*2}],{\alpha }_R=[{\gamma }_{00}{\gamma }_{10}+{\gamma }_{20}^T{\alpha }_{ϵ1|c}{\gamma }_{30}^T{\alpha }_{ϵ1|c}],B_{1\mathit{\text{ij}}}^T={\gamma }_{20}^T+{\gamma }_{30}^T{ϵ}_{\mathit{\text{ij}}}^{*},\text{cov}(R_{\mathit{\text{ij}}}^{*},A_{1\mathit{\text{ij}}}^{*}|b_j,{ϵ}_{\mathit{\text{ij}}}^{*})=B_{1\mathit{\text{ij}}}^T{\Sigma }_{1|c}$ and $C_{\mathit{\text{ij}}}^{*}$ now varies between but not within clusters to imply the provisional model (14) given ${ϵ}_j^{*}=({ϵ}_{1j}^{*},\dots ,{ϵ}_{n_{j}j}^{*})$ for $b_j={[u_{0j}^{*}{\mathit{\nu }}_j^T]}^T$. We compute $h(Y_j)$ and $E(S_{\text{CD}j}|Y_j)$ by (30) $$h(Y_j)E(S_{\text{CD}j}|Y_j)=\int E(S_{\text{CD}j}|{ϵ}_j^{*},Y_j)h(Y_j|{ϵ}_j^{*})\varphi ({ϵ}_j^{*};0,I_{n_j}{\sigma }_{\mathit{\text{cc}}})d{ϵ}_j^{*}$$(30) numerically for $h(Y_j|{ϵ}_j^{*})$ from the provisional model as before. The numerical integral can be computationally intensive, in particular, given large cluster sizes; multivariate Laplace approximation (Pinheiro and Bates 1995; Raudenbush, Yang, and Yosef 2000) and parallel computation of each cluster may result in efficient computation.

5.1 Rules for Choosing Provisionally Known Random Effects

We now provide general rules for selecting PKREs:

1. For an interaction ${ϵ}_{\mathit{\text{ij}}}^{*}{ϵ}_{1\mathit{\text{ij}}}^{*}$ as in (28), we hold ${ϵ}_{\mathit{\text{ij}}}^{*},\le {ϵ}_{1\mathit{\text{ij}}}^{*}$ in dimension, constant. The resulting model quadratic in ${ϵ}_{\mathit{\text{ij}}}^{*}$ will minimize the dimension of AGHQ;

2. For a level-2 interaction ${\nu }_j{\nu }_{1j}$, we hold one with a smaller dimension constant again;

3. For a three-way interaction ${ϵ}_{\mathit{\text{ij}}}^{*}{ϵ}_{1\mathit{\text{ij}}}^{*}{ϵ}_{2\mathit{\text{ij}}}^{*}$ at level 1, hold two terms $\le$ the third one in dimension constant and this applies to a three-way interaction at level 2, too;

4. For cross-level interactions ${\mathit{\nu }}_j{ϵ}_{\mathit{\text{ij}}}^{*}$ as in model (21), hold $u_{1j}^{*}={\gamma }_{11}^T{\mathit{\nu }}_j+u_{1j}$ constant;

5. For the cluster-specific effects $u_{1j}^{*}$ of ${ϵ}_{\mathit{\text{ij}}}^{*}$, we hold $u_{1j}^{*}$ constant;

6. Finally, for a subset of these effects, hold constant the union of the PKREs of the subset.

Our CD model in each case includes a scientific model of interest and induces a provisional joint model (14). Because the models are one-to-one transformations of each other, the scientific model is guaranteed to be compatible with the joint model we estimate.

6 Data Analysis

Rising income inequality in the United States and other nations has recently attracted substantial attention (Piketty 2014). A key question involves the consequence of such inequality for equality of opportunity among children. Following past research, we decompose the association between family income and educational achievement into a contextual component and a child-specific component (Firebaugh 1978; Willms 1986; Lee and Bryk 1989). The contextual component reflects the fact that elementary schools in the United States are quite segregated based on family income. Such segregation reflects and may reinforce residential segregation as a function of family income. Two children having the same family income might differ in educational achievement as a result of their experience in low-income versus high-income schools. The individual component reflects socioeconomic inequality within schools. Children attending the same school who differ with respect to family income may tend to differ with respect to their achievement. However, the magnitude of this within-school disparity may vary from school to school (Raudenbush and Bryk 1986; Lee and Bryk 1989). The contextual effects model (Willms 1986) supports the composition of inequality in achievement that we seek.

First, we decompose family income for child i in school j into between-school and within-school components as in model (19) $$\begin{array}{c}\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}({\text{income}}_{\mathit{\text{ij}}})=C_{\mathit{\text{ij}}}^{*}=\delta +{\nu }_j+{ϵ}_{\mathit{\text{ij}}}^{*},\\ \text{math}{S}_{\mathit{\text{ij}}}=R_{\mathit{\text{ij}}}^{*}=({\gamma }_{00}+{\gamma }_{01}{\nu }_j+u_{0j})\\ +({\gamma }_{10}+{\gamma }_{11}{\nu }_j+u_{1j}){ϵ}_{\mathit{\text{ij}}}^{*}+e_{\mathit{\text{ij}}}^{*}\end{array}$$for the mean of log-income δ, the school-specific deviation from the mean ${\nu }_j\sim N(0,{\sigma }_{\mathit{\text{cc}}})$, and the child-specific component ${ϵ}_{\mathit{\text{ij}}}^{*}\sim N(0,{\sigma }_{\mathit{\text{cc}}})$. A child’s mathematics achievement in spring 1999 (mathS) depends on these components via the model (19) for $\text{math}{S}_{\mathit{\text{ij}}}=R_{\mathit{\text{ij}}}^{*}$. The parameters are ${\theta }_{(19)}^{*}=({\gamma }_{00},{\gamma }_{01},{\gamma }_{10},{\gamma }_{11},{\tau }_{00|\nu },{\tau }_{01|\nu },{\tau }_{11|\nu },{\sigma }^2,\delta ,{\tau }_{\nu \nu },{\sigma }_{\mathit{\text{cc}}})$.

We choose math achievement as our outcome because of its importance in predicting educational attainment and adult earnings (Nomi and Raudenbush 2016; Rivera-Batiz 1982). In model (19), γ₀₁ is the between-school gradient, reflecting the expected difference in $R_{\mathit{\text{ij}}}^{*}$ associated with a unit difference in school mean income; γ₁₀ reflects the average within-school gradient. However, the within-school gradient may depend on school mean income, an interaction effect represented by γ₁₁, and this gradient may also vary randomly over schools as represented by $u_{1j}\sim N(0,{\tau }_{11|\nu })$. School-mean achievement is γ₀₀ and, conditional on income, varies randomly over schools, $u_{0j}\sim N(0,{\tau }_{00|\nu })$. If the within-school gradients were constant (${\gamma }_{11}={\tau }_{11}=0$), the overall linear coefficient for income will be $$E(R_{\mathit{\text{ij}}}^{*}|C_{\mathit{\text{ij}}}^{*}=c+1)-E(R_{\mathit{\text{ij}}}^{*}|C_{\mathit{\text{ij}}}^{*}=c)=\rho {\gamma }_{01}+(1-\rho ){\gamma }_{10}$$where $\rho ={\tau }_{\mathit{\text{cc}}}/({\tau }_{\mathit{\text{cc}}}+{\sigma }_{\mathit{\text{cc}}})$ can be regarded as an index of school segregation as a function of income. Define ${\gamma }_c={\gamma }_{01}-{\gamma }_{10}$ as the “contextual effect” (Willms 1986), the expected difference in math achievement between two students with the same family income who attend two schools that differ by one unit in school mean income. A nation’s income gradient would be $\rho {\gamma }_c+{\gamma }_{10}$, which increases with the within school segregation based on income, ρ, the contextual coefficient, γ_c and within-school gradient γ₁₀. This simple relationship will not hold if the within-school gradients vary over schools, and one purpose of our analysis is to test that proposition.

To do so, we use data from 21,211 children attending kindergarten in 1,018 schools as of fall 1998, a nationally representative sample known as the Early Childhood Longitudinal Study of 1998 (“ECLS”) that is publicly available at https://nces.ed.gov/ecls; see Table 1. Only 8% of the math achievement data are missing. However, family income data are missing for 32% of the sample, a finding that is quite typical in surveys of educational achievement. Fortunately, ECLS (Tourangeau et al. 2009) provides data on auxiliary variables, including the maximum occupational status score of parents (occupation), missing for only 5% of the cases, as well as math achievement in fall 1998 (mathF), which is strongly predictive of math achievement in spring 1999 (mathS). In all, we have four auxiliary variables, correlated with income, outcome or missing patterns.

Table 1 Each variable for analysis with mean (standard deviation (SD), missing %).

Symbol	Name	Mean (SD, missing %)	Symbol	Name	Mean (SD, missing %)
$R_{\mathit{\text{ij}}}^{*}$	mathS	0.00 (1.00, 8)	$A_{1\mathit{\text{ij}}}^{*}$	occupation	0.00 (1.00, 5)
$C_{\mathit{\text{ij}}}^{*}$	income	10.73 (0.71, 32)	$A_{1\mathit{\text{ij}}}^{*}$	age	0.00 (1.00, 6)
$A_{1\mathit{\text{ij}}}^{*}$	mathF	0.00 (1.00, 12)	$A_{2j}^{*}$	enrollment	0.00 (1.00, 17)

In this article, we take convergence to ML to be less than ${10}^{-4}$ in the square root of the summed squared differences between $\widehat{\theta }$ of two consecutive iterations. We estimate the model for covariates ${\mathit{C}}_{\mathit{\text{ij}}}^{*}$ using all observed values by ML via the EM algorithm (Shin and Raudenbush 2007) and a scientific model for $R_{\mathit{\text{ij}}}^{*}$ given the covariates and their sample cluster and overall means by complete case analysis, and transform the estimates to the initial values $\widehat{\theta }$ of the joint model. We carry out complete case analysis by R (RCoreTeam 2017), estimate θ on a Dell XPS laptop with the 11th generation Intel(R) Core(TM) i9-11900H processor at 2.50GHz and 64 GB RAM, and test a hypothesis at a level $\alpha =0.05$.

6.1 Linearly Associated Auxiliary Covariates

Recall from Section 4 that we have two strategies. Following Section 4.1, we model the auxiliary covariates linearly associated with the outcome and income in model (21) where $A_{1\mathit{\text{ij}}}^{*}$ is a vector of mathF, occupation, and age in months at assessment of spring 1999 (age) and $A_{2j}^{*}$ is the square root of kindergarten enrollment (enrollment) by the Box-Cox transformation. Consequently, ${\theta }_{(21)}^{*}=({\gamma }_{00},{\gamma }_{10},{\gamma }_{20},\tau ,{\sigma }^2,\mathit{\delta },{\Sigma }_{ϵ})$ consists of 10 fixed effects (3-by-1 γ₂₀ and 5-by-1 $\mathit{\delta }$) and 39 variances and covariances (7-by-7 τ and 4-by-4 ${\Sigma }_{ϵ}$).

We standardized each variable to have mean 0 and variance 1, except income for interpretation. Estimation of ${\theta }_{(21)}^{*}$ with a provisionally known $u_{1j}^{*}$ and 20 abscissas converged fast to ML in 9 iterations and 13 sec. The transformed estimates ${\widehat{\theta }}_{(19)}^{*}$ and standard errors (SEs), multiplied by 100, are listed under “EM-AGHQ I” in Table 2. The school-mean and within-school components of family income are positively associated with math achievement while their interaction effect is insignificant. The within-school income effects appear to vary at most modestly across schools with the variance estimate ${\widehat{\tau }}_{11|\nu }=0.19$ less than the associated SE 0.25. Based on $\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}({\widehat{\tau }}_{11|\nu })\sim N[\mathrm{\text{ln}}\hspace{0.17em}({\tau }_{11|\nu }),\text{var}({\widehat{\tau }}_{11|\nu })/{\tau }_{11|\nu }^2]$, we find a large-sample 95% confidence interval (CI) for ${\tau }_{11|\nu }$: $(0.01,2.62)$ near zero.

Table 2 Estimates × 100 (standard errors × 100) of model (19) by EM-AGHQ I and II.

Predictor	Parameter	EM-AGHQ I	EM-AGHQ II
1	γ00	0.36 (2.38)	0.20 (2.46)
νj	γ01	78.69 (2.93)**	77.22 (3.01)**
${ϵ}_{\mathit{\text{ij}}}^{*}$	γ10	28.06 (1.36)**	37.38 (1.74)**
${ϵ}_{\mathit{\text{ij}}}^{*}{\nu }_j$	γ11	0.77 (2.03)	14.33 (3.58)**
	${\tau }_{00|\nu }$	9.30	9.64
	${\tau }_{01|\nu }$	0.64	2.35
	${\tau }_{11|\nu }$	0.19 (0.25)	3.08 (0.69)**
	95% CI for ${\tau }_{11|\nu }$	(0.01, 2.62)	(1.98, 4.78)
	${\sigma }^2$	75.21	73.74
	$\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}$ joint ML	–107361.00	–107177.03

${}^{*}$: p-value <0.05; ${}^{**}$: p-value <0.01.

To test this model against the null hypothesis that ${\gamma }_{11}=0$ and ${\tau }_{11|\nu }=0$, we estimated the null model, the multivariate normal distribution of linear associated $(R_{\mathit{\text{ij}}}^{*},C_{\mathit{\text{ij}}}^{*},A_{1\mathit{\text{ij}}}^{*},A_{2j}^{*})$, efficiently by the EM algorithm (Shin and Raudenbush 2007, 2010). The model consisted of 42 parameters comprising 6 fixed intercepts and 6-by-6 level-2 and 5-by-5 level-1 variance covariance matrices, and converged to log ML -107370.60. Compared to the log ML of ${\theta }_{(21)}^{*}$ displayed at the last row of Table 2, the likelihood ratio test statistic to test the two joint models is 19.20 with 7 degrees of freedom to give a conservative p-value < 0.01 (Stram and Lee 1994). Therefore, we infer that the outcome is nonlinearly associated with income. Lastly, estimation using 10 abscissas also produced the same estimates under EM-AGHQ I.

6.2 Nonlinearly Associated Auxiliary Covariates

Preliminary analysis indicates that mathF and occupation may be nonlinearly associated with the outcome and income. To test the hypothesis, following Section 4.2, we model multivariate responses ${\mathit{R}}_{\mathit{\text{ij}}}^{*}=$ (mathS, mathF, occupation) and $A_{1\mathit{\text{ij}}}^{*}=$ age in model (21): (31) $${\mathit{R}}_{\mathit{\text{ij}}}^{*}=({\gamma }_{00}+u_{0j}^{*})+({\gamma }_{10}+u_{1j}^{*}){ϵ}_{\mathit{\text{ij}}}^{*}+{\gamma }_{20}{ϵ}_{1\mathit{\text{ij}}}^{*}+e_{\mathit{\text{ij}}}^{*}$$(31) and ${\mathit{C}}_{\mathit{\text{ij}}}^{*}={[C_{\mathit{\text{ij}}}{A}_{1\mathit{\text{ij}}}^{*}{A}_{2j}^{*}]}^T$ as before for 3-by-1 vectors γ₀₀, γ₁₀ and γ₂₀ of fixed effects and random vectors $u_{0j}^{*},u_{1j}^{*}$ and $e_{\mathit{\text{ij}}}^{*}\sim N(0,{\Sigma }_e)$. Therefore, τ is 9-by-9, Σ_e 3-by-3 and ${\Sigma }_{ϵ}$ 2-by-2; the CD model comprises 12 fixed effects and 54 variances and covariances. With 3-by-1 $u_{1j}^{*}$ provisionally known, we estimated the joint model using 10 abscissae per dimension which converged to ML in 734th iterations and 928 min. The log ML of ${\theta }_{(31)}^{*}$ is shown at the bottom row under EM-AGHQ II. The LRT statistic to test $H_0:$ model (21) versus $H_1:$ model (31) is 367.93. The conservative LRT with 17 degrees of freedom (Stram and Lee 1994) produces a p-value 0 to reject the null in favor of mathF or occupation nonlinearly associated with the outcome and income.

The translated estimates ${\widehat{\theta }}_{(19)}^{*}$ and SEs are listed under EM-AGHQ II in Table 2. Compared to those under EM-AGHQ I, the main effect of within-school income is larger; furthermore, the interaction effect is significant, and so is the random effects of income by the Wald test to produce a 95% CI for τ₁₁ now distant from zero. We conclude that the linearity assumption associated with (21) is violated to attenuate the main, interaction and random effects of within-school income ${ϵ}_{\mathit{\text{ij}}}^{*}$, confounded with the auxiliary covariates. As a result, EM-AGHQ II produces a smaller ${\widehat{\sigma }}^2=73.74$ and, thus, explains more outcome variability within schools than does EM-AGHQ I.

6.3 Known Auxiliary Covariates

Either model (21) or (31) may be extended to control for known auxiliary covariates such as race ethnicity and gender at level 1 and school location and sector at level 2. The joint model may be estimated given the provisionally known $u_{1j}^{*}$ again, and translated to ${\widehat{\theta }}_{(19)}^{*}$ and SEs. See Appendices B and D for detail.

7 Simulation Study

We focus on ML estimation of the scientific HLM (19) after simulating outcome $R_{\mathit{\text{ij}}}^{*}$ and income $C_{\mathit{\text{ij}}}^{*}$ from a joint model conditional on auxiliary covariates within which the HLM is nested. The goal is to compare our estimators (EM-AGHQ) with those by four methods: (a) the benchmark method (BM) given ν_j and ${ϵ}_{\mathit{\text{ij}}}^{*}$; (b) complete-case analysis (CC) given ${\overline{C}}_j^{*}-{\overline{C}}^{*}$ and $C_{\mathit{\text{ij}}}^{*}-{\overline{C}}_j^{*}$ instead; (c) MLE on MI (Shin and Raudenbush 2007, 2010); and (d) the Gibbs sampler (GS) of Enders, Du, and Keller (2020) implemented in software Blimp (Keller and Enders 2021). BM is based on complete data while others are based on data MAR. Therefore, a good method will produce estimates near the BM counterparts. BM and CC estimate the scientific model by the lme4 package (Bates et al. 2015) in R. MLE on MI uses C programs to estimate incompatible MHLM (8) by ML and impute missing values including latent school mean incomes 20 times, more than did past multilevel missing data analyses (Schafer and Yucel 2002; Shin and Raudenbush 2007, 2013; Enders, Du, and Keller 2020), from their predictive distribution given observed data implied by the MHLM at ML (Shin and Raudenbush 2007, 2010); and estimates the HLM given the MI by lme4. GS estimates the joint model, simultaneously generating 20 imputations of missing values excluding latent school mean incomes, by Blimp and, then, the HLM (19) given the MI by Blimp again. EM-AGHQ estimates the joint model by our C program and translates the estimates to the desired MLE in R.

We simulate the ECLS $R_{\mathit{\text{ij}}}^{*}$ and $C_{\mathit{\text{ij}}}^{*}$ closely in terms of sample sizes, correlations and missing rates. Specifically, for n = 20 children in each of J = 1000 schools, we simulate: (i) known auxiliary covariates $X_{21j}\sim \text{Bernoulli}(0.3)$ and $X_{1\mathit{\text{ij}}}\sim \text{Bernoulli}(0.45)$ equal to 1 (0) for a private (public) school and a minority (white) child, respectively, where $X_{1\mathit{\text{ij}}}$ varies within, but not between, schools for simplicity; (ii) random intercepts and slope ${\beta }_{0j}={\gamma }_{00}^T{X}_{2j}+u_{0j}^{*},{\beta }_{1j}={\gamma }_{10}^T{X}_{2j}+u_{1j}^{*}$ and ${\beta }_{\mathit{\text{Cj}}}={\delta }_{00}^T{X}_{2j}+{\nu }_j$ from $N(1+X_{21j},1)$ with covariances 0.8 for $X_{2j}={[1X_{21j}]}^T$; (iii) independent $e_{\mathit{\text{ij}}}^{*},{ϵ}_{\mathit{\text{ij}}}^{*}\sim N(0,10)$ to simulate the CD joint model (32) $$\begin{array}{c}{R}_{\mathit{\text{ij}}}^{*}={\beta }_{0j}+{\beta }_{1j}(C_{\mathit{\text{ij}}}^{*}-{\beta }_{\mathit{\text{Cj}}})+{\gamma }_{20}{X}_{1\mathit{\text{ij}}}+e_{\mathit{\text{ij}}}^{*},\\ C_{\mathit{\text{ij}}}^{*}={\beta }_{\mathit{\text{Cj}}}+{\delta }_{10}{X}_{1\mathit{\text{ij}}}+{ϵ}_{\mathit{\text{ij}}}^{*}.\end{array}$$(32)

The simulated parameters in θ consist of ${\gamma }_{00}^T={\gamma }_{10}^T={\delta }_{00}^T=[11],{\gamma }_{20}={\delta }_{10}=1$, variances ${\tau }_{00}={\tau }_{11}={\tau }_{\nu \nu }=1$ and ${\sigma }_{\mathit{\text{cc}}}={\sigma }^2=10$, and covariances ${\tau }_{01}={\tau }_{0\nu }={\tau }_{1\nu }=0.8$. We marginalize $X_{1\mathit{\text{ij}}}$ and $X_{21j}$ out given their simulated expectations and variances, and translate θ to the scientific model in column two of Table 3 as explained in Appendix D.

Table 3 Scientific model (19) estimated by BM, CC, MLE on MI, GS and EM-AGHQ.

covariate	simulated	BM	CC	MLE on MI	GS	EM-AGHQ
1	γ00=2.34	0.06% (0.03)	–8.32% (0.05)	–8.16% (0.05)	–0.37% (0.07)	0.04% (0.06)
		0.03, 0.96	0.08, 0.11	0.06, 0.09	0.07, 0.94	0.07, 0.92
νj	γ01=1.21	0.12% (0.03)	0.29% (0.03)	–5.81% (0.06)	–0.59% (0.05)	0.07% (0.05)
		0.03, 0.95	0.05, 0.85	0.06, 0.76	0.05, 0.94	0.05, 0.94
${ϵ}_{\mathit{\text{ij}}}^{*}$	γ10=1.32	0.10% (0.02)	0.16% (0.03)	–5.14% (0.03)	–0.26% (0.04)	0.06% (0.03)
		0.02, 0.95	0.04, 0.85	0.03, 0.33	0.04, 0.95	0.04, 0.93
${ϵ}_{\mathit{\text{ij}}}^{*}{\nu }_j$	γ11=0.83	–0.07% (0.02)	–42.05% (0.02)	−54.00% (0.03)	–0.52% (0.03)	–0.04% (0.03)
		0.02, 0.96	0.02, 0.00	0.02, 0.00	0.03, 0.95	0.03, 0.95
	$\sqrt{{\tau }_{00|\nu }}$=0.77	0.00% (-)	60.31% (-)	22.83% (-)	0.36% (-)	0.09% (0.05)
		0.03, -	0.06, -	0.05, -	0.05, -	0.05, 0.96
	${\tau }_{01|\nu }$=0.33	–0.05% (-)	–28.75% (-)	–34.68% (-)	–0.64% (-)	–0.29% (0.04)
		0.02, -	0.06, -	0.03, -	0.04, -	0.04, 0.95
	$\sqrt{{\tau }_{11|\nu }}$=0.61	–0.13% (-)	39.00% (-)	7.31% (-)	2.69% (-)	–0.27% (0.03)
		0.01, -	0.02, -	0.02, -	0.03, -	0.03, 0.95
	$\sqrt{{\sigma }^2}$=3.20	–0.02% (-)	0.07% (-)	25.39% (-)	0.02% (-)	–0.05% (0.02)
		0.02, -	0.02, -	0.05, -	0.02, -	0.02, 0.96

NOTE: Each estimate or cell occupies two rows: % bias (average estimated SE) in the first row and the empirical estimate of true SE over simulated samples and a 95% converge probability of the estimator in the next.

Next, we simulate the ECLS missing rates closely by (33) $$\text{logit}(p_{\mathit{\text{ij}}})={\varphi }_1{X}_{1\mathit{\text{ij}}}+{\varphi }_2(1-X_{1\mathit{\text{ij}}})+z_j,z_j\sim N(0,1)$$(33) given the known level-1 covariate $X_{1\mathit{\text{ij}}}$: missing values drawn from $\text{Bernoulli}(p_{\mathit{\text{ij}}})$ are MAR. Because the mechanism does not provide information about the simulated model (32), the parameter spaces of the missing data mechanism and joint model are also distinct. We simulate higher missing rates for minority than white students by ${\varphi }_1>{\varphi }_2$: ${\varphi }_1=-0.2>{\varphi }_2=-1.2$ for income $C_{\mathit{\text{ij}}}^{*}$ with a 35% missing rate (46% for $X_{1\mathit{\text{ij}}}=1$, 27% for $X_{1\mathit{\text{ij}}}=0$); and ${\varphi }_1=-2>{\varphi }_2=-3$ for response $R_{\mathit{\text{ij}}}^{*}$ with an 11% missing rate (16% for $X_{1\mathit{\text{ij}}}=1$, 7% for $X_{1\mathit{\text{ij}}}=0$) on average.

We repeated simulating data and estimating the scientific model by the approaches 500 times to compute the % bias, average estimated SE (ASE), empirical estimate of the true SE (ESE) over samples and coverage probability (coverage) of each estimator in the next five columns. Each cell or estimate occupies two rows: % bias (ASE) in the first, and ESE and coverage in the next row. The lme4 package is unable to produce ESE and coverage of a variance or covariance estimate. The BM estimates are of course very accurate and precise with $\le 0.13$% bias, small ASE close to ESE, and good coverages near the nominal 0.95 in column three.

The CC estimates in column four, however, are biased despite the large sample sizes. The standard deviations (SDs) $\sqrt{{\tau }_{00|\nu }}$ and $\sqrt{{\tau }_{11|\nu }}$ are 60% and 39% biased upward while the intercept γ₀₀, interaction effect γ₁₁ and covariance ${\tau }_{01|\nu }$ are 8%, 42%, and 29% biased downward, respectively. Only the estimates of γ₀₁, γ₁₀, and ${\sigma }^2$ are comparable in accuracy to those by BM as (A1) reveals in Appendix A. The coverages are low with a zero coverage for γ₁₁. Finally, the uncertainty associated with the estimator of a cluster-level effect γ₀₁ seems underestimated by ASE smaller than ESE.

MLE on MI generates incompatible MI based on MHLM (8) without consideration of the interaction effect and PKRE, thereby producing all estimates biased in column five that do not seem better than the CC estimates. Fixed effects and covariance are biased downward, and SDs upward. SEs are close to but coverages lower than CC counterparts.

On the contrary, GS without consideration of a PKRE produces estimates in column six nearly as accurate as BM counterparts except the SD of the random slope that is biased upward by 2.69% while EM-AGHQ yields all estimates in the last column as accurate or almost as accurate as BM estimates. Overall, both approaches produce estimates slightly less precise than BM estimates; ASE and ESE appear slightly larger than those by BM to reflect extra uncertainty due to latent covariates and missing values.

Computation

We used 20 abscissas to estimate the joint model by EM-AGHQ. Given data MAR, the estimation converged 450 times taking 101.6 iterations on average and 189 iterations at maximum, but did not converge until and was stopped to produce the estimates at the 300th iteration 50 times (10%). This does not appear to be the weakness of our approach as the convergence issue also occurred to each of BM and CC estimations producing 50 or more warnings of a model failing to converge. The convergence issue seems partly due to high missing rates but few covariates to explain missing values and patterns. In our experience thus far, the convergence rates seem positively associated with more covariates or abscissas. For example, this simulation using 10 abscissas resulted in practically identical estimates, but lowered the convergence rate given data MAR.

Blimp estimates 21 models per simulated dataset: joint model (32) and the HLM given each of 20 imputations. We set 20,000 burn-in and 10,000 post burn-in iterations to estimate the joint model and impute MI, and 5000 burn-in and 5000 post burn-in iterations to estimate the HLM given the MI. These settings are based on preliminary analysis of five simulated datasets that produced the potential scale reduction statistics of all estimates of each model lower than or near 1.1 to imply a reasonable convergence to posterior distributions (Gelman and Rubin 1992).

8 Discussion

In this article, we have considered how to estimate a two-level hierarchical linear model (HLM) efficiently where a continuous response $R^{*}$ and continuous covariates $C^{*}$ may be MAR and $C^{*}$ may have interactive, polynomial or randomly varying effects. Nonlinearities of $C^{*}$ imply a nonstandard joint model $h(R^{*},C^{*})=h(Y^{*})$ where $Y^{*}=(Y,Y_{\text{mis}})$ for observed Y and missing $Y_{\text{mis}}$. The key idea is to introduce a unique factorization of the joint model involving “provisionally known” random effects (PKREs) u such that the observed joint model $h(Y|u)=\int h(Y|\nu ,u)g(\nu |u)d\nu$ is an analytically tractable multivariate normal (MN) theory HLM with respect to a high-dimensional random vector ν. We computed the likelihood $h(Y)=\int h(Y|u)g(u)\mathit{\text{du}}$ numerically with respect to a low dimensional u by means of adaptive Gauss-Hermite quadrature (AGHQ). The HLM involved random effects as predictors, reducing bias due to measurement error. The joint model $h(Y^{*}|u)g(u)$ induced by the HLM is guaranteed to be compatible with the HLM. We suggested general rules for selecting the PKREs in a way that minimized the dimension of AGHQ. Although useful for the HLMs considered in this article, they are yet to be extended to other models, for example, for discrete outcomes. We hope that our work will spur research on estimation via PKREs.

The nonlinearities of multiple covariates, multiple outcomes and/or the presence of partially observed discrete variables will increase the dimension of PKREs and, thus, the expense of numerical integration by AGHQ. In that case, integration via multivariate Laplace approximation may contribute to efficient computation (Pinheiro and Bates 1995; Raudenbush, Yang, and Yosef 2000).

Further research may address the problem of highly correlated random effects at the cluster level. One strategy would introduce shared random effects to cope with the “curse of dimensionality” by AGHQ as well as the multicollinearity (Miyazaki and Frank 2006; Sun et al. 2023). In addition, parallel computation of numerical integrals for groups of or single clusters will reduce per-iteration computation time while application of the parameter-extended EM algorithm (Liu, Rubin, and Wu 1998) may reduce the number of iterations to converge.

In related research, Rockwood (2020) estimated a multilevel structural equations model by ML, integrating linear random effects conditional on nonlinear random effects analytically and, then, nonlinear effects numerically by Gaussian quadrature. Our analysis encountered both outcome and predictors quite severely missing. To simulate the analysis closely, we simulated a MAR mechanism due to a known auxiliary predictor. We leave the important extension to a MAR (Grund, Lüdtke, and Robitzsch 2021) or other mechanism due to the fully observed or missing values of the outcome to near future.

It is possible to extend and automate our program that enables a user to specify an analytic HLM and determines PKREs based on a set of rules given the HLM. To that end, we need to develop a more general set of rules, for example, involving discrete covariates MAR.

Often, MI of a binary predictor MAR under multivariate normality is efficiently analyzed (Schafer 1997; Grund, Lüdtke, and Robitzsch 2018). The MI will be, however, incompatible with a HLM having the nonlinear effect of the predictor and, thus, unable to always guarantee unbiased estimation of the HLM. We are currently extending our ML approach via the PKRE idea that will ensure compatibility with and, thus, produce unbiased estimation of a HLM having the nonlinear effects of categorical predictors.

Extension of our approach to MI via Bayesian methods may increase the robustness of findings and is straightforward. In particular, our MN joint model $h(Y^{*},\nu |u;\theta )=h(Y_{\text{mis}},Y,\nu |u;\theta )$ given the PKRE u implies estimation of θ by the Gibbs sampler. The sampler will impute $(Y_{\text{mis}},\nu ,u,\theta )$ from their posteriors compatible with the joint model $h(Y^{*},\nu |u,\theta )g(u|\theta )p(\theta )$ for a reasonably assumed prior $p(\theta )$ by drawing: (i) $Y_{\text{mis}}$ and ν from MN $h(Y_{\text{mis}},\nu |Y,u,\theta )$; (ii) u from nonstandard $g(u|Y^{*},\nu ,\theta )=h(Y^{*},\nu |u,\theta )g(u|\theta )/h(Y^{*},\nu |\theta )$, for example, by importance sampling via Markov chain Monte Carlo integration of $h(Y^{*},\nu |\theta )=E[h(Y^{*},\nu |u,\theta )]$ that samples u from a normal prior $g(u|\theta )$; and (iii) θ from a standard posterior $p(\theta |Y^{*},u,\nu )$ (Schafer and Yucel 2002). A potential virtue of a PKRE is to minimize the dimension of sampling the PKRE from a nonstandard posterior by importance sampling. We find it important to solve the measurement error problem by including level-2 random effects as latent covariates. To that end, we also explain how the Gibbs sampler without consideration of a PKRE (Goldstein et al. 2014; Enders, Du, and Keller 2020) may be modified to be compatible with our scientific model conditional on a latent covariate ν in Appendix E. A valuable future study is to compare the proposed Gibbs sampler estimators with existing estimators of a more sophisticated HLM, for example, involving multiple nonlinear effects or outcomes.

Supplementary Materials

The supplementary files that contain R and C executable codes, data sets and initial values to reproduce Tables 2 and 3 may be downloaded. See README_supplementarymaterials.txt for instruction.

Acknowledgments

We thank two anonymous reviewers and an associate editor for their helpful comments, Craig Enders and Brian Keller for providing Blimp simulation codes, and Dongho Shin for helping Blimp simulation in R environment.

Disclosure Statement

The authors report there are no competing interests to declare.

Additional information

Funding

The research reported here was supported by the Institute of Education Sciences, U.S. Department of Education, through Grant R305D210022. The opinions expressed are those of the authors and do not represent views of the Institute or the U.S. Department of Education.

Appendix A

Problem of Bias in Estimating Model (19)

Let $R_{\mathit{\text{ij}}}^{*}=R_{\mathit{\text{ij}}}$ and $C_{\mathit{\text{ij}}}^{*}=C_{\mathit{\text{ij}}}$ fully observed, δ = 0 to simplify notation, ${\beta }_{\mathit{\text{kj}}}={\gamma }_{k0}+u_{\mathit{\text{kj}}}^{*}\sim N({\gamma }_{k0},{\tau }_{\mathit{\text{kk}}})$ and $\text{cov}(u_{0j}^{*},u_{1j}^{*})={\tau }_{01}$ for k = 0, 1 in model (19). Because ${\overline{C}}_{·j}$ and ${\beta }_j=({\beta }_{0j},{\beta }_{1j},{\nu }_j)$ are independent of $C_{\mathit{\text{ij}}}-{\overline{C}}_{·j}$, and ${ϵ}_{\mathit{\text{ij}}}^{*}|C_{\mathit{\text{ij}}}-{\overline{C}}_{·j}\sim N(C_{\mathit{\text{ij}}}-{\overline{C}}_{·j},{\sigma }_{\mathit{\text{cc}}}/n_j)$, $\left[\begin{array}{c}{R}_{\mathit{\text{ij}}}\\ {\overline{C}}_{·j}\end{array}\right]|{\beta }_j,C_{\mathit{\text{ij}}}-{\overline{C}}_{·j}\sim N(\left[\begin{array}{c}{\beta }_{0j}+{\beta }_{1j}(C_{\mathit{\text{ij}}}-{\overline{C}}_{·j})\\ {\nu }_j\end{array}\right],\left[\begin{array}{cc}{\beta }_{1j}^2{\sigma }_{\mathit{\text{cc}}}/n_j+{\sigma }^2& {\beta }_{1j}{\sigma }_{\mathit{\text{cc}}}/n_j\\ {\beta }_{1j}{\sigma }_{\mathit{\text{cc}}}/n_j& {\sigma }_{\mathit{\text{cc}}}/n_j\end{array}\right])$ implying a mixed model $R_{\mathit{\text{ij}}}|C_{\mathit{\text{ij}}}-{\overline{C}}_{·j}\sim N\left[{\gamma }_{00}+{\gamma }_{10}(C_{\mathit{\text{ij}}}-{\overline{C}}_{·j}),\text{var}(R_{\mathit{\text{ij}}}|C_{\mathit{\text{ij}}}-{\overline{C}}_{·j})\right]$ for $\text{var}(R_{\mathit{\text{ij}}}|C_{\mathit{\text{ij}}}-{\overline{C}}_{·j})={\tau }_{00}+2{\tau }_{01}(C_{\mathit{\text{ij}}}-{\overline{C}}_{·j})+{\tau }_{11}{(C_{\mathit{\text{ij}}}-{\overline{C}}_{·j})}^2+({\tau }_{11}+{\gamma }_{10}^2){\sigma }_{\mathit{\text{cc}}}/n_j+{\sigma }^2$, $\text{cov}(R_{\mathit{\text{ij}}},{\overline{C}}_{·j}|C_{\mathit{\text{ij}}}-{\overline{C}}_{·j})={\gamma }_{01}{\tau }_{\nu \nu }+{\gamma }_{11}{\tau }_{\nu \nu }(C_{\mathit{\text{ij}}}-{\overline{C}}_{·j})+{\gamma }_{10}{\sigma }_{\mathit{\text{cc}}}/n_j$. Let ${\lambda }_j={\tau }_{\nu \nu }/({\tau }_{\nu \nu }+{\sigma }_{\mathit{\text{cc}}}/n_j)$ be the reliability of ${\overline{C}}_{·j}$ as an error-prone measure of ν_j (Raudenbush and Bryk 2002). The implied $R_{\mathit{\text{ij}}}|C_{\mathit{\text{ij}}}-{\overline{C}}_{·j},{\overline{C}}_{·j}\sim N({\mu }_{\mathit{\text{ij}}},V_{\mathit{\text{ij}}})$ has (A1) $$\begin{array}{c}{\mu }_{\mathit{\text{ij}}}={\gamma }_{00}+\left[{\gamma }_{01}-(1-{\lambda }_j)({\gamma }_{01}-{\gamma }_{10})\right]{\overline{C}}_{·j}+{\gamma }_{10}(C_{\mathit{\text{ij}}}-{\overline{C}}_{·j})\\ +{\lambda }_j{\gamma }_{11}{\overline{C}}_{·j}(C_{\mathit{\text{ij}}}-{\overline{C}}_{·j})\\ V_{\mathit{\text{ij}}}={\sigma }^2+\left[{\tau }_{00|\nu }+(1-{\lambda }_j){({\gamma }_{01}-{\gamma }_{10})}^2{\tau }_{\nu \nu }+({\tau }_{11|\nu }+{\gamma }_{11}^2{\tau }_{\nu \nu }){\sigma }_{\mathit{\text{cc}}}/n_j\right]\\ +2\left[{\tau }_{01|\nu }+(1-{\lambda }_j)({\gamma }_{01}-{\gamma }_{10}){\gamma }_{11}{\tau }_{\nu \nu }\right](C_{\mathit{\text{ij}}}-{\overline{C}}_{·j})\\ +\left[{\tau }_{11|\nu }+(1-{\lambda }_j){\gamma }_{11}^2{\tau }_{\nu \nu }\right]{(C_{\mathit{\text{ij}}}-{\overline{C}}_{·j})}^2.\end{array}$$(A1)

The bias terms are complicated functions of cluster sizes n_j and parameters, but revealing in the balanced case of n_j = n where ${\lambda }_j=\lambda$. The interaction effect $\lambda {\gamma }_{11}$ of ${\overline{C}}_{·j}(C_{\mathit{\text{ij}}}-{\overline{C}}_{·j})$ has a downward bias term $-(1-\lambda ){\gamma }_{11}$ that introduces bias $(1-\lambda )({\gamma }_{01}-{\gamma }_{10}){\gamma }_{11}{\tau }_{\nu \nu }$ and $(1-\lambda ){\gamma }_{11}^2{\tau }_{\nu \nu }$ in estimation of ${\tau }_{01|\nu }$ and ${\tau }_{11|\nu }$, respectively. Likewise, the main effect of ${\overline{C}}_{·j}$ has a bias term $-(1-\lambda )({\gamma }_{01}-{\gamma }_{10})$ which propagates bias $(1-{\lambda }_j){({\gamma }_{01}-{\gamma }_{10})}^2{\tau }_{\nu \nu }$ and $(1-\lambda )({\gamma }_{01}-{\gamma }_{10}){\gamma }_{11}{\tau }_{\nu \nu }$ in estimation of ${\tau }_{00|\nu }$ and ${\tau }_{01|\nu }$, respectively. Estimation of ${\tau }_{00|\nu }$ results in an additional upward bias term $({\tau }_{11|\nu }+{\gamma }_{11}^2{\tau }_{\nu \nu }){\sigma }_{\mathit{\text{cc}}}/n$ from the error-prone measure ${\overline{C}}_{·j}$ of ν_j. Consequently, this approach results in biased estimation of $({\gamma }_{01},{\gamma }_{11},{\tau }_{00|\nu },{\tau }_{01|\nu },{\tau }_{11|\nu })$. In particular, the estimate of γ₁₁ is biased downward, but those of ${\tau }_{00|\nu }$ and ${\tau }_{11|\nu }$ upward.

Two special cases are of interest. When ${\gamma }_{11}=0,\text{cov}({\beta }_{0j},{\beta }_{1j}|{\overline{C}}_{·j})={\tau }_{01|\nu }$ and $\text{var}({\beta }_{1j}|{\overline{C}}_{·j})={\tau }_{11|\nu }$ become unbiased, and the estimator of ${\tau }_{00|\nu }$ becomes less biased. As cluster sizes $n_j\to \infty ,{\lambda }_j\to 1,{\overline{C}}_{·j}\to {\nu }_j$ by the laws of large numbers, and all bias terms tend to zero.

Appendix B

The E Step for estimation of model (14)

Let ${\mathit{A}}_{1\mathit{\text{ij}}}^{*}={[C_{\mathit{\text{ij}}}^{*}{A}_{1\mathit{\text{ij}}}^{*T}]}^T,Y_{\mathit{\text{ij}}}^{*}={[R_{\mathit{\text{ij}}}^{*}{\mathit{A}}_{1\mathit{\text{ij}}}^{*T}]}^T p_1$-by-1, $A_{2j}^{*} p_2$-by-1, and ${\mathit{\nu }}_{1j}={[{\nu }_j{\nu }_{1j}^T]}^T$. A reasonably general CD model given known covariates $X_{1\mathit{\text{ij}}}$ at level 1 and $X_{2j}$ at level 2 is (B1) $$\begin{array}{c}{R}_{\mathit{\text{ij}}}^{*}={\gamma }_{00}^T{X}_{2j}+u_{0j}^{*}+B_{1j}^T({\mathit{A}}_{1\mathit{\text{ij}}}^{*}-{\Delta }_{00}{X}_{2j}-{\mathit{\nu }}_{1j})+{\gamma }_{30}^T{X}_{1\mathit{\text{ij}}}+e_{\mathit{\text{ij}}}^{*}\\ {\mathit{A}}_{1\mathit{\text{ij}}}^{*}={\Delta }_{00}{X}_{2j}+{\Delta }_{10}{X}_{1\mathit{\text{ij}}}+{\mathit{\nu }}_{1j}+{\mathit{ϵ}}_{\mathit{\text{ij}}}^{*},A_{2j}^{*}={\Delta }_2{X}_{2j}+{\nu }_{2j}\end{array}$$(B1) for $B_{1j}^T=[{\gamma }_{10}^T{X}_{2j}+u_{1j}^{*}{\gamma }_{20}^T],{\Delta }_{00}={[{\delta }_{001}{\Delta }_{002}^T]}^T$ and ${\Delta }_{10}={[{\delta }_{101}{\Delta }_{102}^T]}^T$. Denote ${\nu }_{0j}={[u_{0j}^{*}{\mathit{\nu }}_{1j}^T]}^T$ to separate all other random effects ${\nu }_j^{*}={[{\nu }_{0j}^T{\nu }_{2j}^T]}^T$ from $u_{1j}^{*}$ at level 2, and let $\text{var}({\nu }_j^{*})=\left[\begin{array}{cc}{T}_{00}& T_{02}\\ T_{20}& T_{22}\end{array}\right],\text{cov}({\nu }_j^{*},u_{1j}^{*})=\left[\begin{array}{c}{T}_{01}\\ T_{21}\end{array}\right]$ and $\text{var}({\nu }_j^{*}|u_{1j}^{*})=\left[\begin{array}{cc}{T}_{00|1}& T_{02|1}\\ T_{20|1}& T_{22|1}\end{array}\right]$. Define matrix $O_{2j}$ that selects observed values $A_{2j}=O_{2j}{A}_{2j}^{*}$ in $A_{2j}^{*}$ such that $\text{var}(A_{2j}|u_{1j}^{*})=O_{2j}{T}_{22|1}{O}_{2j}^T=T_{22|1j}$, and $\text{cov}({\nu }_{\mathit{\text{aj}}},{\nu }_{\mathit{\text{bj}}}|u_{1j}^{*},A_{2j})={\Omega }_{\mathit{\text{abj}}}=T_{\mathit{\text{ab}}|1}-T_{a2|1}{O}_{2j}^T{T}_{22|1j}^{-1}{O}_{2j}{T}_{2b|1}$ for $a,b=0,2$. The likelihood $L(\theta )=\prod _j\int h(Y_j|u_{1j}^{*})\varphi (u_{1j}^{*}|0,{\tau }_{11})d{u}_{1j}^{*}$ has a key component (B2) $$\begin{array}{c}h(Y_j|u_{1j}^{*})\propto {\left(\left|{\Omega }_{00j}{|}^{-1}\right|{\Delta }_j{|}^{-1}\left|T_{22|1j}{|}^{-1}\prod _i\right|{\Sigma }_{\mathit{\text{ij}}}{|}^{-1}\right)}^{1/2}\hspace{0.17em}\mathrm{\text{exp}}\hspace{0.17em}\{-\frac{1}{2}\\ [\sum _i{e}_{o1\mathit{\text{ij}}}^T{\Sigma }_{\mathit{\text{ij}}}^{-1}{e}_{o1\mathit{\text{ij}}}-\sum _i{e}_{o1\mathit{\text{ij}}}^T{\Sigma }_{\mathit{\text{ij}}}^{-1}{O}_{\mathit{\text{ij}}}{\Delta }_j^{-1}(\sum _i{O}_{\mathit{\text{ij}}}^T{\Sigma }_{\mathit{\text{ij}}}^{-1}{e}_{o1\mathit{\text{ij}}}\\ +2{\Omega }_{00j}^{-1}{T}_{02|1}{O}_{2j}^T{T}_{22|1j}^{-1}{e}_{o2j})+e_{o2j}^T(T_{22|1j}^{-1}{O}_{2j}{T}_{20|1}\\ ({\Omega }_{00j}^{-1}-{\Omega }_{00j}^{-1}{\Delta }_j^{-1}{\Omega }_{00j}^{-1})T_{02|1}{O}_{2j}^T{T}_{22|1j}^{-1}+T_{22|1j}^{-1})e_{o2j}]\}\end{array}$$(B2) for $e_{o1\mathit{\text{ij}}}=O_{\mathit{\text{ij}}}(d_{\mathit{\text{ij}}}^{*}-T_{01}{\tau }_{11}^{-1}{u}_{1j}^{*}),e_{o2j}=O_{2j}(A_{2j}^{*}-{\Delta }_2{X}_{2j}-T_{21}{\tau }_{11}^{-1}{u}_{1j}^{*})$ and ${\Delta }_j={\sum }_{i=1}^{n_j}{A}_{\mathit{\text{OOij}}}+{\Omega }_{00j}^{-1}$ where $A_{\mathit{\text{OOij}}}=O_{\mathit{\text{ij}}}^T{\Sigma }_{\mathit{\text{ij}}}^{-1}{O}_{\mathit{\text{ij}}}$ and $d_{\mathit{\text{ij}}}^{*}=\left[\begin{array}{c}{R}_{\mathit{\text{ij}}}^{*}\\ {\mathit{A}}_{1\mathit{\text{ij}}}\end{array}\right]-\left[\begin{array}{c}{\gamma }_{00}^T\\ {\Delta }_{00}\end{array}\right]X_{2j}-\left[\begin{array}{cc}1& B_{1j}^T\\ 0& I_{p_1-1}\end{array}\right]\left[\begin{array}{c}{\gamma }_{30}^T\\ {\Delta }_{10}\end{array}\right]X_{1\mathit{\text{ij}}}$.

Define $\mathcal{E}(A)=E(A|u_{1j}^{*},Y_j),\mathcal{V}(A)=\text{var}(A|u_{1j}^{*},Y_j)$ and $\mathcal{C}(A,B)=\text{cov}(A,B|u_{1j}^{*},Y_j)$. We have multivariate normal $f({\nu }_{0j},{\nu }_{2j},r_{\mathit{\text{ij}}}^{*}|u_{1j}^{*},Y_j)$ for $$\begin{array}{c}\mathcal{E}({\nu }_{0j})=T_{01}{\tau }_{11}^{-1}{u}_{1j}^{*}+{\Delta }_j^{-1}(\sum _i{O}_{\mathit{\text{ij}}}^T{\Sigma }_{\mathit{\text{ij}}}^{-1}{e}_{o1\mathit{\text{ij}}}\\ +{\Omega }_{00j}^{-1}{T}_{02|1}{O}_{2j}^T{T}_{22|1j}^{-1}{e}_{o2j}),\\ \mathcal{E}({\nu }_{2j})=T_{21}{\tau }_{11}^{-1}{u}_{1j}^{*}+T_{22|1}{O}_{2j}^T{T}_{22|1j}^{-1}{e}_{o2j}\\ +{\Omega }_{20j}{\Omega }_{00j}^{-1}\left[\mathcal{E}({\nu }_{0j})-T_{01}{\tau }_{11}^{-1}{u}_{1j}^{*}-T_{02|1}{O}_{2j}^T{T}_{22|1j}^{-1}{e}_{o2j}\right],\\ \mathcal{E}(r_{\mathit{\text{ij}}}^{*})={\Sigma }_j^{*}{A}_{\mathit{\text{OOij}}}[d_{\mathit{\text{ij}}}^{*}-\mathcal{E}({\nu }_{0j})],\\ \mathcal{V}(r_{\mathit{\text{ij}}}^{*})={\Sigma }_j^{*}-{\Sigma }_j^{*}(A_{\mathit{\text{OOij}}}-A_{\mathit{\text{OOij}}}{\Delta }_j^{-1}{A}_{\mathit{\text{OOij}}}){\Sigma }_j^{*},\\ \mathcal{V}({\nu }_{0j})={\Delta }_j^{-1},\\ \mathcal{C}({\nu }_{0j},{\nu }_{2j})={\Delta }_j^{-1}{\Omega }_{00j}^{-1}{\Omega }_{02j},\\ \mathcal{C}({\nu }_{0j},r_{\mathit{\text{ij}}}^{*})=-{\Delta }_j^{-1}{A}_{\mathit{\text{OOij}}}{\Sigma }_j^{*},\\ \mathcal{V}({\nu }_{2j})={\Omega }_{22j}-{\Omega }_{20j}({\Omega }_{00j}^{-1}-{\Omega }_{00j}^{-1}{\Delta }_j^{-1}{\Omega }_{00j}^{-1}){\Omega }_{02j},\\ \mathcal{C}({\nu }_{2j},r_{\mathit{\text{ij}}}^{*})={\Omega }_{20j}{\Omega }_{00j}^{-1}\mathcal{C}({\nu }_{0j},r_{\mathit{\text{ij}}}^{*}).\end{array}$$

Let ${\tilde{A}}_{1\mathit{\text{ij}}}^{*}=A_{1\mathit{\text{ij}}}^{*}-{\Delta }_{002}{X}_{2j}-{\nu }_{1j},{\delta }_{10}=\mathit{\text{vec}}({\Delta }_{10}^T)$ and ${\beta }_j=(I_{p_1+1+p_2}\otimes X_{2j}^T){\gamma }_{\beta }+u_j^{*}$ for ${\gamma }_{\beta }=\mathit{\text{vec}}[{\gamma }_{00}{\gamma }_{10}{\Delta }_{00}^T{\Delta }_2^T]$ and $u_j^{*}={[u_{0j}^{*}{u}_{1j}^{*}{\mathit{\nu }}_j^T]}^T$. The expected CD MLEs are $$\begin{array}{c}{\widehat{\gamma }}_{20}={\gamma }_{20}+{\left(\sum _{j}E\left[\sum _i\mathcal{E}({\tilde{A}}_{1\mathit{\text{ij}}}^{*}{\tilde{A}}_{1\mathit{\text{ij}}}^{*T})|Y_j\right]\right)}^{-1}\sum _{j}E\left[\sum _i\mathcal{E}(e_{\mathit{\text{ij}}}^{*}{\tilde{A}}_{1\mathit{\text{ij}}}^{*})|Y_j\right],\\ {\widehat{\gamma }}_{30}={\gamma }_{30}+{\left(\sum _j\sum _i{X}_{1\mathit{\text{ij}}}{X}_{1\mathit{\text{ij}}}^T\right)}^{-1}\sum _j\sum _i{X}_{1\mathit{\text{ij}}}E\left[\mathcal{E}(e_{\mathit{\text{ij}}}^{*})|Y_j\right],\\ {\widehat{\delta }}_{10}={\delta }_{10}+\mathit{\text{vec}}\left[\left(\sum _j\sum _{i}E\left[\mathcal{E}({\mathit{ϵ}}_{\mathit{\text{ij}}}^{*})|Y_j\right]X_{1\mathit{\text{ij}}}^T\right){\left(\sum _j\sum _i{X}_{1\mathit{\text{ij}}}{X}_{1\mathit{\text{ij}}}^T\right)}^{-1}\right]\\ {\widehat{\gamma }}_{\beta }={\gamma }_{\beta }+\mathit{\text{vec}}\left[\left(\sum _{j}E[\mathcal{E}(u_j^{*})|Y_j]X_{2j}^T\right){\left(\sum _j{X}_{2j}{X}_{2j}^T\right)}^{-1}\right]\\ {\widehat{\sigma }}^2=\sum _{j}E\left[\sum _i\mathcal{E}(e_{\mathit{\text{ij}}}^{*2})|Y_j\right]/N,{\widehat{\Sigma }}_{ϵ}=\sum _{j}E\left[\sum _i\mathcal{E}({\mathit{ϵ}}_{\mathit{\text{ij}}}^{*}{\mathit{ϵ}}_{\mathit{\text{ij}}}^{*T})|Y_j\right]/N,\\ \widehat{\tau }=\sum _{j}E[\mathcal{E}(u_j^{*}{u}_j^{*T})|Y_j]/J\end{array}$$given θ where $\mathcal{E}(e_{\mathit{\text{ij}}}^{*})=B_{1j}^{*T}\mathcal{E}(r_{\mathit{\text{ij}}}^{*})$ and $\mathcal{E}(e_{\mathit{\text{ij}}}^{*2})=B_{1j}^{*T}\mathcal{E}(r_{\mathit{\text{ij}}}^{*}{r}_{\mathit{\text{ij}}}^{*T})B_{1j}^{*}$ for $B_{1j}^{*T}=[1-B_{1j}^T]$.

Appendix C

Numerical Integration by AGHQ

Let $f(u_{1j}^{*})=h(Y_j|u_{1j}^{*})\varphi (u_{1j}^{*};0,{\tau }_{11})$ be a function of $u_{1j}^{*}$. Given ${\tilde{u}}_{1j}^{*}=E(u_{1j}^{*}|Y_j),V_{u1j}=\text{var}(u_{1j}^{*}|Y_j)=L_{u1j}^2/2$, Q-point weights $(w_1,\dots ,w_Q)$ and abscissas $(a_1,\dots ,a_Q)$, (C1) $$\begin{array}{c}h(Y_j)=\int \frac{\varphi (u_{1j}^{*};{\tilde{u}}_{1j}^{*},V_{u1j})}{\varphi (u_{1j}^{*};{\tilde{u}}_{1j}^{*},V_{u1j})}f(u_{1j}^{*})d{u}_{1j}^{*}\\ \approx L_{u1j}\sum _{k=1}^Q{w}_k{e}^{a_k^2}f(z_{\mathit{\text{kj}}}),\end{array}$$(C1) for $z_{\mathit{\text{kj}}}=L_{u1j}{a}_k+{\tilde{u}}_{1j}^{*}$. The $g(u_{1j}^{*}|Y_j)=f(u_{1j}^{*})/h(Y_j)$ is approximately $\varphi (u_{1j}^{*};{\tilde{u}}_{1j}^{*},V_{u1j})$ for large cluster sizes n_j by the Bayesian central limit theorem such that $f(u_{1j}^{*})\propto \varphi (u_{1j}^{*};{\tilde{u}}_{1j}^{*},V_{u1j})$ produces well approximated $h(Y_j)$ by a low degree polynomial. The approximation is exact if $f(u_{1j}^{*})$ is a $2Q-1$ degree polynomial in $u_{1j}^{*}$ (Pinheiro and Bates 1995; Rabe-Hesketh, Skrondal, and Pickles 2002; Carlin and Louis 2009). Likewise, for $E(S_{\text{CD}j}|u_{1j}^{*},Y_j)$ closed-form, (C2) $$\begin{array}{c}E\left(S_{\text{CD}j}|Y_j\right)=\int E(S_{\text{CD}j}|u_{1j}^{*},Y_j)g(u_{1j}^{*}|Y_j)d{u}_{1j}^{*}\\ \approx \frac{L_{u1j}}{h(Y_j)}\sum _{k=1}^{Q}E(S_{\text{CD}j}|z_{\mathit{\text{kj}}},Y_j)w_k{e}^{a_k^2}f(z_{\mathit{\text{kj}}})\end{array}$$(C2)

Let $Y_j^{*}=(Y_j,Y_{\text{mis}j})$, and ${\varphi }_{\tau }$ and ${\varphi }_{\Sigma }$ be vectors of distinct elements of τ and ${\Sigma }_{ϵ}$, respectively. The log-likelihood $l=\sum _j{l}_j$ and score $S=\sum _j{S}_j$ have summands $$\begin{array}{c}{l}_j=\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}h(Y_j)=\hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}\int g_{j}d{Y}_{\text{mis}j}d{\nu }_{0j}d{u}_{1j}^{*},\\ S_j=\frac{\partial l_j}{\partial \theta }=E\left[\mathcal{E}\left(\frac{\partial \hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}{g}_j}{\partial \theta }\right)|Y_j\right].\end{array}$$

for $g_j=\prod _{i}f(R_{\mathit{\text{ij}}}^{*}|{\mathbf{A}}_{1\mathit{\text{ij}}}^{*},u_j^{*};{\gamma }_{20},{\gamma }_{30},{\sigma }^2)f({\mathbf{A}}_{1\mathit{\text{ij}}}^{*}|u_j^{*};{\delta }_{10},{\varphi }_{\Sigma })\varphi (u_j^{*};{\gamma }_{\beta },\tau )$ from (B1). Let $E=\frac{\partial \mathit{\text{vec}}\tau }{\partial {\varphi }_{\tau }^T}$ and $F=\frac{\partial \mathit{\text{vec}}{\Sigma }_{ϵ}}{\partial {\varphi }_{\Sigma }^T}$. The $\mathcal{E}\left(\frac{\partial \hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}{g}_j}{\partial \theta }\right)$ stacks $$\begin{array}{c}\mathcal{E}\left(\frac{\partial \hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}{g}_j}{\partial {\gamma }_{20}}\right)={\sigma }^{-2}\sum _i\mathcal{E}(e_{\mathit{\text{ij}}}^{*}{\tilde{A}}_{1\mathit{\text{ij}}}^{*}),\\ \mathcal{E}\left(\frac{\partial \hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}{g}_j}{\partial {\gamma }_{30}}\right)={\sigma }^{-2}\sum _i\mathcal{E}(e_{\mathit{\text{ij}}}^{*})X_{1\mathit{\text{ij}}},\\ \mathcal{E}\left(\frac{\partial \hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}{g}_j}{\partial {\delta }_{10}}\right)=\mathit{\text{vec}}\left(X_{1\mathit{\text{ij}}}\sum _i\mathcal{E}({\mathit{ϵ}}_{\mathit{\text{ij}}}^{*T}){\Sigma }_{ϵ}^{-1}\right),\\ \mathcal{E}\left(\frac{\partial \hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}{g}_j}{\partial {\sigma }^2}\right)=\frac{1}{2}{A}_{\sigma j},\\ \mathcal{E}\left(\frac{\partial \hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}{g}_j}{\partial {\varphi }_{\Sigma }}\right)=\frac{1}{2}{F}^T{A}_{\Sigma j},\\ \mathcal{E}\left(\frac{\partial \hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}{g}_j}{\partial {\gamma }_{\beta }}\right)=\mathit{\text{vec}}\left(X_{2j}\mathcal{E}(u_j^{*T}){\tau }^{-1}\right),\\ \mathcal{E}\left(\frac{\partial \hspace{0.17em}\mathrm{\text{ln}}\hspace{0.17em}{g}_j}{\partial {\varphi }_{\tau }}\right)=\frac{1}{2}{E}^T\mathit{\text{vec}}[{\tau }^{-1}\mathcal{E}(u_j^{*}{u}_j^{*T}){\tau }^{-1}-{\tau }^{-1}]\end{array}$$for $A_{\sigma j}={\sigma }^{-4}\sum _i\mathcal{E}(e_{\mathit{\text{ij}}}^{*2})-n_j{\sigma }^{-2}$ and $A_{\Sigma j}=\mathit{\text{vec}}({\Sigma }_{ϵ}^{-1}\sum _i\mathcal{E}({\mathit{ϵ}}_{\mathit{\text{ij}}}^{*}{\mathit{ϵ}}_{\mathit{\text{ij}}}^{*T}){\Sigma }_{ϵ}^{-1}-n_j{\Sigma }_{ϵ}^{-1})$. We compute S_j also by AGHQ for $\text{var}(\widehat{\theta })\approx {(\sum _j{S}_j{S}_j^T)}^{-1}$ (Hedeker and Gibbons 1994; Raudenbush, Yang, and Yosef 2000; Olsen and Schafer 2001). Section 7 shows good approximation for the sample sizes analyzed in this article.

Appendix D

Translating model (B1) to ${\widehat{\theta }}_{(19)}^{*}$

Define ${\beta }_{\mathit{\text{kj}}}={\gamma }_{k0}^T{X}_{2j}+u_{\mathit{\text{kj}}}^{*},{\beta }_{\mathit{\text{Cj}}}={\delta }_{001}^T{X}_{2j}+{\nu }_j,\text{cov}(u_{\mathit{\text{kj}}}^{*},u_{k\prime j}^{*}|X_{2j})={\tau }_{\mathit{\text{kk}}\prime },\text{cov}(u_{\mathit{\text{kj}}}^{*},{\nu }_j|X_{2j})={\tau }_{k\nu }$ and $\text{var}({\nu }_j|X_{2j})={\tau }_{\nu \nu }$ for $k,k\prime =0,1$. With $X_{2j}$ marginalized out, ${\beta }_j={[{\beta }_{0j}{\beta }_{1j}{\beta }_{\mathit{\text{Cj}}}]}^T\sim N(\left[\begin{array}{c}{\gamma }_{00}^T\\ {\gamma }_{10}^T\\ {\delta }_{121}^T\end{array}\right]E{X}_2,\left[\begin{array}{ccc}{t}_{00}& t_{01}& t_{0\nu }\\ & t_{11}& t_{1\nu }\\ & & t_{\nu \nu }\end{array}\right]=\left[\begin{array}{ccc}{\tau }_{00}& {\tau }_{01}& {\tau }_{0\nu }\\ & {\tau }_{11}& {\tau }_{1\nu }\\ & & {\tau }_{\nu \nu }\end{array}\right] +\left[\begin{array}{c}{\gamma }_{00}^T\\ {\gamma }_{10}^T\\ {\delta }_{001}^T\end{array}\right]V{X}_2[{\gamma }_{00}{\gamma }_{10}{\delta }_{001}])$ for $E(X_{2j})=E{X}_2$ and $\text{var}(X_{2j})=V{X}_2$. Let ${\tilde{\beta }}_{\mathit{\text{Cj}}}={\beta }_{\mathit{\text{Cj}}}-{\delta }_{001}^{T}E{X}_2\sim N(0,t_{\nu \nu })$ to find (D1) $${\beta }_{\mathit{\text{kj}}}|{\tilde{\beta }}_{\mathit{\text{Cj}}}\sim \left({\gamma }_{k0}^{*}+{\gamma }_{k1}^{*}{\tilde{\beta }}_{\mathit{\text{Cj}}},{\tau }_{\mathit{\text{kk}}}^{*}\right)$$(D1) for ${\gamma }_{k0}^{*}={\gamma }_{k0}^{T}E{X}_2,{\gamma }_{k1}^{*}=t_{k\nu }/t_{\nu \nu }$ and $\text{cov}({\beta }_{\mathit{\text{kj}}},{\beta }_{k\prime j}|{\tilde{\beta }}_{\mathit{\text{Cj}}})={\tau }_{\mathit{\text{kk}}\prime }^{*}=t_{\mathit{\text{kk}}\prime }-{\gamma }_{k1}^{*}{t}_{\nu \nu }{\gamma }_{k\prime 11}^{*}$.

Within cluster j given $u_j^{*}$, denote ${\tilde{\mathit{A}}}_{1\mathit{\text{ij}}}^{*}={\mathit{A}}_{1\mathit{\text{ij}}}^{*}-{\Delta }_{00}{X}_{2j}-{\mathit{\nu }}_{1j}={\Delta }_{10}{X}_{1\mathit{\text{ij}}}+{\mathit{ϵ}}_{\mathit{\text{ij}}}$. Marginalizing $X_{1\mathit{\text{ij}}}$ out using $E{X}_1=E(X_{1\mathit{\text{ij}}})$ and $V{X}_1=\text{var}(X_{1\mathit{\text{ij}}})$, we find $f(Y_{\mathit{\text{ij}}}^{*}|u_j^{*})$ in $$\tilde{{\mathit{A}}_{1\mathit{\text{ij}}}^{*}}={\Delta }_{10}E{X}_1+{\widetilde{\mathit{ϵ}}}_{\mathit{\text{ij}}}^{*},R_{\mathit{\text{ij}}}={\beta }_{0j}+(B_{1j}^T{\Delta }_{10}+{\gamma }_{30}^T)E{X}_1+B_{1j}^T{\widetilde{\mathit{ϵ}}}_{\mathit{\text{ij}}}^{*}+{\tilde{e}}_{\mathit{\text{ij}}}^{*}$$for $\left[\begin{array}{c}{\tilde{e}}_{\mathit{\text{ij}}}^{*}\\ {\widetilde{\mathit{ϵ}}}_{\mathit{\text{ij}}}^{*}\end{array}\right]\sim N(0,\left[\begin{array}{cc}{\sigma }_{\mathit{\text{ee}}}^{*}& {\Sigma }_{eϵ}^{*}\\ {\Sigma }_{ϵe}^{*}& {\Sigma }_{ϵϵ}^{*}\end{array}\right])$ where ${\sigma }_{\mathit{\text{ee}}}^{*}={\gamma }_{30}^{T}V{X}_1{\gamma }_{30}+{\sigma }^2,{\Sigma }_{ϵe}^{*}=\left[\begin{array}{c}{\sigma }_{\mathit{\text{ce}}}\\ {\Sigma }_{1e}\end{array}\right]=\left[\begin{array}{c}{\delta }_{101}^T\\ {\Delta }_{102}\end{array}\right]V{X}_1{\gamma }_{30}$ and ${\Sigma }_{ϵϵ}^{*}=\left[\begin{array}{cc}{\sigma }_{\mathit{\text{cc}}}^{*}& {\Sigma }_{c1}^{*}\\ {\Sigma }_{1c}^{*}& {\Sigma }_{11}^{*}\end{array}\right]=\left[\begin{array}{cc}{\delta }_{101}^{T}V{X}_1{\delta }_{101}+{\sigma }_{\mathit{\text{cc}}}& {\delta }_{101}^{T}V{X}_1{\Delta }_{102}^T+{\Sigma }_{c1}\\ {\Delta }_{102}V{X}_1{\delta }_{101}+{\Sigma }_{1c}& {\Delta }_{102}V{X}_i{\Delta }_{102}^T+{\Sigma }_{11}\end{array}\right]$. Consequently, $R_{\mathit{\text{ij}}}|{\beta }_j,{ϵ}_{\mathit{\text{ij}}}\sim N[E(R_{\mathit{\text{ij}}}^{*}|{\beta }_j,{ϵ}_{\mathit{\text{ij}}}),{\sigma }^{\prime 2}]$ for (D2) $$\begin{array}{c}E(R_{\mathit{\text{ij}}}^{*}|{\beta }_j,{ϵ}_{\mathit{\text{ij}}})={\beta }_{0j}+(B_{1j}^T{\Delta }_{10}+{\gamma }_{30}^T)E{X}_1\\ +({\beta }_{1j}+{\gamma }_{20}^T{\Sigma }_{1c}^{*}/{\sigma }_{\mathit{\text{cc}}}^{*}+{\sigma }_{\mathit{\text{ec}}}^{*}/{\sigma }_{\mathit{\text{cc}}}^{*}){ϵ}_{\mathit{\text{ij}}},\\ {\sigma }^{\prime 2}={\gamma }_{20}^T({\Sigma }_{11}^{*}-{\Sigma }_{1c}^{*}{\Sigma }_{c1}^{*}/{\sigma }_{\mathit{\text{cc}}}^{*}){\gamma }_{20}+2{\gamma }_{20}^T({\Sigma }_{1e}^{*}-{\Sigma }_{1c}^{*}{\sigma }_{\mathit{\text{ce}}}^{*}/{\sigma }_{\mathit{\text{cc}}}^{*})\\ +({\sigma }_{\mathit{\text{ee}}}^{*}-{\sigma }_{\mathit{\text{ce}}}^{*2}/{\sigma }_{\mathit{\text{cc}}}^{*}).\end{array}$$(D2)

As explained by Section 4.1, Equations (D1)–(D2) result in a mixed model (19) for (D3) $$\begin{array}{c}{\gamma }_{00}={\gamma }_{00}^{*}+({\gamma }_{10}^{*}{\delta }_{101}^T+{\gamma }_{20}^T{\Delta }_{102}+{\gamma }_{30}^T)E{X}_1,\\ {\gamma }_{01}={\gamma }_{01}^{*}+{\gamma }_{11}^{*}{\delta }_{101}^{T}E{X}_1,\\ {\gamma }_{10}={\gamma }_{10}^{*}+{\gamma }_{20}^T{\Sigma }_{1c}^{*}/{\sigma }_{\mathit{\text{cc}}}^{*}+{\sigma }_{\mathit{\text{ec}}}^{*}/{\sigma }_{\mathit{\text{cc}}}^{*},\\ {\gamma }_{11}={\gamma }_{11}^{*},\\ {\sigma }^2={\sigma }^{\prime 2}\\ {\tau }_{00|\nu }={\tau }_{00}^{*}+2{\tau }_{01}^{*}{\delta }_{101}^{T}E{X}_1+{\tau }_{11}^{*}{({\delta }_{101}^{T}E{X}_1)}^2,\\ {\tau }_{01|\nu }={\tau }_{01}^{*}+{\tau }_{11}^{*}{\delta }_{101}^{T}E{X}_1,\\ {\tau }_{11|\nu }={\tau }_{11}^{*}.\end{array}$$(D3)

Estimating $(E{X}_1,V{X}_1,E{X}_2,V{X}_2)$ from sample, we find ${\widehat{\theta }}_{(19)}^{*}$ above and compute $\text{var}({\widehat{\theta }}_{(19)}^{*})$ by the delta method. Translations (D3) simplify if $E{X}_1=0$ (e.g., $E{X}_1=E(X_{1\mathit{\text{ij}}}-{\overline{X}}_{1j})$). In Section 6, we used the translations for $X_{2j}=1$ and $X_{1\mathit{\text{ij}}}=0$.

Appendix E

Compatible Gibbs Sampler without a PKRE

Without the merit of a PKRE, a Bayesian joint distribution based on (19) is $f(R_{\mathit{\text{ij}}}^{*}|C_{\mathit{\text{ij}}}^{*},{\nu }_j,u_{0j},u_{1j},\theta )f(u_{0j},u_{1j}|\theta )f(C_{\mathit{\text{ij}}}^{*}|{\nu }_j,\theta )\varphi ({\nu }_j;0,{\tau }_{\nu \nu })p(\theta )$ for a prior $p(\theta )$ and $\theta ={\theta }_{(19)}^{*}$. Our scientific model is an analytic integral $f_1(R_{\mathit{\text{ij}}}^{*}|C_{\mathit{\text{ij}}}^{*},{\nu }_j,\theta )=\int \int f(R_{\mathit{\text{ij}}}^{*}|C_{\mathit{\text{ij}}}^{*},{\nu }_j,u_{0j},u_{1j},\theta )f(u_{0j},u_{1j}|\theta )d{u}_{0j}d{u}_{1j}$. The Gibbs sampler of Enders, Du, and Keller (2020) may be modified to be compatible with our scientific model conditional on a latent covariate ν_j by sampling (i) ν_j from a compatible posterior $$p({\nu }_j|·)=\frac{\prod _{i}f(R_{\mathit{\text{ij}}}^{*}|C_{\mathit{\text{ij}}}^{*},{\nu }_j,u_{0j},u_{1j},\theta )f(C_{\mathit{\text{ij}}}^{*}|{\nu }_j,\theta )\varphi ({\nu }_j;0,{\tau }_{\nu \nu })}{\int \prod _{i}f(R_{\mathit{\text{ij}}}^{*}|C_{\mathit{\text{ij}}}^{*},{\nu }_j,u_{0j},u_{1j},\theta )f(C_{\mathit{\text{ij}}}^{*}|{\nu }_j,\theta )\varphi ({\nu }_j;0,{\tau }_{\nu \nu })d{\nu }_j}$$for the denominator approximated by the MCMC integration, and (ii) a missing $C_{\mathit{\text{ij}}}^{*}$ from a compatible normal posterior $p(C_{\mathit{\text{ij}}}^{*}|·)\propto f(R_{\mathit{\text{ij}}}^{*}|C_{\mathit{\text{ij}}}^{*},{\nu }_j,u_{0j},u_{1j},\theta )f(C_{\mathit{\text{ij}}}^{*}|{\nu }_j,\theta )$ with $$E(C_{\mathit{\text{ij}}}^{*}|·)=\delta +{\nu }_j+\frac{{\beta }_{1j}{\sigma }_{\mathit{\text{cc}}}}{{\beta }_{1j}^2{\sigma }_{\mathit{\text{cc}}}+{\sigma }^2}(R_{\mathit{\text{ij}}}^{*}-{\beta }_{0j}),\text{var}(C_{\mathit{\text{ij}}}^{*}|·)=\frac{{\sigma }_{\mathit{\text{cc}}}{\sigma }^2}{{\beta }_{1j}^2{\sigma }_{\mathit{\text{cc}}}+{\sigma }^2}$$ by (15) for ${\beta }_{0j}={\gamma }_{00}+{\gamma }_{01}{\nu }_j+u_{0j}$ and ${\beta }_{1j}={\gamma }_{10}+{\gamma }_{11}{\nu }_j+u_{1j}$.