What is a good strategy for performing permutation test for coefficients of Generalized Linear Models (GLM) such as logistic regression? Specially when the response is binomial.

Following the questions here and here, and @Geln_b's informative answer (the former link), I was wondering if one could apply the same concept to logistic regression? If not, is there any other way to perform permutation test?

To provide more information about the problem I'm dealing with:

I have case-control cohorts of healthy/diseased individuals with respect to different diseases. The response variable $Y$ is binomial represting whether individual is healthy, and the predictors are genetic mutations.

For example, here is a model one would try to fit that includes pairwise interactions:

$Y = \beta_0 + \beta_1 \; mutation_1 + \beta_2 \; mutation_2 \; + \beta_3 \; mutation_1 * mutation_2 + \beta_4 \; age \; + \beta_5 \; gender \; + \beta_6 \; PC_1 + ... + \beta_9 \; PC_4 \;$

$PC_n$ are first four/five principal components to correct for population sub-structure e.g. ancestry (there are better ways to handle that, but doing this for now).

I'm interested in the interaction coefficient, and find it hard to set up a permutation strategy as I will be testing thousands of pairs of mutations separately (thousands of separate models - the matrix is fat $M << N$). And, I often hear that there is no solid/valid permutation strategy for interactions.

I was hoping to do something like what @GregSnow suggested here for a set of pairwise interactions (permutation distribution of multiple interactions in separate models as a multivariate gaussian distribution):

One possibility is to compute the mean vector from your permutation distribution (should be close to a vector of 0's) and the variance/covariance matrix. Then use that information to compute the Mahalanobis distance for all the points, the p-value would then be the proportion of points with a Mahalanobis distance greater than the original coefficients.

additional info:

Since I'm dealing with diseases, I could use disease prevalence in population, and form the response as a probability instead of a binary variable $Y$. For example, disease A occurs in 0.1% of population, can we somehow use this as the response (like a probability of getting disease) instead of binomial variable, and would it make sense?


1 Answer 1


Testing mutations separately has in large part been a complete failure in genomics. Comprehensive models are recommended, e.g., specifying mutation effects as random effects in a hierarchical (Bayesian or otherwise) model. The two-at-a-time approach you are using will likely result in small non-biologically useful effects upon validation.

If you insist on two-at-a-time testing, this is a setting for the use of Rao's efficient score test and doesn't require expensive permutation tests unless you meant the permutation test to adjust for simultaneous testing of all mutations not just the two. For the score test you fit a single model with all the adjustment variables, then use the score equations to get an instant 3 d.f. test (main effects of two mutations + their interaction) without refitting the adjustment effects. There is a paper in the genomics literature somewhere that I've seen. I wish I had the reference.

  • 1
    $\begingroup$ Thanks very much for your response. I realize the problem you mentioned with the separate testing. I was hoping I could form the problem as an anomaly detection for interaction coefficient. Specifically, by performing permutation, I could compute the covariance matrix for interaction coefficients, even though they are in separate models. This way I could test interaction between sets of mutations... $\endgroup$
    – NULL
    Dec 31, 2017 at 18:27
  • $\begingroup$ ...which could be an alternative to existing methods for gene-set analysis, etc. Existing models are for gene-set analysis of marginal effects, and for interactions we don't have good methods available to my knowledge. $\endgroup$
    – NULL
    Dec 31, 2017 at 18:31
  • $\begingroup$ By putting them in separate models I'm not sure how you'd model the dependence between models. $\endgroup$ Jan 1, 2018 at 18:03
  • 2
    $\begingroup$ The continual creation of somewhat ad hoc frequentist approaches as opposed to spending time specifying full comprehensive models has not served us well IMHO. $\endgroup$ Jan 2, 2018 at 12:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.