Fitting hierarchical logistic regression to large dataset

I am a molecular biologist and I am investigating the effect of genetic variants on a disease. My data consist of 2000 patients, 1000 with the disease in question (asthma) and 1000 who are healthy. A binary outcome. For each patient I have the data of 100000 variants (called SNPs), which takes on 0, 1 or 2 depending on the number of copies the individual posses of each SNP. I want to estimate the effect of each SNP, so I assign a coefficient b to each SNP using a multivariate logistic regression but each SNP is drawn from an adaptive prior, making it a hierarchical model. I am using a standard normal as a prior for b, since I expect most SNPs to have no effect on the disease. My problem is that I quickly run out of memory while fitting the model to this large dataset. I wonder if it is because my model is not constructed well enough, or if it is simply my computer which is not powerful enough (in which case I cannot do anything).

My computer specifications are: Intel® Core™ i5-7200U CPU @ 2.50GHz, with 8 GB RAM.

I am using 500 iterations and using rstans options(mc.cores = 4), as I only have 4 cores.

Do you have any suggestions of what I could do, in order to be able to fit this model to a the large dataset described above? My model is described below:

  int<lower=1> N; // Rows
  int<lower=1> M; // Columns
  int<lower=0, upper=1> y[N]; // Outcome variables
  matrix<lower=0, upper=2>[N, M] x; // Predictor variables

  // Hyper priors
  real mu;
  real<lower=0> sigma;
  // Priors
  real a;
  vector[M] b;

  // Hyper-priors
  mu ~ std_normal();
  sigma ~ cauchy(0, 5);
  // Priors
  a ~ normal(0, 5);
  b ~ normal(mu, sigma);
  // Likelihood
  y ~ bernoulli_logit(a + x * b);

generated quantities{
  vector[M] OR;
  OR = exp(b); // effect size of each SNP

If the problem is running out of memory, the easiest thing to do is to run one chain at a time. This will slow down the overall time to get the solution, but it will give you a better chance of fitting such a large dataset in your memory.

Another thing is to compute the odds ratios outside of stan, so you’ll save up a bit more memory.

You could then look at using bernoulli_logit_glm which should speed up the sampling a bit.

The model itself looks to me as simple as it gets, but perhaps somebody else will be able to offer some advice on that.


Thinking a bit more outside the box, have you looked into ways of reducing the number of SNPs, for example by LD clumping? Or by removing some of the least associated SNPs from large meta-analyses (such as, assuming the data you are using is not part of that)?

Hey! So, N is 2000 and M is 100000? Yeah, that’s kinda large. Did you try excluding the generated quantities block? This only eats up memory and you can easily compute the OR afterwards.

Also try a single chain with as few iterations as sensible (not less, though) as @mcol suggested. You can also try if the save_warmup = FALSE option in the stan call helps.

Good luck! :)

If memory is the concern, then cmdstan (or one of the interfances like cmdstanR) could work, since the output is written straight to the disk rather than being stored in memory.

1 Like

Would cmdstanr also do the trick?

Have you tried running this with fake data. Like 100 patients total and maybe 1000 SNPs? And see how long that takes to run and how much memory it uses up. We have the same issue with microbiome data in our lab.

I would think so.