 Research article
 Open Access
 Published:
Permutationbased variance component test in generalized linear mixed model with application to multilocus genetic association study
BMC Medical Research Methodology volume 15, Article number: 37 (2015)
Abstract
Background
In many medical studies the likelihood ratio test (LRT) has been widely applied to examine whether the random effects variance component is zero within the mixed effects models framework; whereas little work about likelihoodratio based variance component test has been done in the generalized linear mixed models (GLMM), where the response is discrete and the loglikelihood cannot be computed exactly. Before applying the LRT for variance component in GLMM, several difficulties need to be overcome, including the computation of the loglikelihood, the parameter estimation and the derivation of the null distribution for the LRT statistic.
Methods
To overcome these problems, in this paper we make use of the penalized quasilikelihood algorithm and calculate the LRT statistic based on the resulting working response and the quasilikelihood. The permutation procedure is used to obtain the null distribution of the LRT statistic. We evaluate the permutationbased LRT via simulations and compare it with the scorebased variance component test and the tests based on the mixture of chisquare distributions. Finally we apply the permutationbased LRT to multilocus association analysis in the case–control study, where the problem can be investigated under the framework of logistic mixed effects model.
Results
The simulations show that the permutationbased LRT can effectively control the type I error rate, while the score test is sometimes slightly conservative and the tests based on mixtures cannot maintain the type I error rate. Our studies also show that the permutationbased LRT has higher power than these existing tests and still maintains a reasonably high power even when the random effects do not follow a normal distribution. The application to GAW17 data also demonstrates that the proposed LRT has a higher probability to identify the association signals than the score test and the tests based on mixtures.
Conclusions
In the present paper the permutationbased LRT was developed for variance component in GLMM. The LRT outperforms existing tests and has a reasonably higher power under various scenarios; additionally, it is conceptually simple and easy to implement.
Background
LRT for variance component in LMM
In many medical researches the main interest often lies on the problem of testing whether a random effects variance component is equal to zero within the mixed effects models framework. The variance component test has long been a statistical challenge and received much attention in the literature; see for example, Self and Liang [1], Stram and Lee [2], Liang and Self [3], more recently Lindquist, et al. [4], Drikvandi, et al. [5], Nobre, et al. [6], and many others. Beyond the direct interest, some other scientific problems under investigation can be also transformed into this similar issue. For example, to test a parametric null model against a nonparametric alternative in the penalized spline regression, Claeskens [7] first constructed a mixed effects model so that the hypothesis testing was reduced to the problem that whether the variance component of random effects was zero, and then performed a restricted likelihood ratio lackoffit test. Analogous transforms and hypothesis testing problems were also studied by Crainiceanu and Ruppert [810], Crainiceanu, et al. [11] and Greven, et al. [12].
The variance component test is nonstandard in the sense that under the null the variance component is on the boundary of the parameter space. In this situation it is not easy to obtain the null distribution of the likelihood ratio statistic. Under some regularity conditions Self and Liang [1], Stram and Lee [2] and Liang and Self [3] proved that the likelihood ratio statistic followed a 0.50:0.50 mixture of \( {\chi}_0^2 \) and \( {\chi}_1^2 \), where \( {\chi}_0^2 \) is a point mass at zero and \( {\chi}_1^2 \) is a chisquare distribution with one degree of freedom. Whereas Crainiceanu and Ruppert [8] argued that the assumed conditions in these papers were often not guaranteed in practice, they showed the mixture proportion parameter is actually dependent on specific contexts and the equalweight mixture can lead to conservative type I error control. To obtain the null distribution of the likelihood ratio test (LRT) statistic, Crainiceanu and Ruppert [8] developed a simulationbased algorithm using the spectral representation of the LRT statistic.
Instead of using the equalweight mixture as done in Self and Liang [1] and Liang and Self [3], Pinheiro and Bates [13] found that a nonequalweight one may be better and suggested the 0.65:0.35 mixture for some specific longitudinal datasets. Fitzmaurice, et al. [14] evaluated the 0.50:0.50 and 0.65:0.35 mixtures via simulations and concluded that the appropriate mixture is not straightforward to derive.
The work aforementioned shows clearly that it is difficult to obtain an analytical expression for the null distribution of the likelihood ratio statistic. On the other hand, when encountering complex hypothesis testing scenarios in practical data analyses, resorting to resamplingbased methods is a very natural and effective strategy. Faraway [15] and Samuh, et al. [16] applied the parametric bootstrap approach for testing the variance component. Lee and Braun [17] and Samuh, et al. [16] used the permutation procedure to resolve this problem. Their results demonstrated that the bootstrap and permutation tests can control the type I error rate correctly and are more powerful compared to the tests that are based on the usual asymptotic mixture.
LRT for variance component in GLMM
At present most of the work concerning the likelihoodratio based variance component test has been mainly investigated under the context of linear mixed models (LMM), in which the response variable is continuous and the closedform loglikelihood function is easily obtained and calculated [18,19]. However, little literature has been published about LRT of variance component in the generalized linear mixed models (GLMM) framework [20,21], where the response variable is discrete, such as the binary or count variable, and the calculation of the loglikelihood function is not straightforward.
The difficulty of performing the likelihood ratio variance component test in GLMM arises in several aspects: (i) unlike in LMM for the continuous data, the closedform loglikelihood function of GLMM can be no longer yielded except in some extremely limited situations, thus the computation of the likelihood ratio statistic is not easy; (ii) the exact parameter estimation in GLMM is often impossible because of the need of highdimensional numerical integration; (iii) the null distribution of the likelihood ratio statistic in GLMM is complicated, and similar to the case in LMM it cannot be obtained analytically.
Lee and Braun [17] mentioned in their Discussion Section that the permutation test can be applied to the problem of testing for variance components in GLMM, but did not formally give any further results to date. Fitzmaurice, et al. [14] investigated the permutation variance component test in GLMM to detect the cluster effect in a multilevel binary dataset and showed that the permutation test can effectively control the type I error rate and has a higher power compared with the tests based on the asymptotic mixture of chisquare distributions.
Lin [22] developed a global score test to examine variance components in GLMM via the manner of Taylor expansion of the integrated quasi loglikelihood, derived the analytical null distribution of the score statistic, and displayed that this score test was a locally asymptotically most powerful test under some assumptions; see also Zhang and Lin [23] for similar presentations where the score test was constructed for the semiparametric generalized additive mixed models [24]. Sinha [25] constructed a scoretype statistic for the variance component test in GLMM via parametric bootstrap procedure and showed that the bootstrap score test was valid.
A motivating application for the variance component test in GLMM
In this paper we attempt to perform the likelihood ratio variance component test in GLMM via the permutation procedure. To present the permutationbased LRT, we now give an example motivating from genetic statistics concerning multilocus association analyses in case–control studies. One of the main objectives in genetic association studies is to identify the causal markers (e.g., single nucleotide polymorphisms, SNPs) that are significantly related to diseases. Two commonly used methods in genomewide association studies are respectively single marker analysis and multiple markers analysis [26]. Compared with the single marker analysis, the multiple markers analysis is generally more powerful since some useful information among markers is taken into account.
More specifically, suppose that in multilocus association analyses K SNPs within a functional gene are grouped into an SNP set. The objective is to test whether these K SNPs are jointly associated with the disease of interest, say y. Here we assume y is a binary variable. If the number of SNPs (i.e., K) is small, we can perform traditional fixed effects model and test to achieve this aim; however, if the number of SNPs is large, the traditional test such as the Wald test or the score test is typically underpowered for the null due to the large number of degrees of freedom, and the fixed effects model may be subject to some numerical issues, e.g., the problem of multicollinearity since the SNPs within a gene may be highly correlated due to the linkage disequilibrium.
As an effective alternative, if the effects of SNPs are assumed to be random with a common variance component [27], then the fixed effects model becomes the mixed effects model and the test is converted into the variance component test under the framework of GLMM (i.e., the logistic mixed effects model). And now the numerical problems encountered in the fixed effects model mentioned above can be avoided because of the use of mixed effects model [2731].
It should be emphasized that our main objective is focused on the hypothesis testing of variance component in GLMM instead of the parameter estimation. In this paper we attempt to explore the application of LRT in logistic mixed effects model based on the welldeveloped theories for GLMM.
The organization of this paper is given as follows. The likelihood ratio statistic is constructed by taking full of the working response generated in the estimation process of GLMM. To circumvent the derivation of the analytical expression for the null distribution of the likelihood ratio statistic, we use the permutation procedure. We implement simulation studies to evaluate the performance of the proposed permutationbased LRT in the context of genetic association studies. Finally we apply the proposed LRT to multilocus association analysis for the binary genetic data.
Methods
Logistic mixed effects model and the variance component test
Let X = [X_{0}, X_{1}, …, X_{ p }] and Z = [Z_{1}, Z_{2}, …, Z_{ K }] be n × (p + 1) and n × K matrices for covariates, respectively, where X_{0} = 1 corresponds to the intercept term, n is the total sample size; and let y = [y_{1}, y_{2}, …, y_{ n }]’ be a ndimensional vector for a binary response variable coded as 0 and 1, such as y = 0 for the control and y = 1 for the case in case–control studies. Assume the relationship between y, X and Z is characterized via the following logistic model
where α = [α_{0}, α_{1}, …, α_{ p }] are considered as fixed effects, while β = [β_{1}, β_{2}, …, β_{ K }] are considered as random effects following a normal distribution with mean zero and variance component τ^{2}. Under these settings model (1) becomes the logistic mixed effects model which has been widely studied in the literature [20,21,32,33]. Our aim is to test βs simultaneously to determine the overall significance of variables Z. Based on these specifications, the test of β = 0 is immediately equivalent to the test of variance component H_{0}: τ^{2} = 0 in model (1).
Note that model (1) is different from the one considered by Fitzmaurice, et al. [14] and Sinha [25], in which only one random effect was included so that the approximated calculation of the loglikelihood function via numerical integration was feasible; whereas in model (1) there are K random effects with the same variance component and the calculation of the loglikelihood function via numerical integration is generally not possible [21].
Definition of the likelihood ratio statistic
A lot of algorithms have been developed for estimating GLMM, including approximate approaches and Monte Carlo methods [20,21,34]. Here we use the penalized quasilikelihood (PQL) algorithm [20,34] since it has the conceptual and computational advantage compared to others and can be implemented via existing software, such as the glmmPQL function in the R package MASS [35].
We build LRT for the null of H_{0}: τ^{2} = 0 in model (1) by using the working response variable, denoted as Y’ to distinguish the original response variable y, generated after the convergence of PQL algorithm
where R is a n × n diagonal matrix with elements being 1/[μ(1μ)]. Following Breslow and Clayton [20], the corresponding quasi loglikelihood up an unrelated constant for the working response variable Y’ is
with V = τ^{2}ZZ′ + R and \( \widehat{\boldsymbol{\alpha}}\kern0.5em =\kern0.5em {\left({\mathbf{X}}^{\mathbf{\prime}}{\mathbf{V}}^{1}\mathbf{X}\right)}^{1}{\mathbf{X}}^{\mathbf{\prime}}{\mathbf{V}}^{1}{\boldsymbol{Y}}^{\boldsymbol{\prime}}. \) By carefully inspecting we can easily find that Formula (3) is actually the loglikelihood function of LMM with residual term ε and random effects β if Y’ is treated as a continuous response variable. Based on the loglikelihood given in (3), the likelihood ratio statistic for H_{0}: τ^{2} = 0 can be defined as
Here we have some remarks about the LRT statistic given in (4): (i) we calculate T via the working response variable Y’ instead of the original response variable y; (ii) for fair comparisons, the loglikelihood of the null model, i.e., the value of L(τ^{2} = 0), is also calculated according to the quasi loglikelihood, although the null model with τ^{2} = 0 reduces to the general logistic regression whose loglikelihood can be computed directly; in fact we find that negative values of the likelihood ratio statistic are obtained if the loglikelihood of the reduced logistic regression is adopted; more importantly, we indeed do not see any problems with the use of Formula (4) to compute the likelihood ratio statistic in our simulations; (iii) here we only consider LRT, although the restricted likelihood ratio test [810] can be also applied; the difference between the two likelihood ratio tests is generally ignorable as the sample size increases.
Resampling algorithm for the null distribution of T
To obtain the null distribution of the likelihood ratio statistic T, we now use the permutation procedure [3638]. The more general resampling algorithm is described as follows.

