Multiple Comparisons Using Composite Likelihood in Clustered Data

1 Introduction

The prevalence of depression in seniors estimated by the World Health Organization varies between 10% to 20% [1]. Understanding the relationship between depression and other health factors can help prevent the disease and alleviate the symptoms. The health and retirement study (HRS) conducted by the University of Michigan is a longitudinal study measuring various aspects of health, retirement and aging, as well as depression status. In this study, seniors were measured every two years from 1994 to 2012. The objective of our analysis is to estimate the effect of several health factors known to be associated with depression status and compare the effect sizes of different factors. Multiple comparisons on the effect sizes will clarify the relative importance of factors on the disease. For example, sleeplessness and smoking both contribute to the occurrence rate of depression. One might question if one factor is more important than the other for development of the disease. Therefore, to fully understand the effects of the health factors, we perform multiple comparisons on their effect sizes. It is important to note that we only observe associations but do not establish causal relationships of health factors with the disease.

The repeated binary measurements of depression status observed in this data set are correlated within individuals. These repeated measurements can be viewed as “clustered” data since they are recorded from the same experimental unit multiple times. Examples of clustered data arise in many other situations, including measurements coming from siblings or same pedigrees, or measurements taken in close proximity to each other in spatial data. Ignoring existing correlations within clusters leads to invalid individual or multiple inferences.

When performing multiple comparisons in clustered data, one should therefore take into account the correlation structure within the clusters. However, full likelihood analyses on such data often encounter computational challenges. Composite likelihood methods are extensions of the likelihood method that project high-dimensional likelihood functions to low-dimensional ones [2, 3]. This dimension reduction is achieved by compounding valid marginal or conditional densities. It has been shown that, under regularity conditions, the composite likelihood estimator has desirable properties, such as consistency and asymptotic normality [2, 3, 4, 5]. This makes it an appealing alternative in inferential procedures. Furthermore, composite likelihood is more computationally convenient than full likelihood at a cost of some loss of efficiency. The magnitude of this loss depends on the dimension of the multivariate vector and its dependency structure. Composite likelihood methodology has been applied to numerous statistical problems [6, 7, 8], however, the potential of composite likelihood in multiple testing has yet to be explored.

Multiple testing procedures have been developed to control the overall type I error rate when the number of tests is greater than one [9, 10]. Family-wise error rate (FWER) is defined as the probability of falsely rejecting at least one true null hypothesis. The classical Bonferroni method is the simplest procedure to adjust the FWER. It assumes maximum negative correlation between test statistics, and the resulting FWER does not exceed the significance level in any setting given the per comparison error is controlled correctly. However, it is well-known to be very conservative. The Dunn-Sidák procedure [11] uses a slightly less conservative p-value threshold for each comparison. This procedure assumes uncorrelated test statistics that explains the power increase compared to the Bonferroni method. The procedure provides the exact control of the FWER under independence, and it is conservative under positive dependence and it is liberal under negative dependence. For example, testing Ti and −Ti simultaneously is not possible with Sidák. Schéffe [12] established a method for testing all possible linear comparisons among a set of normally distributed variables, which tends to be over-conservative for a finite family of multiple comparisons. Scheffé controls all possible linear combinations, which makes a comparison with the other methods quite difficult, as they are testing different sets of hypotheses, where a smaller set of hypotheses of interest will in most cases result in a higher power. Several stage-wise procedures have also been proposed to improve the power. All of these methods present some shortcuts to avoid testing all intersection hypotheses in closed testing. Simes [13] modified the Bonferroni procedure based on ordered p-values. Holm [14] proposed a multi-stage procedure based on closed testing procedure. The method adjusts the FWER in each step using the number of remaining null hypotheses. Hommel [15] suggested a stage-wise rejective multiple test based on Simes inequality. All of these methods are less conservative and therefore more powerful than the Bonferroni method. However, it is difficult to construct simultaneous confidence intervals based on a stage-wise approach. As another alternative, Hothorn et al. [16] proposed to use equi-coordinate quantiles of the multivariate normal and multivariate t distribution to perform multiple comparisons in parametric methods. This corresponds to the control of the FWER based on the maximum of test statistics. The MNQ method allows for any correlation between test statistics, but needs the assumption of an asymptotic multivariate normal distribution for the test statistics, which is most of the time validated by the multivariate central limit theorem. This method offers sharp control of the FWER. The approach has been employed in many parametric and nonparametric settings to provide both multiple inferences and simultaneous intervals [17, 18].

