Bayesian Analysis for Random Effects Models

Abstract

Random effects models have been widely used to analyze correlated data sets, and Bayesian techniques have emerged as a powerful tool to fit the models. However, there has been scarce literature that systematically reviews and summarizes the recent advances of Bayesian analyses of random effects models. This chapter reviews the use of the Dirichlet process mixture (DPM) prior to approximate the distribution of random errors within the general semiparametric random effects models with parametric random effects for longitudinal data setting and failure time setting separately. In a survival setting with clusters, we propose a new class of nonparametric random effects models which is motivated from the accelerated failure models. We employ a beta process prior to tact clustering and estimation simultaneously. We analyze a new data set integrated from Alzheimer’s disease (AD) study to illustrate the presented model and methods.

1. Introduction

Random effects models have been widely used as a powerful tool for analyzing correlated data [1, 2]. The model features a finite number of random terms acting as latent variables to model unobserved factors; see [3] for a comprehensive review. Some authors have further proposed semiparametric mixed effect models by allowing for infinite dimensional random effects [4, 5]. Most of the aforementioned works draw inferences using frequentist approaches, while Bayesian approaches have been largely ignored because of the lack of computational feasibility and expediency. With the advent of the “supercomputer” era, Bayesian analyses have recently sparked much interest in the setting of random effects models for clustered data or longitudinal settings. However, there is scarce literature that has systematically reviewed the Bayesian works in the area.

By extending the traditional random effects models, recent research focus has shifted to study heterogeneous random effects or nonparametric distributions of random effects, which arise because of skewness of data, missing covariates, or unmeasurable subject-specific covariates [6]. The extended random effects models, termed semiparametric random effects models, improve statistical performance with added interpretability. Bayesian techniques, which provide a convenient means to model non-Gaussian distributions, have recently been proposed for semiparametric random effects model in a variety of settings ([7, 8], among others). The discreteness of the Dirichlet process makes it impossible as a prior for estimating a density. However, as a remedy by convolving with a kernel, Dirichlet process mixture plays an important role [9].

For censored outcome data, transformation models, which transform the time-to-event responses using a monotone function and link them to the covariates of interest, have surged as a strong competitor of the Cox model [10]. Moreover, the transformation model framework is fairly general. The Cox model and the proportional odd model [11] can be viewed as nonparametric transformation linear models with some specific error terms; see [12, 13, 14]. For correlated data, the transformation model naturally extends the semiparametric random effects model by directly incorporating random effects to the transformation functions, treating them as realizations of an underlying random function. Bayesian analyses have found much use in this new area. For example, the beta process has been found to be a reasonable candidate for the prior of the monotone transformation function [15, 16, 17].

This chapter focuses on the Bayesian analysis of the transformed linear model with censored data and in a clustered setting. In many biomedical studies, the observations are naturally clustered. For example, patients in observational studies can be grouped in analysis according to a variety of factors, such as age, race, gender, and hospital, in order to reduce the confounding effects. Following Mallick and Walker [18], we explore using a mixture of beta distributions and the beta process as the candidates for the prior distribution of the random transformation function [17, 19, 20].

The rest of this chapter is structured as follows. Section 2 reviews the use of the Bayesian approach to infer parametric random effects models. In the setting of survival analysis, Section 3 proposes a beta process prior to fit random effects model with nonparametric transformation functions, and Section 4 applies the method to study the progression of Alzheimer’s disease (AD). Section 5 concludes the chapter with future research directions.

2. Dirichlet process mixture prior

In parametric random effects models, we considered the situation that the distribution form of the random error term is unknown. Dirichlet process mixture (DPM) is used as the prior for the baseline distribution in that error terms used to be continuous random variables in most situations.

2.1 Linear mixed effects model

With a longitudinal data set $\left\{Y_i,x_i,z_i\right\}$ , we posit a mixed effects model with an AR(1) serial correlation structure:

$$\begin{array}{l}{\mathbf{y}}_i={\mathbf{x}}_i\mathit{\beta }+{\mathbf{z}}_i{\mathbf{b}}_i+{\mathbf{w}}_i,i=1,\dots ,m;\\ {\mathbf{w}}_i={\left(w_{i1},\dots ,w_{{\mathit{\text{in}}}_i}\right)}^T;w_{\mathit{ij}}=\rho w_{i,j-1}+{\mathrm{ϵ}}_{\mathit{ij}},j=2,\dots ,n_i,\end{array}$$ E1