Resampling algorithm for the null distribution of T

step 1: Calculate the observed value of T according to (4) using the original data D = [y, X, Z], say T^{obs};

step 2: Generate a new dataset D^{*} by resampling;

step 3: Calculate the new value of T according to (4) using D^{*}, say T^{*};

step 4: Repeat the steps 2–3 above B times, get T^{*b}, b = 1, 2, …, B;

step 5: The resamplingbased null distribution of T can be approximated by T^{*b}, and the Monte Carlo p value of T is estimated by the proportion of T^{*b} equal to or greater than T^{obs}.
In step 2 of the resampling algorithm, if we randomly permute y and X as a whole while keeping Z fixed, then we are performing the permutation test [14,38]. Here we implicitly assume that no correlation between X and Z exists. If we only generate the new response variable y (say y*) from the fitted null model (i.e., the fixed effects logistic model with only covariates X) while letting both X and Z unchanged, then we are performing the parametric bootstrap test [15,16,25,36,37]. The permutation and bootstrap algorithms are almost the same except the production of the new resampling dataset in the second step. In this paper we only discuss the permutation test since it requires fewer assumptions compared with the parametric bootstrap test.
Usually, the number of replicates B in the permutation test is set to 1000 ~ 2000 for a relatively large significance level, say α = 0.05. Whereas some authors [14,16] empirically found that much smaller values of B (e.g., 200 or 500) could also behave satisfactorily.
Simulation study and real application
In this section we implemented simulation studies to evaluate the proposed permutationbased LRT method and then applied it to the real genetic data. The simulation data was generated via the following logistic model
where X_{1} is a binary variable with a probability of 0.5 and X_{2} is a standard normal variable, and they are mutually independent. To correspond to the motivating example of multilocus association studies described before and to mimic the true genetic data, in this paper we use the additive genetic model so that Z = 0, 1 and 2 represents the count of the minor allele [31].
For comparisons, besides the proposed permutationbased LRT, the methods of LRT based on the 0.50:0.50 and 0.65:0.35 mixtures and the scorebased variance component test are also implemented. It should be mentioned that the 0.65:0.35 mixture suggested by Pinheiro and Bates [13] is only reasonable for some specific datasets and may not apply to the present simulation settings; but it has been shown that it provided a useful reference in some cases [8,14], thus it is still considered here. The score variance component tests in mixed effects models have been previously investigated by many authors in a wide range of application areas [22,23,2831,39]; these score tests share a very similar rationale to each other. In this paper we use the score test originally developed by Lin [22], which was later applied and called sequence kernel association test (SKAT) by Wu, et al. [31] under the context of rare variants association in sequencing data and is implemented in the R package SKAT (version 0.91) [40]. All the data analyses are performed within the R (version 2.15.3) statistical environment [41].
The sample size was set to 400, 600, 800 and 1000. The value of B in the permutationbased LRT was 1000. The number of simulations was 2000 for controlling the type I error rate and was 1000 for evaluating the statistical power. The estimated type I error rate and power were calculated as the proportion of the pvalues less than the given significance level α = 0.05.
Simulation 1
We first evaluated the type I error and power of the proposed LRT in the context of multilocus association analyses with different number of SNPs and various random effects variance component. To obtain Z = [Z_{1}, Z_{2}, . . ., Z_{ K }], we first generated G = [G_{1}, G_{2}, . . ., G_{ K }] via the standard normal distribution with pairwise correlation between G_{ k } and G_{k’} to be 0.5^{k  k’}, and then set Z_{ k } = 0, 1 and 2 according to whether G_{ k } < −c, −c ≤ G_{ k } ≤ c or G_{ k } > c, where c is the third quartile of the standard normal distribution. We set K = 20, 30 and 40, and β = 0 in model (5) for the type I error simulation and sampled β from a normal distribution with mean zero and variance τ^{2} = 0.15 and 0.20 for the power evaluation. The results are given in Tables 1 and 2.
Simulation 2
This simulation is mainly designed to see the influence of nonnormally distributed random effects in multilocus association studies. Here the variable Z in model (5) was generated via the software COSI [42] based on a coalescent model for European population. Specifically, we first generated a long region with COSI, from which a continuous subregion was selected so that it contained a reasonable number of SNPs; next this subregion was treated as a gene, whose genotypes were specified to Z, i.e., the SNPs included in this subregion are modeled by variables Z. The average number of SNPs included in the simulation is given in Table 3. Usually, not all of the SNPs (i.e., Z) within a gene are causal [31], thus some of them were assumed to have zero effects; i.e., only m out of K SNPs had an impact on the response variable y in model (5). Here we set m = 0 K, 0.1 K, 0.3 K and 0.5 K. The effects β were generated from the uniform distribution U[0.60, 1.10] at random. The results are given in Table 3.
Data analysis
We applied the permutationbased LRT to the wellknown Genetic Analysis Workshop 17 (GAW17) miniexome sequencing data [43], which was constructed based on the 1000 Genomes Project [44]. The GAW17 data contain a binary response variable called affection status (say y here, and y = 1 representing affection and y = 0 representing nonaffection) and five covariates including age (X_{1}, a continuous variable), smoking status (X_{2}, an indicator variable representing smoking or nonsmoking), and three continuous variables Q_{1}, Q_{2} and Q_{4}. The number of individuals is 697. The details of this data were described in Almasy, et al. [43]. We implemented multilocus association analysis for four genes ELAVL4, PRKCB1, PTK2B and SOS2 across all the 200 replicates of the GAW17 data.
The logistic mixed effects model for the GAW17 data is written as
here Z is the genotype of the selected gene, and α_{0} α_{5} are set to be fixed effects while β are set to be random effects. Our interest is to test H_{0}: τ^{2} = 0 in model (6) to determine whether a given gene is associated with affection status. Since in this data we know which SNPs are causal, the proportion of pvalues less than the given significance level provides a measurement of the empirical power for these methods. The results are given in Table 4.
To reduce computation time of the permutationbased LRT, the simulations and data analyses were conducted with the high performance computing at the School of Public Health of Nanjing Medical University. More specifically, we use the function %dopar% in the R package doMC (version 1.3.3) [45] with the function registerDoMC having its argument cores = 30 for parallel execution when implementing the permutation procedure. Our experience showed that this parallel execution could substantially reduce the computation compared to the use of the simple R forloop.
Results
From Table 1 it is observed that the type I error rate of the permutationbased LRT is very close to the nominal significance level 0.05 at various situations. On the other hand, the score test is sometimes slightly conservative, especially when the sample size is relatively small (e.g., 400 and 600); the two tests respectively based on the 0.50:0.50 and 0.65:0.35 mixtures cannot control the type I error rate correctly regardless of the sample size and the number of random effects, and are much more conservative than the score test.
Table 2 shows that the permutationbased LRT is the most powerful method among these tests and the score test has higher power compared to the tests based on mixtures. Typically, as expected the powers of these tests improve as the sample size increases and the number of random effects becomes larger. In Simulation 1 as the sample size improves, the powers of the LRT and the score test approach each other. Here an interesting finding in Table 2 is that the magnitude of the random effects variance component τ^{2} has a greater influence on the statistical power for these tests than the number of random effects (i.e., K). For example, in Table 2 when K becomes from 20 to 40 and τ^{2} = 0.15, the increasing power for LRT is 0.137; when K = 20 and τ^{2} increases from 0.15 to 0.20, the corresponding change of power of LRT is 0.210.
In Table 3 for Simulation 2 where the random effects are not normally distributed, as before the permutation test holds the type I error rate effectively while the score test is slightly conservative and the tests based mixtures give incorrect control. For instance, the type I error rate estimated from the 0.50:0.50 mixture is conservative and the type I error rate from the 0.65:0.35 mixture is liberal, and the incorrect control of the two mixtures does not improve with the increasing sample size. The power of the permutationbased LRT is higher than that obtained from the 0.50:0.50 mixture and lower than that obtained from the 0.65:0.35 mixture. But here only the results of the LRT are valid and believable because the two mixtures cannot correctly control the type I error rate as demonstrated.
In general, as expected the powers of these tests improve as the proportion of causal markers increases. In Table 3 we find that the permutationbased LRT always outperforms the score test even when the sample size is large (e.g., 1000 and 800). In conclusion in this simulation we see that the proposed LRT can maintain a reasonably higher power than the score test even if the distribution of random effects is not normal.
From Table 4 we can see that the permutationbased LRT has higher chance to detect the causal genes compared the other three methods. Here it is somewhat surprising that the score test is extremely underpowered and has the lowest probability to identify the associated genes. For example, when the significance level is 0.001 the chance that these four genes are deemed to be statistically related to the disease status by the score test is less than 5%; relatively, the LRT has a much larger probability than the score test.
As a reviewer pointed out that, in fact, the GAW17 data are not purely real data but rather 200 replicates, which were simulated independently based on the true genotypephenotype association model with the same genotypes but different phenotypes. Therefore, it may be limited for evaluating the performance of these tests. Nevertheless, according to the results presented in Table 4, it seems to be reasonable to conclude that the proposed permutationbased LRT has an empirical advantage to identify association signals compared with other competing tests.
Discussion
In this paper we have applied the LRT to GLMM by taking full use of the PQL algorithm and the resulting working response variable. The main difficulties of using LRT in GLMM are the computation of the loglikelihood function for the alternative model (e.g., the logistic mixed effects model) and the derivation of the null distribution of the likelihood ratio statistic. To avoid the direct computation we use the quasi loglikelihood of the working response, which can be deemed to be linear. By doing this the LRT in GLMM becomes computationally feasible. To obtain the null distribution of LRT we use the permutation procedure which has been extensively employed in practice and has proved to be valid [14,17,38].
In the present paper we only consider the permutation test for the likelihood inference, but extending to other resamplingbased methods is straightforward; for example, the parametric bootstrap method. Still the permutation test has an advantage that it requires few model assumptions compared to other resamplingbased methods [14,17,38]. Via simulations it has been shown that the permutationbased LRT can effectively control the type I error rate under different situations. Similar to the results demonstrated by others [8,14,17], the usuallyused asymptotic mixture distributions such as the 0.50:0.50 and 0.65:0.35 mixtures cannot maintain the type I error rate correctly at the nominal significance level. The simulations show that the score test can control the type I error rate but is sometimes conservative, especially for small sample sizes; the similar result has also been observed by Wu, et al. [31] and Li, et al. [46] under the context of multilocus association studies.
Our studies have shown that the permutationbased LRT has higher power than the score test and the tests based on mixtures under a wide range of scenarios and it still maintains a reasonably high power even when the random effects do not follow a normal distribution. The LRT has the advantage that it is conceptually simple and easy to implement, but the permutationbased LRT is computationally intensive. However, as Lee and Braun [17] pointed out that the computation burden can be reduced significantly via parallel computing; additionally, some approximation methods suggested in linear mixed effects models may be applied [8,12]. In contrast, the tests based on mixtures are relatively rapid since they do not require resampling, and the score test is the most computationally efficient since it only needs to fit the null model and its null distribution can be obtained analytically [2224]. Another advantage of the score test is that its statistic is calculated much more straightforward [2224].
The computation time for the permutationbased LRT depends on the sample size, the number of random effects and the number of replicates B in the permutation algorithm. Although having been performed the parallel computation, we still find that the permutationbased LRT is much slower than other methods, especially than the score test. For example, in Simulation 1 when the sample size is 400 and B = 1000 in the permutationbased LRT, running one simulation requires respectively 975 s, 702 s and 1217 s if K = 20, 30 and 40; whereas the score test only needs about 16 s, 23 s and 30 s if running ten simulations. All the times are obtained via the use of the function %dopar% as mentioned before.
Besides the increasing computational burden of the permutationbased LRT, another disadvantage is that the LRT may be numerically unstable to fit the alternative model when the sample size is limited and the number of the random effects is large. The third shortcoming of the LRT is that its validity and efficiency rely largely on the PQL algorithm, which had been already proved to be seriously biased downward [20,22,47,48]. Although no negative effect on the type I error controlling is observed in our simulations, certainly it can result in power reduction due to the biased estimation. Therefore, designing more efficient and faster algorithms for the proposed LRT is necessary in the future, and developing likelihoodratio based variance component tests in GLMM with the biascorrected PQL estimation is another interesting problem.
Finally, we note in Table 4 for GAW17 data analysis that the overall power for these methods compared in this paper is very low. For example, if we set the significance level at 0.001, none of these methods has an empirical power greater than 0.50. The reasons may be that: (i) the effects of the markers are very weak; (ii) the noncausal SNPs included in the gene lead to the reduction of power; that is, not all the random effects have nonzero coefficients; (iii) the sample size is relatively small for identifying the association signal; (iv) in the GAW17 data most of the SNPs are rare variants, i.e., their minor allele frequencies are very low and typically less than 0.01 [43], which are known to have an extremely limited power to be detected [31,49]. These results also suggest that developing more powerful methods for the lowfrequency variants in case–control sequence data is an urgent demand [31,5052].
Conclusions
In the present paper the permutationbased LRT was developed for variance component in GLMM. Via simulations and real application it has been demonstrated that the developed LRT method outperforms the existing tests and has a reasonably high power under various scenarios; additionally it is conceptually simple and easy to implement.
Abbreviations
 LRT:

