Volume 8 Supplement 1

## Genetic Analysis Workshop 18: Human sequence data in extended pedigrees

# Penalized multivariate linear mixed model for longitudinal genome-wide association studies

- Jin Liu†
^{1}, - Jian Huang†
^{2}and - Shuangge Ma
^{3, 4}Email author

**8(Suppl 1)**:S73

https://doi.org/10.1186/1753-6561-8-S1-S73

© Liu et al.; licensee BioMed Central Ltd. 2014

**Published: **17 June 2014

## Abstract

We consider analysis of Genetic Analysis Workshop 18 data, which involves multiple longitudinal traits and dense genome-wide single-nucleotide polymorphism (SNP) markers. We use a multivariate linear mixed model to account for the covariance of random effects and multivariate residuals. We divide the SNPs into groups according to the genes they belong to and score them using weighted sum statistics. We propose a penalized approach for genetic variant selection at the gene level. The overall modeling and penalized selection method is referred to as the penalized multivariate linear mixed model. Cross-validation is used for tuning parameter selection. A resampling approach is adopted to evaluate the relative stability of the identified genes. Application to the Genetic Analysis Workshop 18 data shows that the proposed approach can effectively select markers associated with phenotypes at gene level.

## Keywords

## Background

The Genetic Analysis Workshop 18 (GAW18) data consists of multiple longitudinal traits and dense genome-wide single-nucleotide polymorphism (SNP) markers. A commonly used approach for identifying markers associated with traits is to conduct single-variant analysis and then adjust for multiple comparisons on each trait. However, for complex polygenic traits, single-variant analysis methods may not be appropriate as they fail to take into account the accumulated and/or joint effects of multiple genetic variants on the traits. In addition, analyzing each trait separately does not take into account the correlation among traits, and thus can be ineffective. To overcome these limitations, we developed a joint analysis approach referred to as the penalized multivariate linear mixed model (PMLMM). This approach takes into account covariance of both random effects and residuals and uses a group minimax concave penalty (MCP) approach [1] for variant selection at the gene level. A resampling approach is adopted to evaluate the relative stability of the identified genes. Our analysis of the GAW18 data indicates that the proposed approach can effectively select markers associated with multiple traits at the gene level.

## Methods

Consider a genetic association study with longitudinal measurements on *N* subjects, *p* genetic variants, and *q* environmental exposure covariates. Here a genetic variant can be a single SNP marker or a score representing a group of SNPs. For subject *i*, suppose that there are ${n}_{i}$ longitudinal measurements on *m* traits. Let ${Y}_{i}$ be the ${n}_{i}\times m$ trait matrix for subject *i*. Let *Y* be the $n\times m$ trait matrix for all the *N* subjects, where $n={\sum}_{i=1}^{N}{n}_{i}$. The transpose of *Y* is $Y\prime =\left({Y}_{1}^{\prime},\cdot \cdot \cdot ,{Y}_{N}^{\prime}\right)$. Let ${X}_{i}$ be the ${n}_{i}\times p$ matrix consisting of the genetic variant scores of subject *i*. Let ${Z}_{i}$ be the ${n}_{i}\times q$ covariate matrix. We center all the measurements to have sample means equal to zero. When $m=1$, this setting simplifies to that in Schelldorfer et al [2].

where *B* is a $p\times m$ matrix representing the effects of *p* genetic variants on *m* traits, and ${C}_{i}$ is a $q\times m$ matrix representing the subject specific effects of the covariates ${Z}_{i}$ for the *i*^{th} subject. We treat ${C}_{i}$ as random effects. Assume that (a) ${E}_{i}~\mathsf{\text{M}}{\mathsf{\text{N}}}_{{n}_{i}\times m}\left(0,\sum _{1},{I}_{{n}_{i}}\right)$, that is, ${E}_{i}$ is row-independent with column covariance matrix $\sum _{1}$, and each ${E}_{i}$ is independent for $i=1,...,N$; (b) ${C}_{i}~\mathsf{\text{M}}{\mathsf{\text{N}}}_{q\times m}\left(0,\sum _{2},D\right)$, where $\sum _{2}$ is the column covariance matrix and *D* is the row covariance matrix, and each ${C}_{i}$ is independent for $i=1,...,N$; (c) each ${E}_{i}$ and ${C}_{i}$ is independent; and (d) $\sum _{1}=\sum _{2}=\sum $.