where ${\mathbf{y}}_i={\left(y_{i1},\dots ,y_{{\mathit{\text{in}}}_i}\right)}^T$ with $y_{\mathit{ij}}$ being the $j$ th response of the $i$ th subject for $i=1,\dots ,m$ , $\beta$ is a $p\times 1$ vector of fixed effect parameters, ${\mathbf{b}}_i$ a $q\times 1$ Gaussian random vector representing the subject-specific random effects, ${\mathbf{x}}_i$ and ${\mathbf{z}}_i$ are $n_i\times p$ and $n_i\times q$ design matrices linking $\beta$ and ${\mathbf{b}}_i$ to ${\mathbf{y}}_i$ , respectively, ${\mathbf{w}}_i={\left(w_{i1},\dots ,w_{{\mathit{\text{in}}}_i}\right)}^T$ is an $n_i\times 1$ vector of model errors, $\rho$ is the autoregressive coefficient, and ${{\mathrm{ϵ}}_{\mathit{ij}}}_{\prime }^{}s$ are i.i.d. noises. When $\left\{{\mathrm{ϵ}}_{\mathit{ij}}\right\}$ is non-normal, we assume a mixture model:

$$f_G\left(\mathrm{ϵ}|{\sigma }^2\right)=\int \phi \left(\mathrm{ϵ}|u,{\sigma }^2\right)\mathit{dG}\left(u\right),$$ E2

where $\phi \left(\cdot |u,{\sigma }^2\right)$ is the probability density function for a normal random variable with mean $u$ and variance ${\sigma }^2$ and $G$ is an unspecified probability distribution of $u$ satisfying $\int \mathit{udG}\left(u\right)=0$ , which ensures that $\mathrm{ϵ}$ comes from a mean-zero mixture distribution.

Replacing the Dirichlet process by an equivalent Pólya urn representation, [8] employed an empirical likelihood approach with the moment constraints and developed a posterior adjusted Gibbs sampler for more precise estimation. The algorithm is computationally feasible.

2.2 Accelerated failure time model

We shift gears to study survival outcomes with a cluster structure. Denote the data set by $\left(T_{\mathit{ij}},X_{\mathit{ij}}\right),i=1,\cdots ,K,j=1,\cdots ,n_i$ , where $T_{\mathit{ij}}$ is the failure time of the $j$ th subject in the $i$ th cluster and $X_{\mathit{ij}}$ is a vector of associated covariates. To accommodate such data, we utilize a general accelerated failure time model:

$$log{T}_{\mathit{ij}}={\mathbf{X}}_{\mathit{ij}}^T\mathit{\beta }+{\epsilon }_{\mathit{ij}},i=1,\cdots ,K\text{and}j=1,\cdots ,n_i,$$ E3

where $\beta$ is a vector of $p$ -dim regression coefficients of interest and ${\epsilon }_{\mathit{ij}}$ are independent random errors following the distribution with density $f_i$ . [7] posed an exponential tilt on the distributions of error terms to incorporate the cluster heterogeneity. That is,

$$\frac{f_i\left(t\right)}{f_1\left(t\right)}=exp\left({\theta }_{0i}+{\mathit{\theta }}_i^{T}q\left(t\right)\right),i=2,\cdots ,K,$$ E4

where $q\left(t\right)$ is a $q$ -dimensional prespecified functions containing potential covariate information and ${\mathit{\theta }}_i$ is the corresponding parameter vector with ${\theta }_{0i}=log\left[{\left(\int exp\left({\mathit{\theta }}_i^{T}q\left(t\right)\right)f_1\left(t\right)\mathit{dt}\right)}^{-1}\right].$ Thus, ${\mathit{\theta }}_i$ represents the parametric random effects in the model. Li et al. [7] place the DPM prior on the baseline density $f_1$ to develop a set of procedures which improves estimation efficiency through information pooling.

3. Beta process prior

We now present a nonparametric random effects model for the clustered survival data with nonparametric monotone link functions. We employ a beta process as the prior for the baseline function.

