Functional principal variance component testing for a genetic association study of HIV progression

TL;DR

This paper introduces a novel functional principal variance component testing framework to analyze the association between genetic markers and HIV progression, capturing nonlinear disease dynamics more effectively than traditional models.

Contribution

The study develops and validates a new FPVC testing method that improves power and computational efficiency in large-scale genetic association studies of HIV progression.

Findings

FPVC testing shows higher power than linear mixed effects models.

The method effectively captures nonlinear disease progression patterns.

Simulation results demonstrate large power gains in detecting genetic associations.

Abstract

HIV-1C is the most prevalent subtype of HIV-1 and accounts for over half of HIV-1 infections worldwide. Host genetic influence of HIV infection has been previously studied in HIV-1B, but little attention has been paid to the more prevalent subtype C. To understand the role of host genetics in HIV-1C disease progression, we perform a study to assess the association between longitudinally collected measures of disease and more than 100,000 genetic markers located on chromosome 6. The most common approach to analyzing longitudinal data in this context is linear mixed effects models, which may be overly simplistic in this case. On the other hand, existing non-parametric methods may suffer from low power due to high degrees of freedom (DF) and may be computationally infeasible at the large scale. We propose a functional principal variance component (FPVC) testing framework which captures the nonlinearity in the CD4 and viral load with potentially low DF and is fast enough to carry out thousands or millions of times. The FPVC testing unfolds in two stages. In the first stage, we summarize the markers of disease progression according to their major patterns of variation via functional principal components analysis (FPCA). In the second stage, we employ a simple working model and variance component testing to examine the association between the summaries of disease progression and a set of single nucleotide polymorphisms. We supplement this analysis with simulation results which indicate that FPVC testing can offer large power gains over the standard linear mixed effects model.

Equations

$$H_0:{\mathbf{y}}_i\perp {\mathbf{z}}_{i\mathcal{S}}\mid {\mathbf{x}}_i$$

$$y_{ir}=Y_i(t_{ir})+{ϵ}_{ir},$$

$$y_{ir}=\mu (t_{ir})+\sum _{k=1}^{\mathcal{K}}{\xi }_{ik}{\varphi }_k(t_{ir})+{ϵ}_{ir}.$$

$$H_0:\{{\xi }_{ik}{\}}_{k=1}^{\mathcal{K}}\perp {\mathbf{z}}_{i\mathcal{S}}\mid {\mathbf{x}}_i.$$

$$Y_i(t)\approx \mu (t)+\sum _{k=1}^K{\xi }_{ik}{\varphi }_k(t),$$

$${\tilde{\xi }}_{ik}={\lambda }_k{\mathbf{\varphi }}_{ik}^{\mathsf{T}}{\Sigma }_{{\mathbf{y}}_i}^{-1}({\mathbf{y}}_i-{\mathbf{\mu }}_i)$$

$$H_0^{†}:\{{\tilde{\xi }}_{ik}{\}}_{k=1}^K\perp {\mathbf{z}}_{i\mathcal{S}}\mid {\mathbf{x}}_i.$$

$$z_{ij}^{\ast }=z_{ij}-{\mu }_{z_j}({\mathbf{x}}_i).$$

$$Q_0={\Vert n^{-\frac{1}{2}}\sum _{i=1}^n{\tilde{\mathbf{\xi }}}_i{\mathbf{z}}_{i\mathcal{S}}^{\ast \mathsf{T}}\Vert }_F^2.$$

$${\hat{\xi }}_{ik}={\hat{\lambda }}_k{\hat{\mathbf{\varphi }}}_{ik}^{\mathsf{T}}{\hat{\Sigma }}_{{\mathbf{y}}_i}^{-1}({\mathbf{y}}_i-{\hat{\mathbf{\mu }}}_i)$$

$${\mu }_{z_j}(\mathbf{x})=g_j({\mathbf{\theta }}_j,\mathbf{x})$$

$$Q=\frac{1}{n}\sum _{j\in \mathcal{S}}\sum _{k=1}^K{\left(\sum _{i=1}^n{\hat{\xi }}_{ik}{\hat{z}}_{ij}^{\ast }\right)}^2={\Vert n^{-\frac{1}{2}}\sum _{i=1}^n{\hat{\mathbf{\xi }}}_i{\hat{\mathbf{z}}}_{i\mathcal{S}}^{\ast \mathsf{T}}\Vert }_F^2.$$

$$y_{ir}=\mu (t_{ir})+\sum _{k=1}^K{\xi }_{ik}{\varphi }_k(t_{ir})+{ϵ}_{ir},$$

$${\mathbf{\xi }}_i=({\xi }_{i1},...,{\xi }_{iK}{)}^{\mathsf{T}}\sim N(\mathbf{B}{\mathbf{z}}_{i\mathcal{S}}^{\ast },\Lambda ),{ϵ}_{ir}\sim N(0,{\sigma }^2)$$

$$H_0:\eta =0.$$

$${\mathbf{y}}_{\mu i}=\sum _{k=1}^K\left(\sum _{j\in \mathcal{S}}\eta {\nu }_{kj}{z}_{ij}^{\ast }+e_{ik}\right){\mathbf{\varphi }}_{ik}+{\mathbf{ϵ}}_i$$

$${\mathbf{y}}_{\mu i}|\mathbf{\nu },\{{\mathbf{z}}_{i\mathcal{S}}^{\ast }{\}}_{i=1}^n\sim N\left(\sum _{j\in \mathcal{S}}\sum _{k=1}^K\eta {\nu }_{kj}{z}_{ij}^{\ast }{\mathbf{\varphi }}_{ik},{\Sigma }_{{\mathbf{y}}_i}\right)$$

$$\log \mathcal{L}(\eta )=-\frac{1}{2}\sum _{i=1}^n\left\{\log \left|{\Sigma }_{{\mathbf{y}}_i}\right|+{\left({\mathbf{y}}_{\mu i}-\eta \sum _{j\in \mathcal{S}}\sum _{k=1}^K{\nu }_{kj}{z}_{ij}^{\ast }{\mathbf{\varphi }}_{ik}\right)}^{\mathsf{T}}{\Sigma }_{{\mathbf{y}}_i}^{-1}\left({\mathbf{y}}_{\mu i}-\eta \sum _{j\in \mathcal{S}}\sum _{k=1}^K{\nu }_{kj}{z}_{ij}^{\ast }{\mathbf{\varphi }}_{ik}\right)\right\}.$$

$$E{\left\{{\frac{\partial \log \mathcal{L}(\eta )}{\partial \eta }|}_{\eta =0}\mid \mathbb{V}\right\}}^2+E\left\{{\frac{{\partial }^2\log \mathcal{L}(\eta )}{\partial {\eta }^2}|}_{\eta =0}\mid \mathbb{V}\right\}$$

$$=E{\left(\sum _{i=1}^n{\mathbf{y}}_{\mu i}^{\mathsf{T}}{\Sigma }_{{\mathbf{y}}_i}^{-1}\sum _{j\in \mathcal{S}}\sum _{k=1}^K{\nu }_{kj}{z}_{ij}^{\ast }{\mathbf{\varphi }}_{ik}\mid \mathbb{V}\right)}^2-E\left(\sum _{i=1}^n\sum _{j,j^{\prime }\in \mathcal{S}}\sum _{k,k^{\prime }=1}^K{\nu }_{kj}{z}_{ij}^{\ast }{\mathbf{\varphi }}_{ik}^{\mathsf{T}}{\Sigma }_{{\mathbf{y}}_i}^{-1}{\nu }_{k^{\prime }{j}^{\prime }}{z}_{i{j}^{\prime }}^{\ast }{\mathbf{\varphi }}_{i{k}^{\prime }}\mid \mathbb{V}\right)$$