In this paper, we propose a new procedure where multiple testing methods are combined with composite likelihood inference. This is motivated by the multiple testing problems for which full likelihood inference can be rather computational complex. We show that the composite likelihood test statistics for multiple hypotheses offer strong control of the FWER with the use of closed testing procedures. We explore in detail different multivariate models for correlated clustered data including the multivariate normal, multivariate probit, and quadratic exponential models to illustrate our multiple comparisons approach. Among these methods, the multivariate normal quantile threshold appears to have the best control of the FWER in most simulation settings.

The structure of this paper is as follows: In Section 2, we develop our composite likelihood based test statistics for multiple inferences and establish their asymptotic properties. In Section 3, we provide details on how to apply the general approach for several multivariate models. In Section 4, we conduct simulation studies to evaluate empirical performance of the proposed method. Finally, we analyze the depression data set to demonstrate the practical utility of the method.

2 Multiple comparisons procedures based on composite likelihood

Let {f(Y;θ),θ∈Θ}, where θ=(θ1,…,θp)T, be a parametric statistical model with parameter space Θ⊂Rp. Let Y=(y1T,⋯,ynT) denote the response variables, where yi=(yi1,⋯,yimi)T is the vector of observations from cluster i, i=1,⋯,n from a study population. It is assumed that observations across different clusters are independent, whereas observations within the same cluster may be dependent. Note that overall sample size is ∑i=1nmi. We assume that the cluster size, mi, is uniformly bounded.

Let

denote the contrast matrix. A family of c linear combinations of the parameters can then be specified by Cθ. Let H0l denote the hypothesis that C(l)θ=0, for l=1,…,c. We focus here on jointly testing the family of hypotheses H0l,l=1,…,c. The multiple testing procedure has weak control of the FWER if the FWER ≤α when all of the null hypotheses are true, whereas one has strong control of the FWER if the FWER ≤α under any combinations of null hypotheses and alternative hypotheses. Composite likelihood was proposed as a pseudo likelihood inference method that has attracted much attention in recent years [3, 4, 5, 19]. Two popular examples of composite likelihood are the “univariate marginal likelihood” and “univariate conditional likelihood”. In the former, we have the composite likelihood functionCL(θ;Y)=∏i=1n∏j=1mif(yij),

where any dependence structure is ignored. In the latter, we take CL(θ;Y)=∏i,jf(yij|yi(−j))=∏i,jf(yi)/f(yi(−j)), where yi(−j) denotes the sub-vector of yi with its j th element removed. We consider both types of approaches in our examples. In general, composite likelihood is a compounded form of marginal or conditional likelihoods, which is often easier to maximize than full likelihood. For Ak⊂{(i,j):j=1,…,mi,i=1,…,n}, let YAk={Yij,(i,j)∈Ak} denote a subset of the data, where k=1,…,K. The composite likelihood function is then defined as

where f(yAk;θ) is the density for the subset vector yAk, and wAk are some suitably chosen weights. In practice, the type of composite likelihood should be chosen so that the resulting composite score equation is consistent for the parameters, and the computation complexity is sufficiently manageable.

The maximum composite likelihood estimate (MCLE) is defined as

Xu and Reid [20] give precise conditions under which θˆnc is consistent for θ. Under appropriate assumptions, n(θˆnc−θ) is also asymptotically normally distributed with mean zero and limiting variance given by the inverse of the the Godambe information matrix [3, 21], where

with H(θ)=limnE(−cl(2)(θ;Y))/n and J(θ)=limnvar(cl(1)(θ;Y))/n. Here, cl(1) is the vector of first derivatives and cl(2) is the matrix of second order derivatives of cl(θ;Y)=logCL(θ;Y) with respect to θ. The matrix H(θ) can be estimated as the negative Hessian matrix evaluated at the maximum composite likelihood estimator, whereas the matrix J(θ) can be estimated as the sample covariance matrix of the composite score vectors. Both estimators, which we denote as Hˆn and Jˆn, are consistent [21].