Likelihood ratio test
 GLMM:

Generalized linear mixed model
 LMM:

Linear mixed model
 GWAS:

Genomewide association study
 PQL:

Penalized quasilikelihood
 GAW17:

Genetic analysis workshop 17
 SNP:

Single nucleotide polymorphism
References
 1.
Self SG, Liang KY. Asymptotic properties of maximum likelihood estimators and likelihood ratio tests under nonstandard conditions. J Roy Stat Soc B. 1987;82(398):605–10.
 2.
Stram DO, Lee JW. Variance components testing in the longitudinal mixed effects model. Biometrics. 1994;50(4):1171–7.
 3.
Liang KY, Self SG. On the asymptotic behaviour of the pseudolikelihood ratio test statistic. J Roy Stat Soc B. 1996;58(4):785–96.
 4.
Lindquist MA, Spicer J, Asllani I, Wager TD. Estimating and testing variance components in a multilevel GLM. Neuroimage. 2012;59(1):490–501.
 5.
Drikvandi R, Verbeke G, Khodadadi A, Partovi Nia V. Testing multiple variance components in linear mixedeffects models. Biostatistics. 2013;14(1):144–59.
 6.
Nobre J, Singer J, Sen P. Utests for variance components in linear mixed models. TEST. 2013;22(4):580–605.
 7.