Let $T_{\mathit{ij}}$ denote the failure time of the $j$ th subject in the $i$ th cluster, $X_{\mathit{ij}}$ be the covariate vector for the subject, and $C_{\mathit{ij}}$ be the potential censoring time to the $j$ th subject in the $i$ th cluster. Assume that $C_{\mathit{ij}}$ is independent of the failure time $T_{\mathit{ij}}$ . Let $Z_{\mathit{ij}}=min\left(T_{\mathit{ij}},C_{\mathit{ij}}\right)$ and let ${\delta }_{\mathit{ij}}=I\left(T_{\mathit{ij}}<C_{\mathit{ij}}\right)$ be the censoring indicator. Then the observed data can be described as

$$\left(Z_{\mathit{ij}},{\delta }_{\mathit{ij}},X_{\mathit{ij}}\right),i=1,\cdots ,n;j=1,\cdots ,n_i.$$ E5

Within each cluster, $T_{\mathit{ij}}$ is linked to $X_{\mathit{ij}}$ via the following transformation model:

$$ln{H}_i\left(T_{\mathit{ij}}\right)=X_{\mathit{ij}}^T\mathit{\beta }+ln{\epsilon }_{\mathit{ij}},i=1,2\cdots ,n,$$ E6

where ${\epsilon }_{\mathit{ij}}$ are i.i.d. variables with a known density function $f_{\epsilon }\left(\cdot \right)$ and $H_i\left(t\right)$ are unknown cluster-specific monotone functions, which are i.i.d. realizations of a random function and can be viewed as a nonparametric version of random effects for independent clusters. In a parametric setting, if we set $H_i\left(t\right)=texp\left(-b_i\right)$ with $b_i$ being a cluster-specific random effect, Eq. (6) reduces to a classical random effects model, which has been discussed in Section 2.2. The challenge, however, lies in how to draw inferences in such a nonparametric setting.

To proceed, let the coefficient vector $\mathit{\beta }$ be a $p$ -dim unknown vector of interest. We further assume $H_i$ ′s are differentiable with derivative $h_i\left(t\right)=H_i^{\prime }\left(t\right)$ , and then the likelihood based on the observed data is

$$L\left(\mathit{\beta },H_1,\cdots ,H_n|\mathit{\text{data}}\right)=\prod _{i=1}^n\prod _{j=1}^{n_i}p\left(T_{\mathit{ij}},X_{\mathit{ij}},{\delta }_{\mathit{ij}}|H_i,\beta \right),$$ E7

where

$$p\left(t,x,\delta |H,\beta \right)={\left(f_{\epsilon }\left(H\left(t\right)e^{-x^T\mathit{\beta }}\right)h\left(t\right)e^{-x^T\mathit{\beta }}\right)}^{\delta }{S}_{\epsilon }{\left(H\left(t\right)e^{-x^T\mathit{\beta }}\right)}^{1-\delta }.$$

Here $S_{\epsilon }$ is the survival function of $\mathit{\text{varepsilon}}$ defined by $S_{\epsilon }\left(s\right)=P\left(\epsilon \ge s\right)$ .

We develop a Bayesian inference procedure based on model (6). We assume that the regression coefficient $\beta$ follows a normal prior:

$$\beta \sim N_p\left(0,{\sigma }_{\beta }^2{I}_p\right),$$ E8

where $I_p$ is the $p\times p$ dimensional identity matrix. Since $H_i$ is assumed differentiable, we model it with a kernel convolution:

$$H_i=\int {\mathrm{\Phi }}_{\sigma }\left(\cdot -s\right){\mathit{dB}}_i\left(s\right),$$

where $B$ is an increasing function and ${\mathrm{\Phi }}_{\sigma }$ is the zero-mean normal distribution with variance ${\sigma }^2$ . Hence, the derivative of $H_i$ is

$$h_i=\int {\varphi }_{\sigma }\left(\cdot -s\right){\mathit{dB}}_i\left(s\right)$$

with ${\varphi }_{\sigma }\left(t\right)=\frac{1}{\sigma }\varphi \left(\frac{t}{\sigma }\right).$ This actually mimics the idea of DPM to smooth beta process by convolution.