Consider the hypothesis test on a family of linear combinations of the parameters: {H0:Cθ=0}. Denote by Γ=G−1(θ) the inverse Godambe information matrix, and let Γˆn denote the consistent estimator of Γ, where Γˆn=Hˆn−1JˆnHˆn−1. We propose the following test statistics for our hypothesis test

The limiting distribution of Tn=(T1,n,⋯,Tl,n)T is multivariate normal MVN(0,V), where

Furthermore, since Vi,i=1, the marginal asymptotic distribution of each individual Tl,n is standard normal. In practice, we estimate V by plugging Γˆn as a consistent estimator of G−1(θ) into eq. (3). This results in a consistent estimator of V.

We illustrate a few multiple comparison procedures as follows.

1. The Bonferroni procedure: The global intersection hypothesis ∩l=1cH0l will be rejected if maxl|Tl,n|>Z1−(α/2c), where Z1−(α/2c) denotes the critical value for standard normal random variable with P(Z>Z1−(α/2c))=α/2c. Each individual hypothesis H0l will be rejected if |Tl,n|>Z1−(α/2c).

2. The Holm’s procedure: For each H0l, evaluate the p-value pl=2P(Z>|Tl,n|). Order the p-values from the least to the greatest as p(1),…,p(c) and the corresponding hypotheses are reordered as H(01),…,H(0c). The global intersection hypothesis ∩l=1cH0l will be rejected if p(1)≤α/c. Let k denote the smallest l so that p(l)>α/(c−l+1). If k>1, then the individual hypotheses H(01),…,H(0,k−1) will be rejected.

3. The MNQ procedure: The global intersection hypothesis ∩l=1cH0l will be rejected if maxl|Tl,n|>Q1−α,V, where Q1−α,V denotes the equi-coordinate quantile for multivariate normal random vector Z=(Z1,…,Zc) with Z∼N(0,V), and P(maxl|Zl|>Q1−α,V)=α. Each individual hypothesis H0l will be rejected if |Tl,n|>Q1−α,V.

As the variance estimators are consistent regardless of the configuration of the hypotheses, the multivariate distribution of any subset of the test statistics still follows a multivariate normal distribution. Therefore, if a closed testing procedure is applied on the test statistics, the procedure has a strong control of the FWER [22, 23].

It is worthy to point out that the test statistics we propose here are Wald-type statistics which are not invariant under re-parametrization. Under re-parametrization, the new statistics follow the same type of limiting distributions, but the values of the statistics are not the same. This is a standard limitation associated with Wald-type statistics.

The multivariate distribution of Tn can be approximated by a multivariate t distribution. The denominator C(l)TΓˆnC(l)/n has an asymptotic equivalent distribution as C(l)TH−1JˆnH−1C(l)/n based on Slutsky’ Theorem. Furthermore, Jˆn=(cl(1))Tcl(1) asymptotically follows a Wishart(J,n). This entails that asymptotically C(l)TH−1JˆnH−1C(l) follows σl2χn2, where σl2=C(l)TH−1JH−1C(l). Reformulate the test statistics as

where the numerator is asymptotically a multivariate normal MVN(0,V), and the denominator is asymptotically χn2/n. Therefore, the multivariate distribution of Tn can be approximated as a multivariate t(V,n), where V is the covariance matrix and n is the degrees of freedom.

The MNQ method also facilitates the construction of simultaneous confidence intervals. The simultaneous (1−α)100% confidence interval for Cθ is

In some applications, the collection of effect sizes are nonlinear monotone transformations of the parameters. For example, we obtain odds ratio from log odds ratio by applying the exponential function. Let G(.) denote a nonlinear monotone transformation. Then the simultaneous (1−α)100% confidence interval for G[(C(l))Tθ],l=1,…,c, is

3 Three multivariate models

To demonstrate the application of our methodology, we consider three specific multivariate distributions: multivariate normal, multivariate probit, and quadratic exponential distributions. In this Section, we give the details of the distribution and computation of the test statistics; Section 4 examines the behavior of the test statistics via simulations.

