Combining controls can improve power in two-stage association studies

James Liley

arXiv:1702.02827·stat.ME·July 11, 2017

Combining controls can improve power in two-stage association studies

James Liley

PDF

TL;DR

This paper introduces a shared-control method for two-stage association studies that enhances statistical power and reduces costs by reusing control samples across discovery and replication phases, with careful error control.

Contribution

The paper proposes a novel shared-control approach for two-stage studies, improving power and cost-efficiency while discussing its theoretical properties and practical applications.

Findings

01

Shared-control method increases power in GWAS.

02

Type-1 error rates are controlled both analytically and empirically.

03

Method is especially useful when control samples are hard to recruit.

Abstract

High dimensional case control studies are ubiquitous in the biological sciences, particularly genomics. To maximise power while constraining cost and to minimise type-1 error rates, researchers typically seek to replicate findings in a second experiment on independent cohorts before proceeding with further analyses. This paper presents a method in which control (or case) samples from the discovery cohort are re-used in the replication study. The theoretical implications of this method are discussed and simulations used to compare performance against the standard method in a range of circumstances. In several common study designs, a shared-control method allows a substantial improvement in power while retaining type-1 error rate control. An important area of potential application arises when control samples are difficult to recruit or ascertain; for example in inter-disease…

Tables1

Table 1. Table 1 : Upper bounds on type 1 error rates with aberrance in cohorts, with β > β ∗ > β ⟂ 𝛽 superscript 𝛽 superscript 𝛽 perpendicular-to \beta>\beta^{*}>\beta^{\perp}

	Aberrant
	None	$C_{0}$	$C_{0}^{'}$	$C_{1}$	$C_{1}^{'}$
M. A	$P_{0}$	$β$	$α$	$β$	$α$
M. B	$P_{0}$	1	$α$	$β^{*}$	$α$
M. C	$P_{0}$	1	1	$β^{⟂}$	$α$

Equations172

Σ_{A} = v a r ((z_{d} z_{r} z_{m})^{t}) Σ_{B} = v a r ((z_{d} z_{s} z_{m})^{t}) Σ_{C} = v a r ((z_{c} z_{s} z_{m})^{t})

Σ_{A} = v a r ((z_{d} z_{r} z_{m})^{t}) Σ_{B} = v a r ((z_{d} z_{s} z_{m})^{t}) Σ_{C} = v a r ((z_{c} z_{s} z_{m})^{t})

2 \int_{z_{α}}^{\infty} \int_{z_{β^{*}}}^{\infty} \int_{z_{γ}}^{\infty} N_{Σ_{B}} (z_{d} z_{r} z_{m}) d z_{m} d z_{r} d z_{d}

2 \int_{z_{α}}^{\infty} \int_{z_{β^{*}}}^{\infty} \int_{z_{γ}}^{\infty} N_{Σ_{B}} (z_{d} z_{r} z_{m}) d z_{m} d z_{r} d z_{d}

= 2 \int_{z_{α}}^{\infty} \int_{z_{β}}^{\infty} \int_{z_{γ}}^{\infty} N_{Σ_{A}} (z_{d} z_{r} z_{m}) d z_{m} d z_{r} d z_{d}

= P r (p_{d} < α, p_{r} < β, p_{m} < γ ∣ H_{0}^{\cap})

P r (p_{d} < α, p_{s} < β^{*} ∣ H_{0}^{\cap}) = P r (p_{c} < α, p_{s} < β^{⊥} ∣ H_{0}^{\cap}) = P r (p_{d} < α, p_{r} < β ∣ H_{0}^{\cap}) = α β

P r (p_{d} < α, p_{s} < β^{*} ∣ H_{0}^{\cap}) = P r (p_{c} < α, p_{s} < β^{⊥} ∣ H_{0}^{\cap}) = P r (p_{d} < α, p_{r} < β ∣ H_{0}^{\cap}) = α β

z_{α} \to \infty lim \frac{z _{β^{*}}}{1 - ρ _{d s}^{2} z _{β} + ρ _{d s} z _{α}} = 1 z_{α} \to \infty lim \frac{z _{β^{⊥}}}{1 - ρ _{cs}^{2} z _{β} + ρ _{cs} z _{α}} = 1

z_{α} \to \infty lim \frac{z _{β^{*}}}{1 - ρ _{d s}^{2} z _{β} + ρ _{d s} z _{α}} = 1 z_{α} \to \infty lim \frac{z _{β^{⊥}}}{1 - ρ _{cs}^{2} z _{β} + ρ _{cs} z _{α}} = 1

E (Z) = \frac{2 n _{0} n _{1}}{n _{0} + n _{1}} \frac{μ _{1} - μ _{0}}{μ ˉ ( 1 - μ ˉ )}

E (Z) = \frac{2 n _{0} n _{1}}{n _{0} + n _{1}} \frac{μ _{1} - μ _{0}}{μ ˉ ( 1 - μ ˉ )}

ζ_{d} = \frac{2 n _{0} n _{1}}{n _{0} + n _{1}} \frac{μ _{1} - μ _{0}}{μ ˉ ( 1 - μ ˉ )} ζ_{r} = \frac{2 n _{0}^{'} n _{1}^{'}}{n _{0}^{'} + n _{1}^{'}} \frac{μ _{1}^{'} - μ _{0}^{'}}{μ ˉ ( 1 - μ ˉ )}

ζ_{d} = \frac{2 n _{0} n _{1}}{n _{0} + n _{1}} \frac{μ _{1} - μ _{0}}{μ ˉ ( 1 - μ ˉ )} ζ_{r} = \frac{2 n _{0}^{'} n _{1}^{'}}{n _{0}^{'} + n _{1}^{'}} \frac{μ _{1}^{'} - μ _{0}^{'}}{μ ˉ ( 1 - μ ˉ )}