$$=\sum _{j\in \mathcal{S}}\sum _{k=1}^K{\left(\sum _{i=1}^n{\mathbf{y}}_{\mu i}^{\mathsf{T}}{\Sigma }_{{\mathbf{y}}_i}^{-1}{\mathbf{\varphi }}_{ik}{\lambda }_k{z}_{ij}^{\ast }\right)}^2-\sum _{j\in \mathcal{S}}\sum _{k=1}^K\left\{\sum _{i=1}^n({\lambda }_k{z}_{ij}^{\ast }{)}^2{\mathbf{\varphi }}_{ik}^{\mathsf{T}}{\Sigma }_{{\mathbf{y}}_i}^{-1}{\mathbf{\varphi }}_{ik}\right\}$$

$$n^{-1}\sum _{j\in \mathcal{S}}\sum _{k=1}^K{\left(\sum _{i=1}^n{\mathbf{y}}_{\mu i}^{\mathsf{T}}{\Sigma }_{{\mathbf{y}}_i}^{-1}{\mathbf{\varphi }}_{ik}{\lambda }_k{z}_{ij}^{\ast }\right)}^2=n^{-1}\sum _{j\in \mathcal{S}}\sum _{k=1}^K{\left(\sum _{i=1}^n{\tilde{\xi }}_{ik}{z}_{ij}^{\ast }\right)}^2=Q_0,$$

$$q_{kj}=n^{-\frac{1}{2}}\sum _{i=1}^n{\hat{\xi }}_{ik}{\hat{z}}_{ij}^{\ast }$$

$${\tilde{q}}_{kj}=n^{-\frac{1}{2}}\sum _{i=1}^n{\tilde{\xi }}_{ik}{\hat{z}}_{ij}^{\ast }$$

$$n^{\frac{1}{2}}({\hat{\mathbf{\theta }}}_j-{\mathbf{\theta }}_j)=n^{-\frac{1}{2}}\sum _{i=1}^n{\mathcal{U}}_j({\mathbf{x}}_i)z_{ij}^{\ast }+o_p(1),$$

$$q_{kj}={\tilde{q}}_{kj}+o_p(1)=n^{-\frac{1}{2}}\sum _{i=1}^n{\mathcal{Q}}_{ikj}+o_p(1)$$

$$\bar{Q}=n^{-1}\sum _{j\in \mathcal{S}}\sum _{k=1}^K{\left[\sum _{i=1}^n{\bar{\xi }}_{ik}{\hat{z}}_{ij}^{\ast }\right]}^2$$

$$q_{kj}^{†}=n^{-\frac{1}{2}}\sum _{i=1}^n{\bar{\xi }}_{ik}{\hat{z}}_{ij}^{\ast }={\tilde{q}}_{kj}+o_p(1).$$

$$H_0:\{{\mathbf{y}}^{(m)}{\}}_{m=1}^M\perp {\mathbf{z}}_{\mathcal{S}}\mid \mathbf{x}$$

$$Q=\Vert n^{-\frac{1}{2}}\sum _{i=1}^n{\hat{\mathbf{\xi }}}_i{\hat{\mathbf{z}}}_{i\mathcal{S}}^{\ast \mathsf{T}}{\Vert }_F^2=\sum _{m=1}^M\Vert n^{-\frac{1}{2}}\sum _{i=1}^n{\hat{\mathbf{\xi }}}_i^{(m)}{\hat{\mathbf{z}}}_{i\mathcal{S}}^{\ast \mathsf{T}}{\Vert }_F^2.$$

$$y_{ir}^{(m)}$$

$$=\sin (t_{ir})+(-1{)}^{m-1}\gamma \left\{\sin (t_{ir}/3)+\cos (t_{ir})\right\}+(1-\gamma )\left\{b_{i0}+0.5b_{i1}^{(m)}\cos (t_{ir}/4)\right\}+$$

Full text

60K35,

60K35,

60K35,

Genomic association studies,

HIV disease progression,

Functional principal component analysis,

Longitudinal data,

Mixed effects models,

Variance component testing,

1 Introduction

An important goal of large-scale genomic association studies is to explore susceptibility to complex diseases. These studies have led to identification of many genomic regions as putatively harboring disease susceptibility alleles for a wide range of disorders. For patients with a particular disease, association studies have also been performed to identify genetic variants associated with progression of disease. The disease progression is often monitored by longitudinally measured biological markers. Such longitudinal measures allow researchers to more clearly characterize clinical outcomes that cannot necessarily be captured in one or even a few measurements.

We are motivated by a large-scale association study of HIV-1 Subtype-C (HIV-1C) progression in sub-Saharan African individuals. HIV-1C is the most prevalent subtype of HIV-1 and accounts for over half of HIV-1 infections worldwide [11]. Sub-Saharan Africa, where HIV-1C dominates, was home to an estimated 69% of people living with HIV in 2012 [29]. While a number of human leukoucyte antigen (HLA) alleles [10, 30, 19, e.g.] and other loci have been identified to be associated with AIDS progression in European males infected by HIV-1B [10, 22], comparatively little research has focused on host genetic influence in this African population and subtype.

In this study, we seek to relate the longitudinal progression of these two markers – log10 CD4 (lCD4) count and log10 viral load (lVL) – to a set of approximately 100,000 single nucleotide polymorphisms (SNPs) located on chromosome 6 in two independent cohorts of treatment-naive individuals in Botswana. We focus on chromosome 6 because it houses the HLA region of genes, which are known to impact immune function. Throughout this paper, we will use ${\bf y}=(y_{1},...,y_{r})^{{\sf\scriptscriptstyle{T}}}$ to denote the longitudinal outcome (in the context of this study, either lCD4 or lVL), measured at times ${\bf t}=(t_{1},...,t_{r})^{{\sf\scriptscriptstyle{T}}}$. We will furthermore let a set of genetic markers of interest be denoted ${\bf z}$ and any potential covariates be ${\bf x}$.

The most common approach to analyzing longitudinal data of this kind is to use linear mixed effects (LME) models [16], which relate ${\bf y}$ linearly to ${\bf z}$, ${\bf x}$, and ${\bf t}$ with both fixed and random effects. However, both in the case of HIV progression measured by lCD4 and lVL and in general, such a linear relationship is likely to be overly simplistic. To incorporate nonlinear trajectories, many generalizations of LMEs have been proposed. Typically, methods assume that ${\bf y}$ is a noisy realization of a smooth underlying function $Y(\cdot)$. These include nonlinear or nonparametric LMEs using a fixed spline basis expansion [18, 25, 12] for ${\bf t}$ and adjusting for ${\bf z}$. Functional regression methods have been proposed for densely sampled trajectories [6] for a single $z$ as well as for irreguarly spaced longitudinal data where ${\bf z}$ is also allowed to change over time [33]. In these methods, estimation proceeds via kernel smoothing for both the population mean function and the covariance process of $Y(\cdot)$. [15] proposed an iterative procedure for fitting functional regression models which accounts for within-subject covariance but does not estimate random effects. Similarly, [23] proposed a ridge-based estimator for the case when ${\bf y}$ is measured on a common, fine grid of points for all subjects, but requires $z$ to be univariate and also ignores random effects. In a similarly dense setting [20] proposed wavelet-based mixed effects models, and inference procedures for the random effects were developed in [2].