Claeskens G. Restricted likelihood ratio lackoffit tests using mixed spline models. J Roy Stat Soc B. 2004;66(4):909–26.
 8.
Crainiceanu CM, Ruppert D. Likelihood ratio tests in linear mixed models with one variance component. J Roy Stat Soc B. 2004;66(1):165–85.
 9.
Crainiceanu CM, Ruppert D. Restricted likelihood ratio tests in nonparametric longitudinal models. Stat Sinica. 2004;14(3):713–30.
 10.
Crainiceanu CM, Ruppert D. Likelihood ratio tests for goodnessoffit of a nonlinear regression model. J Multivariate Anal. 2004;91(1):35–52.
 11.
Crainiceanu C, Ruppert D, Claeskens G, Wand MP. Exact likelihood ratio tests for penalised splines. Biometrika. 2005;92(1):91–103.
 12.
Greven S, Crainiceanu CM, Küchenhoff H, Peters A. Restricted likelihood ratio testing for zero variance components in linear mixed models. J Comput Graph Stat. 2008;17(4):870–91.
 13.
Pinheiro JC, Bates D. Mixedeffects models in S and SPLUS. 2nd ed. New York: Springer; 2009.
 14.
Fitzmaurice GM, Lipsitz SR, Ibrahim JG. A note on permutation tests for variance components in multilevel generalized linear mixed models. Biometrics. 2007;63(3):942–6.
 15.
