Likelihood Ratio Tests in Behavioral Genetics: Problems and Solutions

The likelihood ratio test of nested models for family data plays an important role in the assessment of genetic and environmental influences on the variation in traits. The test is routinely based on the assumption that the test statistic follows a chi-square distribution under the null, with the number of restricted parameters as degrees of freedom. However, tests of variance components constrained to be non-negative correspond to tests of parameters on the boundary of the parameter space. In this situation the standard test procedure provides too large p-values and the use of the Akaike Information Criterion (AIC) or the Bayesian Information Criterion (BIC) for model selection is problematic. Focusing on the classical ACE twin model for univariate traits, we adapt existing theory to show that the asymptotic distribution for the likelihood ratio statistic is a mixture of chi-square distributions, and we derive the mixing probabilities. We conclude that when testing the AE or the CE model against the ACE model, the p-values obtained from using the χ²(1 df) as the reference distribution should be halved. When the E model is tested against the ACE model, a mixture of χ²(0 df), χ²(1 df) and χ²(2 df) should be used as the reference distribution, and we provide a simple formula to compute the mixing probabilities. Similar results for tests of the AE, DE and E models against the ADE model are also derived. Failing to use the appropriate reference distribution can lead to invalid conclusions.

APPENDIX A

The test of the E model against the ACE model corresponds to the test of two parameters on the boundary of the parameter space. The asymptotic LRT distribution for such a test has been shown to be a mixture of χ²(0 df), χ²(1 df) and χ²(2 df), with mixing probabilities $({1/2}-p)$, ${1/2}$, and p (Self and Liang, 1987). The mixing probability p for testing H _E against H _ACE is obtained from

$$ p = {1\over 2\pi}\hbox{arccos}\left({{I^{\rm AC}}\over{\sqrt{I^{\rm AA}I^{\rm CC}}}}\right), $$

(A.1)

where I ^AA, I ^CC and I ^AC are the components of the inverse of the asymptotic covariance matrix for the maximum likelihood estimates of λ ²_A and λ ²_C , evaluated under the null model H _E. The asymptotic covariance matrix for a parameter vector $\varvec{\theta}$ is obtained as the inverse of the expected Fisher information matrix, equal to minus one times the expectation of the Hessian (matrix of second derivatives) of the log-likelihood function. Assuming that the response vector y _i =(y _i1,y _i2) for each twin pair i follow a multivariate normal distribution, the log-likelihood function based on n _MZ MZ twin pairs and n _DZ DZ twin pairs is

$$\eqalign{ \ell({\varvec{\theta}}) &= -(n_{\rm MZ} + n_{\rm DZ}) \log(2\pi) -{{n_{\rm MZ}}\over{2}}\log \vert {\varvec{\Sigma}_{\rm MZ}} \vert - \sum_{\hbox{MZ pairs}}\left({1\over2}({\varvec{y}}_{i} - \varvec{\mu})^{\prime}\varvec{\Sigma}_{\rm MZ}^{-1}({\varvec{y}}_{i} - \varvec{\mu})\right) \cr & -{n_{DZ}\over2}\log \vert \varvec{\Sigma}_{\rm DZ} \vert - \sum_{\hbox{DZ pairs}}\left({1\over2}({\varvec{y}}_{i} - \varvec{\mu})^{\prime}\varvec{\Sigma}_{\rm DZ}^{-1}({\varvec{y}}_{i} - \varvec{\mu})\right),} $$

where $\varvec{\mu}$ is the mean vector and $\varvec{\Sigma}_{\rm MZ}$ and $\varvec{\Sigma}_{\rm DZ}$ are the covariance matrices for MZ and DZ twins given by

$$ {\varvec{\Sigma}_{\rm MZ}} = \left( \begin{array}{ll} \lambda_{\rm A}^{2} +\lambda_{\rm C}^{2} + \lambda_{\rm E}^{2} & \lambda_{\rm A}^{2} + \lambda_{\rm C}^{2} \\ \lambda_{\rm A}^{2} + \lambda_{\rm C}^{2} &\lambda_{\rm A}^{2} + \lambda_{\rm C}^{2} + \lambda_{\rm E}^{2} \end{array} \right) \hbox{ and } \varvec{\Sigma}_{\rm DZ} = \left( \begin{array}{ll} \lambda_{\rm A}^{2} + \lambda_{\rm C}^{2} + \lambda_{\rm E}^{2} & {1\over2}\lambda_{\rm A}^{2} + \lambda_{\rm C}^{2} \\ {1\over2}\lambda_{\rm A}^{2} + \lambda_{\rm C}^{2} &\lambda_{\rm A}^{2} + \lambda_{\rm C}^{2} + \lambda_{\rm E}^{2} \end{array} \right). $$

The components of the expected Fisher information matrix for the variance components ${\varvec{\theta}}=(\lambda_{\rm A}^{2},\lambda_{\rm C}^{2},\lambda_{\rm E}^{2})$ are given by

