I would like to perform a permutation test on a dataset with several batches and would like to take those batches into account during testing. Similar to a linear model variable ~ explanatory_variable + batch.

To illustrate what I mean, let's consider the following toy example:

There are two main types of lung cancer: Adenocarcionoma (LUAD) and squamous cell carcinoma (LUSC). The gene NKX2-1 is a known marker for LUAD tumor cells. Let's assume it is expressed 1.5x as much in LUAD compared to LUSC. Let's also assume our dataset consists of two batches that (for whatever reason) happen to have different baseline expressions of NKX2-1. In each dataset, the LUAD and LUSC patients are highly imbalanced:

batch # LUSC # LUAD NKX2-1 expr in LUSC NKX2-1 expr in LUAD
1 30 3 4 6
2 3 30 2 3

Ignoring batch effects, simply comparing the mean between datasets will incorrectly lead to the conclusion NKX2-1 were expressed at a higher level in LUSC. That's why usually I would use a linear model NKX2-1 ~ tumor_type + batch.

Unfortunately, my real data are not gene expression (for which established statistical models exist) and violate the normality assumption, which is why I would like to perform a permutation test instead. Is there an established way to perform a permutation test that takes batches (or other covariates) into account?

Here's some Python code playing with above example:

import numpy as np
import pandas as pd
import scipy.stats
import statsmodels.formula.api as smf

# for the sake of the example, use a normal distribution. 
# The real data violates normality assumption
lusc1 = np.random.normal(loc=4, scale=1, size=30)
luad1 = np.random.normal(loc=6, scale=1, size=3)
lusc2 = np.random.normal(loc=2, scale=1, size=3)
luad2 = np.random.normal(loc=3, scale=1, size=30)

df = pd.DataFrame().assign(
    expr=np.hstack([lusc1, luad1, lusc2, luad2]),
    tumor_type=["LUSC"] * 30 + ["LUAD"] * 3 + ["LUSC"] * 3 + ["LUAD"] * 30,
    dataset=["batch1"] * 33 + ["batch2"] * 33,

# Compute fold-change of LUSC compared to LUAD
np.mean(df.loc[lambda x: x["tumor_type"] == "LUAD", "expr"]) / np.mean(
    df.loc[lambda x: x["tumor_type"] == "LUSC", "expr"]
# Huh, so the expression is actually less in LUAD!? 
model = smf.ols("expr ~ C(tumor_type, Treatment('LUSC')) + dataset", data=df)
res = model.fit()

Using a linear model, I correctly get a positive coefficient for tumor_type. But how can I achieve the same with a permutation test?

OLS result


1 Answer 1


It's possible to use randomization tests with regressions instead of relying on asymptotic distributions. There was recently a paper on this in the Quarterly Journal of Economics (link to ungated version) and the author provides some Stata code. However, Simonsohn subsequently showed that randomization tests did not outperform the appropriate robust method for computing standard errors.

  • 1
    $\begingroup$ thanks, this got me on the right track! I think the robust linear model is what I actually want! $\endgroup$ Mar 28 at 18:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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