# The hidden factor: accounting for covariate effects in power and sample size
computation for a binary trait

Ziang Zhang
Department of Statistical Science
University of Toronto
Toronto
<EMAIL_ADDRESS>
&Lei Sun
Department of Statistical Science
University of Toronto
Toronto
<EMAIL_ADDRESS>

###### Abstract

Motivation: Accurate power and sample size estimation is crucial to the design
and analysis of genetic association studies. When analyzing a binary trait via
logistic regression, important covariates such as age and sex are typically
included in the model. However, their effects are rarely properly considered
in power or sample size computation during study planning. Unlike when
analyzing a continuous trait, the power of association testing between a
binary trait and a genetic variant depends, explicitly, on covariate effects,
even under the assumption of gene-environment independence. Earlier work
recognizes this hidden factor but implemented methods are not flexible.
Method: We thus propose and implement a generalized method for estimating
power and sample size for (discovery or replication) association studies of
binary traits that a) accommodates different types of non-genetic covariates
E, b) deals with different types of G-E relationships, and c) is
computationally efficient.
Results: Extensive simulation studies show that the proposed method is
accurate and computationally efficient for both prospective and retrospective
sampling designs with various covariate structures. A proof-of-principle
application focused on the understudied African sample in the UK Biobank data.
Results show that, in contrast to studying the continuous blood pressure
trait, when analyzing the binary hypertension trait ignoring covariate effects
of age and sex leads to overestimated power and underestimated replication
sample size.
Availability: The simulated datasets can be found on the online web-page of
this manuscript, and the UK Biobank application data can be accessed at
https://www.ukbiobank.ac.uk. The R package SPCompute that implements the
proposed method is available at CRAN.

_K_ eywords GWAS $\cdot$ Power Computation $\cdot$ Binary Trait $\cdot$
Covariate Effect $\cdot$ Replication Study

##  1 Introduction

Accurate power and sample size estimation is crucial to the design of many
scientific studies, including the ubiquitous genome-wide association studies
(GWAS) of complex and heritable human diseases and traits (Hong and Park,
2012). It is well known that replication studies with underestimated sample
sizes can result in false negatives, missing single nucleotide polymorphisms
(SNPs; $G$’s) that are truly associated with the phenotype of interest ($Y$)
(Patil et al., 2016). Additionally, recent work (Turley et al., 2018) has
shown that failure to correctly estimate power can also result in increased
false positives in pleiotropy studies, where different traits are jointly
analyzed and their GWAS summary statistics are aggregated.

The power and sample size calculation for a continuous trait is well
established, as the phenotype-genotype association analysis is through the
ordinary linear regression, regressing $Y$ on $G$ and important non-genetic
covariates $E$’s. It is then straightforward to show that the power of the
corresponding genetic association test only depends on the effect size and
minor allele frequency (MAF) of the SNP, sample size, and the unexplained
phenotypic variance (Korte and Farlow, 2013). That is, when analyzing a
continuous trait, the sample size for a replication study with sufficient
power is determined by the proportion of phenotypic variance explained by
genetic variants, which is also called narrow-sense heritability (Yang et al.,
2017; J Mayhew and Meyre, 2017).

In contrast, the power calculation for a binary disease outcome requires
additional considerations, as the association analysis typically uses the
logistic or probit regression model (Robinson and Jewell, 1991; Sjölander and
Greenland, 2013). Most heritability estimation methods were rigorously
developed for continuous traits only (Weissbrod et al., 2018; Yang et al.,
2010), and their applications to binary traits have been questioned (Golan et
al., 2014). At the same time, when analyzing a binary outcome $Y$, power of
analyzing a SNP $G$ is affected, explicitly, by the effect size of a non-
genetic covariate $E$, even if $E$ is independent of $G$ and/or there is no
$G$x$E$ interaction effect (Robinson and Jewell, 1991; Pirinen et al., 2012).
Therefore, accurate power and sample estimation for a binary trait-genetic
association analysis must explicitly consider the presence of non-genetic
covariates.

There have been several attempts in the literature to consider the general
problem of power and sample size computation for logistic regression.
Whittemore (1981) derived an approximation method, assuming that the disease
prevalence is small and the covariates have a joint distribution of
multivariate exponential. The approach of Whittemore (1981) was similarly
considered by Hsieh (1989), Hsieh et al. (1998), and Novikov et al. (2010).
Based on the asymptotic power approximation of the score or likelihood ratio
test under local alternatives, Self et al. (1992) and Self and Mauritsen
(1988) proposed an alternative approach that accommodates several categorical
covariates with finite configurations, which was then extended by Shieh (2000)
to allow for one categorical covariate with infinite configurations.

For genetic association studies, Quanto is the most commonly used software in
practice, implementing the method of Gauderman (2002a, b). The method uses the
expected value of a likelihood ratio test (LRT) statistic and accommodates
both continuous and categorical $E$’s for power analysis of $G$x$E$
interaction. However, the approach of Gauderman (2002b) implicitly assumes
that $G$ and $E$ are independent of each other, which may not hold in practice
for complex diseases (Plomin et al., 1977; Scarr and McCartney, 1983; Knafo
and Jaffee, 2013; Zhu et al., 2018; Namjou et al., 2019). Further, the
implemented software Quanto does not accommodate the presence of $E$ unless
the power computation is for $G$x$E$ interaction analysis. That is, $E$ cannot
be included when the power analysis is for the main effect of $G$.

Demidenko (2007), on the other hand, advocated the use of the Wald test to do
the power and sample size computation for logistic regression, and proposed a
method that allows $E$ and $G$ to be dependent through a second-stage logistic
regression model. However, the implemented web-tool (Demidenko, 2008) only
allows for one binary covariate, as otherwise the computation does not admit a
closed-form expression.

Lyles et al. (2007) proposed a different approach to power computation for
generalized linear models, based on the use of an expanded representative
dataset. The idea of expanded representative dataset provides accurate
approximation with good computational efficiency when sample size is small to
medium, but the computation becomes cumbersome when the sample size is large,
while large sample size (and small genetic effect size) is a feature of many
GWAS.

In this paper, we propose and implement a generalized method for estimating
power and sample size for genetic association studies of binary traits that a)
takes into account different types of non-genetic covariates $E$, b) allows
for different types of $G$-$E$ relationship, and c) has good computational
efficiency for large-scale studies. The utility of the proposed method is
illustrated and compared with the existing methods through extensive
simulation studies and an application study of the UK Biobank data (Sudlow et
al., 2015; Bycroft et al., 2017). The proposed method has been implemented as
a R package, SPCompute, available at CRAN.

##  2 Preliminary

###  2.1 Models

For simplicity of the notation, we assume without the loss of generality that
there is only one non-genetic covariate $E$; the method implementation and
application allows for multiple $E$’s. To study the relationship between a
trait $Y$ and a SNP $G$ of interest, conditional on the non-genetic covariate
$E$, we consider the following generalized linear model (glm; McCullagh and
Nelder (2019)),

$g(\mathbb{E}(Y|X))=g(\mu)=\beta_{0}+\beta_{G}G+\beta_{E}E=\eta,$ (1)

where $g(\cdot)$ is a link function, connecting the linear predictor $\eta$
with the mean function of $Y$. This glm model accommodates the analyses of
both continuous and binary traits. Here we focus on binary traits, for which
the logistic regression is the most commonly used model with
$g(\mu)=\log(\frac{\mu}{1-\mu})$.

Let ${X}$ be the design matrix, which has rows
$\\{(1,G_{i},E_{i})\\}_{i=1}^{n}$, where $n$ is the sample size. To ease
notation, we also use $X$ to denote the observed data, and we use
$\beta=(\beta_{0},\beta_{G},\beta_{E})^{T}$ to denote the vector of either all
regression parameters or their true values. The linear predictor $\eta$ is
then expressed as

$\eta=X\beta:=\beta_{0}+\beta_{G}G+\beta_{E}E.$

Following the convention in genetic association studies, the SNP genotypes,
$aa$, $Aa$ and $AA$, are assumed to follow the Hardy–Weinberg equilibrium
(HWE; Mayo (2008)), where $A$ is the minor allele with MAF of $p$. Also by
convention, $G$ is coded additively, tracking the number of allele $A$. Thus,
$\mathbb{P}(G=0)=(1-p)^{2}$, $\mathbb{P}(G=1)=2p(1-p)$ and
$\mathbb{P}(G=2)=p^{2}$.

###  2.2 The Wald test

To test the association between SNP $G$ and trait $Y$, i.e.
$H_{0}:\beta_{G}=0$ vs. $H_{1}:\beta_{G}\neq 0$, one can consider different
tests such as the LRT, Score or Wald tests. These tests have similar
asymptotic behaviors under the null hypothesis, and they are locally
equivalent (Rao et al., 1973; Serfling, 2009). However, as noted by Demidenko
(2007), these three likelihood-based tests differ globally. As Wald tests are
routinely used in GWAS, following the argument of Demidenko (2007), we carry
out our power and sample size computation based on a Wald test.

The Wald test statistic in our setting is expressed as

$T=\frac{{\hat{\beta}_{G}}^{2}}{{I_{{X}}^{-1}(\hat{\beta})_{[2,2]}}},$

where $I_{{X}}^{-1}(\hat{\beta})_{[2,2]}$ denotes the second diagonal element
of the matrix $I_{{X}}^{-1}(\hat{\beta})$, and $\hat{\beta}$ is the maximum
likelihood estimate (MLE) of $\beta$. $I_{{X}}(\beta)$ is the observed or
conditional Fisher information matrix, conditional on the observed $X$,
defined as

$I_{{X}}(\beta)={X}^{T}W(\beta){X},$

where $W(\beta)$ is a $n\times n$ diagonal matrix, with the $i^{th}$ element
as

$w_{i}=(\frac{\partial u_{i}}{\partial\eta_{i}})^{2}/\text{Var}(Y_{i}|X_{i}).$

Under the null hypothesis, $T$ is asymptotically $\chi^{2}_{1}$ distributed.

Using the above expression of $w_{i}$, it is easy to see that when analyzing a
continuous trait with residual variance $\sigma^{2}$ via linear regression
(i.e. using the identity link function), $\frac{\partial
u_{i}}{\partial\eta_{i}}=1$ and $w_{i}=1/\sigma^{2}$, which means
$I_{{X}}(\beta)$ depends on $\sigma^{2}$ but not on any regression
coefficients explicitly. In contrast, when analyzing a binary trait using
logistic regression,

$w_{i}=\frac{\exp(-\eta_{i})}{(1+\exp(-\eta_{i}))^{2}},$

which is a function of both $\beta_{G}$ and $\beta_{E}$. Thus, the size of the
non-genetic covariate effect $\beta_{E}$ explicitly influences the Fisher
information matrix, hence the power analysis of $\beta_{G}$.

###  2.3 The hidden factor in power and sample size computation

Assume the significance level of the test is $\alpha$, and the sample size is
large enough so that the asymptotic distribution of the Wald test statistic
can be used. Let