While some of these methods can be adapted to test for association, none of them are appealing for our study for the following reasons. First, some require restrictive assumptions about the density of measurements [20, e,g] which are clearly not met here. Further, all of these methods were developed with estimation and regression in mind. While many of them could in principle be used to derive testing procedures, the validity of their inference procedures often rely on the model assumptions on the distribution of ${\bf y}$ given the SNPs and the covariates. Under model mis-specification, the resulting test may fail to maintain type I error. Additionally, these procedures would require fitting a complex iterative or smoothing-based model thousands or millions of times for genome-wide studies and hence become computationally infeasible. We, on the other hand, propose a novel testing procedure which is valid regardless of the true distribution of ${\bf y}$, is fast to fit and has a simple limiting distribution, despite the fact that we account for nonlinearity in ${\bf y}$ in a manner akin to previous functional regression methods. To do this, we first capture the nonlinearity in ${\bf y}$ using functional principal components analysis (FPCA) [5, 24, 32, 13, 15], Using a eigenfunction decomposition of the smoothed covariance function of $Y(\cdot)$, we approximate each patient’s $Y(\cdot)$ by a weighted average of the estimated eigenfunctions, with weights corresponding to functional principal component loadings or scores. Then borrowing a variance component test framework and using these scores as pseudo outcomes, we construct a Functional Principal Variance Component (FPVC) test that can capture the nonlinear trajectories without requiring a normality assumption or fitting individual functional regression models. The test statistics can be approximated by a mixture of chi-square and the small number of eigenfunctions needed to approximate the trajectories can result in a low degree of freedom in the resulting test statistic.

Since the data of interest is sampled at irregular time intervals, we use the best linear unbiased predictor (BLUP) to estimate the scores. The BLUP was also the basis for FPCA with sparse longitudinal data in the principal analysis via conditional expectation (PACE) method [32] under a normality assumption. Here, we use BLUP to motivate our testing procedure but do not require normality for the validity of the FPVC test. The test statistic can be derived through the variance component testing framework and viewed as a summary measure of the overall covariance between the estimated subject-specific scores, which characterize the person’s trajectory, and the genetic markers. Similar variance component tests have previously been proposed for standard linear and logistic regressions with observed single outcomes [31].

The primary virtues of FPVC testing are threefold. First, we separate the procedure into two stages of distinct complexity in order to make it feasible at large scale. In the first stage, we model ${\bf y}$ flexibly using FPCA and obtain a succinct summary of disease progression for each patient, once and for all. In the second stage, we perform a rather simple model at large scale. Thus, we segregate the computationally complex stage (which need occur only once) from the large-scale stage (which could require the same computation on the order of millions of times in, e.g., genome-wide association studies). Second, the summary of ${\bf y}$ that we obtain from FPCA is the most succinct summary possible, as the eigenfunctions identified by FPCA are the functions that explain the most variability in ${\bf y}$. Third, our theoretical results suggest that the null distribution of the FPVC test statistic reduces to a simple mixture of $\chi^{2}$ distributions. The variability due to estimating the eigenfunctions does not contribute to the null distribution of the test statistic asymptotically at the first order (see (11) and the derivation of the asymptotic null distribution).

The structure of the rest of the paper is as follows. In section 2, we describe FPCA and introduce FPVC testing and our main theoretical results. In section 3, we give details about the association study for HIV progression. In section 4, we discuss simulation results, and in section 5 we discuss further implications of our procedure.

2 Functional principal variance component testing

2.1 The test statistic

In this section, we propose a testing procedure for assessing the association between a set of genetic markers ${\bf z}$ and a longitudinally measured outcome ${\bf y}$, adjusting for covariates ${\bf x}$. Let the data for analysis consist of $n$ independent random vectors $\mathbb{V}=\left\{{\bf V}_{i}=({\bf y}_{i}^{{\sf\scriptscriptstyle{T}}},{\bf t}_{i}^{{\sf\scriptscriptstyle{T}}},{\bf z}_{i}^{{\sf\scriptscriptstyle{T}}},{\bf x}_{i}^{{\sf\scriptscriptstyle{T}}})^{{\sf\scriptscriptstyle{T}}}\right\}_{i=1}^{n}$, where ${\bf y}_{i}=(y_{i1},...,y_{ir_{i}})^{{\sf\scriptscriptstyle{T}}}$ is a vector of outcome measurements taken at times ${\bf t}_{i}=(t_{i1},...,t_{ir_{i}})^{{\sf\scriptscriptstyle{T}}}\in\mathcal{T}^{r_{i}}$, $\mathcal{T}$ is a closed and bounded interval, ${\bf z}_{i}=(z_{i1},...,z_{ip})^{{\sf\scriptscriptstyle{T}}}$ is a vector of genetic markers of interest, and ${\bf x}_{i}=(1,x_{i1},...,x_{iq})^{{\sf\scriptscriptstyle{T}}}$ is a vector of additional covariates that are potentially related to the outcome, all measured on person $i$. For each $i$, we take $({\bf z}_{i}^{{\sf\scriptscriptstyle{T}}},{\bf x}_{i}^{{\sf\scriptscriptstyle{T}}})^{{\sf\scriptscriptstyle{T}}}$ to be distributed as $({\bf z}^{{\sf\scriptscriptstyle{T}}},{\bf x}^{{\sf\scriptscriptstyle{T}}})^{{\sf\scriptscriptstyle{T}}}$.

Our goal is to test the null hypothesis

[TABLE]

where ${\bf z}_{i{\cal S}}=(z_{ij_{1}},...,z_{ij_{s}})^{{\sf\scriptscriptstyle{T}}}$ is a set of genetic factors to test, identified by the index set ${\cal S}=\{j_{1},...,j_{s}\}\subset\{1,...,p\}$. Special cases include marginal testing, as in traditional genome-wide association studies, where ${\cal S}=\{j\}$ for some $j\in\{1,...,p\}$, or set-based testing where ${\cal S}$ are the indices of the SNPs in a gene or some other related set.

To model the longitudinal trajectory, we assume that $y_{ir}$ is a noisy sample of a smooth underlying function $Y_{i}(\cdot)$, evaluated at the point $t_{ir}$,

[TABLE]

which following the logic of [32] can be written as a linear combination of its population mean $E\{Y(\cdot)\}$ and a set of ${\cal K}$ eigenfunctions $\{\phi_{k}(\cdot)\}$

[TABLE]

where $\xi_{ik}$ is the FPCA score associated with the $k$th eigenfunction, $E(\xi_{ik})=0,\text{Var}(\xi_{ik})=\lambda_{k}$, and the eigenfunctions are ordered such that the $k$th explains the $k$th most variance in $Y(\cdot)$. Thus, the relationship between ${\bf y}_{i}$ and ${\bf x}_{i}$ and ${\bf z}_{i}$ must be captured by the only random quantity in (2), the random coefficients $\{\xi_{ik}\}_{k=1}^{{\cal K}}$. Therefore, testing (1) is equivalent to testing

[TABLE]

However, direct assessment of the association between $\{\xi_{ik}\}_{k=1}^{{\cal K}}$ and ${\bf z}_{i{\cal S}}$ is difficult since $\{\xi_{ik}\}_{k=1}^{{\cal K}}$ are unobservable and ${\cal K}$ could be infinity.

Most methods of estimating a function like $Y_{i}(\cdot)$ require an explicit tuning of the smoothness of the resulting estimator. Here that corresponds to choosing a (typically small) number $K$ of eigenfunctions as an approximation