Faraway JJ. Extending the linear model with R: generalized linear, mixed effects and nonparametric regression models. New York: Chapman & Hall/CRC; 2005.
 16.
Samuh MH, Grilli L, Rampichini C, Salmaso L, Lunardon N. The use of permutation tests for variance components in linear mixed models. Commun StatTheor M. 2012;41(16–17):3020–9.
 17.
Lee OE, Braun TM. Permutation tests for random effects in linear mixed models. Biometrics. 2012;68(2):486–93.
 18.
Verbeke G, Molenberghs G. Linear mixed models for longitudinal data. New York: Springer; 2009.
 19.
Laird NM, Ware JH. Randomeffects models for longitudinal data. Biometrics. 1982;38(4):963–74.
 20.
Breslow N, Clayton D. Approximate inference in generalized linear mixed models. J Am Stat Assoc. 1993;88(421):9–25.
 21.
Diggle P, Heagerty P, Liang KY, Zeger S. Analysis of longitudinal data. 2nd ed. New York: Oxford University Press; 2002.
 22.
Lin X. Variance component testing in generalised linear models with random effects. Biometrika. 1997;84(2):309–26.
 23.
Zhang D, Lin X. Hypothesis testing in semiparametric additive mixed models. Biostatistics. 2003;4(1):57–74.
 24.