${V}_{G,{X}}:=I_{{X}}^{-1}(\beta)_{[2,2]}$

be the variance of $\hat{\beta}_{G}$, the power of the Wald test can be
computed as

$\Phi(-Z_{1-\alpha/2}+\frac{\beta_{G}}{\sqrt{V_{G,{X}}}})+\Phi(-Z_{1-\alpha/2}-\frac{\beta_{G}}{\sqrt{V_{G,{X}}}}),$

where $Z_{1-\alpha/2}$ denotes the $1-\alpha/2$ quantile of the standard
normal distribution. Worth reemphasizing is the fact that ${V}_{G,{X}}$, thus
the power of logistic regression, explicitly depends on both $\beta_{G}$ and
$\beta_{E}$, as discussed in Section 2.2 above.

The above power computation is for conditional power, conditional on the
observed ${X}$. However, for sample size determination for a successful
replication study, the corresponding power analysis is performed prior to
observing any data. In that case, the power is referred to as the
unconditional power (Lyles et al., 2007).

To compute the unconditional power, naturally we replace the conditional
Fisher information matrix $I_{{X}}(\beta)$ above with its unconditional
version $I_{n}(\beta)$. Let $I_{\bm{x}}(\beta)$ be the unconditional Fisher
information for a single observation $\bm{x}=(1,G,E)^{\prime}$. For the
logistic regression considered,