[TABLE]

where $K<\infty$ could be chosen such that the first $K$ directions capture a proportion of the variation at least as large as $\wp\in(0,1]$. Our simulation results (see section 4) suggest that the performance of FPVC testing is not very sensitive to the choice of $K$ provided that $\wp$ is close to 1. Since $\{\xi_{ik}\}_{k=1}^{K}$ are not observable or generally estimable due to sparse sampling of measurement times, we instead infer about the association between $Y_{i}(\cdot)$ and ${\bf z}_{i}$ based on the best linear unbiased predictor (BLUP),

[TABLE]

where $\boldsymbol{\phi}_{ik}=\{\phi_{k}(t_{i1}),...,\phi_{k}(t_{ir_{i}})\}^{{\sf\scriptscriptstyle{T}}},\boldsymbol{\mu}_{i}=\{\mu(t_{i1}),...,\mu(t_{ir_{i}})\}^{{\sf\scriptscriptstyle{T}}}$, and $\Sigma_{{\bf y}_{i}}=\text{Cov}({\bf y}_{i},{\bf y}_{i})$ such that $\left(\Sigma_{{\bf y}_{i}}\right)_{rl}=G(t_{ir},t_{il})+\sigma^{2}\delta_{rl}$, $\text{Cov}\{Y(s),Y(t)\}=G(s,t)$, and $\delta_{rl}=I_{\{r=l\}}$. In the PACE method of [32], ${\widetilde{\xi}}_{ik}$ was obtained as $E(\xi_{ik}|{\bf y}_{i})$ under the assumption that $\xi_{ik}$ and $\epsilon_{ir}$ are jointly normal, but we don’t require normality here. We simply take ${\widetilde{\xi}}_{ik}$ as an observable and reasonable approximation to $\xi_{ik}$ even if normality does not hold, as has been argued in [26] and [14].

Thus, we propose to test (1) by testing

[TABLE]

Taking note that the association we seek to test is conditional on ${\bf x}$, one may construct a test for $H_{0}^{{\dagger}}$ by regressing ${\widetilde{\boldsymbol{\xi}}}_{i}=({\widetilde{\xi}}_{i1},...,{\widetilde{\xi}}_{iK})^{{\sf\scriptscriptstyle{T}}}$ onto $({\bf x}_{i},{\bf z}_{i{\cal S}})$. However, this is only valid if the effect of ${\bf x}_{i}$ on $Y_{i}(\cdot)$ is captured fully based on the model relating ${\bf x}_{i}$ and ${\widetilde{\boldsymbol{\xi}}}_{i}$, which may not be true in general. To remove the effect of ${\bf x}_{i}$ without imposing a strong assumption on how ${\bf x}_{i}$ affects $Y_{i}(\cdot)$, we instead choose to model the conditional expectation of $z_{ij}$ given ${\bf x}_{i}$, $\mu_{z_{j}}({\bf x}_{i})=E(z_{ij}|{\bf x}_{i})$, and center ${\bf z}_{i{\cal S}}$ as ${\bf z}^{*}_{i{\cal S}}=(z^{*}_{ij_{1}},...,z^{*}_{ij_{s}})$ where for any $j$

[TABLE]

To form the test statistic for $H_{0}$, we propose to summarize the overall association between $Y(\cdot)$ and ${\bf z}_{\cal S}$ based on the Frobenius norm of the standardized covariance between ${\widetilde{\boldsymbol{\xi}}}_{i}$ and ${\bf z}_{i{\cal S}}^{*}$

[TABLE]

Though $Q_{0}$ takes a simple form and can be motivated naturally as an estimated covariance (and can thus be considered model-free), it can also be viewed as a variance component score test statistic similar to those considered previously for other regression models [8, 17]. Details on the derivation of the variance component score test statistic are given in section 2.2.

Both ${\widetilde{\boldsymbol{\xi}}}_{i}$ and ${\bf z}^{*}_{i{\cal S}}$ involve various nuisance parameters that remain to be estimated. First, under mild regularity conditions which are outlined in the supplementary article [1], we can use FPCA to estimate the relevant quantities via local linear smoothing as in [13] and [32]. Subsequently, we can estimate ${\widetilde{\xi}}_{ik}$ by

[TABLE]

for ${\widehat{\boldsymbol{\phi}}}_{ik}=\{{\widehat{\phi}}_{k}(t_{i1}),...,{\widehat{\phi}}_{k}(t_{ir_{i}})\}^{{\sf\scriptscriptstyle{T}}},{\widehat{\boldsymbol{\mu}}}_{i}=\{{\widehat{\mu}}(t_{i1}),...,{\widehat{\mu}}(t_{ir_{i}})\}^{{\sf\scriptscriptstyle{T}}}$, and $\left({\widehat{\Sigma}}_{{\bf y}_{i}}\right)_{rl}={\widehat{G}}(t_{ir},t_{il})+{\widehat{\sigma}}^{2}\delta_{rl}$. To estimate $\mu_{z_{j}}({\bf x}_{i})$, various approaches can be taken depending on the nature of ${\bf x}$. For example, when ${\bf x}$ is discrete, $\mu_{z_{j}}({\bf x}_{i})$ can be estimated empirically. With continuous ${\bf x}$, we may impose a parametric model with

[TABLE]

and obtain $\bar{z}_{j}({\bf x})$ as $g_{j}({\widehat{\boldsymbol{\theta}}}_{j},{\bf x})$, where ${\widehat{\boldsymbol{\theta}}}_{j}$ is an estimate of a finite-dimensional parameter $\boldsymbol{\theta}_{j}$. To take two examples that commonly come up in genomics, if $z_{j}$ takes values in $\{0,1\}$, e.g. under the dominant model, then $z_{j}$ can be modeled using logistic regression, and if $z_{j}$ takes values in $\{0,1,2\}$, then we may use a binomial generalized linear model or a proportional odds model. There are two reasons to prefer to remove the effect of ${\bf x}$ from ${\bf z}$ rather than from $\boldsymbol{\xi}$: first, it may in general be easier to specify a model for ${\bf z}$ rather than $\boldsymbol{\xi}$ because of the limited range of ${\bf z}$, and, second, this formulation facilitates asymptotic analysis without the need to derive the asymptotic distributions for the estimated FPCA scores. Finally, based on $\{{\widehat{\xi}}_{ik}\}_{k=1}^{K}$ and $\bar{z}_{j}({\bf x}_{i})$, our proposed test statistic is

[TABLE]

where ${\widehat{\bf z}}_{i{\cal S}}^{*}=({\widehat{z}}_{ij_{1}}^{*},...,{\widehat{z}}_{ij_{s}}^{*})^{{\sf\scriptscriptstyle{T}}}$ and ${\widehat{z}}_{ij}^{*}=z_{ij}-\bar{z}_{j}({\bf x}_{i})$.

2.2 Connection to mixed effects models

In this section, we demonstrate that one can arrive at the quantity (4) via a more familiar mixed effects model. Consider the model

[TABLE]

where ${\bf B}$ is a $K\times s$ matrix with $(k,j)$th entry $\beta_{kj}$ and $\Lambda=\text{diag}(\lambda_{1},...,\lambda_{K})$. We can obtain $Q_{0}$ as the variance component score test statistic for $H_{0}:{\bf B}=0$. Specifically, let $\beta_{kj}=\eta\nu_{kj}$ and we consider a working model such that $\{\nu_{kj}\}$ are independently distributed with $E(\nu_{kj})=0$ and $\text{Var}(\nu_{kj})=\lambda_{k}^{2}$. Under this working model, $H_{0}:{\bf B}=0$ is equivalent to