For the first two distributions, the composite likelihood is constructed as a product of univariate marginal likelihoods. We choose this type of composite likelihood as the univariate marginal likelihood of these two distributions follows the univariate exponential family for which the estimation of parameters is not difficult. For the quadratic exponential distribution, the composite likelihood is constructed as a conditional likelihood. For the quadratic exponential model the normalizing constant is computationally intensive to obtain, and in this composite likelihood the normalizing constant is canceled out, greatly simplifying the computational burden of estimation.

In order to include covariates into our modelling scheme, let Xi denote an mi×p matrix containing the values of p covariates for the mi individuals in the ith cluster and β=(β1,…,βp)T denote the vector of regression coefficients. Let x⃗ij denote the jth row of the matrix Xi (this is the vector of covariates for individual j in cluster i).

4 Simulation results

We evaluate the validity of our proposed approach on three different multivariate models from Section 3 using simulations. We randomly simulate covariates from independent normal distributions. We consider different values of the regression parameters and correlation structures. We simulate multivariate clustered data under varying cluster sizes and varying number of clusters. We perform multiple comparisons on the regression parameters using the proposed composite likelihood method. For the multivariate normal case, we also compare the proposed method with the likelihood approach.

We test two different global null hypotheses on the regression coefficients β1,⋯,βp: many-to-one comparisons, H0:∩i=2p{β1=βi}; and all pairwise comparisons H0:∩1≤i≠j≤p{βi=βj}. The results for many-to-one comparisons are summarized here while the results for all pairwise comparisons are provided in the supplementary material. We choose a collection of multiple testing methods including one-step Bonferroni method and the MNQ method based on the multivariate distribution of the test statistics. The equi-coordinate critical values for multivariate normal and multivariate t distributions are obtained using the R package mvtnorm [25]. We also examine the performance of Dunn-Sidák method, Holm’s stage-wise procedure, and Scheffé’s projection method and the results are provided in the supplementary material.

Part of our goal is to show practitioners what happens if the correlation structure in the clustered data is ignored. To this end, we also include a “misspecified” scenario, where independence is erroneously assumed within the clusters. Due to the specific composite likelihood methods we use (univariate marginals and univariate conditionals), such a misspecification is equivalent to H(θ)=J(β) in eq. (1). This results in an estimate of Γˆn=Hˆn−1. This misspecified scenario is included for comparison, and we consider it only with the MNQ multiple comparison method (that is, the MNQ cutoff is calculated based on V estimated by plugging in Γˆn=Hˆn−1). Throughout, it is referred to as the “naive” approach.

In our simulations, we study the three models described in the previous section. For each model, a different sample size is needed for our asymptotic approximations to be valid. We determine this sample size with an initial simulation, before we proceed with our more in-depth investigations. For each simulation setting, 10000 simulated data sets are generated and the FWER is set to 0.05. The standard deviation for the observed FWER is hence approximately 0.002. These preliminary simulation results are given in Table 1. We observe that n=200,500 and 700 are required for the multivariate normal, multivariate probit and quadratic exponential models to maintain FWER within two standard deviations away from 0.05, respectively. These are the sample sizes used for the simulation results which follow.

To evaluate the power of each methods, we consider two different alternative scenarios: one alternative configuration a1 with only one non-zero parameter having a large effect size, and a second alternative configuration a2 with five true non-zero parameters but having small effect sizes for all. We are interested in the ability of the test to reject the global null hypothesis, but also in the ability of the test to reject the individual null hypotheses. Under the alternative scenario a1, we calculate the power to reject the global hypothesis (denoted as “a1” in the tables) and for the alternative configuration a2, we calculate both the power to reject the global null hypothesis (denoted as “a2” in the tables) and the average individual powers for the five individual true alternatives (denoted as “ind a2” in the tables).

5 Analysis of depression data