$\displaystyle I_{1}(\beta)=\mathbb{E}_{\mathbb{F}_{X}}[w\bm{x}^{T}\bm{x}]$
(2)
$\displaystyle=\mathbb{E}_{\mathbb{F}_{X}}\bigg{[}\frac{\exp(-(\beta_{0}+\beta_{G}G+\beta_{E}E)}{(1+\exp(-(\beta_{0}+\beta_{G}G+\beta_{E}E)))^{2}}\begin{bmatrix}1&G&E\\\
G&G^{2}&GE\\\ E&GE&E^{2}\end{bmatrix}\bigg{]},$

where the expectation is taken over, $\mathbb{F}_{X}$, the distribution of the
covariate space $X$ (i.e. both $G$ and $E$).

Once $I_{1}(\beta)$ has been computed for a given $\mathbb{F}_{X}$, the
unconditional Fisher information matrix for a random sample of size $n$ is

$I_{n}(\beta):=nI_{1}(\beta).$

The unconditional power is then

$\displaystyle\Phi(-Z_{1-\alpha/2}+\frac{\beta_{G}}{\sqrt{V_{G,n}}})+\Phi(-Z_{1-\alpha/2}-\frac{\beta_{G}}{\sqrt{V_{G,n}}}),$
(3)

where

$V_{G,n}:=I_{n}^{-1}(\beta)_{[2,2]}$

is based on the unconditional Fisher information matrix, $I_{n}(\beta)$.

To plan a successful replication study at the $\alpha$ level, the sample size
$n$ required to achieve a desirable power can be computed by simply inverting
the power function, which is monotonic with respective to $n$. Although the
sample size computation is for a specific genetic effect $\beta_{G}$, it is
clear that, similar to the conditional Fisher information, the unconditional
Fisher information in (2), therefore $V_{G,n}$ in (3), also depends
$\beta_{E}$. Thus, this hidden factor must be explicitly accounted for when
performing sample size calculation for a binary trait.

##  3 Methods

###  3.1 Designing the covariate space

To compute the unconditional Fisher information matrix $I_{n}(\beta)$, one
needs to compute the moments and covariance of a random sample pair
$(G_{i},E_{i})$ from the corresponding covariate space $\mathbb{F}_{X}$. An
appropriately designed covariate space $\mathbb{F}_{X}$ should be flexible
enough to accommodate potential complex dependence structure between $G$ and
$E$, while conceptually simple enough so that that practitioners can make use
of their domain knowledge.

In the work of Gauderman (2002b), the author implicitly assumed independence
between $G$ and $E$, requiring only the marginal distributions of $G$ and $E$.
Although this makes the method easy-to-implement, the assumption may not hold
in practice (Plomin et al., 1977; Scarr and McCartney, 1983; Knafo and Jaffee,
2013; Zhu et al., 2018; Namjou et al., 2019). Furthermore, the implemented
software Quanto only allows users to specify $E$ when the target analysis is
the $G\times E$ interaction effect, not the main effect of $G$.

The work of Demidenko (2007), on the other hands, allows $\mathbb{F}_{X}$ to
accommodate dependence between a binary $G$ and a binary $E$, by introducing a
second-stage logistic regression,

$\displaystyle\log\bigg{(}\frac{\mathbb{P}(G=1|E)}{\mathbb{P}(G=0|E)}\bigg{)}=\gamma_{0}^{*}+\gamma_{E}^{*}E,$
(4)

where $\gamma_{0}^{*}$ is determined by user-specified marginal probabilities
of $G$ and $E$. As a result, users will only need to additionally input the
knowledge about $\gamma_{E}^{*}$ to fully specify $\mathbb{F}_{X}$. However,
the method of Demidenko (2007) is designed for a binary $G$ (hence the typical
GWAS additive coding of $G$ not applicable) and a binary $E$, and its
generalization to different types of $G$ and $E$ is nontrivial.

Here, we utilize the idea of a second-stage regression of Demidenko (2007) but
extend it to a more general setting. Instead of treating $E$ as a covariate in
the second-stage regression, we consider it to be the response variable such
that,

$\displaystyle g_{2}(\mathbb{E}(E|G))=\gamma_{0}+\gamma_{G}G,$ (5)

where $g_{2}$ is the link function, being identity when $E$ is continuous and
logit when $E$ is binary.

Compared with the second-stage regression model in (4), the proposed method
can accommodate different types of $E$ in a unified framework. When $E$ is
continuous, the regression model also requires $\text{Var}(E|G)$ in order to
be fully specified. The value of $\text{Var}(E|G)$ can be computed based on
user-provided information such as $\mu_{E},\sigma_{E}$ and $p$.

We note that the proposed method can account for multiple $E$’s by specifying
the corresponding $g_{2}$ function for each $E$ considered. Later, we will
demonstrate the utility of our approach in a UK Biobank data application of
hypertension, for which both age and sex are important covariates to consider
for power and sample size computation. In the rest of this section, we assume
there is only one covariate $E$ for clarity of the presentation.

###  3.2 Proposed method 1: Semi-Simulation (P1.SS)

The estimation of the unconditional power heavily depends on the computation
of $I_{n}(\beta):=nI_{1}(\beta)$. Unfortunately, unless in some special cases
such as when both $G$ and $E$ are binary, $I_{1}(\beta)$ in (2) does not have
a closed-form expression for a general $\mathbb{F}_{X}$ (Demidenko, 2007).
Thus, to estimate $I_{n}(\beta)$, we propose to use a sample estimate.

Specifically, for a large integer $B$, we simulate independent observations
$\\{G_{i},E_{i}\\}_{i=1}^{B}$ from the covariate space $\mathbb{F}_{X}$, and
for each $\bm{x}_{i}=(1,G_{i},E_{i})^{\prime}$ we compute the corresponding
conditional Fisher information matrix,

$I_{\bm{x}_{i}}(\beta)={\bm{x}_{i}}^{T}w_{i}(\beta){\bm{x}_{i}},$

where

$w_{i}(\beta)=\frac{\exp(-(\beta_{0}+\beta_{G}G_{i}+\beta_{E}E_{i})}{(1+\exp(-(\beta_{0}+\beta_{G}G_{i}+\beta_{E}E_{i})))^{2}}.$

By a simple application of the law of large number, the sample estimate,

$\displaystyle\tilde{I}_{n}(\beta):=\frac{\sum_{i=1}^{B}nI_{\bm{x}_{i}}(\beta)}{B},$
(6)

converges almost surely to the true expected matrix $I_{n}(\beta)$ as $B$
grows.

As we will illustrate later in the simulation studies, $\tilde{I}_{n}(\beta)$
exhibits little variation for large $B$ (e.g. >10,000). Furthermore, the
proposed semi-simulation method is scalable, as for each
$I_{\bm{x}_{i}}(\beta)$ we only compute the observed Fisher information matrix
for one single observation. Thus, the computational load depends on $B$ but is
independent of the target sample size $n$. Once $I_{n}(\beta)$ is replaced by
$\tilde{I}_{n}(\beta)$, the power computation can proceed using equation (3),
and sample size estimation by inverting the power function.

###  3.3 Proposed method 2: Representative Dataset (P2.RD)

An alternative method that does not rely on plugging in the sample estimate of
$I_{n}(\beta)$ is through the use of a representative dataset, an idea that
was originally suggested by O’brien (1986) and later extended by Lyles et al.
(2007).

In our setting, given a sample size $n$, assume there exists a representative
covariate sample
$\\{\bm{x}_{i}\\}_{i}^{n}=\\{(1,G_{i},E_{i})^{\prime}\\}_{i=1}^{n}$ from the
covariate space $\mathbb{F}_{X}$, which we define later. We then expand
$\\{\bm{x}_{i}\\}_{i}^{n}$ to consider both possible outcomes of the binary
trait, so that each observation $\bm{x}_{i}$ splits into
$\\{\bm{x}_{i},y_{i}=0\\}$ and $\\{\bm{x}_{i},y_{i}=1\\}$. Additionally, each
expanded observation is given a weight, so that
$\delta_{i}^{0}+\delta_{i}^{1}=1$, where

$\delta_{i}^{l}=\mathbb{P}(y_{i}=l|\bm{x}_{i})=\frac{\exp(\beta_{0}+\beta_{G}G_{i}+\beta_{E}E_{i})^{l}}{1+\exp(\beta_{0}+\beta_{G}G_{i}+\beta_{E}E_{i})},$
(7)

for $l=$ 0 and 1.

Thus, the original representative dataset $\\{\bm{x}_{i}\\}_{i=1}^{n}$ is now
expanded into the following representative dataset,

$\displaystyle RD:=\begin{Bmatrix}\bm{x}_{i},&y_{i}=0,&\delta_{i}^{0}\\\
\bm{x}_{i},&y_{i}=1,&\delta_{i}^{1}\\\ \end{Bmatrix}_{i=1}^{n},$ (8)

which has $2n$ (weighted) observations. Standard MLE of $V_{G,n}$, derived
from the corresponding weighted log-likelihood, yields $\hat{V}_{G,n}$, which
can be directly plugged into equation (3) to complete the power computation
(Lyles et al., 2007).

It remains to be discussed what is a representative
$\\{\bm{x}_{i}\\}_{i=1}^{n}$ and how the expanded representative dataset $RD$
can be obtained in our study setting. In the case of conditional power
analysis where covariates are already observed, the observed
$\\{\bm{x}_{i}\\}_{i=1}^{n}$ can be directly used in (7) to establish the
representative dataset of (8).

For the unconditional power analysis, $\\{\bm{x}_{i}\\}_{i=1}^{n}$ can be
obtained by using user-provided $\mathbb{F}_{X}$. Lyles et al. (2007) provided
examples on how to define the notion of representative dataset for different
types of $\mathbb{F}_{X}$. We follow the procedures of Lyles et al. (2007) for
the types of $\mathbb{F}_{X}$ considered in Section 3.1.

When $E$ is binary and the link function in (5) is logistic, we can compute
the expected counts for category $\\{(G=i,E=j)\\}$, $i=0,1$, and $2$, and
$j=0$ and $1$, as

$n_{i,j}=n\mathbb{P}(G=i,E=j),$

using the available information such as MAF and the inheritance mode, with
appropriate rounding to ensure that $n_{i,j}$’s are integers and sum to $n$.

When $E$ is continuous and the link function is identity with
$\text{Var}(E|G)=\sigma_{E}^{2}$, we first categorize the dataset based on $G$
such that $n_{i}=n\mathbb{P}(G=i)$, $i=0,1$, and 2. Then for each of the
$j=1,\ldots,n_{i}$ observations of $\\{G_{j}=i\\}_{j=1}^{n_{i}}$,

$E_{j}=\gamma_{0}+\gamma_{G}i+\sigma_{E}\Phi^{-1}[(j-0.375)/(n_{i}+0.25)],$

where $\Phi^{-1}$ is the inverse of the cumulative distribution function of
the standard normal.

##  4 Simulation Studies

###  4.1 Overview of the simulation design

We compared the power and sample size computed using the proposed P1.SS and
P2.RD methods (implemented as the R package SPCompute available at CRAN) with
those computed using Quanto of Gauderman (2002b) (version 1.2.4 downloaded
from http://hydra.usc.edu/GxE), and the method of Demidenko (2007) using its
web-platform (dartmouth.edu/$\sim$eugened/power-samplesize.php).

We considered three different scenarios for $\mathbb{F}_{X}$, including no
covariate $E$ (as Quanto does not allow for $E$), $E$ being binary (as the
method of Demidenko (2007) only allows for binary $E$), and $E$ being
continuous. Although we only considered one $E$ in the simulation studies for
method comparison, our implemented SPCompute R package allows for multiple
$E$’s, as demonstrated in our UK Biobank application study in Section 5.
Finally, for the simulation studies we also considered three study designs,
where S1 is case-control retrospective (Quanto only allows for the case-
control study design), while S2 and S3 are prospective to reflect the design
of the emerging biobank-sized data such as the UK Biobank data used in our
application.

The accuracy of each method was assessed by comparing the computed power (and
sample size) with the empirical values obtained through 1,000 independent
replications. Given the large number of replications, the empirical values
were treated as the oracle values and used to benchmark. We calculated the
average and the maximum of the absolute error (AE) of each computed power (and
sample size) as compared with the oracle values. The more accurate method is
expected to have smaller mean AE and max AE. By convention, for replication
sample size computation, the desirable power was set to be 80% at the
significance level of $0.05$, and for consistency, the same significance level
was used for power computation.

####  4.1.1 Scenario 1: No covariate $E$ with a case-control retrospective
study design

The choice of no covariate effect was to accommodate the implementation of
Quanto of Gauderman (2002b). Without loss of generality, the disease
prevalence was assumed to be 20%, and the observations were obtained
independently with a retrospective sampling design and the standard case-to-
control ratio of 1-to-1.

The associated SNP has a MAF of $0.1$, with a dominant effect $\beta_{G}$
ranging from $\log(1.1)$ to $\log(2.5)$. That is, the odds ratio (OR) of the
associated SNP ranged from 1.1 to 2.5. The choice of a dominant genetic effect
was to accommodate the implementation of Demidenko (2007) method, which only
allows for a binary $G$. Finally, we used $\beta_{0}=-2$, though we note that
the intercept parameter does not affect the power of a case-control study.

####  4.1.2 Scenario 2: Binary $E$ with a prospective study design

Similar to S1 above, in the second scenario the disease prevalence is also
20%, and the associated SNP with MAF of $0.1$ has a dominant effect
$\beta_{G}$ ranging from $\log(1.1)$ to $\log(2.5)$. However, the observations
were obtained independently with a prospective sampling design as in the UK
Biobank data. Additionally, the non-genetic covariate $E$ has a population
exposure rate of $\mathbb{P}(E=1)=0.3$ with effect $\beta_{E}=\log(2.5)$.
Finally, $\gamma_{G}=\log(0.2)$, quantifying the dependency between $G$ and
$E$ as defined in (5).

To implement Quanto of Gauderman (2002b), which only allows for case-control
study design, we used a case-to-control ratio of 1-to-4 to approximate the
result for a disease with prevalence of 20%. Additionally, when the power
analysis is about $G$ main effect (as opposed to $G\times E$ interaction
effect), Quanto does not consider the presence of $E$. Thus, we only input the
information about $G$ in the implementation of Quanto. The method of Demidenko
(2007) accommodates the presence of one binary covariate $E$ for the power
(and sample size) computation; $G$ must be binary, hence the dominant genetic
model was assumed for method implementation.

####  4.1.3 Scenario 3: Continuous $E$ with a prospective study design

For this last scenario, without loss of generality, the covariate $E$ was
assumed to follow the standard normal distribution conditional on $G$. The
dependency between $G$ and $E$ was set to $\gamma_{G}=\log(0.5)$. All other
model specifications are the same as in S2 above, including the disease
prevalence (20%), MAF (0.1), the genetic effect (ranging from $\log(1.1)$ to
$\log(2.5)$), and the non-genetic covariate effect ($\log(2.5)$.

As in the previous scenario, we ignored the information on $E$ for the
implementation of Quanto. For the method of Demidenko (2007), which only
allows for a binary $E$, we considered two approaches. We first omitted the
continuous covariate $E$ (corresponding results∗), and we then dichotomized
$E$ by defining $\tilde{E}:=\mathbb{I}(E>0)$. This corresponds to creating two
misspecified models,

$\displaystyle\log\bigg{(}\frac{\mathbb{P}(Y=1|G,\tilde{E})}{\mathbb{P}(Y=0|G,\tilde{E})}\bigg{)}$
$\displaystyle=\beta_{0}+\beta_{G}G+\tilde{\beta}_{E}\tilde{E},$ (9)

and

$\displaystyle\log\bigg{(}\frac{\mathbb{P}(\tilde{E}=1|G)}{\mathbb{P}(\tilde{E}=0|G)}\bigg{)}$
$\displaystyle=\tilde{\gamma}_{0}+\tilde{\gamma}_{G}G.$ (10)

As the parameter values specified for the true models (1) and (5) cannot be
directly used for the two misspecified models, we used estimated
$\tilde{\gamma}_{G}$ and $\tilde{\beta}_{E}$. We first simulated a large a
number of observations $\\{G_{i},E_{i},Y_{i}\\}_{i}^{3\times 10^{5}}$ using
the true model. We then dichotomized the continuous $E$ to obtain $\tilde{E}$
as specified above. Finally, we regressed $Y$ on $G$ and $\tilde{E}$, and
$\tilde{E}$ on $G$ to obtain sample estimates of $\tilde{\beta}_{E}$ and
$\tilde{\gamma}_{G}$ for the second implementation of the method of Demidenko
(2007).

####  4.1.4 Methods comparison across the three scenarios

Figure 1 shows the computed power and sample size curves, across the three
scenarios, and the empirical results are consistent with our analytical
expectation: Ignoring covariate effect, at the (replication) study planning
stage, can lead to overestimated power and underestimated replication sample
size for studying a binary trait.

In the left panel of Figure 1, the OR of the associated SNP was fixed at 1.5
in (a,c) and 1.3 in (e), and Figures 1(a), 1(c) and 1(e) show the estimated
power for different sample size at the significance level of 0.05. Results
clearly show that, if there are no covariates, all methods perform well
(Figure 1(a)). In the presence of a binary covariate, Quanto tends to
overestimate the power of an association study, while both Demidenko and the
proposed two methods provide power estimates close to the Oracle values
(Figure 1(c)). However, if the influential covariate is continuous (e.g. age
as in our UK Biobank application for hypertension), only the proposed two
methods (P1.SS and P2.DD) perform well (Figure 1(e)). The superior
performances of the proposed two methods are also shown in Table 1, as P1.SS
and P2.RD have smaller average and smaller maximum absolute error, as compared
with the true power.

In the right panel of Figure 1, Figures 1(b), 1(d) and 1(f) show the estimated
sample sizes, necessary to achieve 80% power at the 0.05 level, to
successfully replicate an associated SNP with OR ranges from 1.1 to 2.5.
Similar conclusion can be drawn here, as the existing methods tend to
underestimate the necessary sample size for a successful replication study in
the presence of influential covariate, while the proposed P1.SS and P2.RD
methods are accurate.

###  4.2 Choice between the proposed P1.SS and P2.RD methods from the
computational perspective

To select between the two proposed methods, P1.SS and P2.RD as respectively
described in Section 3.2 and Section 3.3, here we study factors influencing
the computational efficiencies of the two methods and make recommendations to
practitioners.

Conceptually, the computational efficiency of the semi-simulation P1.SS method
depends on $B$, the number of independent observations drawn from
$\mathbb{F}_{X}$ in order to obtain $\tilde{I}_{n}(\beta)$ in (6). As
$\tilde{I}_{n}(\beta)$ is based on $B$ replicates of one-sample
$I_{\bm{x}_{i}}(\beta)$, the targeted sample size $n$ does not have a direct
impact on computational time. In contrast, the computational time of the P2.RD
method depends on $n$, as the method first creates a representative dataset of
size $n$ from $\mathbb{F}_{X}$ then expand it to weighted $2n$ observations as
in (8).

To numerically demonstrate the computational properties of the two methods,
without loss of generality, we considered simulation scenario S2 in Section
4.1.2 and used $\beta_{G}=\log(1.5)$ for illustration. Results in Figure 2(a)
confirm our analytical expectation: The computational time of P2.RD grows
linearly with respect to $n$, while that of P1.SS is independent of $n$.

However, the accuracy of P1.SS depends on large $B$; we used $B=$ 10,000
($\log_{10}(B)=4$) in Figure 2(a). Figure 2(b) shows the stability of P1.SS
with respect to $\log_{10}(B)$. For each choice of $\log_{10}(B)$ from 3.0 to
4.3, the $\log_{10}$ sample standard error of the power (SE) computed by
P1.SS, obtained from 1,000 independently simulated replicates, was shown in
Figure 2(b). Results clearly show that $B\geq$ 10,000 ($\log_{10}(B)\geq 4$)
leads to negligible SE of less than 0.01 ($\log_{10}(\text{SE})<-2$). Thus,
the default value for $B$ in our method implementation is 10,000.
Interestingly, the relationship between $B$ and SE appears to be approximately
log-log linear.

We note that Y-axis in Figure 2(a) measures the run time in seconds for
computation of one set of parameter values (i.e. per computation). In
practice, it is often necessary to run power and sample size analysis for a
fine grid of a large number of possible parameter values. Thus, the run time
difference aggregates and can differ significantly between P1.SS and P2.RD. In
general, when $n$ is larger than 25,000, P1.SS is preferred over P2.RD. Thus,
although the implemented software SPCompute includes both methods, it sets
P1.SS as the default method.

##  5 Application to the UK Biobank data

To illustrate the practical utility of the proposed power and sample size
computation methods, we applied them to the UK Biobank data (Sudlow et al.,
2015; Bycroft et al., 2017), focusing on the understudied participants with
African background. Without loss of generality, we chose hypertension as the
binary trait of interest. For completeness, we also analyzed (diastolic) blood
pressure, a continuous trait to contrast with the binary trait.

###  5.1 Sample and SNP data quality control

We started with the 3,460 participants with self-reported ancestry being
African. First, we followed the standard practice (Marees et al., 2018) to
filter out individuals with genotype missingness higher than 20 percent. To
remove related individuals, we then filtered out individuals with kinship
coefficient larger than 0.25, which ended up with 3,182 (approximately)
unrelated self-reported Africans.

To account for reporting error and other ancestry biases, we then performed
principle component analysis (PCA) using the overall principle components
(PCs) provided by UKB (Data-Field 22009). Figure 3(a) shows the first two PCs
of all UKB samples, stratified by self-reported Africans and Others, which
suggests reporting error. We then applied a K-mean algorithm with $K=4$
(Hartigan and Wong, 1979) to the 3,182 unrelated self-reported Africans
(Figure 3(b-d)) using the overall PCs provided. Cluster 1 contains 2,601
individuals, 75% of all self-reported African participants (Figure 3(b-d)).
Following the common practice, we also computed new PCs using only the 3,182
individuals (Supplementary Figure S1). Among the 2,601 individuals in Cluster
1, we then removed 91 individuals whose new PCs were four standard deviations
away from the mean. Thus, the final GWAS sample consists of $n=$ 2,510 (1232
females and 1278 males) unrelated individuals with PC-confirmed African
ancestry.

For the genetic data, we started with 784,256 genotyped autosomal SNPs (Data-
Field 22418), and we then filtered out SNPs based on the thresholds of HWE
p-value$<$1e-10, MAF$<0.01$ and missing rate $>$0.2. The X chromosome was not
included in our analysis due to the recent report of previously unrecognized
data quality issue of the X chromosome (Wang et al., 2022). In total, 379,003
common, good quality autosomal SNPs were selected for the subsequent analyses.

###  5.2 GWAS of hypertension and blood pressure

We considered two phenotypes, one binary (hypertension) and the other
continuous (diastolic blood pressure). In this application, we only considered
their measurements at the initial assessment, as longitudinal data analysis is
beyond the scope of this work. There are two measurements of blood pressure
during the initial assessment, and we used the average. Additionally, for
blood pressure we considered the automated reading instead of the manual
reading. Among the $n=$ 2,510 analyzed individuals, the prevalence rate of
hypertension is 39.48 percent, and the diastolic blood pressure has mean 85.35
and standard deviation 10.75.

The GWAS of both traits included age and sex as covariates. The analysis of
the binary hypertension trait used logistic regression, and the analysis of
continuous trait used linear regression as in convention. The two GWAS results
are displayed in Figure 4(a) and (b), respectively, for hypertension and blood
pressure. Given the small sample size, it is not surprising that none of the
SNPs reached the genome-wide significance of 5e-8 (Dudbridge and Gusnanto,
2008).

###  5.3 Power and sample size computation

To illustrate the importance of including covariates in power and sample size
computation for a binary trait, as a proof-of-principle, we focused on the top
five ranked SNPs from each GWAS. The effect size estimates of both SNPs and
covariates were used for the corresponding power and sample size computation,
though we recognize the potential issue of winner’s curse (Sun et al., 2011);
this issue does not change characteristically the conclusion drawn from the
methods comparison.

To illustrate the important role of age and sex in study planning of the
binary trait of hypertension, we computed powers and required replication
sample sizes twice. We first ignore the two covariates as commonly done
(without E), which is equivalent to using Quanto; the method of Demidenko
(2007) is not applicable as there are two covariates. We then accounted for
covariate effects (with E) using the proposed method P1.SS, the default method
implemented in SPcompute. Finally, in addition to $\alpha=0.05$, the standard
significance level used for planning a successful replication study, we also
considered 5e-8, the genome-wide significance level to demonstrate the power
and sample size needed for the discovery GWAS to achieve 80% power. For
comparison, we also analyzed the continuous trait of blood pressure, where the
covariate effects are implicitly incorporated through the specification of the
residual variance.

For each trait analyzed, the computed powers and for the top SNPs are shown in
Figure 5, which are consistent with our simulation results. For the binary
hypertension trait (the left panel of Figure 5), the higher blue bars in
Figure (a) show that ignoring covariate effects leads to overestimated power
of our discovery study (at $\alpha=$ 5e-8); the power for $\alpha=$ 0.05 is
close to 100% as expected, thus uninteresting and not shown. The shorter blue
bars in Figure (c) show that ignoring covariate effects leads to
underestimated discovery sample size (for 80% power at $\alpha=$ 5e-8).
Similarly, the shorter blue bars in Figure (e) show that ignoring covariate
effects leads to underestimated replication sample size (for 80% power at
$\alpha=$ 0.05) when planning a replication study of a binary trait. In
contrast, when studying the continuous blood pressure trait, age and sex
effects are already implicitly incorporated through the residual variance.
Thus, covariate effects $\beta_{E}$’s do not have to be explicitly included in
the power and sample size computation for a continuous trait.

##  6 Discussion

Adjusting for covariate effect is a standard practice in GWAS, but it is
rarely done at the study planning stage and in replication studies. When
analyzing a continuous trait through linear regression, covariate effects are
implicitly accounted for through residual variance. However, when analyzing a
binary trait through logistic regression, covariate effects ($\beta_{E}$’s)
must be explicitly specified and included in power and sample size
computation, in addition to the genetic effect of interest ($\beta_{G}$).

This phenomenon is known in the statistics literature, but tools available to
practitioners are limited. For example, the most well-known software Qaunto
does not consider $\beta_{E}$ unless the power analysis is for $G$x$E$
interaction analysis (Gauderman, 2002b, a), while the method of Demidenko
(2007) only allows for one binary $E$. In this work, we developed and
implemented a flexible software SPCompute for accurate and efficient power and
sample size computation for a binary trait. We applied the proposed method to
the UK Biobank data, analyzing the binary hypertension trait and
simultaneously accounting for age and sex covariate effects in power and
sample size computation. We also conducted extensive simulation studies to
demonstrate the accuracy and efficiency of the proposed method.

However, there are still several limitations of the proposed method that
require future work to address. For example, winner’s curse where the effect
size estimates of significant SNPs are biased upward is known to be a common
problem in GWAS (Sun and Bull, 2005; Zöllner and Pritchard, 2007; Zhong and
Prentice, 2010; Sun et al., 2011). Therefore, it would be of interest to
investigate how SPCompute accounts for the winner’s curse. Another direction
of extension would be to account for mis-classification (particularly the
control data), which affects the power of an association study (Rekaya et al.,
2016; Zhang and Yi, 2020; Lin et al., 2020). Additionally, the proposed
framework can be further generalized to accommodate the simultaneous analysis
of multiple rare variants (Derkach et al., 2014). Finally, the proposed method
assumes a random sample of unrelated individuals. Power and sample size
computation for related individuals are worthy future work.

The proposed method, and the default setting of the implemented software
SPCompute, assumes all the parameters are specified based on a prospective
model as in the UKB application. However, it can also be applied to data
collected through the case-control design by modifying the disease prevalence
parameter to reflect the case-control ratio used, as shown in Demidenko
(2007). Additionally, our extensive simulation studies in Section 4
demonstrated that the proposed method is accurate when parameters were
specified based on a prospective model while data were collected through case-
control design. However, it should be noted that unlike the regression
parameters $\beta_{G}$ and $\beta_{E}$, the covariate space $\mathbb{F}_{X}$
can be very different once being conditioned on the case-control ratio
(Prentice and Pyke, 1979; Self and Mauritsen, 1988).

Our UKB application in Section 5 only serves as a proof-of-principle and
highlights the practical utility of SPCompute, such as its ability to handle
both binary and continuous covariates. We made some simplifying assumptions to
make the example easier to be understood. For example, we accounted for the
covariate effects of age and sex simultaneously by introducing two separate
models in the second stage regression of (5). This implicitly assumed the two
covariates are conditionally independent given the SNP $G$, an assumption that
might be unrealistic in more complex settings. Finally, the framework of the
proposed method can be generalized to incorporate the gene-gene and gene-
environment interaction effects, which we will provide as future software
updates.

## Acknowledgements

ZZ is a trainee of the CANSSI-ONTARIO STAGE (Strategic Training for Advanced
Genetic Epidemiology) training program at the University of Toronto.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any
commercial or financial relationships that could be construed as a potential
conflict of interest.

## Author Contributions

ZZ developed the method, performed the analyses, summarized the results, and
drafted the manuscript. LS conceptualized and supervised the study. Both
authors read and approved the final manuscript.

## Funding

This research was funded by the Natural Sciences and Engineering Research
Council of Canada (NSERC, RGPIN-04934), the Center for Addition and Mental
Health (CAMH) Discovery Fund Seed Funding, and the University of Toronto Data
Sciences Institute (DSI) Catalyst Grant.

## Data Availability Statement

The UK Biobank data were used under the license for this study (application
number 64875). Data are available at www.ukbiobank.ac.uk/ with the permission
of UK Biobank. The simulated datasets can be found at the online repository of
this manuscript.

## References

  * Bycroft et al. (2017) Bycroft, C., Freeman, C., Petkova, D., Band, G., Elliott, L. T., Sharp, K., Motyer, A., Vukcevic, D., Delaneau, O., O’Connell, J., et al. (2017).  Genome-wide genetic data on~ 500,000 uk biobank participants.  BioRxiv, page 166298.
  * Demidenko (2007) Demidenko, E. (2007).  Sample size determination for logistic regression revisited.  Statistics in medicine, 26(18), 3385–3397.
  * Demidenko (2008) Demidenko, E. (2008).  Sample size and optimal design for logistic regression with binary interaction.  Statistics in medicine, 27(1), 36–46.
  * Demidenko (2020) Demidenko, E. (2020).  Approximations of the power functions for wald, likelihood ratio, and score tests and their applications to linear and logistic regressions.  Model Assisted Statistics and Applications, 15(4), 335–349.
  * Derkach et al. (2014) Derkach, A., Lawless, J. F., and Sun, L. (2014).  Pooled association tests for rare genetic variants: a review and some new results.  Statistical Science, 29(2), 302–321.
  * Dudbridge and Gusnanto (2008) Dudbridge, F. and Gusnanto, A. (2008).  Estimation of significance thresholds for genomewide association scans.  Genetic Epidemiology: The Official Publication of the International Genetic Epidemiology Society, 32(3), 227–234.
  * Forstmeier and Schielzeth (2011) Forstmeier, W. and Schielzeth, H. (2011).  Cryptic multiple hypotheses testing in linear models: overestimated effect sizes and the winner’s curse.  Behavioral ecology and sociobiology, 65(1), 47–55.
  * Garner (2007) Garner, C. (2007).  Upward bias in odds ratio estimates from genome-wide association studies.  Genetic Epidemiology: The Official Publication of the International Genetic Epidemiology Society, 31(4), 288–295.
  * Gauderman (2002a) Gauderman, W. J. (2002a).  Sample size requirements for association studies of gene-gene interaction.  American journal of epidemiology, 155(5), 478–484.
  * Gauderman (2002b) Gauderman, W. J. (2002b).  Sample size requirements for matched case-control studies of gene–environment interaction.  Statistics in medicine, 21(1), 35–50.
  * Golan et al. (2014) Golan, D., Lander, E. S., and Rosset, S. (2014).  Measuring missing heritability: inferring the contribution of common variants.  Proceedings of the National Academy of Sciences, 111(49), E5272–E5281.
  * Hartigan and Wong (1979) Hartigan, J. A. and Wong, M. A. (1979).  Algorithm as 136: A k-means clustering algorithm.  Journal of the royal statistical society. series c (applied statistics), 28(1), 100–108.
  * Hong and Park (2012) Hong, E. P. and Park, J. W. (2012).  Sample size and statistical power calculation in genetic association studies.  Genomics & informatics, 10(2), 117.
  * Hsieh (1989) Hsieh, F. (1989).  Sample size tables for logistic regression.  Statistics in medicine, 8(7), 795–802.
  * Hsieh et al. (1998) Hsieh, F. Y., Bloch, D. A., and Larsen, M. D. (1998).  A simple method of sample size calculation for linear and logistic regression.  Statistics in medicine, 17(14), 1623–1634.
  * J Mayhew and Meyre (2017) J Mayhew, A. and Meyre, D. (2017).  Assessing the heritability of complex traits in humans: methodological challenges and opportunities.  Current genomics, 18(4), 332–340.
  * Knafo and Jaffee (2013) Knafo, A. and Jaffee, S. R. (2013).  Gene–environment correlation in developmental psychopathology.  Development and psychopathology, 25(1), 1–6.
  * Korte and Farlow (2013) Korte, A. and Farlow, A. (2013).  The advantages and limitations of trait analysis with gwas: a review.  Plant methods, 9(1), 1–9.
  * Lin et al. (2020) Lin, Y.-C., Brooks, J. D., Bull, S. B., Gagnon, F., Greenwood, C. M., Hung, R. J., Lawless, J., Paterson, A. D., Sun, L., and Strug, L. J. (2020).  Statistical power in covid-19 case-control host genomic study design.  Genome medicine, 12(1), 1–8.
  * Lyles et al. (2007) Lyles, R. H., Lin, H.-M., and Williamson, J. M. (2007).  A practical approach to computing power for generalized linear models with nominal, count, or ordinal responses.  Statistics in Medicine, 26(7), 1632–1648.
  * Marees et al. (2018) Marees, A. T., de Kluiver, H., Stringer, S., Vorspan, F., Curis, E., Marie-Claire, C., and Derks, E. M. (2018).  A tutorial on conducting genome-wide association studies: Quality control and statistical analysis.  International journal of methods in psychiatric research, 27(2), e1608.
  * Mayo (2008) Mayo, O. (2008).  A century of hardy–weinberg equilibrium.  Twin Research and Human Genetics, 11(3), 249–256.
  * McCullagh and Nelder (2019) McCullagh, P. and Nelder, J. A. (2019).  Generalized linear models.  Routledge.
  * Namjou et al. (2019) Namjou, B., Lingren, T., Huang, Y., Parameswaran, S., Cobb, B. L., Stanaway, I. B., Connolly, J. J., Mentch, F. D., Benoit, B., Niu, X., et al. (2019).  Gwas and enrichment analyses of non-alcoholic fatty liver disease identify new trait-associated genes and pathways across emerge network.  BMC medicine, 17(1), 1–19.
  * Novikov et al. (2010) Novikov, I., Fund, N., and Freedman, L. (2010).  A modified approach to estimating sample size for simple logistic regression with one continuous covariate.  Statistics in medicine, 29(1), 97–107.
  * O’brien (1986) O’brien, R. (1986).  Using the sas system to perform power analyses for log-linear models.  In Proceedings of the eleventh annual SAS users group international conference, pages 778–782. SAS Institute Inc., Cary, North Carolina.
  * Patil et al. (2016) Patil, P., Peng, R. D., and Leek, J. T. (2016).  What should researchers expect when they replicate studies? a statistical view of replicability in psychological science.  Perspectives on Psychological Science, 11(4), 539–544.
  * Pirinen et al. (2012) Pirinen, M., Donnelly, P., and Spencer, C. C. (2012).  Including known covariates can reduce power to detect genetic effects in case-control studies.  Nature genetics, 44(8), 848–851.
  * Plomin et al. (1977) Plomin, R., DeFries, J. C., and Loehlin, J. C. (1977).  Genotype-environment interaction and correlation in the analysis of human behavior.  Psychological bulletin, 84(2), 309.
  * Prentice and Pyke (1979) Prentice, R. L. and Pyke, R. (1979).  Logistic disease incidence models and case-control studies.  Biometrika, 66(3), 403–411.
  * Rao et al. (1973) Rao, C. R., Rao, C. R., Statistiker, M., Rao, C. R., and Rao, C. R. (1973).  Linear statistical inference and its applications, volume 2.  Wiley New York.
  * Rekaya et al. (2016) Rekaya, R., Smith, S., El Hamidi Hay, N. F., and Aggrey, S. E. (2016).  Analysis of binary responses with outcome-specific misclassification probability in genome-wide association studies.  The application of clinical genetics, 9, 169.
  * Robinson and Jewell (1991) Robinson, L. D. and Jewell, N. P. (1991).  Some surprising results about covariate adjustment in logistic regression models.  International Statistical Review/Revue Internationale de Statistique, pages 227–240.
  * Scarr and McCartney (1983) Scarr, S. and McCartney, K. (1983).  How people make their own environments: A theory of genotype→ environment effects.  Child development, pages 424–435.
  * Self and Mauritsen (1988) Self, S. G. and Mauritsen, R. H. (1988).  Power/sample size calculations for generalized linear models.  Biometrics, pages 79–86.
  * Self et al. (1992) Self, S. G., Mauritsen, R. H., and Ohara, J. (1992).  Power calculations for likelihood ratio tests in generalized linear models.  Biometrics, pages 31–39.
  * Serfling (2009) Serfling, R. J. (2009).  Approximation theorems of mathematical statistics, volume 162\.  John Wiley & Sons.
  * Shieh (2000) Shieh, G. (2000).  On power and sample size calculations for likelihood ratio tests in generalized linear models.  Biometrics, 56(4), 1192–1196.
  * Sjölander and Greenland (2013) Sjölander, A. and Greenland, S. (2013).  Ignoring the matching variables in cohort studies–when is it valid and why?  Statistics in medicine, 32(27), 4696–4708.
  * Sudlow et al. (2015) Sudlow, C., Gallacher, J., Allen, N., Beral, V., Burton, P., Danesh, J., Downey, P., Elliott, P., Green, J., Landray, M., et al. (2015).  Uk biobank: an open access resource for identifying the causes of a wide range of complex diseases of middle and old age.  PLoS medicine, 12(3), e1001779.
  * Sun and Bull (2005) Sun, L. and Bull, S. B. (2005).  Reduction of selection bias in genomewide studies by resampling.  Genetic Epidemiology: The Official Publication of the International Genetic Epidemiology Society, 28(4), 352–367.
  * Sun et al. (2011) Sun, L., Dimitromanolakis, A., Faye, L. L., Paterson, A. D., Waggott, D., and Bull, S. B. (2011).  Br-squared: a practical solution to the winner’s curse in genome-wide scans.  Human genetics, 129(5), 545–552.
  * Turley et al. (2018) Turley, P., Walters, R. K., Maghzian, O., Okbay, A., Lee, J. J., Fontana, M. A., Nguyen-Viet, T. A., Wedow, R., Zacher, M., Furlotte, N. A., et al. (2018).  Multi-trait analysis of genome-wide association summary statistics using mtag.  Nature genetics, 50(2), 229–237.
  * Veall and Zimmermann (1994) Veall, M. R. and Zimmermann, K. F. (1994).  Evaluating pseudo-r 2’s for binary probit models.  Quality and Quantity, 28(2), 151–164.
  * Wang et al. (2022) Wang, Z., Sun, L., and Paterson, A. D. (2022).  Major sex differences in allele frequencies for x chromosomal variants in both the 1000 genomes project and gnomad.  PLoS genetics, 18(5), e1010231.
  * Weissbrod et al. (2018) Weissbrod, O., Flint, J., and Rosset, S. (2018).  Estimating snp-based heritability and genetic correlation in case-control studies directly and with summary statistics.  The American Journal of Human Genetics, 103(1), 89–99.
  * Whittemore (1981) Whittemore, A. S. (1981).  Sample size for logistic regression with small response probability.  Journal of the American Statistical Association, 76(373), 27–32.
  * Yang et al. (2010) Yang, J., Benyamin, B., McEvoy, B. P., Gordon, S., Henders, A. K., Nyholt, D. R., Madden, P. A., Heath, A. C., Martin, N. G., Montgomery, G. W., et al. (2010).  Common snps explain a large proportion of the heritability for human height.  Nature genetics, 42(7), 565–569.
  * Yang et al. (2017) Yang, J., Zeng, J., Goddard, M. E., Wray, N. R., and Visscher, P. M. (2017).  Concepts, estimation and interpretation of snp-based heritability.  Nature genetics, 49(9), 1304–1310.
  * Zhang and Yi (2020) Zhang, Q. and Yi, G. Y. (2020).  Genetic association studies with bivariate mixed responses subject to measurement error and misclassification.  Statistics in Medicine, 39(26), 3700–3719.
  * Zhong and Prentice (2010) Zhong, H. and Prentice, R. L. (2010).  Correcting “winner’s curse” in odds ratios from genomewide association findings for major complex human diseases.  Genetic Epidemiology: The Official Publication of the International Genetic Epidemiology Society, 34(1), 78–91.
  * Zhu et al. (2018) Zhu, Z., Zheng, Z., Zhang, F., Wu, Y., Trzaskowski, M., Maier, R., Robinson, M. R., McGrath, J. J., Visscher, P. M., Wray, N. R., et al. (2018).  Causal associations between risk factors and common diseases inferred from gwas summary data.  Nature communications, 9(1), 1–12.
  * Zöllner and Pritchard (2007) Zöllner, S. and Pritchard, J. K. (2007).  Overcoming the winner’s curse: estimating penetrance parameters from case-control data.  The American Journal of Human Genetics, 80(4), 605–615.

(a) S1: Retrospective without covariate (b) S1: Retrospective without
covariate (c) S2: Prospective, binary covariate (d) S2: Prospective, binary
covariate

(e) S3: Prospective, continuous covariate (f) S3: Prospective, continuous
covariate

Figure 1: Simulation results for the three scenarios considered in Section
4.1. Scenario 1 (S1) is the retrospective case-control sampling design without
$E$. Scenario 2 (S2) and Scenario 3 (S3) are the prospective sampling design
with, respectively, a binary and continuous covariate $E$. Figures (a), (c)
and (e) on the left panel compare the power computation when $\beta_{G}$ is
fixed at log(1.5) (i.e. OR of 1.5, for (a-c)) or log(1.3) (for (e)), and
Figures (b), (d) and (f) on the right panel compare the sample size
computation to achieve power of 80% at the significance level of 0.05 across
different $\exp(\beta_{G})$. The red curves are for the ‘semi-simulation’
method (P1.SS in Section 3.2), blue curves for the ‘representative dataset’
method (P2.RD in Section 3.3), purple curves for Quanto of Gauderman (2002b),
and green and pink curves for the method of Demidenko (2007); in S3, the
method of Demidenko (2007) was implemented by dichotomizing $E$ or without
considering $E$). The black cures represent the oracle power and replication
sample size estimated empirically.

| Scenario 1 (S1) | Scenario 2 (S2) | Scenario 3 (S3)
---|---|---|---
Methods |  Average AE |  Maximum AE |  Average AE |  Maximum AE |  Average AE |  Maximum AE
P1.SS |  0.008 |  0.022 |  0.010 |  0.021 |  0.015 |  0.032
P2.RD |  0.009 |  0.019 |  0.012 |  0.023 |  0.012 |  0.027
Demidenko |  0.009 |  0.019 |  0.012 |  0.024 |  0.124 |  0.194
|  |  |  |  |  0.097* |  0.150*
Quanto |  0.009 |  0.022 |  0.052 |  0.084 |  0.112 |  0.173

Table 1: The average and maximum absolute error (AE), across different sample
sizes, between the oracle and computed power using the different methods for
the three scenarios considered. P1.SS is the proposed ‘semi-simulation’ method
in Section 3.2, and P2.RD is the proposed ‘representative dataset’ method in
Section 3.3. In Scenario 3, the method of Demidenko (2007) was implemented by
dichotomizing $E$ or without considering $E$ (results*). See legend to Figure
1 for additional details.

(a) Runtime comparison (b) Stability of P1:SS

Figure 2: Figure (a) displays the relationship between the run time in seconds
per computation for each method across different sample sizes, using
simulation Scenario 2; results for the other two scenarios are
characteristically similar. The Proposed 1 is the ‘semi-simulation’ method
(P1:SS) in Section 3.2, The proposed 2 is the ‘representative dataset’ method
(P2:RD) in Section 3.3. Figure (b) shows that there is a linear relationship
between the $\log_{10}$ standard error (SE) of estimated power and the
$\log_{10}$ number of replicates (B) used for the proposed ‘semi-simulation’
method 1.

(a) All UKB individuals: PC2 vs PC1 (b) Self-reported Africans: PC2 vs PC1 (c)
Self-reported Africans: PC3 vs PC1 (d) Self-reported Africans: PC3 vs PC2

Figure 3: Population principle components plots for (a) the whole UK Biobank
sample, stratified by $n=$ 3,460 self-reported African sample vs. Others, and
(b)–(d) the self-reported African sample. In Figures (b)–(d), the four
clusters were identified by a K-mean algorithm as discussed in Section 5. The
GWAS shown in Figure 4 used the $n=$ 2,510 individuals in Cluster 1,
identified based on this PCA analysis.

(a) The (binary) hypertension trait

(b) The (continuous) blood pressure trait

Figure 4: GWAS Manhattan results for (a) the binary hypertension trait and (b)
the continuous diastolic blood pressure trait, using the African sample ($n=$
2,510) identified through a PCA analysis of the self-identified African sample
of the UK Biobank data as discussed in Section 5. The association analyses
included age and sex as important covariates. No SNPs reached genome-wide
significance level of 5e-8.

(a) Overestimated power for the (binary) hypertension trait if not accounting
for E (b) Same power for the (continuous) blood pressure trait (c)
Underestimated (Discovery) sample size for the (binary) hypertension trait if
not accounting for E (d) Same (Discovery) sample size for the (continuous)
blood pressure trait

(e) Underestimated (Replication) sample size for the (binary) hypertension
trait if not accounting for E (f) Same (Replication) sample size for the
(continuous) blood pressure trait

Figure 5: Powers (Figures (a) and (b)) and sample sizes (Figures (c)-(e))
estimation for study planning of the top five-ranked SNPs identified in GWAS
of the binary hypertension trait (Figures (a), (c), (e)) and the continuous
diastolic blood pressure trait (Figures (b), (d) and (f)), using the African
sample ($n=$ 2,510) identified through a PCA analysis of the self-identified
African sample of the UK Biobank data as discussed in Section 5. The genetic
effects of these SNPs used for power and sample size computations are based on
a standard GWAS, shown in Figure 4, where age and sex were included as
important covariates. For (replication) study planning, the red bars are the
computed power or sample size with adjustment for age and sex, and the blues
bars are the values without explicitly considering age and sex. The two
approaches do not have difference in power and sample size planning for the
continuous blood pressure trait, as age and sex effects are incorporated
through residual variance. In contrast, when analyzing a binary trait, the
higher blue bars in Figure (a) show that ignoring covariate effects leads to
overestimated power of our discovery study (at $\alpha=$ 5e-8); power for
$\alpha=$ 0.05 is close to 100% as expected, thus not shown. The shorter blue
bars in Figure (c) show that ignoring covariate effects leads to
underestimated discovery sample size (for 80% power at $\alpha=$ 5e-8).
Similarly, the shorter blue bars in Figure (e) show that ignoring covariate
effects leads to underestimated replication sample size (for 80% power at
$\alpha=$ 0.05) when planning a replication study of a binary trait.


<EMAIL_ADDRESS>
# GRAVITY acquisition camera: characterization results

Narsireddy Anugu*a  Paulo Garciaa  Antonio Amorimb  Erich Wiezorrekc  Ekkehard
Wieprechtc  Frank Eisenhauerc  Thomas Ottc  Oliver Pfuhlc  Paulo Gordob  Guy
Perrind  Wolfgang Brandnere  Christian Straubmeierf and Karine Perrautg
aUniversidade do Porto  Faculdade de Engenharia  Unidade FCT CENTRA  Rua Dr.
Roberto Frias  s/n  P-4200-465 Porto  Portugal; bCENTRA-SIM and FCUL -
Edificio C8  Campo Grande  P1749016 Lisboa  Portugal;
cMax Planck Institute for extraterrestrial Physics  PO Box 1312
Giessenbachstr  85741 Garching  Germany;
dObservatoire de Paris Meudon  Paris  France;
eMaxPlanck-Institut fur Astronomie  Heidelberg  Germany;
fUniversität zu Köln Germany;
gInstitut de Planétologie et d´Astrophysique de Grenoble (IPAG)  Grenoble
France

###### Abstract

GRAVITY acquisition camera implements four optical functions to track multiple
beams of Very Large Telescope Interferometer (VLTI): a) pupil tracker: a
$2\times 2$ lenslet images four pupil reference lasers mounted on the spiders
of telescope secondary mirror; b) field tracker: images science object; c)
pupil imager: reimages telescope pupil; d) aberration tracker: images a Shack-
Hartmann. The estimation of beam stabilization parameters from the acquisition
camera detector image is carried out, for every 0.7 s, with a dedicated data
reduction software. The measured parameters are used in: a) alignment of
GRAVITY with the VLTI; b) active pupil and field stabilization; c) defocus
correction and engineering purposes. The instrument is now successfully
operational on-sky in closed loop. The relevant data reduction and on-sky
characterization results are reported.