[TABLE]

This formulation follows the logic of variance component score tests that have been proposed previously [31] and recalls, for example, the likelihood ratio test proposed in [9]. To obtain the variance component test statistic, rewrite the model as

[TABLE]

for centered outcome ${\bf y}_{\mu i}=\{y_{i1}-\mu(t_{i1}),...,y_{ir_{i}}-\mu(t_{ir_{i}})\}^{{\sf\scriptscriptstyle{T}}}$, error vector $\boldsymbol{\epsilon}_{i}=(\epsilon_{i1},...,\epsilon_{ir_{i}})^{{\sf\scriptscriptstyle{T}}}$, and random effects ${\bf e}_{i}=(e_{i1},...,e_{iK})^{{\sf\scriptscriptstyle{T}}}\sim N(0,\Lambda)$. Then

[TABLE]

where $\Sigma_{{\bf y}_{i}}=\sum_{k=1}^{K}\lambda_{k}\boldsymbol{\phi}_{ik}\boldsymbol{\phi}_{ik}^{{\sf\scriptscriptstyle{T}}}+\sigma^{2}I_{r_{i}}$ and $I_{r_{i}}$ is the $r_{i}\times r_{i}$ identity matrix.

The log-likelihood for ${\bf y}_{\mu i}$ can then be written

[TABLE]

Because the target of inference is $\eta$, we marginalize over the nuisance parameter $\boldsymbol{\nu}$ conditional on the observed data to obtain ${\cal L}^{*}(\eta)=E\{{\cal L}(\eta)|\mathbb{V}\}$ where the expectation is taken over the distribution of $\boldsymbol{\nu}$. We follow the argument in [8] and note that the score at the null value is 0: $\lim_{\eta\to 0}\partial\log{\cal L}^{*}(\eta)/\partial\eta=E\left(\sum_{i=1}^{n}{\bf y}_{\mu i}^{{\sf\scriptscriptstyle{T}}}\Sigma_{{\bf y}_{i}}^{-1}\sum_{j\in{\cal S}}\sum_{k=1}^{K}\nu_{kj}z^{*}_{ij}\boldsymbol{\phi}_{ik}\mid\mathbb{V}\right)=0$. So we instead consider the score with respect to $\eta^{2},\lim_{\eta\to 0}\partial\log{\cal L}^{*}(\eta)/\partial(\eta^{2})$, and we show in the supplementary article [1] that it can be written

[TABLE]

up to a scaling constant.

To finally obtain $Q_{0}$, we standardize by $n^{-1}$ and drop the second term because it converges to a constant, yielding the score statistic

[TABLE]

taking note of the form of ${\widetilde{\xi}}_{ik}$ from (3). Thus, our proposed test statistic can be obtained as a variance component test under a normal mixed model framework. On the other hand, we can also view $Q_{0}$ as a simple summary of the overall covariance between the scores of the FPCA and the genetic markers. We next derive the null distribution of the FPVC test statistic without requiring the normal mixed model to hold.

2.3 Estimating the null distribution of the test statistic

In order to obtain p-values for FPVC testing, we must identify the null distribution of $Q$. To this end, we show in [1] that the key quantity in $Q$

[TABLE]

is asymptotically equivalent to

[TABLE]

under $H_{0}$, i.e. $q_{kj}-{\widetilde{q}}_{kj}=o_{p}(1)$ for each $j$ and $k$. The key idea for deriving the null distribution of $q_{kj}$ is that, since ${\widehat{z}}_{ij}^{*}$ is approximately mean 0 conditional on ${\bf x}_{i}$, the variability due to approximating ${\widetilde{\xi}}_{ik}$ by ${\widehat{\xi}}_{ik}$ does not contribute any additional noise to $q_{kj}$ (compared to ${\widetilde{q}}_{kj}$) at the first order under $H_{0}$. Thus, we can obtain the limiting distribution of $Q$ by analyzing the quantity ${\widetilde{Q}}=\sum_{j\in{\cal S}}\sum_{k=1}^{K}{\widetilde{q}}_{kj}^{2}.$

To characterize the null distribution of ${\widetilde{Q}}$, we need to account for the variability in the estimated model parameters for $\mu_{z_{j}}({\bf x}_{i})=g_{j}(\boldsymbol{\theta}_{j},{\bf x}_{i})$ in ${\widehat{z}}_{ij}^{*}$. Without loss of generality, we assume that for each $j$

[TABLE]

where $\mathcal{U}(\cdot)$ is some ($q+1$)-dimensional function of ${\bf x}_{i}$ with $E\{\mathcal{U}({\bf x}_{i})^{2}\}<\infty$. It follows that

[TABLE]

where ${\cal Q}_{ikj}=\left\{{\widetilde{\xi}}_{ik}-\mathbb{A}_{kj}\mathcal{U}({\bf x}_{i})\right\}z^{*}_{ij}$, \mathbb{A}_{kj}=E\{{\widetilde{\xi}}_{ik}\dot{{\mbox{\boldmathg}}}_{j}(\boldsymbol{\theta}_{j},{\bf x}_{i})^{{\sf\scriptscriptstyle{T}}}\}, and \dot{{\mbox{\boldmathg}}}_{j}(\boldsymbol{\theta}_{j},{\bf x}_{i})=\partial g_{j}(\boldsymbol{\theta}_{j},{\bf x}_{i})/\partial\boldsymbol{\theta}_{j}. We show in [1] that the limiting null distribution of $Q$ is a mixture of $\chi^{2}_{1}$ random variables, $Q\sim\sum_{l=1}^{sK}a_{l}\chi^{2}_{1}$, with mixing coefficients determined by the eigenvalues of the covariance matrix of $\{{\cal Q}_{ikj}\}_{j\in{\cal S},1\leq k\leq K}$. So finally we obtain a p-value for the association between the set ${\bf z}_{\cal S}$ and $Y(\cdot)$ as $P(\sum_{l=1}^{sK}{\widehat{a}}_{l}\chi^{2}_{1}>Q\mid\mathbb{V})$, where ${\widehat{a}}_{l}$ is an empirical estimate of $a_{l}$.

By a similar argument, one could construct an asymptotically equivalent test statistic by estimating ${\widetilde{\boldsymbol{\xi}}}_{i}$ in two stages. Instead of obtaining an estimator directly from FPCA via equation (5), FPCA can be used to estimate only $\mu(\cdot)$ and $\{\phi_{k}(\cdot)\}_{k=1}^{K}$. By plugging the estimated ${\widehat{\mu}}(\cdot)$ and $\{{\widehat{\phi}}_{k}(\cdot)\}_{k=1}^{K}$ into the mixed model (8), one can obtain what we will call the re-fitted test statistic

[TABLE]

where $\bar{\boldsymbol{\xi}}_{i}=(\bar{\xi}_{i1},...,\bar{\xi}_{iK})^{{\sf\scriptscriptstyle{T}}}$ is the BLUP from the model $y_{ir}-{\widehat{\mu}}(t_{ir})=\sum_{k=1}^{K}\xi_{ik}{\widehat{\phi}}_{k}(t_{ir})+\epsilon_{ir}$ with $\text{Cov}(\boldsymbol{\xi}_{i})=D$, for some unspecified positive definite matrix $D$. By the same argument above, estimation of ${\widetilde{\xi}}_{ik}$ by $\bar{\xi}_{ik}$ contributes no additional variability to the test statistic at the first order. It follows that