Lin X, Zhang D. Inference in generalized additive mixed models by using smoothing splines. J Roy Stat Soc B. 1999;61(2):381–400.
 25.
Sinha SK. Bootstrap tests for variance components in generalized linear mixed models. Can J Stat. 2009;37(2):219–34.
 26.
Balding D. A tutorial on statistical methods for population association studies. Nat Rev Genet. 2006;7(10):781–91.
 27.
Zeng P, Zhao Y, Liu J, Liu L, Zhang L, Wang T, et al. Likelihood ratio tests in rare variant detection for continuous phenotypes. Ann Hum Genet. 2014;78(5):320–32.
 28.
Tzeng J, Zhang D. Haplotypebased association analysis via variancecomponents score test. Am J Hum Genet. 2007;81(5):927–38.
 29.
Kwee L, Liu D, Lin X, Ghosh D, Epstein M. A powerful and flexible multilocus association test for quantitative traits. Am J Hum Genet. 2008;82(2):386–97.
 30.
Wu MC, Kraft P, Epstein MP, Taylor DM, Chanock SJ, Hunter DJ, et al. Powerful SNPset analysis for case–control genomewide association studies. Am J Hum Genet. 2010;86(6):929–42.
 31.
Wu MC, Lee S, Cai T, Li Y, Boehnke M, Lin X. Rarevariant association testing for sequencing data with the sequence kernel association test. Am J Hum Genet. 2011;89(1):82–93.
 32.
