<jupyter_start><jupyter_text>In this notebook, we are going to show how pycombat works. This package implements Combat, a technique for data harmonisation based on a linear mixed model in which location and scale random effects across batches are adjusted using a bayesian approach (Johnson, 2007):$$ Y_{ijk} = \alpha_k + X\beta_k +  \gamma_{ik} + \delta_{ik}\epsilon_{ijk}$$Original Combat tecnique also allowed to keep the baseline effects $\alpha_k$ and the effects of interest $\beta_k$ by reintroducing these after harmonisation:$$Y \longrightarrow Y^{combat}_{ijk} = \frac{Y_{ijk} - \hat{\alpha}_k - X\hat{\beta}_k - \hat{\gamma}_{ik}}{\hat{\delta}_{ik}} +  \hat{\alpha}_k + X\hat{\beta}_k $$One extension of this python package is the possibility of removing unwanted variables' effect by no reintroducing them again. Using the same linear mixed model of the begining, we now separte source of covariation $C$ from sources of effects of interest $X$:$$ Y_{ijk} = \alpha_k + X\beta_k^{x} + C\beta_c^{c}  +  \gamma_{ik} + \delta_{ik}\epsilon_{ijk}$$And then in this case, combat adjustment will be given by:$$Y \longrightarrow Y^{combat}_{ijk} = \frac{Y_{ijk} - \hat{\alpha}_k - X\hat{\beta}_k^{x} - C\hat{\beta}_c^{c} - \hat{\gamma}_{ik}}{\hat{\delta}_{ik}} +  \hat{\alpha}_k + X\hat{\beta}_k^{x}  $$Such modification to combat has been recently proposed and applied by some authors (Wachinger, 2020).In this notebook, we are going to show how this package works by making use of a data on gene expression measurements from a bladder cancer study (Dyrskjot, 2004). We will then compare our results with those from neuroCombat (Fortin, 2017), a known python implementation of Combat (though this does not include the modification we have implemented).*References*:- <NAME>, <NAME>, <NAME>, Adjusting batch effects in microarray expression data using empirical Bayes methods, Biostatistics, Volume 8, Issue 1, January 2007, Pages 118–127, https://doi.org/10.1093/biostatistics/kxj037- <NAME>, <NAME>, <NAME>, <NAME>, <NAME>,<NAME>, and <NAME>. Gene expression in the urinary bladder: acommon carcinoma in situ gene expression signature exists disregardinghistopathological classification.Cancer Res., 64:4040–4048, Jun 2004.- <NAME>, <NAME>, <NAME>. Detect and Correct Bias in Multi-Site Neuroimaging Datasets. arXiv:2002.05049 - <NAME>., <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME>, <NAME> and <NAME> (2017). "Harmonization of cortical thickness measurements across scanners and sites." Neuroimage 167: 104-120.<jupyter_code>import numpy as np
import pandas as pd
import sys
from pycombat import Combat<jupyter_output><empty_output><jupyter_text>Following the spirit of scikit-learn, Combat is a class that includes a method called **fit**, which finds the fitted values of the linear mixed model, and **transform**, a method that used the previously learning paramters to adjust the data. There is also a method called **fit_transform**, which concatenates both methods. So the first thing that you need to do is to define a instance of this class<jupyter_code>combat = Combat()<jupyter_output><empty_output><jupyter_text>When you define the instance, you can pass it the following parameters:    - method: which is either "p" for paramteric or "np" for non-parametric (not implemented yet!!)    - conv: the criterion to decide when to stop the EB optimization step (default value = 0.0001) Now, you have to call the method **fit**, passsing it the data. These data will consist of the following ingredients:    - Y: The matrix of response variables, with dimensions [observations x features]    - b: The array of batch label for each observation. In principle these could be labelled as numbers or strings.    - X: The matrix of effects of interest to keep, with dimensions [observations x features_interest]    - C: The matrix of covariates to remove, with dimensions [observations x features_covariates]    ***Important:*** If you have effects of interest or covariates that involve categorical features, make sure that you drop the first level of these categories when building the independent matrices, otherwise they would be singular. You can easily accomplished this using the pandas and **pd.get_dummies** with the option *drop_first* checked. Let's then see how this package works<jupyter_code># Y is the matrix of response variables
Y = np.load('bladder-expr.npy')

print("the matrix of response variables has %d observations and %d outcome variables" % (Y.shape[0], Y.shape[1]))
# This loads the set of independent variables, including the batch labels
pheno = pd.read_csv('bladder-pheno.txt', delimiter='\t')

pheno.head()
b = pheno.batch.values

print("We have %d different batches" % len(np.unique(b)))<jupyter_output>We have 5 different batches<jupyter_text>We also have information about the type of cancer and age in this data<jupyter_code>pheno.cancer.value_counts()
pheno.describe()['age']<jupyter_output><empty_output><jupyter_text>Say we want to keep these two effects after combat harmonisation. We can then build our matrix X from these two variables<jupyter_code>X = np.column_stack((pd.get_dummies(pheno.cancer.values, drop_first=True).values,
                    pheno.age.values))

print(X[:10,:])<jupyter_output>[[0 1 1]
 [0 1 2]
 [0 1 3]
 [0 1 4]
 [0 1 5]
 [0 1 6]
 [0 1 7]
 [0 1 1]
 [1 0 2]
 [1 0 3]]<jupyter_text>And now we **fit** the combat model with these pieces<jupyter_code>combat.fit(Y, b, X=X, C=None) # X and C are None by default, so no need here to write C=None<jupyter_output><empty_output><jupyter_text>And then, we can adjust this dataset calling the **transform** method<jupyter_code>Y_combat = combat.transform(Y=Y, b=b, X=X)<jupyter_output><empty_output><jupyter_text>We can check that this gives the same result as that applying neuroCombat<jupyter_code>from neuroCombat import neuroCombat
discrete_cols = ['cancer']
continuous_cols = ['age']
batch_col = 'batch'

Y_neurocombat = neuroCombat(data=Y,
                            covars=pheno,
                            batch_col=batch_col,
                            discrete_cols=discrete_cols,
                            continuous_cols=continuous_cols)
np.corrcoef(Y_combat.flat, Y_neurocombat.flat)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10441830896376351, 'fraction_numerical': 0.01795931341385887, 'mean_word_length': 4.789328426862926, 'pattern_counts': {'":': 0, '<': 59, '<?xml version=': 0, '>': 59, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '8', 'dir': 'jupyter-structured-clean-dedup', 'id': '160603', 'license': 'MIT', 'n_tokens_mistral': 1962, 'n_tokens_neox': 1800, 'n_words': 745, 'path': 'test/test.ipynb', 'repo_name': 'CoAxLab/pycombat', 'size': '6465'}