[TABLE]

and hence $\bar{Q}$ has the same limiting null distribution as $Q$. Not surprisingly, simulation results suggest that the performance of $\bar{Q}$ is quite similar to the performance of $Q$. This equivalence indicates that effectively our proposed testing procedure uses FPCA to estimate potentially nonlinear bases and assesses the effect of genetic markers by fitting a mixed model with these basis functions. On the other hand, the test statistics can also be viewed as a simple summary of covariances, and – since we estimate the null distribution without relying on the normality assumption required by the mixed models – our testing procedure remains valid regardless of the adequacy of the mixed model.

2.4 Combining multiple sources of outcome information

In the HIV progression study, we seek to test the overall association between SNPs and both lCD4 and lVL simultaneously because more and distinct information about HIV progression is captured in both measures than in either one alone. FPVC testing, as outlined above, can be easily adapted to perform a test for the overall association between ${\bf z}_{\cal S}$ and all outcomes of interest. To use information in multiple outcomes, $\{{\bf y}^{(m)}\}_{m=1}^{M}$, we simply perform FPCA separately on each ${\bf y}^{(m)}$ and obtain FPCA scores for each person and each outcome. Subject $i$’s scores for ${\bf y}_{i}^{(m)}$ would be ${\widehat{\boldsymbol{\xi}}}_{i}^{(m)}=({\widehat{\xi}}_{i1}^{(m)},...,{\widehat{\xi}}^{(m)}_{iK_{m}})^{{\sf\scriptscriptstyle{T}}}$, as in (5), and the full set of scores for person $i$ would be ${\widehat{\boldsymbol{\xi}}}_{i}=({\widehat{\boldsymbol{\xi}}}_{i}^{(1){\sf\scriptscriptstyle{T}}},...,{\widehat{\boldsymbol{\xi}}}_{i}^{(m){\sf\scriptscriptstyle{T}}})^{{\sf\scriptscriptstyle{T}}}$. Then we simply proceed by testing

[TABLE]

as before based on

[TABLE]

Since each outcome may be measured on a different scale, one may use scaling or weighting to allow scores from each outcome to contribute similarly to the test statistic. See section 5 for further discussion of scaling/weighting.

3 Association study for HIV progression

In this study, two independent cohorts were recruited in Botswana to detect sets of SNPs related to HIV disease progression as measured by lCD4 and lVL. The first cohort, which we will denote BHP010, was a natural history observational prospective cohort study recruited from clinics in Gaborone. This cohort included HIV-1C-infected individuals with CD4 cell counts above 400 cells per $\mu$l and not yet qualified for the Botswana highly active antiretroviral treatment (HAART) program. Patients were not enrolled if they were younger than 18, had an active AIDS-defining illness requiring the initiation of HAART, presented with an AIDS-related malignancy, or had previously been exposed to HAART during pregnancy or breast feeding.

Follow-up visits occurred at approximately three-month intervals with an additional visit one month after enrollment. VL was generally collected at six-month intervals, and most patients in this cohort do not have VL measurements after two years of follow-up. Follow-up began in 2005 and lasted for up to 255 weeks. The mean follow-up time was 41 months. At least two CD4 measurements were required for measuring disease progression, and 449 patients satisfied this criterion. Of these, the median age at baseline was 34 years old with an interquartile range (IQR) of (28, 39). There were 366 women. In 2008, 143 patients were genotyped on an Illumina LCG BeadChip, the chip used in the second cohort. After making exclusions for quality control – call rate greater than 0.99, gender matching listed gender – 137 individuals were included in the association study.

The second cohort we will denote BHP011. This cohort came from a randomized, multifactorial, double-blind placebo-controlled trial conducted between December 2004 and July 2009 [3]. The purpose of the trial was to determine the efficacy of micronutrient supplementation (supplementation of multivitamin, selenium, or both) in improving immune function in HIV-1C-infected individuals. It was composed of 878 treatment-naive patients with CD4 higher than 350 cells/$\mu$l, as well as body mass index (BMI) greater than 18 for women and 18.5 for men (calculated as weight in kilograms divided by height in meters squared), age of 18 years or older, no current AIDS-defining conditions or history of AIDS-defining conditions, and no history of endocrine or psychiatric disorders.

Patients were followed up for a maximum of 169 weeks. They returned to clinics approximately every 3 months to measure CD4 and approximately every 6 months to measure VL. The mean follow-up time was 696 days. Of these, 838 had at least two CD4 measurements, of whom 613 were women. The median age at baseline was 33 years old with an interquartile range (IQR) of (28, 39). In this cohort, 326 individuals were genotyped on Illumina LCG BeadChips, with 320 entered into the association study after quality control exclusions.

FPCA was performed on each cohort separately for both lCD4 and lVL. Three eigenfunctions were chosen for lCD4 and two for lVL in each cohort, which corresponds to 99% of proportion of variance explained for each. The form of the eigenfunctions look similar for both lCD4 and lVL in each cohort and lend themselves to reasonable interpretations. The first eigenfunction tends to serve as a mean shift or an intercept; the second eigenfunction acts something like a slope; and the third eigenfunction behaves approximately as a quadratic term. The vector of estimated ”re-fitted” scores (refer to (12)) for each individual to be used in testing can be written ${\widehat{\boldsymbol{\xi}}}_{i}=({\widehat{\xi}}_{ik}^{(m)})_{k=1,2,3;m=1,2}$ where ${\widehat{\xi}}_{ik}^{(m)}$ is the estimated score corresponding to the $k$th estimated eigenfunction of lCD4 when $m=1$ and lVL when $m=2$.

A total of 155,007 SNPs on chromosome 6 were genotyped. After requiring less than 5% missingness and at least five individuals with any minor alleles in each cohort, $p=108,665$ SNPs, ${\bf z}_{i}=(z_{i1},...,z_{ip})$ remained for association testing, where the ordering in ${\bf z}_{i}$ corresponds to position on the chromosome. The dominant model was used for analysis such that $z_{ij}=1$ if any minor alleles are present and $z_{ij}$ = 0 if none are present. Missing values were imputed as the minor allele frequency for that SNP. In order to gain power by pooling information in nearby SNPs, sets of 10 contiguous SNPs were constructed as ${\bf z}_{i{\cal S}_{j}}=(z_{ij},...,z_{ij+9})$ for $j=1,...,108,656$. Here the choice of 10 merely serves as example for illustration and sets could in principle be constructed with more SNPs, but due to the small sample size in each cohort, we kept the size of the sets modest. Tests were performed on each cohort separately, and p-values were combined using the Fisher method. The false discovery rate was controlled at 0.1 using the Benjamini-Hochberg procedure [4], which is expected to remain valid since although the moving window construction of sets induces high correlations for nearby regions, SNP sets in distant regions are not expected to be correlated [28].

Tests were adjusted for age and gender to remove any possible confounding, so that we are testing for the effect of SNP sets on disease progression conditional on age and gender. Logistic regression was used to remove their effect. The method appears to be robust to this specification, as results were not markedly changed either when no adjustment was made or by specifying a probit model (results not shown).

The manhattan plot for the 108,656 tests is given in figure 1.