Molenberghs G, Verbeke G. Models for discrete longitudinal data. New York: Springer; 2005.
 33.
Fitzmaurice GM, Laird NM, Ware JH. Applied longitudinal analysis. New York: John Wiley & Sons; 2004.
 34.
Wolfinger R, O’Connell M. Generalized linear mixed models a pseudolikelihood approach. J Stat Comput Sim. 1993;48(3–4):233–43.
 35.
Venables WN, Ripley BD. Modern applied statistics with S. 4th ed. New York: Springer; 2002.
 36.
Davison AC, Hinkley DV. Bootstrap methods and their application. Cambridge: Cambridge University Press; 1997.
 37.
Efron B, Tibshirani R. An introduction to the bootstrap. New York: Chapman & Hall/CRC; 1993.
 38.
Good P. Permutation, parametric, and bootstrap tests of hypotheses. 3rd ed. New York: Springer; 2005.
 39.
Goeman JJ, van de Geer SA, de Kort F, van Houwelingen HC. A global test for groups of genes: testing association with a clinical outcome. Bioinformatics. 2004;20(1):93–9.
 40.
Lee S, Miropolsky L, Wu M. SKAT: SNPset (Sequence) kernel association test. R package version 0.91. 2013. URL http://CRAN.Rproject.org/package=SKAT.
 41.