We apply our proposed method to the health and retirement study (HRS) dataset. Information about health, financial situation, family structure, and health factors were collected by the RAND center at the University of Michigan. Depression status is recorded as a binary response variable. Seven health factors include “age” (in months), “smoking”, “restless sleep”, “diabetes”, “high blood pressure”, “activity”, and “difficulty in walking” are considered as predictors and each predictor is highly significant with unadjusted p-values all less than 10−4. The effect size estimates and unadjusted p-values are reported in Table 6. As we are interested in comparing the importance of different predictors, we perform multiple comparisons on the effects of these seven health factors. For each individual we include only the years for which all of the factors were recorded. In total, there are 33,636 people included in the analysis and the number of repeated measurements vary across individuals. As the response variable is binary and the cluster sizes vary, we propose to use a quadratic exponential model to model this data set. The w parameter in the quadratic exponential model allows us to account for the interaction effect among the repeated measurements for the same individuals. The full likelihood approach is computationally challenging for this model, and hence we use the proposed composite likelihood method to perform inference.

To compare the effect sizes of all the seven health factors, we perform both all pairwise comparisons and many-to-one comparisons on the seven parameters. For all pairwise comparisons, the MNQ approach rejects 15 of the 21 hypothesis. The results are given in Table 7. Based on the estimates of the effect sizes (Table 6), we note that “restless sleep” and “smoking” are the two health factors with the largest effect sizes. The pairwise comparisons between restless sleep with all other health factors and smoking with all other factors are rejected. This shows that “restless sleep” and “smoking” are the two leading health factors for the occurrence of depression. “High blood pressure”, “diabetes”, and “difficulty in walking” have less influence on the occurrence rate of depression. When we examine the three pairwise comparisons among these three factors, the three null hypotheses are not rejected, indicating that these three health factors have similar effects on the disease. Furthermore, when we compare “high blood pressure” with “age” and “activity”, the test rejects the two comparisons, indicating that “high blood pressure” is more important than “age” and “activity” with regard to the disease development.

To show how different the results will be if the within-patient correlations are ignored, we also compare the result of the MNQ with the naive MNQ method. Both the MNQ and the naive MNQ reject the global null hypothesis that all pairs of health factors have equal effects on the depression status. The MNQ method rejects 15 hypotheses, whereas the naive MNQ method rejects 18 out of the total 21 hypotheses. MNQ and naive MNQ are in agreement in all the aforementioned comparisons. However, when we compare the effect sizes between age and diabetes, diabetes and activity, age and activity, the MNQ method cannot reject these three null hypotheses while the naive method rejects all three. The difference between the two methods is due to the correlation among the repeated measurements, which is estimated as wˆ=0.2877. By ignoring this correlation, as in the naive method, the standard errors are underestimated leading to more rejections.

We also conduct many-to-one comparisons on the seven health factors using “diabetes” as the baseline factor. Table 8 shows that the comparisons are consistent with our findings in Table 7. “Restless sleep” and “smoking” are two factors which are much more influential than “diabetes” in terms of increasing the risk of depression, whereas “activity” is one factor which is more important than “diabetes” in terms of decreasing the risk of depression. The “age”, “high blood pressure” and “difficulty in walking” seem to be have similar effect sizes as “diabetes”. Naive MNQ method rejects one more comparison than MNQ method indicating that naive method is also more liberal in many-to-one comparisons.

6 Discussion

For many correlated multivariate datasets, it is often computationally convenient to perform multiple comparisons based on the composite likelihood. Theory is developed based on the asymptotic properties of the composite likelihood test statistics. Sample sizes greater than 50 will be sufficient to apply the proposed multiple testing procedures in various models examined in the simulations. We establish the strong control of the FWER of our proposed composite likelihood test statistics with closed testing procedures. The method is illustrated for three different models: multivariate normal, multivariate probit and quadratic exponential. The equi-coordinate quantile of a multivariate normal distribution is used as a threshold for test statistics compared to some well-known traditional methods. This MNQ method, which is based on composite likelihood test statistics and uses multivariate normal quantiles to derive cut-off values for the test statistics, shows a better control of the FWER in most simulation settings, compared to the other test procedures.

7 Supplementary Files

We provide a supplementary file containing additional simulation results and technical derivations. An R package named CLMC is developed and can be downloaded from github account “m-azad”.