We are in a position to select an appropriate stochastic process used as the prior of $B_i$ . Beta process, as studied by [16, 17], is an ideal candidate for the prior of a monotone function. Specifically, beta process $\mathrm{BP}\left(\gamma ,B_0\right)$ with concentration parameter $\gamma$ and a base measure $B_0$ is an increasing Lévy process with independent increments of the form

$$\mathit{dB}\left(t\right)\sim \text{Beta}\left(\gamma {\mathit{dB}}_0\left(t\right),\gamma \left(1-{\mathit{dB}}_0\left(t\right)\right)\right).$$

Teh et al. [20] showed that a sample from $\mathrm{BP}\left(\gamma ,B_0\right)$ could be represented as

$$B_i\left(y\right)=\sum _{l=1}^{\infty }{p}_{\mathit{il}}I\left({\theta }_{\mathit{il}}\le y\right),$$ E9

where $p_{\mathit{il}}={\prod }_{j=1}^l{\nu }_{\mathit{il}}$ and $\left({\theta }_{\mathit{il}},{\nu }_{\mathit{il}}\right)$ follows

$${\theta }_{\mathit{il}}\sim B_0\left(\theta \right),{\nu }_{\mathit{il}}\sim \text{Beta}\left(\gamma ,1\right)l=1,2,\cdots .$$

In practice, we need to approximate samples of $\mathrm{BP}\left(\gamma ,B_0\right)$ with a finite dimensional form. Since beta process $\mathrm{BP}\left(\gamma B_0\right)$ can be represented by a stick-breaking process defined in Eq. (9), a natural approximation is obtained by retaining its first $L$ components. That is,

$$B_i^{\ast }=\sum _{l=1}^L{p}_{\mathit{il}}{\delta }_{{\theta }_{\mathit{il}}},$$

with $p_{\mathit{il}}={\prod }_{j=1}^l{\nu }_{\mathit{il}},l=1,\cdots ,L$ . Denote ${\mathit{\xi }}_i={\left({\nu }_{i1},\cdots ,{\nu }_{\mathit{iL}},{\theta }_{i1},\cdots ,{\theta }_{\mathit{iL}}\right)}^T$ and define

$$H_{\sigma }^{\ast }\left(z,{\mathit{\xi }}_i\right)=\sum _{l=1}^L{p}_{\mathit{il}}{\mathrm{\Phi }}_{\sigma }\left(z-{\theta }_{\mathit{il}}\right),h_{\sigma }^{\ast }\left(z,{\mathit{\xi }}_i\right)=\sum _{l=1}^L{p}_{\mathit{il}}{\varphi }_{\sigma }\left(z-{\theta }_{\mathit{il}}\right).$$

The approximated posterior based on the truncated DP is

$$\pi \left(\beta \right)\prod _{i=1}^n\left[{\pi }^{\xi }\left({\mathit{\xi }}_i\right)\prod _{j=1}^{n_i}f(Z_{\mathit{ij}},X_{\mathit{ij}},{\delta }_{\mathit{ij}}|\beta ,{\mathit{\xi }}_i)\right],$$ E10

where

$$\begin{array}{ll}f\left(z,x,\delta |\beta ,\mathit{\xi }\right)& ={\left(p_{\epsilon }\left(H_{\sigma }^{\ast }\left(z,\mathit{\xi }\right)exp\left(-x^T\mathit{\beta }\right)\right)h_{\sigma }^{\ast }\left(z,\mathit{\xi }\right)exp\left(-x^T\mathit{\beta }\right)\right)}^{\mathit{\delta }}\\ & \times {\left(P_{\epsilon }\left(H_{\sigma }^{\ast }\left(z,\mathit{\xi }\right)exp\left(-x^T\mathit{\beta }\right)\right)\right)}^{1-\delta }.\end{array}$$