###### keywords:

optical and infrared interferometry, GRAVITY, wavefront sensing, adaptive
optics, pupil tracking

##  1 INTRODUCTION

GRAVITY is a high precision narrow angle astrometry and interferometric
imaging instrument [1]. It has been built for the Very Large Telescope
Interferometer (VLTI) of the European Southern Observatory and was born with
the goal of monitoring the stellar sources in the vicinity of the Galactic
Center super massive black hole and the actual emission of the black hole
close environment environment. It combines four beams (in the K-band)
coherently of either Unit Telescopes (UTs) or Auxiliary Telescopes (ATs),
delivering an astrometry of $\sim 10$ micro-arcseconds ($\mu$as) and a angular
resolution of around $\sim 4$ milli-arcseconds (mas).

The main error sources for the GRAVITY astrometric measurements are the
atmospheric turbulence and the error in the baseline length ($B$) between any
two telescopes. The GRAVITY-Coudé Infrared Adaptive Optics (CIAO) measures and
corrects wavefront aberrations induced by the atmospheric turbulence in the
light path from the sky to the Coudé laboratory [2]. However the CIAO
corrections do not include the VLTI tunnel seeing induced fast moving tip-
tilts. These are tracked by launching external lasers at the telescopes (at
the star separator) and imaging them in the GRAVITY beam combiner [3]. The
measured tip-tilts are corrected by dedicated tip-tilt-piston mirrors. But the
above corrections involve tip-tilt residuals of more than 10 mas. These tip-
tilt errors limit the performance of GRAVITY in two ways [4]: a) efficiency of
star light injection into the single mode fibers which are used to transport
the light to the integrated optics [5] to make the beam combination; b)
unwanted tip-tilt error causes additional astrometric error.

