Linear mixed models with marginally symmetric nonparametric random effects

Abstract

Linear mixed models (LMMs) are used as an important tool in the data analysis of repeated measures and longitudinal studies. The most common form of LMMs utilizes a normal distribution to model the random effects. Such assumptions can often lead to misspecification errors when the random effects are not normal. One approach to remedy the misspecification errors is to utilize a point-mass distribution to model the random effects; this is known as the nonparametric maximum likelihood-fitted (NPML) model. The NPML model is flexible but requires a large number of parameters to characterize the random-effects distribution. It is often natural to assume that the random-effects distribution be at least marginally symmetric. The marginally symmetric NPML (MSNPML) random-effects model is introduced, which assumes a marginally symmetric point-mass distribution for the random effects. Under the symmetry assumption, the MSNPML model utilizes half the number of parameters to characterize the same number of point masses as the NPML model; thus the model confers an advantage in economy and parsimony. An EM-type algorithm is presented for the maximum likelihood (ML) estimation of LMMs with MSNPML random effects; the algorithm is shown to monotonically increase the log-likelihood and is proven to be convergent to a stationary point of the log-likelihood function in the case of convergence. Furthermore, it is shown that the ML estimator is consistent and asymptotically normal under certain conditions, and the estimation of quantities such as the random-effects covariance matrix and individual a posteriori expectations is demonstrated. A simulation study is used to illustrate the gains in efficiency of the MSNPML model over the NPML model under the assumption of symmetry. A pair of real data applications are then used to demonstrate the manner in which the MSNPML model can be used to draw useful statistical inference.

Introduction

Linear mixed models (LMMs) are used as a fundamental tool for the statistical analysis of longitudinal data and data with repeated measurements; see McCulloch and Searle (2001, Ch. 6), Pinheiro and Bates (2000, Ch. 1), and Verbeke and Molenberghs (2000) for introductions on the topic. In the style of Laird and Ware (1982), an LMM can be characterized as follows.

Let $Y_j={(Y_{j1},\dots ,Y_{j{n}_j})}^T$ be a vector of $n_j$ responses belonging to individual $j$, for $j=1,\dots ,n$, where $n$ is the total number of individuals. Further, let each measurement $Y_{jk}$, for $k=1,\dots ,n_j$, be dependent upon covariate vectors $x_{jk}\in {\mathbb{R}}^p$ and $z_{jk}\in {\mathbb{R}}^q$, and for each $j$, let $B_j\in {\mathbb{R}}^q$ be a vector of individual random effects arising from the a priori probability distribution $F_B\left(b\right)$ with density $f_B\left(b\right)$. We say that the data arises from an LMM if for each $j$ and $k$,

$$Y_{jk}|(B_j=b_j)=x_{jk}^T\beta +z_{jk}^T{b}_j+E_{jk}\text{,}$$

where $\beta \in {\mathbb{R}}^p$ is a vector of fixed effects and $E_{jk}$ is a random error with probability density $f_E\left(e\right)$. Here, a superscript $T$ indicates matrix transposition.

The main difficulty that arises in the use of LMMs is the evaluation and manipulation of the marginal density of $y_j$, which has the general form

$$f_{Y_j}\left(y_j\right)={\int }_{{\mathbb{R}}^q}\left[\prod _{k=1}^{n_k}{f}_E(y_{jk}-x_{jk}^T\beta -z_{jk}^T{b}_j)\right]f_B\left(b_j\right)\text{d}{b}_j\text{,}$$

and can often be quite complex due to the integration form.

The traditional approach in dealing with the complexities of (2) is to set the error density as

$$f_E\left(e\right)=\varphi (e;0,{\sigma }^2)\text{,}$$

where

$$\varphi (e;\mu ,{\sigma }^2)=\frac{1}{\sqrt{2\pi {\sigma }^2}}exp[-\frac{{(e-\mu )}^2}{2{\sigma }^2}]$$

is the normal density function with mean $\mu$ and variance ${\sigma }^2$, and to let $f_B\left(b\right)$ be a multivariate normal density. This approach results in $f_{Y_j}\left(y_j\right)$ having the form of a multivariate normal density function, and allows for simple inference by maximum likelihood (ML) estimation via an expectation–maximization (EM) algorithm; see Laird and Ware (1982, Sec. 4) and McLachlan and Krishnan (2008, Sec. 5.9) for details.

It is well known that the estimation of the fixed effects $\beta$ is robust to the specification of the random-effects distribution. However, this robustness does not extend to the characterization of the random effects in the case of misspecification. The robustness as well as the effects of misspecification are explored in Agresti et al. (2004), Butler and Louis (1992), and McCulloch and Neuhaus (2011). For instance, in all three articles, the authors note that the estimates for the fixed effects tended not to be influenced greatly by the choice of the random-effects model. However, Agresti et al. (2004) note that the usual normal model can be highly inefficient when the true random-effects model is polarizing, such as in the case of binary or mixture random-effects distributions. It is further remarked in McCulloch and Neuhaus (2011) that the estimation of random intercept coefficients can be biased by making incorrect assumptions regarding the shape of the underlying random-effects distribution.