R Core Team. R: a language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2013. URL http://www.Rproject.org/.
 42.
Schaffner SF, Foo C, Gabriel S, Reich D, Daly MJ, Altshuler D. Calibrating a coalescent simulation of human genome sequence variation. Genome Res. 2005;15(11):1576–83.
 43.
Almasy L, Dyer T, Peralta J, Kent J, Charlesworth J, Curran J, et al. Genetic analysis workshop 17 miniexome simulation. BMC Proc. 2011;5 Suppl 9:S2.
 44.
The 1000 Genomes Project Consortium. A map of human genome variation from populationscale sequencing. Nature. 2010;467(7319):1061–73.
 45.
Revolution Analytics. doMC: Foreach parallel adaptor for the multicore package. R package version 1.3.3. 2014. URL http://CRAN.Rproject.org/package=doMC.
 46.
Li M, He Z, Zhang M, Zhan X, Wei C, Elston RC, et al. A generalized genetic random field method for the genetic association analysis of sequencing data. Genet Epidemiol. 2014;38(3):242–53.
 47.
Breslow NE, Lin X. Bias correction in generalised linear mixed models with a single component of dispersion. Biometrika. 1995;82(1):81–91.
 48.
Lin X, Breslow NE. Bias correction in generalized linear mixed models with multiple components of dispersion. J Am Stat Assoc. 1996;91(435):1007–16.
 49.
Bansal V, Libiger O, Torkamani A, Schork NJ. Statistical analysis strategies for association studies involving rare variants. Nat Rev Genet. 2010;11(11):773–85.
 50.
Lee S, Abecasis Gonçalo R, Boehnke M, Lin X. Rarevariant association analysis: study designs and statistical tests. Am J Hum Genet. 2014;95(1):5–23.
 51.
Moutsianas L, Morris AP. Methodology for the analysis of rare genetic variation in genomewide association and resequencing studies of complex human traits. Brief Funct Genomics. 2014. doi: 10.1093/bfgp/elu012.
 52.
Derkach A, Lawless JF, Sun L. Pooled association tests for rare genetic variants: a review and some new results. Stat Sci. 2014;29(2):302–21.
Acknowledgements
We would like to thank the reviewers for their careful reading of the paper and suggestive comments which substantially improve our original manuscript; we would also like to thank the editors and the associate editor for their supports. We gratefully thank the GAW17 committee for providing the GAW17 data. Preparation of the Genetic Analysis Workshop 17 Simulated Exome Data Set was supported in part by NIH grant R01 MH059490 and used sequencing data from the 1000 Genomes Project (www.1000genomes.org). The Genetic Analysis Workshop is supported by NIH grant R01 GM031575. Our work in this paper was financially supported by the National Natural Science Foundation of P. R. China [grant numbers: 81473070, 81373102 and 81402765], the Statistical Science Research Project from National Bureau of Statistics of P. R. China [grant number: 2014LY112], the College Philosophy and Social Science Foundation from Education Department of Jiangsu Province of P. R. China [grant numbers: 2013SJB790059, 2013SJD790032], and the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD). The work was also supported by the Cultivation Foundation of Excellent Doctoral Dissertation from the School of Public Health of Nanjing Medical University and the Foundation of Young Medical Talents Training Program of Pudong Health Bureau of Shanghai [grant number: PWRq201320].
Author information
Additional information
Competing interests
The authors declare that they have no competing interests.
Authors’ contributions
PZ and FC conceived and conducted the study, and PZ drafted the manuscript. YZ and TW participated in the design and data generation. PZ and HL conducted the simulation and interpretation of the analyses. All the authors read and approved the final manuscript.
Rights and permissions
About this article
Cite this article
Zeng, P., Zhao, Y., Li, H. et al. Permutationbased variance component test in generalized linear mixed model with application to multilocus genetic association study. BMC Med Res Methodol 15, 37 (2015). https://doi.org/10.1186/s1287401500301
Received:
Accepted:
Published:
Keywords
 Likelihood ratio test
 Permutation procedure
 Variance component
 Multilocus association analysis
 Working response
 Penalized quasilikelihood algorithm
 Generalized linear mixed model
 Logistic mixed effects model