The factor limiting the precision in measuring the accurate baseline (RMS
error of a few mm for a $\sim 100$ m baseline) length between two telescopes
is pupil position errors. These errors are mainly due to the misalignments and
temperature gradients. In between the telescope pupil and the fiber-fed beam
combiner, there exists several optical pupils that are under motion due to the
optical train vibrations while tracking the object of interest. The typical
lateral and longitudinal pupil position errors experienced by the VLTI are
around 5% in the pupil diameter and 1 m, respectively, for the 80 mm beam. Eq.
1 gives the astrometric error [6] associated with the given tip-tilt error
($\Delta\alpha$) and pupil lateral position ($\Delta L$) errors. For example,
for 10 mas tip-tilt error and 0.4 m (5% of 8 m pupil) lateral pupil position
error, the astrometric error becomes 40 $\mu$as for a single beam combiner

$\sigma=\dfrac{\Delta\alpha\times\Delta L}{B}.$ (1)

The CIAO measures and corrects the incoming distorted wavefronts with respect
to a flat wavefront generated in the Coudé laboratory. However, there are many
optics and tunnel seeing in between the CIAO and the GRAVITY beam combiner.
They will also introduce additional wavefront aberrations to the incoming
beam. The quasi-static wavefront aberrations and the non-common path errors
that exist between the CIAO and the GRAVITY beam combiner also contribute to
the astrometric error [4].