In all, 126 tests passed the FDR threshold, corresponding to 4 broad regions of the chromosome. Six contiguous tests rejected in the region between positions 6,784,416 and 6,793,116 (region 1), which fall between the LY86 and BTF3P7 genes; 117 tests rejected between positions 31,022,266 and 31,080,899 (region 2), including SNPs on the HCG22 and C6orf15 genes; two tests rejected between 122,990,817 and 123,014,708 (region 3) on the PKIB gene; and one test rejected representing SNPs in the region between 162,250,522 and 162,254,546 (region 4) including SNPs on the PARK2 gene. Notably, the C6orf15 gene has been reported to be associated with susceptibility to follicular lymphoma [27], and genes in linkage disequilibrium with HCG22 and C6orf15 have demonstrated associations to total white blood cell counts [21] and multiple myeloma [7].

Furthermore, regions 1, 3, and 4 – whose -log10 p-values are depicted in figure 2 (a), (c), and (d), respectively – only have strong signals in one of the two cohorts. Region 1 demonstrates association largely in BHP010, as can be seen in the figure where the large triangles in the figure (representing the set-based p-value in BHP010) tend to lie above the large dots (representing the combined set-based p-value), while the diamonds for BHP011 tend to be very low. Conversely, in regions 3 and 4 the association is apparent only in BHP011. Whereas in region 2, associations tend to be strong in both cohorts, and the combined p-values tend to be lower (higher in the figure) than either of the component p-values.

To get a better sense of what the test is picking up on, we looked at average disease progression within groups of patients with similar minor allele burden. To do this, we identified the SNP set with the smallest p-value, which lay in the HCG22 portion of region 2 and included the following SNPs: rs2535308, rs2535307, rs2535306, rs2535305, rs3130955, rs2535304, rs12527394, rs2535303, kgp9442190, and rs3130959. Patients were grouped according to the number of loci among these 10 at which they had any minor alleles. Within these groups, we averaged the estimated mean lCD4 and lVL in each cohort over time ${\widehat{Y}}(\cdot)={\widehat{\mu}}(\cdot)+\sum_{k=1}^{K}{\widehat{\xi}}_{ik}{\widehat{\phi}}_{k}(\cdot)$. As a demonstration, the results for those with 2, 3, 8, and 9 loci with minor alleles are depicted in figure 3. We selected these groups to demonstrate burden extremes (very few individuals had 0 or 1 loci affected, so they were not shown).

Those with only 2 and 3 affected loci looked to generally be the healthiest group, having higher and more stable CD4 and lower and gently increasing VL throughout the study period in each cohort. Those with 3 affected loci tended to have a more negative CD4 slope and higher and increasing VL over the study period. Those with 9 affected loci in general had the worst progression: low and declining CD4 counts in both cohorts, and high and relatively stable VL in both cohorts. Those with 8 affected loci tend to fall in the middle. Smoothing the raw data directly in each of these groups yielded similar results and nearly identical conclusions (results not shown).

4 Simulation results

We have performed simulation studies to assess the finite sample performance of our proposed testing procedure and compare its power to the standard linear-mixed-model-based procedures. For simplicity, we focused on a single marker $z$ in the absence of covariates and two potential functional outcomes generated from

[TABLE]

where $b_{ij}^{(m)}\sim N(0,0.25),j=0,1$ are independent and identically distributed (iid) random effects and $\epsilon_{ir}^{(m)}\sim N(0,0.25)$ are iid errors, for $m=1,2$. For each subject $i$, we generate the number of observations from a Poisson distribution $r_{i}\sim$ Poisson($\lambda$) + 2, and we generate $t_{ir}$ uniformly over the time interval $(0,2\pi)$. The parameter $\beta$ controls the magnitude of the genetic effect. The parameter $\alpha$ controls how linear the genetic affect is – when $\alpha=0$ the genetic effect is entirely linear, and when $\alpha=1$ the effect is entirely nonlinear. On the other hand, $\gamma$ controls the complexity of the mean process and the amount of inter-subject variability – when $\gamma=0$, the mean process is relatively simple but the inter-subject variability is high, and when $\gamma=1$ the mean process is complex and the inter-subject variability is low. The genetic factor $z_{i}$ is generated according to a binomial(2, 0.1), with 0.1 the minor allele frequency.

We examined the performance of the FPVC test statistic $Q$ (defined in (13), here denoted by “FPCA”) and its asymptotically equivalent counterpart $\bar{Q}$ (defined in the context of a single outcome in (12), here denoted ”Re-fitted”). For the purposes of comparison, we also examined the performance of a similar test statistic that does not use FPVC but instead employs a pre-specified basis. Consider the test statistic $Q_{\text{lin}}=\frac{1}{n}\sum_{m=1}^{2}\sum_{k=1}^{2}\left[\sum_{i=1}^{n}\xi^{(m){\dagger}}_{ik}{\widehat{z}}^{*}_{i}\right]^{2}$, where $\xi^{(m){\dagger}}_{ik}$ is the BLUP from the linear mixed model $y_{ir}=\beta_{0}+\beta_{1}t_{ir}+\xi_{i1}+\xi_{i2}t_{ir}+\epsilon_{ir}$. In the following, we denote results for $Q_{\text{lin}}$ by ”Linear”.

The number of FPCA scores for the $m$th outcome, $K_{m}$, was selected as the smallest $K$ such that the fraction of variation explained (FVE), $\sum_{k=1}^{K}{\widehat{\lambda}}_{k}/(\sum_{k}{\widehat{\lambda}}_{k})$, was at least $\wp=0.99$. To ensure that the scores for each outcome contributed comparably to the test statistics, we centered and scaled each outcome as $y_{ir}^{*(m)}=(y_{ir}^{(m)}-\bar{y}^{(m)})/{\widehat{\sigma}}_{y}^{(m)}$, prior to obtaining ${\widehat{\xi}}^{(m)}_{ik}$ and $\xi_{ik}^{(m){\dagger}}$, where ${\widehat{\sigma}}^{(m)}_{y}=\sqrt{(n-1)^{-1}\sum_{i,r}(y^{(m)}_{ir}-\bar{y}^{(m)})^{2}}$ and $\bar{y}^{(m)}=n^{-1}\sum_{i,r}y^{(m)}_{ir}$.

In the following we report power as the proportion of 1000 simulations for which the testing procedure produced a p-value below 0.05 to demonstrate the relative performance of the various testing procedures. To ensure that the asymptotic null distribution of the test statistic yields a valid testing procedure, we evaluate the entire distribution of p-values under the null hypothesis, including at levels much lower than 0.05.

4.1 Type I error

In the following we take $\lambda=6$. The empirical type I error rates for testing at the 0.05 level ranged from 0.040 ($\gamma$ = 0.25) to 0.048 ($\gamma$ = 0) for $Q$; from 0.036 ($\gamma$ = 1) to 0.047 ($\gamma$ = 0) for $\bar{Q}$; and from 0.040 ($\gamma$ = 0.75) to 0.059 ($\gamma$ = 1) for $Q_{\text{lin}}$. However, levels much smaller than 0.05 are necessary to control error rates in large-scale testing. Thus, in order to establish the validity of our testing procedure for performing many tests, we establish that the resulting p-values are approximately uniform under the null hypothesis. We performed $10^{6}$ simulations under the null, with $n=200,\gamma=1,$ and $\alpha=0$, and we obtained the type I error of FPVC testing at each of the following levels: $1\times 10^{-6},...,9\times 10^{-6},1\times 10^{-5},...,9\times 10^{-5},1\times 10^{-4},...,9\times 10^{-4},1\times 10^{-3},...,9\times 10^{-3}$. Results are depicted in Figure 4. The Figure shows that the level is preserved at all levels of testing. Further simulations would provide better approximations of type I error rates at smaller levels, but these results suggest that the asymptotic null distribution fits quite well in small samples.

