Code
reticulate::use_python('/usr/local/bin/python3')Logistic SER coverage under polygenic model
We are interested in knowing if the inferences under the logistic SER are biased when phenotypes are simulated under a polygenic model. In these simulations we simulate one causal variant in the locus, and model unlinked genetic effects with a random intercept. The SER fits the ordinary/marginal regression model at each step, without accounting for this extra variability in liability. Nonetheless, we see good coverage of the credible sets.
Simulation
We simulate binary \(X\) under a markov model. The first column of \(X\) are iid \(Bin(1, p)\) to generate \(x_{k+1} | x_{k}\) we sample \(\pi\) indices and flip the bits. By setting \(\pi\) to a small value we get highly correlated \(X\), or by setting \(\pi =0.5\) we get uncorrelated \(X\).
To simulate our phenotypes, we specify \(h^2\) of liability for the trait, and \(\rho\) the fraction of \(h^2\) explained by the causal variant in the locus. We fix \(h^2 = 0.2\) and vary \(\rho\) from \(0.02\) to \(1.0\) (corresponding to odds ratios of \(\approx 1.1 - 2\)). For each simulation we set the intercept to achieve the desired fraction of cases \(k=0.1\). Note the effect sizes in these simulations are not random, they are computed deterministically.
These simulation settings are meant to be a reasonable approximation to what we might see in GWAS. But admitidely, the way we generate \(X\) is a little odd.
We can improve these simulations by: 1. Simulating more realistic genotype type data 1. Being a bit more thoughtful about how we pick effect sizes e.g. have an MAF dependent prior OR.
Results
We see good coverage of the \(95%\) CSs across the different simulation scenarios. This is in contrast to earlier simulations, which I’ve been having trouble replicating, where we see a dramatic loss of coverage. It could be that choosing a more realistic genotype would change the results here, and we should follow up on that.
Code
import numpy as np
import pandas as pd
from susiepy.logistic_susie import logistic_ser
import itertoolsCode
import numpy as np
def flip_bits(x, flip):
"""
Flips the values of x based on the flip array.
Parameters:
x (ndarray): Input array
flip (ndarray): Array of 0s and 1s indicating which values to flip
Returns:
ndarray: Flipped array
"""
return np.where(flip == 1, 1 - x, x)
def sim_binary_X(n: int=1000, p: int=100, pi1: float = 0.1, rho=0.1):
"""
n : number of rows
p: number of columns
pi1: fraction of ones in first column
rho: probability of flipping bits
"""
X = np.zeros((n, p))
X[:, 0] = np.random.choice([0, 1], size=n, p=[1 -pi1, pi1])
for j in range(1, p):
flip = np.random.binomial(1, rho, size=n)
X[:, j] = flip_bits(X[:, j-1], flip)
return X
def augment(X, y):
yaug = np.append(y, [1., 1., 0., 0.])
Xaug = np.column_stack([np.append(X[:, i], [1., 0., 1., 0.]) for i in range(X.shape[1])])
Xaug = (Xaug - Xaug.mean(0)[None])/Xaug.std(0)[None]
return Xaug, yaug
def fit_ser_aug(X, y):
Xaug, yaug = augment(X, y)
ser = logistic_ser(Xaug, yaug)
return ser
def compute_cs(alpha, target=0.95):
sorted_indices = np.argsort(alpha)[::-1]
cumulative_sum = np.cumsum(alpha[sorted_indices])
top_indices = sorted_indices[:(cumulative_sum < target).sum() + 1]
return top_indices
def simulate_case_control(x, k=0.1, h2=0.2, rho=0.1):
nu = np.pi**2/3 * h2/(1-h2)
sigma = np.sqrt(rho * nu)
s = np.sqrt(nu - sigma**2)
xnorm = (x - x.mean()) / x.std()
b = sigma
# b = np.random.normal(scale=sigma)
logit = xnorm * b + np.random.normal(scale=s, size=n)
b0 = np.quantile(logit + np.random.logistic(size=n), k)
logit = b0 + logit
y = np.random.binomial(1, 1/(1 + np.exp(-logit)))
sim = dict(k=k, h2=h2, rho=rho, nu=nu, sigma=sigma, s=s, b0=b0, b=b, logit=logit, y=y)
return sim