Then ${Z}_{i}{C}_{i}+{E}_{i}~\mathsf{\text{M}}{\mathsf{\text{N}}}_{{n}_{i}\times m}\left(0,\sum ,{Z}_{i}D\underset{i}{\overset{\prime}{Z}}+{I}_{{n}_{i}}\right)$ and $Y~\mathsf{\text{M}}{\mathsf{\text{N}}}_{n\times m}\left(XB,\sum ,V\right)$ where $V=Diag\left({V}_{1},...,{V}_{N}\right)$ and ${V}_{i}={Z}_{i}D{Z}_{i}^{\prime}+{\mathsf{\text{I}}}_{{n}_{i}}$, where ${\mathsf{\text{I}}}_{{n}_{i}}$ is an ${n}_{i}\times {n}_{i}$ identity matrix. A more detailed description of this model can be found in Liu et al [3].

*D*by using the estimates from

*m*univariate linear mixed models and subsequently get the estimate $\widehat{V}$ of

*V*as ${\widehat{V}}_{i}={Z}_{i}\widehat{D}{Z}_{i}^{\prime}+{\mathsf{\text{I}}}_{{n}_{i}}$. Given $\widehat{\sum}$ and $\widehat{V}$, we can transform the negative likelihood function into a weighted least squares criterion for estimating

*B*, which is $\mathsf{\text{tr}}\left({{\displaystyle \widehat{\sum}}}^{-1}\left(y\mathsf{\text{-}}\phantom{\rule{0.3em}{0ex}}XB\right)\prime {\widehat{V}}^{-1}\left(y\mathsf{\text{-}}\phantom{\rule{0.3em}{0ex}}XB\right)\right)$. For variant selection, we adopt the group MCP approach [6]. The overall penalized objective function is

where ${B}_{j}$ is the *j*^{th} row of B and $\rho \left(t;\lambda ,\gamma \right)=\lambda {\int}_{0}^{\left|t\right|}{\left(1-x/\left(\gamma \lambda \right)\right)}_{+}dx$ is the MCP with tuning parameter $\lambda $ and regularization parameter $\gamma $[7].

For computation, we use a group coordinate descent algorithm [1]. The group MCP involves a regularization parameter and a tuning parameter. Generally speaking, smaller values of $\gamma $ are better at retaining the unbiasedness of the MCP penalty for large coefficients, but they have the risk of creating objective functions that have problems with nonconvexity [8], are difficult to optimize, and yield solutions that are discontinuous with respect to $\lambda $. Simulation studies by Breheny and Huang [8] suggest that $\gamma =6$ is a reasonable choice. Therefore, we fix it to be 6 in our analyses. We search for the optimal value of $\lambda $ using 5-fold cross-validation.

## Results

The GAW18 data set consists of dense genome-wide markers with longitudinal measurements on systolic and diastolic blood pressure (SBP and DBP) and other covariates. Other measurements include gender, age, year of examination, use of antihypertensive medications, and tobacco smoking at up to 4 time points. In this study, we analyze the 157 unrelated individuals using SBP and DBP as traits and other medical and demographic covariates as random effects. Gene annotations for SNP data are obtained from http://www.scandb.org. SNPs in each gene are scored using weighted sum statistics to generate gene-level measurements [9]. After quality control, we have the genetic scores of 10,400 genes for further analysis. SBP, DBP, and genetic scores are standardized to have zero means and unit variances. This procedure removes the estimation of intercepts and makes the genes comparable.

Genes identified by PMLMM: estimates for SBP and DBP, and OOI

Gene | SBP | DBP | OOI | Gene | SBP | DBP | OOI |
---|---|---|---|---|---|---|---|

| 0.002 | −0.002 | 0.333 |
| 0.027 | 0.033 | 0.403 |

| 0.085 | 0.060 | 0.697 |
| 0.024 | 0.006 | 0.247 |

| 0.071 | −0.032 | 0.323 |
| 0.025 | −0.007 | 0.507 |