Multiple strategies have been considered for remedying random-effects misspecification. For example, Pinheiro et al. (2001) and Song et al. (2007) considered $t$ distributed random-effects and noise models, whereas Arellano-Valle et al. (2005), Lachos et al. (2010), and Ho and Lin (2010) considered the use of skew normal and $t$ distributed random and noise models. Although a rich class, the use of such distributions often do not allow for simplification of the marginal density (2) and do not allow for enough flexibility to model random-effects distributions with multiple modes or deviations from bell-shaped curves.

Based on the nonparametric maximum likelihood (NPML) principle of Laird (1978), Aitkin (1999) and Butler and Louis (1992) suggested using the point-mass density

$$f_B\left(b\right)=\sum _{i=1}^g{\pi }_i\delta (b-{\lambda }_i)$$

for the random effects, where $\delta \left(x\right)$ is the Dirac delta function, $g\ge 1$ is the number of point masses, and ${\lambda }_i\in {\mathbb{R}}^q$ and ${\pi }_k>0$ are the point-mass locations and weights, for $i=1,..,g$, respectively. To ensure that the total probability of the point masses adds up to unity and that the mean of $f_B\left(b\right)$ is zero, we also require the restrictions that ${\sum }_{i=1}^g{\pi }_i=1$ (which implies that ${\pi }_g=1-{\sum }_{i=1}^{g-1}{\pi }_i$) and ${\sum }_{i=1}^g{\pi }_i{\lambda }_i=0$, where $0$ is a zero vector of appropriate size. We shall refer to densities of the form (4) as NPML-fitted (NPML) densities.

In Agresti et al. (2004) it was shown that NPML densities offered improvements in efficiency for estimating the characteristics of the random-effects density such as the covariance and individual a posteriori estimates of the random effects (i.e. $\mathbb{E}\left(B_j\right|Y_j=y_j)$) when compared to the use of parametric random-effects densities in situations, where the true random-effect densities deviated from the assumed parametric form. However, this improvement comes at a cost of modeling complexity, since for any given $g$ and $q$, the number of parameters required for the specification and estimation of the NPML density is $(g-1)(q+1)$. This number can grow quickly if $g$ or $q$ are large.

We note that although unimodality or bell shape cannot be assumed in general, it may still be acceptable to assume symmetry in the random-effects distribution. Thus, we introduce the marginally symmetric NPML (MSMPML) random-effects density

$$f_B\left(b\right)=\sum _{i=1}^g\frac{{\pi }_i}{2}[\delta (b-{\lambda }_i)+\delta (b+{\lambda }_i)]\text{,}$$

where $g$ and ${\pi }_i$ are as specified for (4) and no restrictions are made on ${\lambda }_i$ as density (5) has zero mean by definition. For any $g$, density (5) comprises $2g$ point masses at ${\lambda }_i$ and $-{\lambda }_i$, each with weights ${\pi }_i/2$. Therefore, in situations where marginal symmetry can be assumed, the MSNPML model doubles the number of point masses when compared to the MPML density, for the same number of parameters.

Aside from the aforementioned articles, our work shares commonalities with the mixture random-effects densities of Verbeke and Lesaffre (1996) and, by extension, the mixture of mixed effects models of Celeux et al. (2005) and Ng et al. (2006); in the aforementioned articles, mixtures of normal densities are used to model the random effects, rather than point-mass densities. When $n_j=N$ for all $j$, the MSNPML model corresponds to a symmetry-restricted and homoscedastic version of the repeated-measures linear mixture regression model of Grun and Leisch (2008). Furthermore, the use of parameter constraints to enforce marginal symmetry is similar to those used in Benaglia et al. (2009) and Chauveau and Hunter (2013). A kernel density approach to mixture modeling with symmetric components is presented in Chee and Wang (2013).

In this article, we shall consider the estimation of LMMs with the MSNPML random-effect density under the assumption of normal random error, using ML estimation via a multicycle expectation–conditional maximization (ECM) algorithm (Meng and Rubin, 1993). This algorithm is known to monotonically increase the log-likelihood and to lead to a stationary point of the said function in the case of convergence. Furthermore, we demonstrate how the covariance matrix and a posteriori expectations of the random effects can be estimated under model (5), and also consider various statistical matters such as consistency and the selection of the number of point masses. We then use numerical simulations to study the efficiency gains due to MSMPML over MPML, in the style of Agresti et al. (2004). A pair of example applications chosen from Pinheiro and Bates (2000) are subsequently used to demonstrate the inferential process of using MSMPML random effects.

The article shall proceed as follows. In Section  2, we shall discuss matters regarding ML estimation. In Section  3, statistical matters and computation of statistics are discussed. Results from simulations are then presented in Section  4, followed by example applications in Section  5. Conclusions are then drawn in Section  6.