4.2 Power

In figure 5, we display results for $n=100$ and all levels of $\gamma$ and $\alpha$. There, the figure demonstrates that, despite the fact that the true effect was linear when $\alpha=0$, both FPC-based tests dominate the linear-based tests, and the advantage of using FPVC, in terms of power, did not disappear, even when the true effect was linear. Notably, as $\gamma$ varied, we saw power gains by using the FPVC-based $Q$ and $\bar{Q}$, with the gains increasing as $\gamma$ approached 1 and the functional form of $Y_{i}^{(m)}(\cdot)$ became more complex, and the need to flexibly model it increased.

In all of our simulations, the FPVC methods dominated the linear method in terms of power while maintaining desirable type I error rates. We wanted to ensure that the improvement we were seeing was not simply due to the fact that the linear model used only two scores, a random intercept $\xi^{(m){\dagger}}_{i1}$ and a random slope $\xi^{(m){\dagger}}_{i2}$, for each outcome whereas the FPVC-based methods used $K_{m}$ scores, where $K_{m}$ was often selected larger than 2. Thus, we also considered the performance of scores based on fixed-basis expansions of ${\bf t}$, using either polynomial or spline bases.

Specifically, we fit models with $K=2,3,...,6$ degrees of freedom. For polynomial bases corresponded to the model $y_{ir}=\sum_{k=1}^{K}(\beta_{k}+\xi_{ik}){\widetilde{t}}_{ir}^{k-1}+\epsilon_{ir}$ for ${\widetilde{t}}_{ir}$ a centered and scaled version of $t_{ir}$. And for the spline basis, we used cubic B-splines constructed with the specified degrees of freedom with the bs function in the splines R package. Because, in some sense, FPCA does model selection by choosing the basis that explains the most variability in ${\bf y}$, we also perform model selection on the pre-specified bases in order to make the comparison fair. We select the model with the lowest AIC and use the ${\widehat{\xi}}_{ik}$s from that model in the testing procedure. We will call the test statistic based on B-splines $Q_{B}$ and the model based on polynomial bases $Q_{p}$.

Results are found in figure 6. We found that there were some situations when using the pre-specified B-spline basis could outperform the FPVC tests, particularly when $\gamma$ was near 0 (low mean complexity) and $\alpha$ was near 1 (linear genetic effect). However, the polynomial basis never outperformed FPVC. And as the complexity of the trajectory $\gamma$ increased, the desirability of FPVC testing always increased. This suggests that in simple problems, using a pre-specified basis may be preferable, but for complex effects and complex trajectories, FPVC will likely be preferred. And in general, if the complexity of the trajectory is unknown, FPVC testing offers a generally powerful method for all settings that is insensitive to tuning parameter selection.

5 Discussion

We have proposed functional principal variance component testing, a FPCA-based testing procedure for assessing the association between a set of genetic variants ${\bf z}_{\cal S}$ and a complexly varying longitudinal outcome ${\bf y}$ that is feasible on the genome-wide scale, allowing adjustment for other covariates. Unlike the standard mixed-model-based approaches, we do not model the trajectories $\{Y_{i}(\cdot)\}_{i=1}^{n}$ parametrically but use the data to identify the most parsimonious summaries of the trajectory patterns via FPCA. We subsequently test the association between the random coefficients $\boldsymbol{\xi}_{i}$ and the markers of interest using a test statistic motivated by variance component testing. Our procedure could potentially be much more powerful than procedures based on pre-specified bases, which might suffer power loss due to either high degrees of freedom or inability to capture the complexity in the trajectories. Furthermore, our FPVC testing is computationally efficient as we are able to perform thousands or even millions of tests quickly by separating the time-intensive FPCA from the testing. This makes our method feasible on the genome-wide scale where millions of marginal tests may be necessary. As an example, computing test statistics and p-values for FPVC testing typically takes less than 0.1 seconds for a set of 10 SNPs and both lCD4 and lVL combined. On the other hand, fitting a single linear mixed effects model for only lCD4 with a random effect for a small pre-specified B-spline basis takes more than 2 seconds. At the genome-wide scale we would observe a speed-up on the order of hours. Code for FPVC testing is available at https://github.com/denisagniel/fpvc.

It is important to note that while we make mild assumptions on the longitudinal outcome ${\bf y}$ to obtain the form of our proposed test statistic, the validity of FPVC testing requires no assumption about the relationship between ${\bf y}$ and ${\bf z}_{\cal S}$. FPVC testing remains valid even if the working mixed model (8) fails to hold. Additionally, while one can motivate the quantity ${\widetilde{\xi}}_{ik}$ as the conditional expectation of $\xi_{ik}$ under a normality assumption on $\xi_{ik}$ and $\epsilon_{ir}$, even when this normality fails to hold, testing based on $Q$ remains valid since the estimated eigenvalues and eigenfunctions from functional PCA are uniformly convergent to their limits [13]. In fact, one can consider FPVC model-free in that the test statistic $Q$ could be motivated simply as an estimated covariance. Furthermore, we assume that the errors $\epsilon_{ir}$ are iid with mean 0 and variance $\sigma^{2}$, but some relaxation of this assumption is possible for some ”degree of weak dependence and in cases of non-identical distribution” [13], while still maintaining the validity of our procedure.

FPVC testing can also simultaneously consider multiple sources of outcome information to better characterize complex phenotypes. With multiple longitudinal outcomes, one might wish to ensure that scores for all outcomes are roughly on the same scale, so that each outcome contributes comparably to the test statistic. To this end, one may consider a weighted version of (13) as

[TABLE]

where $\omega_{m}$ are nonnegative outcome-specific weights that can be pre-specified or data adaptive. Alternatively, in the absence of relevant weights, one can simply scale each ${\bf y}^{(m)}$ so that the magnitude of ${\widehat{\boldsymbol{\xi}}}_{i}^{(m)}$ is comparable across different values of $m$. Let $y^{*(m)}_{ir}=y^{(m)}_{ir}/{\widehat{\sigma}}^{(m)}_{y}$ where ${\widehat{\sigma}}^{(m)}_{y}=\sqrt{(n-1)^{-1}\sum_{i,r}(y^{(m)}_{ir}-\bar{y}^{(m)})^{2}}$ and $\bar{y}^{(m)}=n^{-1}\sum_{i,r}y^{(m)}_{ir}$. Then obtain ${\widehat{\boldsymbol{\xi}}}_{i}^{*(m)}$ via FPCA on $\{{\bf y}_{i}^{*(m)}\}_{i=1}^{n}$ and construct the test statistic $\sum_{m=1}^{M}\|n^{-\frac{1}{2}}\sum_{i=1}^{n}{\widehat{\boldsymbol{\xi}}}_{i}^{*(m)}{\widehat{\bf z}}_{i{\cal S}}^{*}\|_{F}^{2}$. Such a strategy appears to work well in simulation studies.

While we use FPCA to summarize the longitudinal trajectories for the purpose of testing with low degrees of freedom, in principle another suitably parsimonious nonparametric method could be used instead. For example, if observations were measured on a common, fine grid of points, then one could imagine using the methods in [20] to first regress ${\bf y}$ on ${\bf t}$, obtain the random effects estimates (similar to the $\boldsymbol{\xi}_{i}$ employed here), and use these in testing. However, no available approaches are as widely applicable as our FPCA-based approach, which can be used even when data are sparsely observed; other approaches may not have as small effective degrees of freedom as an FPCA-based method, and the resulting testing procedure may be more sensitive to correct tuning.