Therefore a beam stabilization system is required to address the above issues.
The GRAVITY acquisition camera has been built to image and analyze the VLTI
beams and to: a) enable the stabilization of star light injection into the
single mode fibers; b) minimize the astrometric errors induced by the field
and the pupil errors.

##  2 Acquisition camera

Field lensPupil trackerAbr. sensorPupil imagerField imager$2\times
2$lenslet$9\times 9$lensletPupil lensAcquisition cameraFiber fed TTPPMCK-band
beam H-band beam Ref. laser Starbeam combinerVCMM2M1M4M3 Figure 1:  The
acquisition camera working principle and optical layout (for clarity only one
telescope beam case is shown). The astrophysical target’s K-band beam is used
for science measurements and the H-band beam is used for beam stabilization.
The H-band beam feeds the three imaging modes of the acquisition camera:
aberration sensor, pupil imager, field imager. The pupil tracker images
external pupil reference laser light
($1.2\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$). The box in red color
encloses the acquisition camera optics.

Figure 1 presents the optical layout of the acquisition camera for one
telescope beam. The acquisition camera provides four optical imaging modes: a)
field imager – images and tracks the field (tip-tilt) of input beams; b) pupil
tracker – tracks the telescope pupil motions by imaging four pupil beacons
that are installed on the spiders of secondary mirror (M2) of the telescope;
c) aberration sensor – implements four $9\times 9$ lenslets to sense the
incoming distorted wavefronts; d) pupil imager – images telescope pupil, which
is used for monitoring the pupil visually during the observations. Pupil
tracker images four pupil guiding laser beacons of
$1.2\text{\,}\mathrm{\SIUnitSymbolMicro}\mathrm{m}$ wavelength with a $2\times
2$ lenslet. The other three imaging modes image astrophysical targets in the
H-band. The acquisition camera uses a $2048\times 2048$ pixel Hawaii-2RG
detector to image all four input telescope beams. The optical design of the
camera [7] and its optical alignment and integration [8] are presented in the
previous SPIE proceedings.