ζ_{s} = \frac{2 ( n _{0} + n _{0}^{'} ) n _{1}^{'}}{n _{0} + n _{0}^{'} + n _{1}^{'}} \frac{μ _{1}^{'} - \frac{μ _{0} n _{0} + μ _{0}^{'} n _{0}^{'}}{n _{0} + n _{0}^{'}}}{μ ˉ ( 1 - μ ˉ )} ζ_{c} = \frac{2 ( n _{0} + n _{0}^{'} ) n _{1}}{n _{0} + n _{0}^{'} + n _{1}} \frac{μ _{1} - \frac{μ _{0} n _{0} + μ _{0}^{'} n _{0}^{'}}{n _{0} + n _{0}^{'}}}{μ ˉ ( 1 - μ ˉ )}

ζ_{m} = \frac{2 ( n _{0} + n _{0}^{'} ) ( n _{1} + n _{1}^{'} )}{n _{0} + n _{0}^{'} + n _{1} + n _{1}^{'}} \frac{\frac{μ _{1} n _{1} + μ _{1}^{'} n _{1}^{'}}{n _{1} + n _{1}^{'}} - \frac{μ _{0} n _{0} + μ _{0}^{'} n _{0}^{'}}{n _{0} + n _{0}^{'}}}{μ ˉ ( 1 - μ ˉ )}

\int_{z_{α} - ζ_{d}}^{\infty} \int_{z_{β} - ζ_{r}}^{\infty} \int_{z_{γ} - ζ_{m}}^{\infty} N_{Σ_{A}} ((z_{d} z_{r} z_{m})^{t}) d z_{m} d z_{r} d z_{d}

\int_{z_{α} - ζ_{d}}^{\infty} \int_{z_{β} - ζ_{r}}^{\infty} \int_{z_{γ} - ζ_{m}}^{\infty} N_{Σ_{A}} ((z_{d} z_{r} z_{m})^{t}) d z_{m} d z_{r} d z_{d}

+ \int_{z_{α} + ζ_{d}}^{\infty} \int_{z_{β} + ζ_{r}}^{\infty} \int_{z_{γ} + ζ_{m}}^{\infty} N_{Σ_{A}} ((z_{d} z_{r} z_{m})^{t}) d z_{m} d z_{r} d z_{d}

\int_{z_{α} - ζ_{d}}^{\infty} \int_{z_{β} - ζ_{s}}^{\infty} \int_{z_{γ} - ζ_{m}}^{\infty} N_{Σ_{B}} ((z_{d} z_{s} z_{m})^{t}) d z_{m} d z_{s} d z_{d}

\int_{z_{α} - ζ_{d}}^{\infty} \int_{z_{β} - ζ_{s}}^{\infty} \int_{z_{γ} - ζ_{m}}^{\infty} N_{Σ_{B}} ((z_{d} z_{s} z_{m})^{t}) d z_{m} d z_{s} d z_{d}

+ \int_{z_{α} + ζ_{d}}^{\infty} \int_{z_{β} + ζ_{s}}^{\infty} \int_{z_{γ} + ζ_{m}}^{\infty} N_{Σ_{B}} ((z_{d} z_{s} z_{m})^{t}) d z_{m} d z_{s} d z_{d}

\frac{α}{2 2 π} (\frac{k}{1 - ρ ^{2}} - 1) z_{β} ≪ \frac{α}{2}

\frac{α}{2 2 π} (\frac{k}{1 - ρ ^{2}} - 1) z_{β} ≪ \frac{α}{2}

k = \frac{ζ _{s}}{ζ _{r}} \approx \frac{( n _{0} + n _{0}^{'} ) ( n _{0}^{'} + n _{1}^{'} )}{n _{0}^{'} ( n _{0} + n _{0}^{'} + n _{1}^{'} )}

k = \frac{ζ _{s}}{ζ _{r}} \approx \frac{( n _{0} + n _{0}^{'} ) ( n _{0}^{'} + n _{1}^{'} )}{n _{0}^{'} ( n _{0} + n _{0}^{'} + n _{1}^{'} )}

k^{'} = \frac{ζ _{s}}{ζ _{r}} \approx \frac{n _{0}^{'} ( n _{0}^{'} + n _{1}^{'} )}{( n _{0} + n _{0}^{'} ) ( n _{0} + n _{0}^{'} + n _{1}^{'} )}

k^{'} = \frac{ζ _{s}}{ζ _{r}} \approx \frac{n _{0}^{'} ( n _{0}^{'} + n _{1}^{'} )}{( n _{0} + n _{0}^{'} ) ( n _{0} + n _{0}^{'} + n _{1}^{'} )}

cor (m_{1 i} - m_{0 i}, m_{1 j} - m_{0 j}) = \frac{n _{0 i} n _{0 j} n _{1 ij} + n _{1 i} n _{1 j} n _{0 ij}}{n _{0 i} n _{0 j} n _{1 i} n _{1 j} \frac{1}{n _{0 i}} + \frac{1}{n _{0 j}} \frac{1}{n _{0 j}} + \frac{1}{n _{1 j}}}

cor (m_{1 i} - m_{0 i}, m_{1 j} - m_{0 j}) = \frac{n _{0 i} n _{0 j} n _{1 ij} + n _{1 i} n _{1 j} n _{0 ij}}{n _{0 i} n _{0 j} n _{1 i} n _{1 j} \frac{1}{n _{0 i}} + \frac{1}{n _{0 j}} \frac{1}{n _{0 j}} + \frac{1}{n _{1 j}}}

d e t (Σ_{A})

d e t (Σ_{A})

= \frac{( n _{0} n _{1}^{'} - n _{0}^{'} n _{1} ) ^{2}}{( n _{0} + n _{0}^{'} ) ( n _{1} + n _{1}^{'} ) ( n _{0} + n _{1} ) ( n _{0}^{'} + n _{1}^{'} )}

d e t (Σ_{B})

= \frac{n _{0}^{'} n _{1}^{2}}{( n _{0} + n _{1} ) ( n _{0} + n _{0}^{'} + n _{1}^{'} ) ( n _{1} + n _{1}^{'} )}

z_{i} = r \in 1.. s \sum α_{i r} (m_{1 i}^{r} - m_{0 i}^{r})

z_{i} = r \in 1.. s \sum α_{i r} (m_{1 i}^{r} - m_{0 i}^{r})

α_{i r} \propto \frac{n _{0 i}^{r} n _{1 i}^{r}}{n _{0 i}^{r} + n _{1 i}^{r}}

α_{i r} \propto \frac{n _{0 i}^{r} n _{1 i}^{r}}{n _{0 i}^{r} + n _{1 i}^{r}}

cor (z_{i}, z_{j})

cor (z_{i}, z_{j})

\approx \frac{\sum _{r \in 1.. s} α ( \frac{n _{0 ij}^{r}}{n _{0 i}^{r} n _{0 j}^{r}} + \frac{n _{1 ij}^{r}}{n _{1 i}^{r} n _{1 j}^{r}} )}{( \sum α _{i r}^{2} ( \frac{1}{n _{0 i}^{r}} + \frac{1}{n _{1 i}^{r}} ) ) ( \sum α _{j r}^{2} ( \frac{1}{n _{0 j}^{r}} + \frac{1}{n _{1 j}^{r}} ) )}

z \propto \frac{1}{∣ C _{1} ∣} k \in C_{1} \sum c_{1}^{k} g_{1}^{k} - \frac{1}{∣ C _{0} ∣} k \in C_{0} \sum c_{0}^{k} g_{0}^{k}

z \propto \frac{1}{∣ C _{1} ∣} k \in C_{1} \sum c_{1}^{k} g_{1}^{k} - \frac{1}{∣ C _{0} ∣} k \in C_{0} \sum c_{0}^{k} g_{0}^{k}

E (z)

E (z)

∣ C_{0} ∣, ∣ C_{1} ∣ \to \infty lim E (z)

\propto \int (c_{1} (x) f_{1} (x) - c_{0} (x) f_{0} (x)) E (g ∣ x) d x

c_{1} (x) f_{1} (x) - c_{0} (x) f_{0} (x) \equiv 0 ⟹ c_{1} (x) \propto \frac{f ( x )}{f _{1} ( x )}, c_{0} (x) \propto \frac{f ( x )}{f _{0} ( x )}

c_{1} (x) f_{1} (x) - c_{0} (x) f_{0} (x) \equiv 0 ⟹ c_{1} (x) \propto \frac{f ( x )}{f _{1} ( x )}, c_{0} (x) \propto \frac{f ( x )}{f _{0} ( x )}

co v (z_{i}, z_{j})

co v (z_{i}, z_{j})

\to \frac{\frac{n _{0 s}}{n _{0 i} n _{0 j}} \int f _{0 s} ( x ) \frac{f _{i} ( x ) f _{j} ( x )}{f _{0 i} ( x ) f _{0 j} ( x )} d x + \frac{n _{1 s}}{n _{1 i} n _{1 j}} \int f _{1 s} ( x ) \frac{f _{i} ( x ) f _{j} ( x )}{f _{1 i} ( x ) f _{1 j} ( x ) d x}}{\frac{1}{n _{0 i}} \int \frac{f _{i} ( x ) ^{2}}{f _{0 i} ( x )} d x + \frac{1}{n _{1 i}} \int \frac{f _{i} ( x ) ^{2}}{f _{1 i} ( x )} d x \frac{1}{n _{0 j}} \int \frac{f _{j} ( x ) ^{2}}{f _{0 j} ( x )} d x + \frac{1}{n _{1 j}} \int \frac{f _{j} ( x ) ^{2}}{f _{1 j} ( x )} d x}

Σ_{d m} = v a r ((z_{d} z_{m})^{t} ∣ H_{0}^{\cup}) = (1 ρ_{d m} ρ_{d m} 1)

Σ_{d m} = v a r ((z_{d} z_{m})^{t} ∣ H_{0}^{\cup}) = (1 ρ_{d m} ρ_{d m} 1)

z_{γ}^{'} \to \infty/ z_{α}^{'} \to \infty lim \frac{∣ Σ _{A} ∣ ( ( ρ _{d s} ρ _{d m} - ρ _{s m} ) z _{γ}^{'} + ( ρ _{d m} ρ _{s m} - ρ _{d s} ) z _{α}^{'} + ∣ Σ _{d m} ∣ z _{β^{*}} )}{∣ Σ _{B} ∣ ( - ρ _{r m} z _{γ}^{'} + ρ _{d m} ρ _{r m} z _{α}^{'} + ∣ Σ _{d m} ∣ z _{β} )} = 1

z_{γ}^{'} \to \infty/ z_{α}^{'} \to \infty lim \frac{∣ Σ _{A} ∣ ( ( ρ _{d s} ρ _{d m} - ρ _{s m} ) z _{γ}^{'} + ( ρ _{d m} ρ _{s m} - ρ _{d s} ) z _{α}^{'} + ∣ Σ _{d m} ∣ z _{β^{*}} )}{∣ Σ _{B} ∣ ( - ρ _{r m} z _{γ}^{'} + ρ _{d m} ρ _{r m} z _{α}^{'} + ∣ Σ _{d m} ∣ z _{β} )} = 1

α \to 0 lim \frac{z _{β^{*}}}{1 - ρ _{d s}^{2} z _{β} + ρ _{d s} z _{α}} = 1

α \to 0 lim \frac{z _{β^{*}}}{1 - ρ _{d s}^{2} z _{β} + ρ _{d s} z _{α}} = 1

\int_{R} e^{- \frac{1}{2} x^{t} Σ^{- 1} x} (Φ (A_{1}^{t} x + C) - Φ (A_{0}^{t} x + C_{0})) d x_{1} d x_{2} ... d x_{N} = 0

Peer Reviews

No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.

Videos

No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.

Full text

Combining controls can improve power in two-stage association studies

James Liley

Department of Medicine, University of Cambridge, Addenbrooke’s Hospital, Cambridge, CB2 0SP, UK. E: [email protected]

Abstract

High dimensional case control studies are ubiquitous in the biological sciences, particularly genomics. To maximise power while constraining cost and to minimise type-1 error rates, researchers typically seek to replicate findings in a second experiment on independent cohorts before proceeding with further analyses.

This paper presents a method in which control (or case) samples from the discovery cohort are re-used in the replication study. The theoretical implications of this method are discussed and simulations used to compare performance against the standard method in a range of circumstances. In several common study designs, a shared-control method allows a substantial improvement in power while retaining type-1 error rate control.

An important area of potential application arises when control samples are difficult to recruit or ascertain; for example in inter-disease comparisons, or studies on degenerative diseases. Using similar methods, a procedure is proposed for ‘partial replication’ using a new independent cohort consisting of only controls. This methods can be used to provide some validation of findings when a full replication procedure is not possible.

The new method has differing sensitivity to confounding in study cohorts compared to the standard procedure, which must be considered in its application. Type-1 error rates in these scenarios are analytically and empirically derived, and an online tool for comparing power and error rates is provided.

Although careful consideration must be made of all necessary assumptions, this method can enable more efficient use of data in genome-wide association studies (GWAS) and other applications.

Introduction

High-dimensional case-control studies have become a mainstay of investigation of pathophysiology in complex diseases and traits. An important part of their analysis is the process of replication [1], in which the results of a high-dimensional study are used to inform the design of a second study at a subset of the original variables, with a joint analysis used to determine overall association.

Replicating studies in this way has the advantage of increasing the effective study sample sizes without requiring measurement of all variables in all samples. It also serves to protect against false-positives due to systematic errors in the original datasets, by re-testing association in a second nominally independent dataset.

Replication has a significant cost, and can require large numbers of samples, especially when associated variables have small effects (ie [2]). There is therefore a need to minimise the number of additional samples which need to be analysed. This paper presents a method to perform replication by combining controls in both the original ‘discovery’ and second ‘replication’ datasets, potentially reducing the number of new samples required. Shared-control approaches can improve study efficiency in many related applications in which studies are compared [3, 4, 5, 6, 7, 8].

Results from original and replication datasets for which some or all controls are shared cannot be directly compared due to the correlation between test statistics directly resulting from shared controls even under the null hypothesis [5]; use of the same thresholds in a shared-control design as used in an independent-controls design will lead to higher type-1 error rates. This paper demonstrates a simple adaptation to a standard design to account for the changed correlation structure and retain control of type-1 error rate, only requiring a change to one p-value threshold.

The action of sharing control samples results in a different spectrum of sensitivity to confounding in study groups. It necessitates a sacrifice of type-1 error rate control in variables affected by confounding in the discovery-phase control cohort, but improves type-1 error rate control in variables affected by confounding in the replication-phase control cohort. Performance is largely equivalent to an independent-controls design for variables affected by confounding in either case cohort.

The new spectrum of false positive rates can be advantageous in circumstances where control samples in the replication cohort are less well-ascertained than those in the discovery cohort. This may be the case in studies on degenerative disease, where control ascertainment is generally uncertain, and population-sourced controls may be used for replication. The shared-control design can reduce power losses from mis-specified controls in the replication cohort, as well as reducing false-positive rates caused by confounding in the cohort.

When used with shared cases instead of controls, this method can be adapted to a ‘partial replication’ procedure where only a new control set is used. Although not equivalent to a full replication in an independent dataset, the procedure enables improvement in type-1 error rates and control over confounding. This is applicable in studies on rare traits, where all available samples need to be included in the discovery analysis for adequate power.

Throughout this paper, GWAS terminology will be used (SNPs, allele frequency, variants etc) although the method is applicable to any high-dimensional case control study. ‘Controls’ will be considered to generally be samples unaffected by a disease or trait of interest, although the method can be applied with case/control labels swapped, or applied to comparisons between subgroups of a case group. Asymptotic analytical results are established where possible, but all type 1/type 2 error rates are readily tractable empirically to good accuracy given study sizes and proposed p-value thresholds, and a tool is provided to do this at https://wallacegroup-liley.shinyapps.io/replication_shared/.

Results

Overview of method

We assume a GWAS dataset of a set of cases $C_{1}$ and controls $C_{0}$ used in a ‘discovery’ phase of a GWAS or similar study, and corresponding sets of cases and controls $C_{1}^{\prime}$ , $C_{0}^{\prime}$ in the replication phase. We assume that $C_{0}$ and $C_{1}$ are genotyped at a set of SNPs $S$ and $C_{0}^{\prime}$ , $C_{1}^{\prime}$ at a set $S^{\prime}\subseteq S$ .

For each SNP we designate $\mu_{1}$ , $\mu_{0}$ , $\mu_{1}^{\prime}$ , $\mu_{0}^{\prime}$ as the population minor allele frequency in the corresponding group, and $m_{1}$ , $m_{0}$ , $m_{1}^{\prime}$ , $m_{0}^{\prime}$ as the observed allele frequency (so $E(m_{i})=\mu_{i}$ ). We designate two null hypotheses; $H_{0}^{\cup}:(\mu_{1}=\mu_{0})\cup(\mu_{1}^{\prime}=\mu_{0}^{\prime})$ and $H_{0}^{\cap}:(\mu_{1}=\mu_{0}=\mu_{1}^{\prime}=\mu_{0}^{\prime})$ , noting that $H_{0}^{\cup}\supseteq H_{0}^{\cap}$ . In a typical conservative GWAS approach, we seek to test against $H_{0}^{\cup}$ , since $\mu_{1}\neq\mu_{0}$ or $\mu_{1}^{\prime}\neq\mu_{0}^{\prime}$ may hold at non-disease associated SNPs due to confounding in the original or replication studies respectively.

A typical two-stage genetic testing procedure [1], which we will refer to as method A, begins by comparing genotypes of $C_{1}$ and $C_{0}$ at SNPs $S$ generating p-values $p_{d}$ (discovery). A subset $S^{\prime}$ of SNPs reaching putative significance level $p_{d}<\alpha$ are genotyped in $C_{0}^{\prime}$ and $C_{1}^{\prime}$ , with genotypes compared to generate p-values $p_{r}$ (replication stage). Finally, genotypes are compared between $C_{0}\cup C_{0}^{\prime}$ and $C_{1}\cup C_{1}^{\prime}$ at SNPs $S^{\prime}$ to generate p-values $p_{m}$ (meta-analytic stage). SNPs are designated as ‘hits’ if $p_{d}<\alpha,p_{r}<\beta,p_{m}<\gamma$ for some $\beta$ , $\gamma$ , and all effects have the same direction.

The main modification proposed in this paper, denoted as method B, differs at the replication stage in that $C_{1}^{\prime}$ is compared with $C_{0}\cup C_{0}^{\prime}$ at $S^{\prime}$ instead of just $C_{0}^{\prime}$ . The p-values resulting from the modified replication stage are termed $p_{s}$ , and the criterion to designate a hit changed to $p_{d}<\alpha,p_{s}<\beta^{*},p_{m}<\gamma$ , with all effects in the same direction. The threshold $\beta^{*}$ is chosen to conserve type-1 error rate between methods (see methods section, appendix 1.1).

A second modification, denoted method C, combines $C_{0}$ and $C_{0}^{\prime}$ at both the discovery and replication phase. This is analagous to a situation in which only a single control cohort is available, and a choice must be made to split it between discovery and replication procedures or to use it for both. In this case, $C_{0}\cup C_{0}^{\prime}$ is compared with $C_{1}$ at SNPs $S$ in the discovery phase to produce p-values $p_{c}$ , then $C_{0}\cup C_{0}^{\prime}$ is compared with $C_{1}^{\prime}$ at SNPs $S^{\prime}$ at the replication phase and compared with $C_{1}\cup C_{1}^{\prime}$ at the meta-analytic stage to produce p-values $p_{s}$ and $p_{m}$ as before. A hit is determined by $p_{d}<\alpha,p_{s}<\beta^{\perp},p_{m}<\gamma$ , with all effects in the same direction. Again, $\beta^{\perp}$ is chosen to maintain the type-1 error rate between methods.

General properties

For SNPs in $H_{0}^{\cap}$ , the overall type-1 error rate is conserved between methods by the definition of $\beta^{*}$ , $\beta^{\perp}$ (equation 2) at a level $P_{0}$ . It is shown in appendix 1.2.2 that $\beta>\beta^{*}>\beta^{\perp}$ . For SNPs in $H_{0}^{\cup}\setminus H_{0}^{\cap}$ the type-1 error rates differ between methods. Such SNPs may be characterised by the group(s) amongst $C_{0}$ , $C_{1}$ , $C_{0}^{\prime}$ , $C_{1}^{\prime}$ in which their expected MAF is aberrant from the expected MAF in the population which the group ostensibly represents. ‘Aberrance’ is taken to mean an incorrect expected value from systematic measurement error or uncorrected confounding, rather than random deviance around a correct expected value.

Bounds on type-1 error rates with aberrance in each group are shown in table 1. Methods B and C necessitate sacrificing bounds on error rates with aberrance in $C_{0}$ and $C_{0}$ , $C_{0}^{\prime}$ respectively. The bound on error with aberrance in $C_{1}^{\prime}$ improves through methods A-C. In the methods section, it is shown that the type-1 error with aberrance in $C_{0}^{\prime}$ decreases from methods A to B, and the error with aberrance in $C_{1}^{\prime}$ increases from A through C, although the upper bound is the same for both.

Simulation

The power difference between methods B and A was analysed systematically across a range of values of $(n_{0},n_{1},n_{0}^{\prime},n_{1}^{\prime})$ . Average power difference and maximum power difference were compared (see methods section). Figure 3 shows power difference at various study sizes for typical $\alpha$ , $\beta$ , $\gamma$ values ( $\alpha=5\times 10^{-6}$ , $\beta=5\times 10^{-4}$ , $\gamma=5\times 10^{-8}$ ) and minor allele frequency 0.1. The difference is typically highest when the ratio of controls to cases is high in the discovery cohort and low or equal in the replication cohort, and the number of cases in the discovery cohort is larger than the number in the replication cohort. Power to detect SNPs in $H_{1}$ is typically highest in method C, second-highest in method B, and lowest in method A.

Recommended applications

To demonstrate areas where this approach is applicable, several examples are constructed or sourced from the GWAS field in which the procedure of sharing controls or cases will improve power or type-1 error profile of the two-stage testing procedure or enable some form of orthogonal replication to be performed.

Assumptions

In order to use method B or C, it must be assumed that cohort $C_{0}$ and $C_{0}^{\prime}$ are sampled from similar enough populations to be comparable to $C_{1}$ and $C_{1}^{\prime}$ (possibly with the inclusion of strata or covariates in the genetic risk model). An important caveat of methods B and C is sacrifice of control over errors arising from aberrance in $C_{0}$ (method B) or $C_{0}\cup C_{0}^{\prime}$ (method C), so an assumption must be made that variables affected by confounding or measurement error in these cohorts are understood to be distinguishable from true associations by quality-control measures only.

Post-hoc assessment of all putative hits should be performed to check for genotyping errors [9] and assess whether the hit could have arisen from aberrance in $C_{0}$ .

Conventional GWAS

Method B is applicable in several cases in large conventional GWAS, particularly when then ratio of controls to cases in the discovery cohort is larger than that in the replication cohort. In a relatively recent GWAS on rheumatoid arthritis [10] with comparable sample populations for discovery and replication cohorts, method B could be used to attain greater power than method A for a fixed type-1 error rate. Assuming that summary statistics are well-approximated by binomial tests of allelic differences (so covariates and strata used in computation of summary statistics have only small effects), the improvement in power is around 4% for SNPs with an odds-ratio of $1.3$ , MAF 0.1, and is positive across all odds ratios. More than 2000 additional controls in $C_{0}^{\prime}$ would be needed to increase power by this amount (figure 4(a)).

Small power advantages such as this may make minimal difference in a single study, although since they require no extra cost, are worth attaining if possible. The power of method B is generally considerably higher than method A when $n_{0}>n_{1}$ and $n_{0}^{\prime}\approx n_{1}^{\prime}$ , Power advantages may be more substantial in some cases; for example, a study with $(n_{0},n_{1},n_{0}^{\prime},n_{1}^{\prime})=(15000,5000,5000,5000)$ , method B enables a power increase of up to 8% (Figure 4(b)). To achieve comparable performance with method A, around 2000 additional controls would be necessary in the replication cohort. Method B with $(n_{0},n_{0}^{\prime})=(15000,5000)$ is also more powerful than method A would be if controls were divided equally between $C_{0}$ and $C_{0}^{\prime}$ (see Figure 4(b)).

Difficult control ascertainment

An important application of the method presented in this paper is in studies for which ‘control’ samples are expensive or difficult to ascertain. This is often the case in comparative studies between disease subtypes. In such studies, sharing controls can improve power substantially, especially if a proportion of samples in the replication cohort are falsely assigned to the control cohort (see methods section).

An international GWAS on fronto-temporal dementia in 2014 [11] is an example in which sharing controls may be beneficial. The study had sample sizes $(n_{0},n_{1},n_{0}^{\prime},n_{1}^{\prime})=(4308,2154,5094,1372)$ . Control samples in the discovery phase were assessed for current neurological disease, and were used in previous studies on Parkinson’s disease, indicating a high degree of reliability. Control samples in the replication phase were collected from the same geographic distribution as cases, but were not explicitly used in previous neurological studies, suggesting better control ascertainment amongst the discovery cohort.

In this study, sharing controls could allow for a more strongly-ascertained control cohort, and reduce the effects of confounders affecting $C_{1}^{\prime}$ . At typical values $\alpha=1\times 10^{-4}$ , $\beta=1\times 10^{-3}$ , $\gamma=5\times 10^{-8}$ , power is nearly equivalent between the two methods assuming all controls are genuine. However, with 10% misascertainment in $C_{1}^{\prime}$ , the power advantage of method B is up to 5%. Given the near-identical distribution of cases in the discovery and validation cohort, cases could alternatively be shared, leading to a power increase of up to 6%.

Prospective study design

Studies may be planned and powered with the assumption that samples may be shared. For certain restrictions on sample numbers, this can provide the potential for greater power than would be attainable by restricting to an independent-controls design. For instance, if we seek to validate hits on a GWAS with 10000 controls and 5000 cases, and can afford to genotype a further 10000 samples, power is higher after recruiting 4000 additional controls and 6000 additional cases and sharing controls than can be achieved from any independent-control study design (Figure 4(d)).

This may be a common scenario if controls are sourced from a known database rather than specifically recruited for the study.

Partial replication

In circumstances where case recruitment is difficult, as in studies of rare diseases, an assessment of replicability may be made by re-testing results from a discovery phase with a new control set only. This can enable the use of control cohorts which only partially match the case cohort.

In a GWAS on pemphigus vulgaris [12], a rare disease primarily affecting individuals of Ashkenazi Jewish ethnicity, the discovery cohorts were sampled from Jewish populations, with age- and population- matched controls. Control cohorts were small ( $(n_{0},n_{1},n_{0}^{\prime},n_{1}^{\prime})=(100,400,59,285)$ ), potentially due to difficulty recruiting both ethnically- and geographically-matched controls.

Method C could be used in this instance to enable a larger control set and greater power. If a control cohort of Ashkenazi individuals could be assembled without requiring geographic matching with the case set, it would be inappropriate to use as a sole control cohort against the existing case cohort, due to the potential for geographic confounding. However, such a cohort could be used as either $C_{0}$ or $C_{0}^{\prime}$ in method C, with the existing ethnically- and geographically- matched controls serving as the other cohort. In this way, the power advantage of the larger cohort could be used while maintaining control over potential aberrance in the larger control group.

Method C enables computation of power and type-1 error rates, and comparison to alternative designs with cases split into smaller independent discovery and validation cohorts (method A). Testing a case cohort against two separate control cohorts is almost always more powerful for a fixed type-1 error rate than splitting the case cohort in two and performing method A (see supplementary figures 7,9).

Discussion

This paper proposes a method to improve efficiency of data use in a replication procedure, adding to the body of methods for comparison of high-dimensional case-control studies. For many common study sizes, the method can reduce the cost of replication, or increase power of discovery. The adapted method is simple to apply, only requiring modification of a single association threshold. A standard replication procedure (or more general comparison of case-control studies) with independent control datasets does not make use of the information that expected values of variables in control datasets are, in principle, the same. In this way, the same dataset can in theory yield more information when controls are shared.

The most important caveat of these methods is the loss of systematic type-1 error rate control for null SNPs which are aberrant in $C_{0}$ . Control of such errors must not be sacrificed entirely, but in some circumstances it may be satisfactory to assess such errors on a SNP-by-SNP basis. Such assessment is important and standard for all proposed GWAS hits under any method [13] in the interests of quality control. In method C, control over aberrance in $C_{0}^{\prime}$ is additionally lost; however, since this method is largely applicable when $C_{0}\cup C_{0}^{\prime}$ is a single homogenous control (or case) cohort, there is no way that aberrance in the cohort can be systematically identified by comparison with other cohorts.

Somewhat better control of the type-1 error rate can often be achieved for SNPs with aberrance in $C_{1}$ or $C_{0}^{\prime}$ . This may incentivise the use of this method when confidence in the representativeness of these cohorts is low compared to that of $C_{0}$ . The type 1 error rate is somewhat increased for SNPs with aberrance in $C_{1}^{\prime}$ , although as it remains bounded by $\alpha$ , this increase is not a major problem.

The two-stage validation procedure is similar to a meta-analysis of the discovery and validation experiments, for which several adaptations to shared-control designs have been proposed [3, 4]. However, there are several important distinctions which necessitate an alternative approach in this case. Firstly, not all variables are measured in the second (replication) study; we are restricted to analysis of variables reaching a given observed effect size. Secondly, the studies to be ‘meta-analysed’ are not complete, in the sense that there may be residual confounding; a strong effect size in the meta-analysis alone is not adequate evidence for association and some level of association (with consistent direction) is additionally required in both constituent studies.

The method is inapplicable when replication is performed on cohorts from completely distinct geographic groups, although there can be some difference in geographic distribution between control sets if this is controlled for in computing summary statistics. The method is most applicable when control groups are sampled from similar populations and genotyped on similar platforms.

The widespread discoveries of the GWAS field have led to corresponding increases in complexity of phenotypic definitions, with ever-finer delineations of disease types of ever-rarer prevalence. The genetic analysis of such complex phenotypes is necessarily comparative; there is little use understanding the genetics of a rare disease subtype except in the context of the genetics of the disease in general. Such analyses necessitate GWAS and other comparative studies between rare phenotypic types [14], with ‘controls’ meaning the better-characterised disease subphenotype in this sense, as well as between cases and controls. Rare disease subtypes are often afflicted with ascertainment difficulties, leading to varying degrees of expected aberrance in disease cohorts. Within this paradigm, the applicability of this method is likely to expand.

Acknowledgments

I would like to thank Dr Chris Wallace and Dr Jenn Asimit for helpful comments and review of this work. JL is funded by the NIHR Cambridge Biomedical Research Centre and is on the Wellcome Trust PhD programme in Mathematical Genomics and Medicine at the University of Cambridge. He is also supported by the Wallace group, funded by the Wellcome Trust (089989,107881, CW) and the MRC (MC_UP_1302/5, CW). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Conflicts of interest

None declared

Methods

Definitions

Denote $z_{x}$ for $x\in\{d,r,s,m,c\}$ as the signed z-score ( $\pm\Phi^{-1}(p_{x}/2)$ ) corresponding to $p_{x}$ , and $z_{x}$ for $x\in{\alpha,\beta,\beta^{*},\gamma}$ as the positive corresponding threshold $-\Phi^{-1}(x/2)$ , where $\Phi,\Phi^{-1}$ are the standard normal CDF and quantile functions. Other than $(z_{d},z_{r})$ , all pairs of z-scores are correlated under $H_{0}^{\cap}$ , with correlation estimable from sample sizes or empirically if covariates are used (Appendix 1.1). Denote $\rho_{xy}$ as the correlation between $z_{x}$ and $z_{y}$ , $(x,y)\in\{d,r,s,m\}^{2}$ , and set

[TABLE]

For $i\in\{d,r,s,m,c\}$ define $\zeta_{i}=E(z_{i})$ , where the expectation is conditional on the SNP in question. For SNPs in $H_{0}^{\cap}$ , $\zeta_{i}\equiv 0$ , but this may not hold for SNPs in $H_{0}^{\cup}\setminus H_{0}^{\cap}$ . In theoretical working, aberrance in groups is characterised by values $\zeta_{i}$ rather than log-odds ratios, noting that the values $\zeta_{i}$ are asymptotically proportional to the corresponding log-odds ratios. Define $R_{A}$ , $R_{B}$ , $R_{C}$ as the false-positive rates for a SNP of interest in methods A, B and C respectively.

General type 1 error rate

The values $\beta^{*}$ , $\beta^{\perp}$ are chosen to satisfy

[TABLE]

thus conserving the type 1 error rate (denoted $P_{0}$ ) against $H_{0}^{\cap}$ between methods (Figure 5). If no threshold is used on $p_{m}$ (ie, $\gamma=1$ ), then $\beta^{*}$ , $\beta^{\perp}$ satisfy

[TABLE]

since $z_{d}\mathrel{\text{\scalebox{1.07}{$ \perp\mkern-10.0mu\perp $}}}z_{r}|H_{0}^{\cap}$ . Definition 2 will be considered a generalisation of definition 3, with results established first for $\beta^{*}$ as per definition 3 and extending where possible to definition 2.

For $\beta^{*}$ defined as per definition 3 we have (see Appendix 1.2.1)

[TABLE]

approaching from above, so $z_{\beta^{*}}>max\left(z_{\beta},\sqrt{1-\rho_{ds}^{2}}z_{\beta}+\rho_{ds}z_{\alpha}\right)$ , $z_{\beta^{\perp}}>max\left(z_{\beta},\sqrt{1-\rho_{cs}^{2}}z_{\beta}+\rho_{cs}z_{\alpha}\right)$ . As defined by equation 3, $z_{\beta^{*}}$ , $z_{\beta^{\perp}}$ are also asymptotically linear in $z_{\alpha}$ , $z_{\gamma}$ , $z_{\beta}$ as the former two tend to $\infty$ , with some constraints (Appendix 1.2.1), although the limit does not necessarily approach from above. For both definitions, $\beta^{\perp}<\beta^{*}<\beta$ (Appendix 1.2.2)

Study sizes, odds ratios and allele frequencies

Consider a study with $n_{0}$ controls and $n_{1}$ cases, with underlying allele frequencies $\mu_{0}$ and $\mu_{1}$ in cases and observed allele frequencies $m_{0}$ , $m_{1}$ . Let $Z$ be a signed Z-score derived from a GWAS p-value against the null hypothesis $\mu_{0}=\mu_{1}$ . To first order,

[TABLE]

where $\bar{\mu}=\frac{n_{0}\mu_{0}+n_{1}\mu_{1}}{n_{0}+n_{1}}$ . Hence

[TABLE]

where $\bar{\mu}$ varies between definitions (though it is sometimes taken to be approximately equal). These formulae allow $\zeta_{i}$ to be estimated in empirical computations. Estimation of $\zeta_{i}$ is more complex if covariates or strata are used in the computation of $z_{i}$ (appendix 1.1).

In all empirical computations, systematic allelic differences (whether due to aberrance or true association) are characterised by odds ratios and minor allele frequency rather than $\zeta_{i}$ . Values $\mu_{0}$ and $\mu_{1}$ (in general terms) are readily computable from odds-ratios $R=\frac{\mu_{1}(1-\mu_{0})}{\mu_{0}(1-\mu_{1})}$ and minor allele frequency $\mu=\frac{\mu_{0}+\mu_{1}}{2}$ . The use of this characterisation of minor allele frequency (rather than the ‘weighted’ characterisation $\frac{n_{0}\mu_{0}+n_{1}\mu_{1}}{n_{0}+n_{1}}$ ) is so that the correspondence between $(R,\mu)$ and $(\mu_{0},\mu_{1})$ is independent of $n_{0},n_{1}$ ; for instance, for a disease-associated variant satisfying $\mu_{1}=\mu_{1}^{\prime}$ , $\mu_{0}=\mu_{0}^{\prime}$ will have identical odds ratios $\frac{\mu_{1}(1-\mu_{0})}{\mu_{0}(1-\mu_{1})}$ and $\frac{\mu_{1}^{\prime}(1-\mu_{0}^{\prime})}{\mu_{0}^{\prime}(1-\mu_{1}^{\prime})}$ between the discovery and validation experiments, but generally different weighted MAFs $\frac{n_{0}\mu_{0}+n_{1}\mu_{1}}{n_{0}+n_{1}}$ and $\frac{n_{0}^{\prime}\mu_{0}^{\prime}+n_{1}^{\prime}\mu_{1}^{\prime}}{n_{0}^{\prime}+n_{1}^{\prime}}$ .

False ascertainment

In general, for a true association, $\mu_{0}=\mu_{0}^{\prime}$ and $\mu_{1}=\mu_{1}^{\prime}$ . If some proportion $\kappa$ of samples in $C_{0}^{\prime}$ are incorrectly assigned and come from the case population, then $\mu_{0}^{\prime}=(1-\kappa)\mu_{0}+\kappa\mu_{1}$ . This lowers the absolute values of $\zeta_{r}$ , $\zeta_{s}$ and $\zeta_{m}$ , reducing the power to detect the SNP.

Empirical computations

Define $N_{\Sigma}(\mathbf{z})$ as the $pdf$ of the multivariate normal with mean [math] and variance $\Sigma$ at $\mathbf{z}$ . Determination of covariance is described in Appendix 1.1. Given $\zeta_{d}$ , $\zeta_{r}$ , $\zeta_{s}$ , $\zeta_{m}$ , the probability of rejecting the null for a given SNP using method A is

[TABLE]

and using method B

[TABLE]

If $\frac{n_{0}}{n_{1}}=\frac{n_{0}^{\prime}}{n_{1}^{\prime}}$ , matrix $\Sigma_{A}$ is singular (Appendix 1.1), in which case $z_{m}=\rho_{dm}z_{d}+\rho_{vm}z_{v}$ and the expression above may be reduced to a two-dimensional integral over a more complex region (Figure 5). Matrix $\Sigma_{C}$ is generally singular, so the formula $z_{m}=\frac{\rho_{cs}\rho_{sm}-\rho_{cm}}{\rho_{cs}^{2}-1}z_{d}+\frac{\rho_{cs}\rho_{cm}-\rho_{sm}}{\rho_{cs}^{2}-1}z_{s}$ is used to reduce the integral in a similar way. A similar formula may be used if $\Sigma_{B}$ is nearly singular.

In Figure 3, mean power difference is determined as the integral of the power difference with respect to the log-odds ratio over the real line.

Type 1 error rates

Aberrance in $C_{1}$

For SNPs aberrant in only $C_{1}$ we have $\zeta_{d}\neq 0$ , $\zeta_{c}\neq 0$ , $\zeta_{m}\neq 0$ , and $\zeta_{r}=\zeta_{s}=0$ .

$R_{A}$ , $R_{B}$ , $R_{C}$ can be considered as functions of $\zeta_{d}$ . As $\zeta_{d}\to 0$ , $R_{A},R_{B},R_{C}\to P_{0}$ (equation 2). As $\zeta_{d}\to\pm\infty$ , $R_{A}\to\frac{\beta}{2}$ , $R_{B}=\frac{\beta^{*}}{2}$ and $R_{C}=\frac{\beta^{\perp}}{2}$ . For positive $\zeta_{d}$ both $R_{A}$ and $R_{B}$ are increasing (and both are symmetric in $\zeta_{d}$ ) so $R_{A}<\frac{\beta}{2}$ , $R_{B}<\frac{\beta^{*}}{2}$ , $R_{C}<\frac{\beta^{\perp}}{2}$ for all $\zeta_{d}$ .

Since $\beta^{\perp}<\beta^{*}<\beta$ (often substantially), methods B and C are generally better at rejecting $H_{0}^{\cap}$ for such SNPs. In the simplified case where $z_{\gamma}=1$ , $R_{A}\geq R_{B}$ universally (Appendix 1.3.1. This typically holds for all $z_{\gamma}$ , except for small deviations in pathological cases.

In general, we consider aberrance which is only still present after any strata or covariates have been accounted for in the computation of $z$ scores. If strata or covariates remove the effective aberrance between groups, the type-1 error rate is equivalent to that under $H_{0}^{\cap}$ .

Aberrance in $C_{1}^{\prime}$

For SNPs aberrant in $C_{1}^{\prime}$ , we have $\zeta_{d}=0$ , $\zeta_{c}=0$ , $\zeta_{r}\neq 0$ , $\zeta_{s}\neq 0$ and $\zeta_{m}\neq 0$ .

Again, $R_{A},R_{B},R_{C}\to P_{0}$ as $\zeta_{r}\to 0$ . As $\zeta_{r}\to\pm\infty$ , $R_{A},R_{B},R_{C}\to\frac{\alpha}{2}$ , and both are bounded by $\frac{\alpha}{2}$ . Although $R_{B}$ and $R_{C}$ are typically higher than $R_{A}$ in this case, since both have the same (typically conservative) upper bound, this is not typically a large sacrifice in type 1 error.

In the simplified case where $\gamma=1$ , an approximate upper bound on $R_{B}-R_{A}$ is given by (Appendix 1.4)

[TABLE]

where

[TABLE]

In practice, there is typically a similarly small difference between $R_{C}$ , $R_{B}$ and $R_{A}$ in the general case.

Aberrance in $C_{0}^{\prime}$

For SNPs aberrant in $C_{0}^{\prime}$ , $\zeta_{d}=0$ , $\zeta_{r}\neq 0$ , $\zeta_{c}\neq 0$ , $\zeta_{s}\neq 0$ and $\zeta_{m}\neq 0$ . As for SNPs with aberrance in $C_{1}^{\prime}$ , $R_{A},R_{B},R_{C}\to P_{0}$ as $\zeta_{r}\to 0$ and as $\zeta_{r}\to\pm\infty$ , $R_{A},R_{B}\to\frac{\alpha}{2}$ , both bounded above by $\frac{\alpha}{2}$ . $R_{C}$ , however, tends to 1 as $\zeta_{d}\to\infty$ .

In method B the cohort $C_{0}$ has a correcting effect on the replication study, meaning $|\zeta_{s}|<|\zeta_{r}|$ and $R_{B}<R_{A}$ .

For the simplified case where $\gamma=1$ , a similar bound to 9 holds for the difference $R_{A}-R_{B}$ (note signs are reversed) with

[TABLE]

in the place of $k$ . The improvement in type-1 error rate for a SNP with aberrance in $C_{0}^{\prime}$ is generally larger than the loss with the same aberrance in $C_{1}^{\prime}$ (see methods), meaning that if aberrances are of similar prevalence and size in $C_{1}^{\prime}$ and $C_{0}^{\prime}$ , method B will typically have a lower type-1 error rate than method A.

Aberrance in $C_{0}$

Aberrance in $C_{0}$ represents a serious problem in case-control study comparison. False-positive rates are generally worse under method B, and tend to 1 as $E(z)\to\infty$ . If aberrances of this type are expected to be very frequent, this may preclude use of methods B or C.

However, aberrances of this type may be best detected retrospectively by examining aberrances between control groups at SNPs declared ‘hits’. This procedure is already a necessary quality-control procedure in method A [13, 9], as method A does not provide any control over differences between $C_{0}$ and $C_{0}^{\prime}$ . The number of SNPs reaching significance in the two-stage procedure is usually small enough that this examination is readily tractable.

Aberrance in two or more cohorts

If SNPs are aberrant in both $C_{1}$ and $C_{1}^{\prime}$ , or in both $C_{0}$ and $C_{0}^{\prime}$ , the effect on $R_{A}$ and $R_{B}$ is similar. If both cohorts are aberrant in the same direction, there is no way to differentiate the SNP from a genuine association on the basis of the genotype data alone. If cohorts are aberrant in different directions, then in both methods, the type-1 error rate is lower than for a null SNP with no aberration or aberration in only one cohort, as effect sizes for the discovery and replication cohorts are biased in opposite directions. The same typically holds if $C_{0}^{\prime}$ and $C_{1}$ , or $C_{0}$ and $C_{1}^{\prime}$ , are biased in the same direction.

If $C_{0}^{\prime}$ and $C_{1}^{\prime}$ or $C_{0}$ and $C_{1}$ are both biased in the same direction, $R_{A}$ is generally lower than $R_{B}$ , as $\zeta_{s}\neq 0$ . Both $R_{A}$ and $R_{B}$ are bounded by $\frac{\alpha}{2}$ in this case. In addition, a systematic bias in both replication groups (or both discovery groups) is likely to be due to a known confounder, the effect of which can be removed by performing a stratified test (as is typically good practice when confounders are known). Aberrance in opposite directions leads to $R_{B}>R_{A}$ in the first case, and a scenario similar to aberrance in $C_{0}$ in the second case.

Aberrance in three or more cohorts corresponds to a chaotic scenario in which neither methods A,B, or C will reliably provide FPR control. Aberrance of this extent is typically detectable and removable using quality control procedures.

Supplementary figures

1 Appendices

1.1 Covariance between Z scores due to shared samples

The matching of type-1 error rates between methods relies on establishing the covariance between triples of z-scores under $H_{0}^{\cap}$ . The covariance can be readily approximated when z-scores are assumed to be derived from tests of equality of binomial random variables $m_{i}$ . Z-scores comparing proportions without using additional information (such as strata or covariates) which are monotonic to allelic difference and asymptotically have $N(0,1)$ distribution must be asymptotically equivalent to those derived from binomial comparisons, so this assumption is reasonable.

If strata or covariates are used, either an assumption must be made that computed z-scores are well-approximated by comparisons of binomial proportions, or correlations must be approximated allowing for the covariate or strata structure. This is sometimes tractable analytically, but can also be estimated empirically either by using known non-associated variants or by simulating variants with the same covariate structure.

The presence of strata or covariates also affects the values $\zeta_{i}$ , and if the effects of covariates are large, the approximations in equations 6 in the main paper may be poor. Values $\zeta_{i}$ can be estimated as functions of allelic differences by simulating variants with the same covariate structure.

1.1.1 No covariates or stratification

Assume study $i$ and $j$ have $n_{0}^{i}$ , $n_{0}^{j}$ controls and $n_{1}^{i}$ , $n_{1}^{j}$ cases respectively, of which $n_{0}^{ij}$ controls and $n_{1}^{ij}$ cases are shared between both studies. Let $m_{0},m_{1},m_{0}^{\prime},m_{1}^{\prime}$ , denote the observed allele frequencies of a SNP in the respective cohort, and $\mu_{0},\mu_{1},\mu_{0}^{\prime},\mu_{1}^{\prime}$ the expected allele frequency.

If no strata or covariates are used in the calculation of summary statistics, $z$ scores $z_{d},z_{r},z_{s},z_{m}$ are asymptotically proportional to the allelic differences $m_{1}-m_{0}$ , $m_{1}^{\prime}-m_{0}^{\prime}$ , $m_{1}^{\prime}-\frac{m_{0}^{\prime}n_{0}^{\prime}+m_{0}n_{0}}{n_{0}+n_{0}^{\prime}}$ , $\frac{m_{1}n_{1}+m_{1}^{\prime}n_{1}^{\prime}}{n_{1}+n_{1}^{\prime}}-\frac{m_{0}^{\prime}n_{0}^{\prime}+m_{0}n_{0}}{n_{0}+n_{0}^{\prime}}$ respectively, since $z$ scores are monotonic with allelic differences and allelic differences are asymptotically normal. Since $m_{0},m_{1},m_{0}^{\prime},m_{1}^{\prime}$ are independent and asymptotically normal the multivariate random variables $(z_{d},z_{r},z_{m})$ and $(z_{d},z_{s},z_{m})$ have multivariate normal distributions.

For studies $i$ on $n_{0i}$ , $n_{1i}$ controls and cases and $j$ on $n_{0j}$ , $n_{1j}$ controls and cases in which $n_{0ij}$ and $n_{1ij}$ controls and cases are shared between studies, the correlation between the observed allelic differences $m_{1i}-m_{0i}$ , $m_{1j}-m_{0j}$ is given by

[TABLE]

This holds under $H_{0}^{\cap}$ and approximately holds in general. Expressions for $\rho_{ds}$ , $\rho_{dm}$ , $\rho_{rm}$ and $\rho_{sm}$ may be derived in terms of $n_{0}$ , $n_{1}$ , $n_{0}^{\prime}$ , and $n_{1}^{\prime}$ . Specifically

[TABLE]

so $\Sigma_{A}$ is singular if $\frac{n_{0}}{n_{1}}=\frac{n_{0}^{\prime}}{n_{1}^{\prime}}$ , and $\Sigma_{B}$ if $n_{0}^{\prime}n_{1}=0$ .

1.1.2 Z scores with stratification

If computation of Z scores is performed with correction for strata or covariates, formula 12 will not asymptotically hold and may be a poor approximation to the true covariance. The true covariance can be computed in some cases.

If samples are divided into strata $1,2,...s$ , and $n_{pq}^{r}$ , $m_{pq}^{r}$ , $\mu_{pq}^{r}$ denote the number of samples and observed and expected minor allele frequencies in cohort $p$ , study $q$ , stratum $r$ respectively, then the $z$ score $z_{i}$ for study $q=i$ is asymptotically given by

[TABLE]

for positive values $\alpha_{ir}$ depending on the values $n_{pi}^{r}$ . If the Cochran-Mantel-Hanszel test is used, then

[TABLE]

Suppose that $n_{0ij}^{r}$ controls and $n_{1ij}^{r}$ cases are shared between studies $i$ and $j$ in stratum $r$ . Since the values $m_{pq}^{r}$ are dependent only within the same values of $p$ and $r$ , the correlation between $z_{i}$ and $z_{j}$ under the null hypothesis $\mu_{0i}^{r}\equiv\mu_{1i}^{r},\mu_{0i}^{r}\equiv\mu_{1i}^{r}$ is given by

[TABLE]

where all sums are over the values of $r\in 1..s$ for which the relevant values of $n_{pq}^{r}$ are positive.

1.1.3 Z scores with covariates

If $z$ scores are computed adjusting for one or more covariates, the estimation of correlation is more difficult.

Assume that in a case population $C_{1}$ and a control population $C_{0}$ the values of some covariate(s) $x$ have different known distributions $f_{1}$ , $f_{0}$ , and that genotypes $g$ at some SNP of interest may vary with $x$ . We will assume the populations are large and that $f_{1}$ , $f_{0}$ , and $E(g|x)$ are continuous functions of $x$ .

Let $g_{p}^{k}$ denote the genotype of individual $k$ in cohort $p$ ( $p\in C_{0},C_{1}$ ) and $x_{p}^{k}$ denote covariate value(s), where $g_{p}^{k}$ is an observation of a random variable $g$ . An idealised $z$ -score testing association of $g$ with case/control status should be monotonic with each $g_{p}^{k}$ and have expectation 0 if $x$ is independent of case/control status, whatever the form of the function $E(g|x)$ . Because individual genotypes are assumed to be independent between individuals, cross-terms of the form $\prod_{i}g_{i}$ should carry no additional information from singleton genotypes. We thus assume that $z$ can thus be decomposed into a weighted linear sum of individual genotypes:

[TABLE]

where the (positive) values $c_{1}^{k},c_{0}^{k}$ depend only on the values $x_{1}^{\cdot},x_{0}^{\cdot}$ ; that is, not on the relationship between $g$ and $x$ , and the constant of proportionality depends on only on the observed allele frequency. Let function $c_{0}(x),c_{1}(x)$ denote the values of $c_{i}$ corresponding to covariate value(s) $x$ in $C_{0}$ , $C_{1}$ .

For a null SNP, $E(g|x)$ is independent of case/control status, but may take any (continuous) form. We have

[TABLE]

From a standard result from the calculus of variations, this implies that

[TABLE]

for some function $f$ , so the values $c_{1}^{k},c_{0}^{k}$ effectively reweight the contribution of individuals to a common density $f(x)$ across $x$ . The procedure of weighting observation $k$ in this way is analogous to a limiting case of stratification, in which weights are defined by the frequency of stratum $r$ (see above). For a constant allelic difference across the range of $x$ , the best common distribution to ‘map to’ does not depend on the relationship between $g$ and $x$ , and hence the best values of $c_{i}$ should be constant for all functions $E(g|x)$ .

Let $z_{q}$ denote a z-score for study $q$ ; $n_{pq}$ , $f_{pq}=f_{pq}(x)$ and $C_{pq}$ denote the number of samples, density function of $x$ , and set of samples in cohort $p$ , study $q$ ; $g_{pq}^{i}$ and $c_{pq}^{k}$ denote the normalised genotype of sample $k$ in cohort $p$ , study $q$ and its coefficient in $z_{q}$ ; $n_{0s}$ , $n_{1s}$ , $f_{0s}$ , $f_{1s}$ and $C_{0}^{s}$ , $C_{1}^{s}$ the number of shared controls/cases between studies, the density of $x$ amongst the shared samples, and the sets of shared samples; and $f_{q}$ the common density function to which cases and controls are weighted in study $q$ (equation 20). Then

[TABLE]

with integrals over the domain of $x$ , and the limit as sample sizes tend to infinity while ratios between them remain bounded.

Logistic regression models with continuous covariates can only model simple (generally linear) relationships between $c_{i}$ and $x_{i}$ , and property 20 may not hold. If the values $c_{pq}^{k}$ are known, the correlation can be determined using equation 21. If not, some methods for estimating correlation are outlined below.

1.1.4 Practical estimation of covariance

Although the asymptotic correlation between $z$ scores may be intractable, as long as the $z$ score permits an expansion of the form 18, the correlation will be nearly invariant with allele frequency and change only minimally for SNPs associated with the covariate.

In practical terms, one method to estimate the correlation between $z$ scores is to simply use the sample correlation at a set of variants presumed to be not associated with the main trait of interest. This approach may be unreliable and have limited power due to the difficulty of identifying such variants

Another option is to permute existing genotypes without permuting covariates, and compute correlation between resultant $z$ scores. This has the disadvantage that it is difficult to permute whilst maintaining potential relationships between genotypes and confounders.

Since the correlation should only depend on the sample sizes and structure of covariate distributions, a more convenient and powerful method is to simply simulate sets of genotypes unassociated with the trait, but potentially associated with covariates in a range of different ways, and compute correlation between the resultant $z$ scores. Given the shortcomings of standard methods such as logistic regression in fully accounting for covariate effects, this is an advisable procedure in any analysis including covariates.

All results in the main paper which require conditions on sample sizes are only approximate when using studies with stratification or covariates, with the approximation worsening with greater differences in covariate values between groups and lower effective sample sizes.

1.2 Properties of $\beta^{*}$

1.2.1 Asymptotic properties of $\beta^{*}$

In this appendix, an asymptotic approximation is established for $\beta^{*}$ and it is shown that $\beta^{*}>\beta$ for all $n_{0}^{i}$ , $n_{0}^{j}$ , $n_{1}^{i}$ , $n_{1}^{j}$ , $z_{\alpha}$ , $z_{\gamma}$ . Define $\Sigma_{A}$ and $\Sigma_{B}$ as per equations 1 in the main paper, and note that $\Sigma_{A}$ and $\Sigma_{B}$ only differ in their middle row/column. Further define

[TABLE]

Let $(z_{\alpha}^{\prime}\,z_{\gamma}^{\prime})$ be the point in $\{z_{d}>z_{\alpha},z_{m}>z_{\gamma}\}$ at minimal Mahalanobis distance from the origin with respect to $\Sigma_{dm}$ (ie, minimal $(z_{d}\,z_{m})\Sigma_{dm}^{-1}(z_{d}\,z_{m})^{t}$ ). Then for $z_{\gamma}^{\prime}-\rho_{dm}z_{\alpha}^{\prime}$ held constant, we have

[TABLE]

Specifically, for $\beta^{*}$ defined as per equation 3, we have

[TABLE]

and $z_{\beta^{*}}>max(\beta,\sqrt{1-\rho_{ds}^{2}}z_{\beta}+\rho z_{\alpha})$ for all $z_{\alpha}$ . Firstly the following lemma and corollary are established:

Lemma 1.

Let $\Sigma$ be a positive definite matrix of dimension $N$ , $\mathbf{x}$ be the vector $(x_{1}\,x_{2}...x_{n})^{t}$ , $\mathbf{A_{1}}$ , $\mathbf{A_{0}}$ , and $\mathbf{Z}=(z_{1}\,z_{2}...z_{n})^{t}$ constant vectors of dimension $N$ with $\mathbf{A_{1}}\neq\mathbf{A_{0}}\neq 0$ , $C_{0}$ a constant, and $R$ the (closed) region $x_{1}\geq z_{1},x_{2}\geq z_{2},,,x_{N}\geq z_{N}$ .

Define $C$ as the (unique) value satisfying

[TABLE]

Denote $\mathbf{y}=(y_{1}\,y_{2}...y_{N})$ as the point in R at minimal Mahalanobis distance $M(\mathbf{y})$ from the origin with respect to $\Sigma$ (usually, $\mathbf{y}=\mathbf{Z}$ ). Consider all regions $R$ for which the corresponding value of $\mathbf{y}$ lies on the hyperplane $\mathbf{A_{0}^{t}y}+C_{0}^{\prime}=0$ , $C_{0}^{\prime}\neq C_{0}$ . Then

[TABLE]

Proof.

The value $C$ is unique since the function $\Phi(\mathbf{A_{1}^{t}x}+C)$ is continuous and monotonically increasing in $C$ for all $\mathbf{x}$ , and hence so is the integrand (and integral).

We proceed from the formal definition of a limit

[TABLE]

Because $\mathbf{A_{0}^{t}y}+C_{0}^{\prime}=0$ , the right-hand side is equivalent to

[TABLE]

We will show that there exists $Y$ such that $M(y)>Y$ implies that when $C$ takes values at the endpoints of the interval in the integral 29, the integral 26 takes different signs. Since the integral is increasing in $C$ and must be 0, $C$ must lie in the interval in 29 for $M(y)>Y$ .

If $C$ takes the upper value, then at $\mathbf{x}=\mathbf{y}$ , the value of the integrand is

[TABLE]

the sign of which depends on the sign of $C_{0}-C_{0}^{\prime}$ . We shall assume it is positive (with analogous arguments if it is negative). Because $\epsilon>0$ , point $\mathbf{y}$ does not lie on the hyperplane $(\mathbf{A_{1}^{t}}-\mathbf{A_{0}^{t}})\mathbf{x}+(1+\epsilon)(C_{0}-C_{0}^{\prime})-\mathbf{A_{1}^{t}y}-C_{0}=0$ (on which the integrand of 26 is 0). The distance from $\mathbf{y}$ to the hyperplane is given by

[TABLE]

which is independent of $\mathbf{y}$ . Consider a hypersphere centred at $\mathbf{y}$ of radius $d<D$ . Each point in the hypersphere can be expressed as $\mathbf{y}+\mathbf{\kappa}$ with $|\mathbf{\kappa}|\leq d$ , so within the hypersphere we have

[TABLE]

Thus $d$ can be chosen independently of $\mathbf{y}$ such that $\Phi(\mathbf{A_{1}^{t}x}+C)-\Phi(\mathbf{A_{0}^{t}x}+C_{0})$ is bounded below in the hypersphere by a constant $X$ also independent of $\mathbf{y}$ . The function $\Phi(\mathbf{A_{1}^{t}x}+C)-\Phi(\mathbf{A_{0}^{t}x}+C_{0})$ is obviously bounded by $\pm 2$ . Let $R^{\prime}$ be the intersection of $R$ and the hypersphere. The integral 26 now satisfies

[TABLE]

Because $d$ (the radius of the hypersphere) does not depend on $\mathbf{y}$ , by the properties of the Gaussian integral a value $M_{+}$ can be chosen such that $M(y)>M_{+}$ implies that the ratio

[TABLE]

is arbitrarily large (namely, $>2/X$ ), and hence integral 33 is positive. In a similar way, a value $M_{-}$ can be chosen such that if $C$ takes the lower value of interval 29, the integral is negative for $M(y)>M_{-}$ . For $M(y)>max(M_{+},M_{-})$ , the value of $C$ satisfying equation 26 lies within the interval 29, and the limit is established.

∎

Corollary 1.

Given $b,c,y\in\mathbb{R}^{+}$ , define $a$ such that

[TABLE]

then

[TABLE]

and $a>by+c\,\,\forall\,y$

Proof.

We note firstly that the function $\Phi(c)-\Phi(a-bx)$ is increasing for all $x$ . If the integral is 0, the (smooth) integrand must cross 0 at some finite $x\in(y,\infty)$ , and hence its value at $x=y$ must be negative. As $\Phi$ is increasing, we have $\Phi(a-by)>\Phi(c)\implies a>by+c$

The proof of the limit proceeds in a similar way to the proof of the lemma above.

∎

Now (recalling definition 2 in the main paper)

[TABLE]

where

[TABLE]

The asymptotic property of $\beta^{*}$ follows from corollary 1.

If $\gamma=1$ , we have from definition 3 in the main paper

[TABLE]

from which the result follows from an application of lemma 1.

1.2.2 Size of $\beta$ , $\beta^{*}$ and $\beta^{\perp}$

To show that $\beta^{*}<\beta$ , we show that if we set $z_{\beta^{*}}=z_{\beta}$ in the integral 37, then the integral is positive. Since it is decreasing with $z_{\beta}^{*}$ (as $\Phi$ is increasing) we must have $z_{\beta}^{*}>z_{\beta}$ if the integral is to be 0. A similar argument can be used to show that $\beta^{\perp}<\beta^{*}$ . Denote by $I(z_{d},z_{m})$ the value of the integrand of 37 with $z_{\beta^{*}}=z_{\beta}$ .

Consider the line $a_{1}z_{d}+b_{1}z_{m}+c_{1}=a_{0}z_{d}+b_{0}z_{m}+c_{0}$ on the $(z_{d},z_{m})$ plane on which the integrand of 37 is 0. The gradient of this line is

[TABLE]

Since $|n_{0}^{\prime}n_{1}-n_{0}n_{1}^{\prime}|\geq(n_{0}^{\prime}n_{1}-n_{0}n_{1}^{\prime})$ the numerator of the second fraction is greater than or equal to

[TABLE]

and similarly the denominator is greater than or equal to

[TABLE]

so the gradient is positive. If $b_{1}-b_{0}>0$ , $I(z_{d},z_{m})$ is positive if $(z_{d},z_{m})$ falls above the line, and negative if below it; if $b_{1}-b_{0}<0$ , the other way around. Assume for the moment that $b_{1}-b_{0}<0$ .

If the point $(z_{\alpha},z_{\gamma})$ lies above the line, then since $I(z_{d},z_{m})$ is negative in the region $(-\infty,z_{\alpha})\times(z_{\gamma},\infty)$ , we have

[TABLE]

If the point lies below the line, let $z_{\gamma}^{\prime}>z_{\gamma}$ be defined such that the point $(z_{\alpha},z_{\gamma}^{\prime})$ lies on the line. Since $I(z_{d},z_{m})$ is positive in the region $(z_{\alpha},\infty)\times(z_{\gamma},z_{\gamma}^{\prime})$ and negative in the region $(-\infty,z_{\alpha})\times(z_{\gamma},\infty)$ , we have

[TABLE]

so it is sufficient to prove that the integral is positive when the range $(z_{\alpha},\infty)$ is replaced with $(-\infty,\infty)$ . Similar arguments can be used when $b_{1}-b_{0}>0$ , in which case it is sufficient to prove positivity when $z_{\gamma}=0$ .

This enables $z_{d}$ (or $z_{m}$ ) to be integrated out, namely reducing to showing that

[TABLE]

The second part of the integrand is monotonically increasing in $z_{s}$ as $\rho_{sm}>\rho_{rm}$ . Thus the integral is minimised as $z_{\beta}\to-\infty$ , at which the value is $\Phi(z_{\gamma})$ , which is positive.

1.3 SNPs with aberrant allele frequency in one group

1.3.1 $R_{B}<R_{A}$ for SNPs with aberrance in $C_{1}$

If SNPs have aberrant MAF in $C_{1}$ only, we have $E(z_{d})=\zeta_{d}\neq 0$ , $E(z_{m})=\zeta_{m}\neq 0$ and $E(z_{s})=E(z_{r})=0$ . As noted in the main text, as $\zeta_{d}\to 0$ , $R_{B},R_{A}\to P_{0}$ (equation 2 in the main paper) and

[TABLE]

and similarly, $R_{A}\to\frac{\beta}{2}$ , $R_{B}\to\frac{\beta^{*}}{2}$ as $\zeta_{d}\to\pm\infty$ , with $\beta^{*}<\beta$ as shown above. For $\beta^{*}$ defined by 3 in the main paper, we show here that $R_{A}>R_{B}$ for all $\zeta_{d}$ . For the more general definition of $\beta^{*}$ (equation 2 in the main paper), the inequality $R_{B}<R_{A}$ may not hold for all $\zeta_{d}$ . However, in practice, the inequality holds for almost all $\zeta_{d}$ and any deviation is small and near $\zeta_{d}=0$ .

Define the shorthand $N_{\rho}(x,y)$ as the value at $(x,y)$ of the bivariate normal PDF with mean $\left(\begin{smallmatrix}0\\ 0\end{smallmatrix}\right)$ and variance $\left(\begin{smallmatrix}1&\rho\\ \rho&1\end{smallmatrix}\right)$ , and $erfc(x)=2\left(1-\Phi(\sqrt{2}x)\right)$ as the complementary error function. In this section, $\rho$ refers to $\rho_{ds}$ .

Consider the value $R_{A}-R_{B}$ as a function of $\zeta_{d}$ . We will show that the partial derivative $\frac{\delta}{\delta\zeta_{d}}\left(R_{A}-R_{B}\right)$ is zero only when $\zeta_{d}=0$ . Because $R_{A}-R_{B}$ passes through the origin and is asymptotically positive in both directions, it is positive for all $\zeta_{d}\neq 0$ . We have

[TABLE]

Showing that $\frac{\delta}{\delta\zeta_{d}}\left(R_{A}-R_{B}\right)>0$ when $\zeta_{d}>0$ is equivalent to showing that $(a-b)-(pa-qb)>0$ where $a=e^{-\frac{1}{2}(\zeta_{d}-z_{\alpha})^{2}}$ , $b=e^{-\frac{1}{2}(\zeta_{d}+z_{\alpha})^{2}}$ , $p=\frac{erfc\left(\frac{z_{\beta^{*}}+\rho(\zeta_{d}-z_{\alpha})}{\sqrt{2(1-\rho^{2})}}\right)}{erfc\left(\frac{z_{\beta}}{2}\right)}$ and $q=\frac{erfc\left(\frac{z_{\beta^{*}}+\rho(\zeta_{d}-z_{\alpha})}{\sqrt{2(1-\rho^{2})}}\right)}{erfc\left(\frac{z_{\beta}}{2}\right)}$ .

Since $(\zeta_{d}-z_{a})^{2}<(\zeta_{d}+z_{a})^{2}$ for $\zeta_{d}>0$ , we have $a>b$ . Because $erfc$ is strictly decreasing, we have $p<q$ . Because $\frac{\delta p}{\delta\zeta_{d}}<0$ , we have

[TABLE]

where the second inequality arises because $z_{\beta^{*}}>\sqrt{1-\rho^{2}}z_{\beta}+\rho z_{\alpha}$ . Thus $pa-qb<pa-pb=p(a-b)<a-b$ , and $\frac{\delta}{\delta\zeta_{d}}\left(R_{A}-R_{B}\right)>0$ as required.

1.4 Upper bound on $R_{B}-R_{A}$ with aberrance in $C_{1}^{\prime}$

For SNPs with aberrance in $C_{1}^{\prime}$ , we have $E(z_{d})=0$ , $E(z_{r})=\zeta_{r}\neq 0$ , $E(z_{s})=\zeta_{s}\neq 0$ and $E(z_{m})=\zeta_{m}\neq 0$ . As above $R_{A},R_{B}\to P_{0}$ as $\zeta_{r}\to 0$ , and by similar arguments to the section above, $R_{A},R_{B}\to\frac{\alpha}{2}$ as $\zeta_{r}\to\pm\infty$ .

For $\beta^{*}$ defined as per equation 3 in the main paper, it is possible to derive an approximate (asymptotically accurate) upper bound on $R_{B}-R_{A}$ , corresponding to the most serious increase in FPR. The approach is not readily applied to the general definition of $\beta^{*}$ , but again the difference is typically small in practice.

To first order

[TABLE]

Now

[TABLE]

Define $z_{r}^{+}$ , $z_{r}^{-}$ such that

[TABLE]

From equation 25 in Appendix 1.2.1, we have $z_{\beta^{*}}-\zeta_{s}\approx\sqrt{1-\rho^{2}}z_{r}^{-}-\rho z_{\alpha}$ and $z_{\beta^{*}}+\zeta_{s}\approx\sqrt{1-\rho^{2}}z_{r}^{+}-\rho z_{\alpha}$ .

Noting that $\int_{a}^{\infty}\int_{b}^{\infty}N_{0}(x,y)dxdy=\Phi(-a)\Phi(-b)$ and $\Phi(x)=1-\Phi(-x)$ we now have

[TABLE]

Applying the identity $\Phi(-z_{\alpha})=\frac{\alpha}{2}$ and approximations $z_{\beta}^{*}\approx\sqrt{1-\rho^{2}}z_{\beta}+\rho z_{\alpha}$ , $\zeta_{s}\approx k\zeta_{r}$ , yields

[TABLE]

Considered as a function of $\zeta_{r}$ , the value $R_{B}-R_{A}$ will be 0 at $\zeta_{r}=0$ and tend to 0 as $\zeta_{r}\to\pm\infty$ . It will be maximised approximately at the points where $\Phi(z_{\beta}-\zeta_{r})$ or $\Phi(z_{\beta}+\zeta_{r})$ are changing most rapidly; that is, $\zeta_{r}=\pm z_{\beta}$ . At $\zeta_{r}=z_{\beta}$ , the contribution to the value $R_{B}-R_{A}$ from the difference $\Phi\left(z_{\beta}+\frac{k}{\sqrt{1-\rho^{2}}}\zeta_{r}\right)-\Phi(z_{\beta}+\zeta_{r})$ is negligible (and similarly for the other difference when $\zeta_{r}=-z_{\beta}$ ). Using the first-order approximation for $\Phi(z_{\beta}-\zeta_{r})$ about $\zeta_{r}=z_{\beta}$ yields

[TABLE]

In general, this value is substantially less than $\alpha$ .

All instances of ‘approximately equal’ are asymptotic limits as $z_{a}\to\infty$ and $n_{0}$ , $n_{0^{\prime}}$ , $n_{1}$ , $n_{1}^{\prime}\to\infty$ with $z_{0}^{\prime}$ held finite.

1.5 Aberrance in $C_{0}^{\prime}$

For SNPs aberrant in $C_{0}^{\prime}$ , again $E(z_{d})=0$ , $E(z_{r})=\zeta_{r}\neq 0$ , $E(z_{s})=\zeta_{s}\neq 0$ and $E(z_{m})=\zeta_{m}\neq 0$ . As above $R_{A},R_{B}\to P_{0}$ as $\zeta_{r}\to 0$ , and $R_{A},R_{B}\to\frac{\alpha}{2}$ as $\zeta_{r}\to\pm\infty$ . In this case, $R_{B}$ is typically less than $R_{A}$ .

1.6 General aberrance in replication cohorts

For $\beta^{*}$ defined according to 3 in the main paper, the increase in FPR $R_{B}-R_{A}$ for method B for a SNP with aberrance in $C_{1}^{\prime}$ is generally smaller than the decrease in FPR $R_{A}-R_{B}$ for a SNP with a similarly-sized aberrance in $C_{0}^{\prime}$ , in that the integral of the difference over the range of $\zeta_{r}$ is generally smaller in the former.

We define $k$ as in the section above and $k_{1}=\frac{\zeta_{s}}{\zeta_{r}}\big{|}C_{0}^{\prime}\textrm{ aberrant}=\sqrt{\frac{n_{0}^{\prime}(n_{0}^{\prime}+n_{1}^{\prime})}{(n_{0}+n_{0}^{\prime})(n_{0}+n_{0}^{\prime}+n_{1}^{\prime})}}$ similarly. Now for $C_{0}^{\prime}$ aberrant

[TABLE]

Since $\int_{0}^{x}\Phi(z)dz=x\Phi(x)+\frac{1}{\sqrt{2\pi}}\left(e^{-\frac{x^{2}}{2}}-1\right)$ , we have

[TABLE]

Thus with aberrant $C_{0}^{\prime}$

[TABLE]

Comparing $R_{A}$ and $R_{B}$ under the two aberrance scenarios with the same $\zeta_{d}$

[TABLE]

For this to be $>1$ , a necessary condition is $\left(1-\frac{\sqrt{1-\rho^{2}}}{k_{2}}\right)>\left(\frac{\sqrt{1-\rho^{2}}}{k_{1}}-1\right)$ From the definitions of $\rho_{ds}$ (Appendix 1.1), $k$ (equation 50) and $k_{1}$ , this is equivalent to

[TABLE]

The final term in this product is of the form $x+\frac{1}{x}$ so is greater than 2. A sufficient condition is thus

[TABLE]

which holds in most study designs.

Bibliography14

The reference list from the paper itself. Each links out to its DOI / PubMed record.

1[1] Wason JM, Dudbridge F (2012) A general framework for two-stage analysis of genome-wide association studies and its application to case-control studies. The American Journal of Human Genetics 90: 760–773.
2[2] Fuchsberger C, Flannick J, Teslovich TM, Mahajan A, Agarwala V, et al. (2016) The genetic architecture of type 2 diabetes. Nature .
3[3] Lin D, Sullivan PF (2009) Meta-analysis of genome-wide association studies with overlapping subjects. American Journal of Human Genetics 85: 862-872.
4[4] Han B, Duong D, Sul JH, de Bakker PI, Eskin E, et al. (2016) A general framework for meta-analyzing dependent studies with overlapping subjects in association mapping. Human Molecular Genetics 90.
5[5] Bhattacharjee S, Rajaraman P, Jacobs KB, Wheeler WA, Melin BS, et al. (2012) A subset-based approach improves power and interpretation for the combined analysis of genetic association studies of heterogeneous traits. The American Journal of Human Genetics 90: 821–835.
6[6] Zaykin DV, Kozbur DO (2010) P-value based analysis for shared controls design in genome-wide association studies. Genetic Epidemiology 34: 725-738.
7[7] Liley J, Wallace C (2015) A pleiotropy-informed bayesian false discovery rate adapted to a shared control design finds new disease associations from gwas summary statistics. PLOS Genetics .
8[8] Fortune MD, Guo H, Burren O, Schofield E, Walker NM, et al. (2015) Statistical colocalization of genetic risk variants for related autoimmune diseases in the context of common controls. Nature Genetics 47: 839-846.

TL;DR

Contribution

Findings

Abstract

Peer Reviews

Videos

Combining controls can improve power in two-stage association studies

Abstract

Introduction

Results

Overview of method

General properties

Simulation

Recommended applications

Assumptions

Conventional GWAS

Difficult control ascertainment

Prospective study design

Partial replication

Discussion

Acknowledgments

Conflicts of interest

Methods

Definitions

General type 1 error rate

Study sizes, odds ratios and allele frequencies

False ascertainment

Empirical computations

Type 1 error rates

Aberrance in C1C_{1}C1​

Aberrance in C1′C_{1}^{\prime}C1′​

Aberrance in C0′C_{0}^{\prime}C0′​

Aberrance in C0C_{0}C0​

Aberrance in two or more cohorts

Supplementary figures

1 Appendices

1.1 Covariance between Z scores due to shared samples

1.1.1 No covariates or stratification

1.1.2 Z scores with stratification

1.1.3 Z scores with covariates

1.1.4 Practical estimation of covariance

1.2 Properties of β∗\beta^{*}β∗

1.2.1 Asymptotic properties of β∗\beta^{*}β∗

Lemma 1**.**

Proof.

Corollary 1**.**

Proof.

1.2.2 Size of β\betaβ, β∗\beta^{*}β∗ and β⊥\beta^{\perp}β⊥

1.3 SNPs with aberrant allele frequency in one group

1.3.1 RB<RAR_{B}<R_{A}RB​<RA​ for SNPs with aberrance in C1C_{1}C1​

1.4 Upper bound on RB−RAR_{B}-R_{A}RB​−RA​ with aberrance in C1′C_{1}^{\prime}C1′​

1.5 Aberrance in C0′C_{0}^{\prime}C0′​

1.6 General aberrance in replication cohorts

Aberrance in $C_{1}$

Aberrance in $C_{1}^{\prime}$

Aberrance in $C_{0}^{\prime}$

Aberrance in $C_{0}$

1.2 Properties of $\beta^{*}$

1.2.1 Asymptotic properties of $\beta^{*}$

Lemma 1.

Corollary 1.

1.2.2 Size of $\beta$ , $\beta^{*}$ and $\beta^{\perp}$

1.3.1 $R_{B}<R_{A}$ for SNPs with aberrance in $C_{1}$

1.4 Upper bound on $R_{B}-R_{A}$ with aberrance in $C_{1}^{\prime}$

1.5 Aberrance in $C_{0}^{\prime}$