The samples for $\mathit{\beta }$ and $\left({\mathit{\xi }}_1,\dots ,{\mathit{\xi }}_n\right)$ based on the posterior can be obtained with Markov chain Monte Carlo (MCMC) [21]. In our simulation, we use the R-package MCMC (https://cran.r-project.org/web/packages/mcmc/index.html) to draw samples for ${\mathit{\xi }}_1,\dots ,{\mathit{\xi }}_n$ and $\beta$ and use the Metropolis algorithm with a normal working distribution.

4. An application to Alzheimer’s disease neuroimaging initiative

Alzheimer’s Disease Neuroimaging Initiative (ADNI) is a multisite cooperative study for the purpose of improving the prevention and treatment of Alzheimer’s disease. The subjects in the study fall into three groups, cognitively normal (CN) individuals, mild cognitive impairment (MCI) patients, and early AD patients. ADNI provides a rich array of patients’ information, including functional magnetic resonance imaging (fMRI), positron emission tomography (PET), longitudinal functional cognitive tests scores, blood samples, genetics data, and censored failure time outcomes. Details of the study can be found at http://adni.loni.usc.edu.

We focus on the MCI group. MCI is recognized as a transitional stage between normal cognition and Alzheimer’s disease. The failure time is defined to be the time that a MCI patient is diagnosed with AD, which will be censored if a MCI patient remains at the MCI stage at the end of the follow-up time. Wide heterogeneities are exhibited among the failure times, which may be due to demographics and a variety of functional clinical biomarkers, such as the brain areas of the hippocampus, ventricles, and entorhinal cortex. The goal of the analysis is to study the impact of risk factors on progression to AD.

Using the same data as analyzed by [14], we demonstrate our methodology by modeling the failure time (the observed time of AD diagnosis from MCI stage in year) of 281 MCI patients on gender (0 = female, 1 = male), years of education, the number of apolipoprotein E alleles (0, 1, or 2), and the baseline hippocampal volume.

As age is a strong confounder but the functional form of its impact has not reached consensus, we elect to model its impact nonparametrically. Specifically, we use age to form two strata (below and above the median age) and use model (6) to estimate the stratum-specific transformation functions and the effects of other covariates. For comparisons, we also fit model (6) with age as a continuous variable and with a common transformation function. That is, we do not assume the data are clustered. For both models, the regression errors $\epsilon$ ’s are assumed to follow an exponential distribution with mean 10. In our calculation, we approximate the BPs by a finite truncation with $L=20$ . We assume the precision parameter $\alpha =1$ and scale parameter ${\sigma }^2\sim 1/{\sigma }^2$ .

Figure 1 illustrates the estimated transformation function $H$ of the failure time without clustering. The posterior means (PM) and standard errors (SE) of the regression coefficients in the model are reported in Table 1 . We run the MCMC for 20,000 iterations with the first 4000 draws discarded as burn-in samples and use Geweke’s statistic to ensure the convergence of the chains.

Figure 1.

Smoothed transformation function without clustering.

	RID	AGE	PTGENDER	PTEDUCAT	APOE4	Hipp.
PM	−0.9635	0.0069	−0.1453	−0.0231	−0.1817	0.2710
SE	1.3288	0.0841	1.2331	0.1835	0.8616	0.5333

Table 1.

Posterior estimates of regression coefficients with standard errors.

Figure 2 illustrates the estimated transformation functions with age-stratified data, and Table 2 summarizes the posterior means and standard errors of the other regression coefficients.

Figure 2.

Smoothed transformation functions with two age-strata: The left curve is the smoothed transformation function for group aged below the average age; the right curve is the smoothed transformation function for the group aged over the average age.

	RID	PTGENDER	PTEDUCAT	APOE4	Hipp.
PM	−0.6399	−0.0706	−0.0072	−0.1349	0.1919
SE	0.9273	0.8491	0.1267	0.6098	0.3716

Table 2.

Posterior estimators of regression coefficients with standard errors.

The left curve is relatively flat, while the right curve has a sharper slope. This is consistent with the recognition that AD is an aging disease: elder people above a certain age threshold tend to progress faster from MCI to AD.

Both Tables 1 and 2 show that none of the biomarkers are significant, whereas they are statistically significant in the analysis of [14]. One possible conjecture is that our nonparametric transformation functions may have well captured the effects of unobserved confounders, which may leave little to be explained by the observed covariates. More thorough investigation is warranted.

5. Future directions

Following [12], we can extend the transformation model (6) by allowing the error function $f_{\epsilon }$ to be unspecified. In this case, we need to specify the regression coefficient $\beta$ to obey some constraints such as ${\beta }_1=1$ or $\parallel \mathit{\beta }\parallel =1$ for identifiability. We will propose to model the error function using a Dirichlet processes mixture model:

$$f_{\epsilon }\left(t\right)=\int \phi \left(t|\mu ,{\sigma }^2\right)\mathit{dG}\left(\mu ,{\sigma }^2\right),G\sim \mathrm{DP}\left(\alpha ,G_0=\mathrm{N}\left(\mu |{\mu }_0,{\sigma }_0^2\right)\times \mathrm{IG}\left({\alpha }_1,{\alpha }_2\right)\right),$$

where $\phi \left(t|\mu ,{\sigma }^2\right)$ is a normal kernel with mean $\mu$ and variance ${\sigma }^2$ and $G$ are samples from a Dirichlet process $\mathrm{DP}\left({\alpha }_1,G_0=\mathrm{N}\left(\mu |{\mu }_0,{\sigma }_0^2\right)\times \mathrm{IG}\left(a,b\right)\right)$ , where ${\alpha }_1$ is the mass parameter and $\mathrm{IG}\left(\cdot |a,b\right)$ is the inverse gamma distribution with shape parameter $a$ and scale parameter $b$ .

In a slightly different context, we may also consider clustering observations by developing a new nested beta-Dirichlet process prior with companion MCMC algorithms. As there are limited works on functional random effects models that accommodate clustering structures observed, for example, from neural studies, we may propose a nested Dirichlet process [19] as the prior of Dirichlet process to cluster cumulative distribution functions successfully. We envision that such a nested Bayesian procedure will provide substantial computational expedience for practitioners and can certainly be applied to studies that cover beyond the neurodegenerative and aging diseases.

Acknowledgments

Shen’s research is partially supported by Beijing Natural Science Foundation 1192006 and National Natural Science Foundation of China; Liu’s research is partially supported by General Research Fund, Research Grants Council, Hong Kong, 15327216, and the Hong Kong Polytechnic University grant YBTR. Data collection and sharing for this project were funded by the Alzheimer’s Disease Neuroimaging Initiative (ADNI) (National Institutes of Health Grant U01 AG024904) and DOD ADNI (Department of Defense award number W81XWH-12-2-0012). ADNI is funded by the National Institute on Aging, the National Institute of Biomedical Imaging and Bioengineering, and through generous contributions from the following: AbbVie, Alzheimer’s Association; Alzheimer’s Drug Discovery Foundation; Araclon Biotech; BioClinica, Inc.; Biogen; Bristol-Myers Squibb Company; CereSpir, Inc.; Cogstate; Eisai Inc.; Elan Pharmaceuticals, Inc.; Eli Lilly and Company; EuroImmun; F. Hoffmann-La Roche Ltd and its affiliated company Genentech, Inc.; Fujirebio; GE Healthcare; IXICO Ltd.; Janssen Alzheimer Immunotherapy Research Development, LLC.; Johnson & Johnson Pharmaceutical Research Development LLC.; Lumosity; Lundbeck; Merck Co., Inc.; Meso Scale Diagnostics, LLC.; NeuroRx Research; Neurotrack Technologies; Novartis Pharmaceuticals Corporation; Pfizer Inc.; Piramal Imaging; Servier; Takeda Pharmaceutical Company; and Transition Therapeutics. The Canadian Institutes of Health Research is providing funds to support ADNI clinical sites in Canada. Private sector contributions are facilitated by the Foundation for the National Institutes of Health (www.fnih.org). The grantee organization is the Northern California Institute for Research and Education, and the study is coordinated by the Alzheimer’s Therapeutic Research Institute at the University of Southern California. ADNI data are disseminated by the Laboratory for Neuro Imaging at the University of Southern California.

Data used in preparation of this article were obtained from the Alzheimers Disease Neuroimaging Initiative (ADNI) database (adni.loni.usc.edu). As such, the investigators within the ADNI contributed to the design and implementation of ADNI and/or provided data but did not participate in analysis or writing of this report. A complete listing of ADNI investigators can be found at: http://adni.loni.usc.edu/wp-content/uploads/how_to_apply/ADNI_Acknowledgement_List.pdf.