To extract the beam stabilization parameters from the detector images, a
dedicated software has been developed based on the image analysis methods
presented in Anugu et al. (2014) [9]. The software was written in C and C++
using CLIP library [10] and integrated in the GRAVITY observational software
[11]. This software works on the instrument workstation and does the image
analysis in two steps. First, it reads the acquisition camera detector image
from the instrument shared memory for every 0.7 s (detector frame rate).
Second, it evaluates the beam’s tip-tilts, pupil shifts and the wavefront
aberrations by analyzing the input detector image. The evaluated parameters
are written to the instrument database. They feed the beam stabilization which
uses: a) the tip-tilt and piston (TTP); b)the pupil motion controller (PMC)
[3] and; c) the Variable Curvature Mirror (VCM) [12]. By stabilizing the
field, it enables the injection of the K-band light of the astrophysical
targets into single mode fibers which transport the beams towards the
interferometric beam combiner. The image acquisition, analysis and beam
correction are carried out on-line in the instrument workstation in parallel
to the observations for all four telescopes.

##  3 Characterization results

GRAVITY passed Preliminary Verification Europe tests (PAE) in early 2015 and
was moved in August 2015 to the Paranal Observatory where it got its first
light in November 2015. Figure 2 presents the acquisition camera detector
image acquired on-sky with the ATs.

meter

---

meter

arc-sec

---

arc-sec

Figure 2:  On-sky measurements of the acquisition camera imaging modes (left
to right): pupil tracker, pupil imager, field tracker and aberration sensor.
The astrophysical target is $\theta^{1}$ Ori A (North up and East on right).
The aberration sensor image is rotated counter clockwise to the field image
due to the mirror reflection.

###  3.1 Characterization for PAE

During the PAE, the acquisition camera was characterized using the GRAVITY
calibration unit [13]. The calibration unit is a GRAVITY subsystem which
provides laboratory beams and artificial stars to test GRAVITY
functionalities. It generates two artificial stars and four pupil guiding
reference beams to test the acquisition camera imaging modes.

####  3.1.1 Field tracking

Characterization of the field tracker is accomplished in three fronts: a) RMS
and absolute tracking accuracy; b) RMS accuracy as a function of target
magnitude; c) target flux injection into the fibers. To characterize the
absolute accuracy of the field tracker, known tip-tilts ($\theta_{\rm i}$) are
applied to the incoming beams by manipulating the TTP device and measured
backm ($\theta_{\rm o}$), using the field tracker function.

Figure 3 (left panel) presents the absolute tracking accuracy and is around
$\sim 2$ mas and the RMS error is $\sim 1$ mas. These measurements are carried
out with a ”star” with an H-band magnitude of 15. Figure 3 (middle panel)
presents the RMS error of the field tracker as a function of ”star” magnitude.
These ”stars” are realized by the Calibration Unit and by varying the voltage
of the lamp. Figure 3 (right panel) presents how the flux injection in the
fiber is reduced. By enabling the field stabilization, the coupling efficiency
of the fiber is maintained around 80%.

$|\theta_{\rm i}$ \- $\theta_{\rm o}|$ (mas)

---

$\theta_{\rm i}$ (arcsec)

RMS error (mas)

Normalized flux in the fiber
---

Target magnitude

Star shift from the fiber center (arcsec)

Figure 3:  Performance of the field tracker in the laboratory. Left: the
absolute field tracking error $|\theta_{\rm i}-\theta_{\rm o}|$ as a function
of true object position $\theta_{\rm i}$. Middle: the RMS error of the field
tracking as a function of target magnitude. Right: Normalized star’s flux
injection into the fiber as function of tip-tilt/star shift from the fiber
center.

####  3.1.2 Lateral pupil tracking

The lateral pupil tracking accuracy tests are implemented in two steps by
manipulating the PMC actuators. First, known lateral pupil shifts
($L_{\rm{x_{0}}}$) are applied to the incoming beams by actuating the PMCs.
Second, the input pupil shifts are sensed (say, $L_{\rm x}$) using the pupil
tracker function. Figure 4 presents the lateral pupil tracking accuracy as a
function of the input pupil shifts. The lateral pupil tracking absolute
accuracy is better than $4\text{\,}\mathrm{mm}$ at the UT beam magnification
(0.05% of the UT diameter) and the RMS error is $2\text{\,}\mathrm{mm}$.

$|L_{\rm{x_{0}}}-L_{\rm x}|$ (mm)

---

$L_{\rm{x_{0}}}$ (m)

Figure 4:  The absolute error in the lateral pupil position measurement
$|L_{\rm x}-L_{\rm{x_{0}}}|$ as a function of input lateral pupil shift
$L_{x_{0}}$. The measurements are in the UT scale and where obtained in the
laboratory.

####  3.1.3 Longitudinal pupil tracking

Longitudinal pupil tracking precision is characterized by manipulating the VCM
positions. In this experiment, known longitudinal pupil shifts are applied to
the beam by moving the VCM and those are measured back by the pupil tracker.
Figure 5 (left panel) presents the longitudinal pupil characterization results
and it can be seen that the accuracy is better than 40 mm for the 80 mm beam.
The middle and right panel present, respectively, the UT pupil before and
after closing the pupil guiding loop.

$|L_{\rm z}-L_{\rm{z_{0}}}|$ (mm)

meter |
---

meter

---

$L_{\rm{z_{0}}}$ (meter)

meter

meter

Figure 5: Left: On-sky absolute longitudinal pupil position error measurement
($|L_{\rm z}-L_{\rm{z_{0}}}|$) at 80 mm as a function of input longitudinal
pupil position $L_{z_{0}}$. Center, right: On-sky pupil images obtained before
and after closing the pupil guiding closed loop.

###  3.2 On-sky beam guiding performance

The acquisition camera beam guiding has been tested for the ATs and for the
UTs. The guiding difference between the ATs and the UTs is that the ATs are
not corrected by CIAO. For this reason, the field imager is used to measure
slow (tunnel) atmospheric tip-tilts. Whereas for the UTs, the field imager is
used to measure the residual tip-tilts of the CIAO. The ATs telescope focus
correction is implemented by tracking the focus using the aberration sensor
and correcting it with their M2 mirrors.

As stated previously, the GRAVITY field and pupil offset corrections are
achieved using two types of actuators: a) the VLTI star separator actuators
[14] and the delay line VCM; b) GRAVITY internal field (TTP) and pupil (PMC)
actuators. The VLTI star separator has a field selector mirror and a pupil
position mirror (M14, laterally and longitudinally movable). Small field and
pupil corrections occurring during interferometric observations are corrected
in closed loop using GRAVITY internal actuators for the speed and accuracy.
Large field and pupil offsets (which usually take place during the initial
alignment of GRAVITY with the VLTI) are corrected using star separator
actuators. While operating with GRAVITY internal actuators in the closed loop,
when the offsets are larger than GRAVITY internal actuators range, they are
offloaded to the VLTI actuators.

During the installation, verification and commissioning phases, the guiding
loops of the field, the pupil were extensively tested. Characterization of the
field guiding is realized by observing several dual-field stars and checking
the stability of injection of star’s light into single mode fibers.
Characterization of the pupil guiding is implemented by applying known pupil
shifts with pupil correcting actuators and checking their stability over an
hour. The pupil tracker experiences high backgrounds from bright astrophysical
targets due to the closeness of: a) the operating wavelengths of the pupil
tracker (operating at 1.2 $\mu$m) and the field tracker (operates at $H$-band
– 1.65 $\mu$m) and; b) no adequate dichroic filter for the pupil tracker. The
operation of the pupil guiding for an interferometric observation depends on
two factors: a) the flux of the pupil guiding laser beacons; b) the magnitude
of the astrophysical target. To remove the astrophysical target background, a
new software mode, BLINK, is implemented. The basic idea of this mode is to
remove the background by switching OFF and ON the pupil beacons. During the
turn OFF period, the background is stored and the is subtracted from the
subsequent frame. The performance of the pupil guiding is improved by 2
magnitudes when using this BLINK mode. Currently it can operate up to
magnitude 1 for the ATs.

The non-common path errors between the CIAO and GRAVITY are calibrated by
inputting known wavefronts to the GRAVITY beams that are generated by
manipulating the MACAO deformable mirror [15] and measuring the input
wavefront aberrations using the aberration sensor. The offsets are then taken
into account by CIAO to compensate the non-common path aberrations.

The following summarizes the current on-sky performance of the acquisition
camera:

  * •

Pupil guiding residuals: The standard deviation of the lateral and the
longitudinal pupil guiding residuals are remain smaller than $\pm~{}0.2$
pixels RMS on the detector (i.e., 12 mm at the UTs) and the 50 mm for 80 beam
size respectively in the presence of faint astrophysical target.

  * •