| 0.018 | 0.032 | 0.563 |
| −0.008 | −0.031 | 0.540 |

| −0.028 | −0.022 | 0.333 |
| −0.042 | −0.007 | 0.623 |

| 0.002 | 0.006 | 0.307 |
| −0.068 | −0.032 | 0.647 |

| 0.006 | 0.010 | 0.373 |
| −0.041 | 0.048 | 0.693 |

| −0.004 | 0.001 | 0.337 |
| −0.004 | −0.005 | 0.297 |

| 2E-04 | 0.045 | 0.490 |
| −0.027 | 0.006 | 0.573 |

| −0.033 | −0.030 | 0.627 |
| 0.023 | 0.028 | 0.353 |

| 0.005 | 0.011 | 0.217 |
| 0.098 | −0.012 | 0.880 |

| −0.013 | −0.028 | 0.493 |
| 0.047 | −0.014 | 0.490 |

| 0.017 | 0.001 | 0.450 |
| 0.026 | 0.019 | 0.283 |

| 0.034 | 0.008 | 0.627 |
| 0.024 | 0.013 | 0.250 |

| −0.014 | −0.014 | 0.417 |
| 0.002 | 0.002 | 0.270 |

| 0.003 | −0.052 | 0.677 |
| 0.015 | −0.004 | 0.170 |

| 0.005 | 0.007 | 0.290 |
| 0.002 | 0.021 | 0.657 |

| 0.019 | −0.029 | 0.560 |
| 0.019 | −0.005 | 0.210 |

| 0.001 | −0.024 | 0.553 |
| 0.003 | −0.001 | 0.187 |

| −0.069 | 0.006 | 0.627 |
| −0.056 | −0.012 | 0.377 |

| 0.123 | 0.038 | 0.777 |
| 0.013 | 0.048 | 0.577 |

| 0.035 | −0.014 | 0.660 |
| −0.001 | −0.063 | 0.550 |

| −0.034 | −0.041 | 0.453 |
| 4E-04 | 0.015 | 0.333 |

| −2E-04 | −3E-04 | 0.197 |
| 0.024 | −0.034 | 0.603 |

| 0.029 | −0.001 | 0.630 |
| 0.018 | 0.005 | 0.127 |

| 0.013 | 0.005 | 0.333 |
| 0.048 | −0.040 | 0.547 |

| −0.015 | −0.040 | 0.480 |
| 0.028 | 0.082 | 0.720 |

| 0.006 | −0.004 | 0.437 |
| 0.020 | −0.014 | 0.233 |

| 0.004 | 0.024 | 0.467 |
| 0.004 | 2E-04 | 0.327 |

| 0.004 | −1E-04 | 0.300 |
| 0.027 | 0.003 | 0.330 |

| 0.057 | −0.008 | 0.710 |
| −0.116 | −0.074 | 0.940 |

| 0.023 | −0.017 | 0.490 |
| 0.002 | 0.003 | 0.233 |

Overlap of selected genes between PMLMM and PLMM

PMLMM | PLMM* | PLMM† | |
---|---|---|---|

PMLMM | 64 | 24 | 16 |

PLMM | 40 | 0 | |

PLMM | 29 |

## Discussion

In this study, we analyze the GAW18 data and develop a PMLMM approach. A multivariate linear mixed model is used to model variance components among traits and longitudinal measurements. A penalization approach is adopted for variant selection. In the estimation procedure, it can be considered heuristic to use $\widehat{\sum}$ and $\widehat{V}$ as proposed. Assumptions (a) to (c) are standard in mixed models, but the assumption that $\sum _{1}=\sum _{2}$ may be restrictive. Because our study is to identify multitrait-associated markers at the gene level, the restriction on variance components does not affect the selection result significantly. We are currently developing a similar approach to update variance components with more relaxed assumptions on $\sum _{1}$ and $\sum _{2}$. An iterative algorithm can be implemented to solve for *B*, $\sum $, and *V*. In variant selection, our method is designed to search for genes associated with all the traits considered. When different sets of genetic variants are suspected to be associated with different phenotypes, the sparse group penalization approach [11] can be applied.

