To evaluate the effective sample size required to have power in LDSR analyses, we consider the relationship between \(N_{eff}\) and the observed standard error (SE) for the LDSR \(h^2_g\) estimate for the 4178 Round 2 phenotypes.

Takeaway: There’s a clear inverse relationship between \(N_{eff}\) and the SE of the \(h^2_g\) estimate. The SE is also very large (e.g. wider than the 0-1 range for \(h^2\)) at small effective sample sizes.

For both visualization and modelling, it’s useful to look at this relationship in terms of the inverse variance (\(1/SE^2\)):

Takeaway: There’s a nearly linear relationship between \(N_{eff}\) and \(1/SE^2\). This relationship is highly heteroskedastic, with much greater variability in the inverse variance when \(N_{eff}\) is large.

We can model this observed relationship in order to establish the “expected” inverse variance as a function of \(N_{eff}\). Given the observed heteroskedasticity, and that our primary interest is in the error variance at small sample sizes, we consider a weighted regression with weights inversely proportional to sample size. Based on the plot above we also allow for quadratic curvature in the relationship with \(N_{eff}\), and also compare regression to a loess fit.

Note: Orange = regression, green = loess.

Of these, the loess fit appears to better fit the overall trend of the data, especially when focusing on smaller sample sizes (as is visible when zooming the above plot). Note the smaller sample sizes are the range relevant to the current question of the minimum viable sample size for ldsc.

We adopt the loess model (green) for predicting SEs from sample size for the rest of this section.

Given the above model for expected SEs as a function of sample size, and assuming the SEs are well calibrated, we can then turn these modelled SEs into a power analysis.

Specifically, we ask for a given true SNP heritability, sample size and p-value threshold what the probability is of the estimated \(h^2_g\) divided by its SE exceeding to corresponding critical value for the Z statistic, assuming the sampling variation of the estimated estimated \(h^2_g\) and it’s estimated SE both match the expected SE based on sample size given by the above model. We focus on here power at \(p < .05\) without any multiple testing correction in the interest of providing results relevant to hypotheses that include look-ups of \(h^2_g\) for a single phenotype (as opposed to scans across all available UKB phenotypes). We’ll return to the question of which results we’d consider significant in the context of looking at all 4178 phenotypes later.

Note that this is a bit more empirical than a formal power analysis. We’re relying on the average observed SE reported for the ldsc \(h^2_g\) estimates rather than some theoretical expectation for the SE. There’s also likely a relationship between the SE and the true \(h^2_g\) (and corresponding genetic architecture), as evidenced by the previous plots of precision vs. sample size colored by \(h^2_g\) estimate, that we do not account for in this power estimate. Nevertheless, we use this modelled estimate of power as a useful rough benchmark for evaluating what range of sample sizes to include in this \(h^2_g\) analysis.

Now that we have rough estimates of power as a function of \(h^2_g\) and sample size, the question is what \(h^2_g\) is relevant to have power to detect.

If we look at the UKB phenotypes with the larges sample sizes (i.e. \(N_{eff}\)>300k, where there’s no worry about power) we find that the vast majority of phenotypes have \(h^2_g\) estimates \(\leq 0.3\), with a somewhat bimodal distribution (overall mean=0.118, median=0.078).

From the observed \(h^2_g\) distribution, is appears the phenotypes with larger estimated SNP heritability tend to have estimates around \(h^2_g=0.2\) or slightly larger. We can then orient our decisions on sample size around having power to detect \(h^2_g \geq 0.2\).

One possible hypothesis is that the change in the shape of attenuation reflects the changing mix of phenotypes in the Round 2 GWAS release. In particular, there’s a large number of additional diet items with \(N_{eff} \approx 50,000\) and mental health items at \(N_{eff} \approx 120,000\). It’s possible these subsets of items are influencing the shape of the attentuation curve by confounding true \(h^2_g\) with \(N_{eff}\). If we restrict to variables observed in nearly all samples:

Note: zoom out for \(h^2_g\) outliers.

Takeaway: the milder attentuation of \(h^2_g\) as a function of \(N_{eff}\) is still observed when restricting to 963 phenotypes with \(N > 300,000\). Similar trends are observed for subsetting on other dimensions (not shown).

Conclusion: The mix of phenotypes doesn’t appear to explain the weaker attentuation effect.

Instead of looking across phenotypes, we can also consider the impact of sample size on \(h^2_g\) estimates within a phenotype. We evaluate this by looking at how the LDSR estimates change when subsetting individuals from UKB for phenotypes with high \(N\)s and strong \(h^2_g\) estimates.

Specifically, for a given downsampling experiment we split UKB individuals into 300 chunks, perform a GWAS within each chunk (controlling for the standard GWAS covariates), and then meta-analyze increasing numbers of chunks. We then run LDSR for each of these meta-analyses to assess how the \(h^2_g\) estimates change with growing meta-analyses.

Note: Above are two specific examples among a broader set of ongoing analyses. Plots by Nikolas Baya.

Takeaway: There are signs of attenutation of the LDSR \(h^2_g\) estimate at low \(N\) in height, but leg impedence shows no such attentuation, and instead may have an upward bias. In both cases, estimates are clearly unstable at low \(N\).

Conclusion: These results are preliminary, and appear fairly unstable across phenotypes, but suggest that attenuation may occur irregularly across phenotypes and than in some instances there may be upward bias rather than downward attenuation at low sample sizes.

The downsampling experiments, while inconclusive, do suggest instability of the LSDR \(h^2_g\) at low \(N\), with variability possibly exceeding the SE estimates. Noteably, there were signs of attenuation at low \(N\) for height, a phenotype where the GWAS is likely affected by population stratification (intercept = \(1.313\), p = \(9.29\times 10^{-10}\), ratio = \(0.084\)), but not leg impedence where stratification is weaker (intercept = \(1.096\), p = \(5.53\times 10^{-5}\), ratio = \(0.053\)).

This leads us to evaluate whether the relationship between sample size and LDSR \(h^2_g\) varies as a function of stratification. We start by noting that the LDSR intercept is expected to estimate \(1+Na\) where \(N\) is sample size and \(a\) is an index of population stratification or other counfounding that is, under a simple stratification model, proportional to \(F_{st}\). It follows that \(a=(intercept-1)/N\) should give an estimate of stratification that is invariant to sample size.

On that basis, we look at LDSR \(h^2_g\) estimates as a function of \(N_{eff}\) split by deciles of \(a=(intercept-1)/N\).