$$ \eqalign{ I_{hk}(\varvec{\theta}) &= - \hbox{E}\left({{\partial^2 \ell(\varvec{\theta})}\over{\partial \theta_{h} \partial \theta_{k}}}\right) \cr &= {n_{\rm MZ}\over2}tr\left(\varvec{\Sigma}_{\rm MZ}^{-1}{{\partial \varvec{\Sigma}_{\rm MZ}}\over{\partial \theta_{h}}}\varvec{\Sigma}_{\rm MZ}^{-1}{{\partial \varvec{\Sigma}_{\rm MZ}}\over{\partial \theta_{k}}}\right)+ {n_{\rm DZ}\over2}tr\left(\varvec{\Sigma}_{\rm DZ}^{-1}{{\partial \varvec{\Sigma}_{\rm DZ}}\over{\partial \theta_{h}}}\varvec{\Sigma}_{\rm DZ}^{-1} {{\partial \varvec{\Sigma}_{\rm DZ}}\over{\partial \theta_{k}}}\right),} $$

(A.2)

where the second equality is obtained from general results for differentiation of matrices (Searle et al., 1992). Using Σ_MZ and Σ_DZ given above, the expected Fisher information matrix obtained from expression (A.2) is evaluated under the E model, which corresponds to the parameter vector ${\varvec{\theta}_{0}}=(0,0,\lambda_{\rm E}^{2})$. The inverse of this matrix, expressed in terms of the ratio $r={n_{\rm MZ}\over n_{\rm DZ}}$, is

$$ I^{-1}({\varvec{\theta}_{0}}) = {4\lambda_{\rm E}^{4} \over n_{\rm MZ}} \left( \begin{array}{lll} r+1 & -r - 1/2 & -1/2 \\ -r - 1/2 & r + 1/4 & 1/4\\ -1/2 & 1/4 & {2r + 1 \over 4r + 4} \end{array} \right). $$

The components I ^AA, I ^CC and I ^AC are obtained as the elements of the inverse of the 2×2 submatrix of $I^{-1}(\varvec{\theta}_{0})$ corresponding to the asymptotic covariance matrix of the maximum likelihood estimates of λ ²_A and λ ²_C ,

$$ \eqalign{ I^{\rm AA} &= {{n_{\rm DZ}} \over {\lambda_{\rm E}^{4}}}\left(r + {1\over 4}\right) \cr I^{\rm CC} &= {{n_{\rm DZ}} \over{\lambda_{\rm E}^{4}}}(r + 1) \cr I^{\rm AC} &= {{n_{\rm DZ}} \over{\lambda_{\rm E}^{4}}}\left(r + {1 \over 2}\right).} $$

Inserting these components into equation (A.1) finally gives the mixing probability

$$ p = {1 \over 2\pi}\hbox{arccos}\left({{r + {1\over2}}\over{\sqrt{(r + {1 \over 4})(r + 1)}}}\right). $$

The mixture probabilities for the likelihood ratio test of the E model against the ADE model can be derived similarly. The only difference is that the covariance matrices for MZ and DZ twins are

$${\varvec{\Sigma}_{\rm MZ}} = \left( \begin{array}{ll} \lambda_{\rm A}^{2} + \lambda_{\rm D}^{2} + \lambda_{\rm E}^{2} & \lambda_{\rm A}^{2} + \lambda_{\rm D}^{2} \\ \lambda_{\rm A}^{2} + \lambda_{\rm D}^{2} &\lambda_{\rm A}^{2} + \lambda_{\rm D}^{2} + \lambda_{\rm E}^{2} \end{array} \right) \hbox{ and } {\varvec{\Sigma}_{\rm DZ}} = \left( \begin{array}{ll} \lambda_{\rm A}^{2} + \lambda_{\rm D}^{2} + \lambda_{\rm E}^{2} & {1\over2}\lambda_{\rm A}^{2} + {1\over4}\lambda_{\rm D}^{2} \\ {1\over2}\lambda_{\rm A}^{2} +{1\over4}\lambda_{\rm D}^{2} & \lambda_{\rm A}^{2} + \lambda_{\rm D}^{2} + \lambda_{\rm E}^{2} \end{array} \right). $$

The components of the expected Fisher information matrix for the variance components ${\varvec{\theta}}=(\lambda_{\rm A}^{2},\lambda_{\rm D}^{2},\lambda_{\rm E}^{2})$ are given by formula (A.2), and the inverse of this matrix evaluated at ${\varvec{\theta}}_{0}=(0,0,\lambda_{\rm E}^{2})$ is

$$ I^{-1}({\varvec{\theta}_{0})} = {{16\lambda_{\rm E}^{4}} \over{n_{MZ}}} \left( \begin{array}{lll} r+1/16 & -r - 1/8 & 1/16 \\ -r - 1/8 & r + 1/4 &-1/8 \\ 1/16 & -1/8 & {2r +1 \over 16r + 16} \end{array} \right). $$

This gives

$$\eqalign{ I^{\rm AA} &= {{n_{\rm DZ}}\over{\lambda_{\rm E}^{4}}}(r + 1/4) \cr I^{\rm DD} &= {{n_{\rm DZ}}\over{\lambda_{\rm E}^{4}}}(r + 1/16) \cr I^{\rm AD} &= {{n_{\rm DZ}}\over{\lambda_{\rm E}^{4}}}(r + 1/8),} $$

and the mixing probability for the test of the E model against the ADE model becomes

$$ p^{*} = {1\over 2\pi}\hbox{arccos}\left({{I^{\rm AD}}\over{\sqrt{I^{\rm AA}I^{\rm DD}}}}\right) = {1\over 2\pi}\hbox{arccos}\left({{r + {1\over8}}\over{\sqrt{(r + {1\over4})(r + {1\over16})}}}\right). $$