## Conclusions

We have presented a penalized multivariate linear mixed model (PMLMM) for detecting pleiotropic genetic associations among multiple traits in the presence of pedigree structures. The proposed approach combines the advantages of mixed models that allow for elegant correction for pedigree-based family data and integrative analysis of multiple traits. Compared with PLMM which considers one trait at a time, the proposed PMLMM can achieve better performance when the pleiotropic effect is appropriately modeled.

## Notes

## Declarations

### Acknowledgements

This study was supported by NIH grants CA142774, CA165923, and CA120988, the VA Cooperative Studies Program of the Department of Veterans Affairs, Office of Research and Development, and 2012LD001 from National Bureau of Statistics of China. The GAW18 whole genome sequence data were provided by the T2D-GENES Consortium, which is supported by NIH grants U01 DK085524, U01 DK085584, U01 DK085501, U01 DK085526, and U01 DK085545. The other genetic and phenotypic data for GAW18 were provided by the San Antonio Family Heart Study and San Antonio Family Diabetes/Gallbladder Study, which are supported by NIH grants P01 HL045222, R01 DK047482, and R01 DK053889. The Genetic Analysis Workshop is supported by NIH grant R01 GM031575.

This article has been published as part of *BMC Proceedings* Volume 8 Supplement 1, 2014: Genetic Analysis Workshop 18. The full contents of the supplement are available online at http://www.biomedcentral.com/bmcproc/supplements/8/S1. Publication charges for this supplement were funded by the Texas Biomedical Research Institute.

## Authors’ Affiliations

## References

- Huang J, Wei F, Ma S: Semiparametric regression pursuit. Stat Sin. 2012, 22: 1403-1426.PubMed CentralPubMedGoogle Scholar
- Schelldorfer J, van de Geer S: Estimation for high-dimensional linear mixed-effects models using L1-penalization. Scand Stat Theory Appl. 2011, 38 (2): 197-214. 10.1111/j.1467-9469.2011.00740.x.View ArticleGoogle Scholar
- Liu J, Yang C, Shi X, Zhao H, Huang J, Ma S: A penalized multiple-trait mixed model for association mapping with population structure correction. Technical Report (arXiv preprint arXiv:1305.4413). 2013Google Scholar
- Dawid A: Some matrix-variate distribution theory: notational considerations and a Bayesian application. Biometrika. 1981, 68: 265-274. 10.1093/biomet/68.1.265.View ArticleGoogle Scholar
- Hastie T, Tibshirani R, Friedman J: The Elements of Statistical Learning: Data Mining, Inference, and Prediction. 2009, New York, Springer-Verlag, 2View ArticleGoogle Scholar
- Liu J, Huang J, Ma S: Analysis of genome-wide association studies with multiple outcomes using penalization. PLoS One. 2012, 7: e51198-10.1371/journal.pone.0051198.PubMed CentralView ArticlePubMedGoogle Scholar
- Zhang CH: Nearly unbiased variable selection under minimax concave penalty. Ann Stat. 2010, 38: 894-942. 10.1214/09-AOS729.View ArticleGoogle Scholar
- Breheny P, Huang J: Coordinate descent algorithms for nonconvex penalized regression, with applications to biological feature selection. Ann Appl Stat. 2011, 5: 232-253. 10.1214/10-AOAS388.PubMed CentralView ArticlePubMedGoogle Scholar
- Madsen B, Browning S: A groupwise association test for rare mutations using a weighted sum statistic. PLoS Genet. 2009, 5: e1000384-10.1371/journal.pgen.1000384.PubMed CentralView ArticlePubMedGoogle Scholar
- Huang J, Ma S: Variable selection in the accelerated failure time model via the bridge method. Lifetime Data Anal. 2010, 16: 176-195. 10.1007/s10985-009-9144-2.PubMed CentralView ArticlePubMedGoogle Scholar
- Liu J, Huang J, Xie Y, Ma S: Sparse group penalized integrative analysis of multiple cancer prognosis datasets. Genet Res (Camb). 2013, 95 (2-3): 68-77. 10.1017/S0016672313000086.View ArticleGoogle Scholar

## Copyright

This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.