Field guiding residuals: The standard deviation of the field guiding residuals
are smaller than $\pm~{}0.36$ pixels (6.4 mas). When the tests of the field
guiding for the ATs (not equipped with the CIAO) took place, the laser based
tunnel seeing tip-tilt tracking system was not installed. That is why the
residuals are large.

  * •

Focus guiding The standard deviation of the focus guiding residuals for the
ATs are within $\pm~{}\lambda/8$ RMS.

##  4 Discussion and conclusions

The GRAVITY acquisition camera was installed at the Very Large Telescope
Interferometer during the last months of 2015. Since then it has been
successfully working in the closed loop to stabilize the slow (tunnel)
atmospheric turbulence field motions and the pupil motions.

The camera is characterized using laboratory generated telescope beams at MPE
and on-sky at the Paranal Observatory. The characterization results revealed
that it is able to analyze the telescope beams with the accuracy required by
$10\,\mu$as astrometry, namely: a) field tracking with 2 mas RMS; b) lateral
pupil tracking with 4 mm RMS (at the UT scale); c) longitudinal pupil tracking
with 50 mm RMS (at 80 mm beam scale); and d) quasi-static higher order
wavefront aberration measurements with 80 nm RMS.

The acquisition camera measured beam parameters are used in the stabilization
of the field, the pupil and the focus correction. With the stabilization of
the field, the star’s light injection into single mode fibers is improved. The
quasi-static higher order aberrations measurements are used for the non-common
path errors corrections between the CIAO and GRAVITY. The pupil imager is used
for visual monitoring. By stabilizing the telescope beams, the astrometric
error induced by the field (accuracy of 2 mas RMS, cf. Section 3.1.1) and the
pupil errors (accuracy of 12 mm RMS) is minimized to $0.34\,\mu$as).

## Acknowledgments

Anugu acknowledges FCT-Portugal grant SFRH/BD/52066/2012. The research leading
to these results was partially supported by FCT-Portugal grant PTDC/CTE-
AST/116561/2010 and the European Community’s Seventh Framework Programme under
Grant Agreement 312430.

## References

  * [1] Eisenhauer, F., Brandner, W., Straubmeier, C., Perraut, K., Gillessen, S., Lima, J., Henning, T., Garcia, P., Kellner, S., Avila, G., Perrin, G., Brandner, W., Straubmeier, C., Perraut, K., Amorim, A., Schöller, M., Gillessen, S., Kervella, P., Benisty, M., Araujo-Hauck, C., Jocou, L., Lima, J., Jakob, G., Haug, M., Clénet, Y., Henning, T., Eckart, A., Berger, J.-P., Garcia, P., Abuter, R., Kellner, S., Paumard, T., Hippler, S., Fischer, S., Moulin, T., Villate, J., Avila, G., Gräter, A., Lacour, S., Huber, A., Wiest, M., Nolot, A., Carvas, P., Dorn, R., Pfuhl, O., Gendron, E., Kendrew, S., Yazici, S., Anton, S., Jung, Y., Thiel, M., Choquet, É., Klein, R., Teixeira, P., Gitton, P., Moch, D., Vincent, F., Kudryavtseva, N., Ströbele, S., Sturm, S., Fédou, P., Lenzen, R., Jolley, P., Kister, C., Lapeyrère, V., Naranjo, V., Lucuix, C., Hofmann, R., Chapron, F., Neumann, U., Mehrgan, L., Hans, O., Rousset, G., Ramos, J., Suarez, M., Lederer, R., Reess, J.-M., Rohloff, R.-R., Haguenauer, P., Bartko, H., Sevin, A., Wagner, K., Lizon, J.-L., Rabien, S., Collin, C., Finger, G., Davies, R., Rouan, D., Wittkowski, M., Dodds-Eden, K., Ziegler, D., Cassaing, F., Bonnet, H., Casali, M., Genzel, R., and Lena, P., “GRAVITY: Observing the Universe in Motion,” The Messenger 143(March), 16–24 (2011).
  * [2] Kendrew, S., Hippler, S., Brandner, W., Clénet, Y., Deen, C., Gendron, E., Huber, A., Klein, R., Laun, W., Lenzen, R., Naranjo, V., Neumann, U., Ramos, J., Rohloff, R.-R., Yang, P., Eisenhauer, F., Fedrigo, E., Suarez-Valles, M., Amorim, A., Perraut, K., Perrin, G., and Straubmeier, C., “The GRAVITY Coude Infrared Adaptive Optics (CIAO) system for the VLT Interferometer,” eprint arXiv:1207.2945 (2012).
  * [3] Pfuhl, O., Haug, M., Eisenhauer, F., Kellner, S., Haussmann, F., Perrin, G., Gillessen, S., Straubmeier, C., Ott, T., Rousselet-Perraut, K., Amorim, A., Lippa, M., Janssen, A., Brandner, W., Kok, Y., Blind, N., Burtscher, L., Sturm, E., Wieprecht, E., Schoeller, M., Weber, J., Hans, O., and Huber, S., “The fiber coupler and beam stabilization system of the GRAVITY interferometer,” in [SPIE Astronomical Telescopes + Instrumentation ], Rajagopal, J. K., Creech-Eakman, M. J., and Malbet, F., eds., 914623, International Society for Optics and Photonics (jul 2014).
  * [4] Lacour, S., Eisenhauer, F., Gillessen, S., Pfuhl, O., Woillez, J., Bonnet, H., Perrin, G., Lazareff, B., Rabien, S., Lapeyrere, V., Clenet, Y., Kervella, P., and Kok, Y., “Reaching micro-arcsecond astrometry with long baseline optical interferometry; application to the GRAVITY instrument,” Astronomy & Astrophysics (apr 2014).
  * [5] Jocou, L., Perraut, K., Nolot, A., Berger, J. P., Moulin, T., Labeye, P., Lacour, S., Perrin, G., Lebouquin, J. B., Bartko, H., Thiel, M., and Eisenhauer, F., “The GRAVITY integrated optics beam combination,” in [Optical and Infrared Interferometry II. Edited by Danchi ], 7734, 773430–773430–12 (2010).
  * [6] Lacour, S. and Bonnet, H., “GRAVITY Astrometric Error Budget,” VLT-TRE-GRA-15882-3201 (2011).
  * [7] Amorim, A., Lima, J., Anugu, N., Eisenhauer, F., Graeter, A., Haug, M., Ott, T., Pfuhl, O., Thiel, M., Wieprecht, E., Carvas, P., Garcia, P., Perrin, G., Brandner, W., Straubmeier, C., and Perraut, K., “The final design of the GRAVITY acquisition camera and associated VLTI beam monitoring strategy,” in [Optical and Infrared Interferometry III. Proceedings of the SPIE, Volume 8445, article id. 844534, 14 pp. (2012). ], Delplancke, F., Rajagopal, J. K., and Malbet, F., eds., 8445, 844534 (sep 2012).
  * [8] Gordo, P., Amorim, A., Abreu, J., Eisenhauer, F., Anugu, N., Garcia, P., Pfuhl, O., Haug, M., Sturm, E., Wieprecht, E., Perrin, G., Brandner, W., Straubmeier, C., Perraut, K., Naia, M. D., and Guimarães, M., “Integration and testing of the GRAVITY infrared camera for multiple telescope optical beam analysis,” in [Proceedings of the SPIE, Volume 9146, id. 91461V 8 pp. (2014). ], Rajagopal, J. K., Creech-Eakman, M. J., and Malbet, F., eds., 9146, 91461V (jul 2014).
  * [9] Anugu, N., Garcia, P. J. V., Wieprecht, E., Amorim, A., Burtscher, L., Ott, T., Gordo, P., Eisenhauer, F., Perrin, G., Brandner, W., Straubmeier, C., and Perraut, K., “The GRAVITY/VLTI acquisition camera software,” in [Proceedings of the SPIE, Volume 9146, id. 91462C 13 pp. (2014). ], Rajagopal, J. K., Creech-Eakman, M. J., and Malbet, F., eds., 9146, 91462C (jul 2014).
  * [10] Ballester, P., Biereichel, P., Kaufer, A., Kiekebusch, M., and Lorch, H., “CLIP: Bridging Pipelines to Instrument Control Software,” Astronomical Data Analysis Software and Systems ASP Conference Series 394 (2008).
  * [11] Burtscher, L., Wieprecht, E., Ott, T., Kok, Y., Yazici, S., Anugu, N., Dembet, R., Fedou, P., Lacour, S., Ott, J., Others, Paumard, T., Lapeyrere, V., Kervella, P., Abuter, R., Pozna, E., Eisenhauer, F., Blind, N., Genzel, R., Gillessen, S., Hans, O., Haug, M., Haussmann, F., Kellner, S., Lippa, M., Pfuhl, O., Sturm, E., Weber, J., Amorim, A., Brandner, W., Rousselet-Perraut, K., Perrin, G. S., Straubmeier, C., and Schöller, M., “The GRAVITY instrument software/high-level software,” in [SPIE Astronomical Telescopes+ Instrumentation ], Rajagopal, J. K., Creech-Eakman, M. J., and Malbet, F., eds., 91462B—-91462B, International Society for Optics and Photonics, International Society for Optics and Photonics (jul 2014).
  * [12] Ferrari, M., Lemaitre, G. R., Mazzanti, S. P., Derie, F., Huxley, A., Lemerrer, J., Lanzoni, P., Dargent, P., and Wallander, A., “Variable curvature mirrors: implementation in the VLTI delay-lines for field compensation,” in [Interferometry for Optical Astronomy II ], 4838, 1155–1162 (2003).
  * [13] Blind, N., Eisenhauer, F., Haug, M., Gillessen, S., Lippa, M., Burtscher, L., Hans, O., Haussmann, F., Huber, S., Janssen, A., Kellner, S., Kok, Y., Ott, T., Pfuhl, O., Sturm, E., Weber, J., Wieprecht, E., Amorim, A., Brandner, W., Perrin, G., Perraut, K., and Straubmeier, C., “GRAVITY: the Calibration Unit,” Optical and Infrared Interferometry IV 9146, 12 (jul 2014).
  * [14] Nijenhuis, J., Visser, H., de Man, H., Dekker, B., Mekking, J., and Kamphues, F., “Simultaneous observation of two stars using the PRIMA Star Separator,” in [Optical and Infrared Interferometry ], 7013, 70133F (jul 2008).
  * [15] Arsenault, R., Alonso, J., Bonnet, H., Brynnel, J., Delabre, B., Donaldson, R., Dupuy, C., Fedrigo, E., Spyromilio, J., Erm, T., Farinato, J., Hubin, N., Ivanescu, L., Kasper, M., Oberti, S., Paufique, J., Rossi, S., Tordo, S., Stroebele, S., Lizon, J.-L., Gigan, P., Pouplard, F., Delplancke, F., Silber, A., Quattri, M., and Reiss, R., “MACAO-VLTI first light: adaptive optics at the service of interferometry,” The Messenger (ISSN 0722-6691) (2003).