From dc0219a33e99c61f7597d6f91997fb6ee96e9a07 Mon Sep 17 00:00:00 2001 From: coldfire79 Date: Mon, 19 Sep 2022 15:55:39 -0700 Subject: [PATCH 1/2] add in np.var to get the unbiased estimator of the variance as in the original paper --- combat/pycombat.py | 6 +++--- combat/test_unit.py | 12 ++++++------ 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/combat/pycombat.py b/combat/pycombat.py index 8f01dda..650f26f 100644 --- a/combat/pycombat.py +++ b/combat/pycombat.py @@ -88,7 +88,7 @@ def compute_prior(prior, gamma_hat, mean_only): if mean_only: return 1 m = np.mean(gamma_hat) - s2 = np.var(gamma_hat) + s2 = np.var(gamma_hat, ddof=1) if prior == 'a': return (2*s2+m*m)/s2 elif prior == 'b': @@ -521,11 +521,11 @@ def fit_model(design, n_batch, s_data, batches, mean_only, par_prior, precision, else: for i in batches: # feed incrementally delta_hat list_map = np.transpose(np.transpose(s_data)[i]).var( - axis=1) # variance for each row + axis=1, ddof=1) # variance for each row delta_hat.append(np.squeeze(np.asarray(list_map))) gamma_bar = list(map(np.mean, gamma_hat)) # vector of means for gamma_hat - t2 = list(map(np.var, gamma_hat)) # vector of variances for gamma_hat + t2 = gamma_hat.var(axis=1, ddof=1).flatten().tolist()[0] # vector of variances for gamma_hat # calculates hyper priors for gamma (additive batch effect) a_prior = list( diff --git a/combat/test_unit.py b/combat/test_unit.py index b0a148c..29c13db 100644 --- a/combat/test_unit.py +++ b/combat/test_unit.py @@ -34,13 +34,13 @@ ########## # import function used for unit testing -from .pycombat import model_matrix, all_1 +from pycombat import model_matrix, all_1 #covariate_model_matrix -from .pycombat import compute_prior, postmean, postvar, it_sol, int_eprior -from .pycombat import check_mean_only, define_batchmod, check_ref_batch, treat_batches, treat_covariates, check_NAs -from .pycombat import calculate_mean_var, calculate_stand_mean -from .pycombat import standardise_data, fit_model, adjust_data -from .pycombat import pycombat +from pycombat import compute_prior, postmean, postvar, it_sol, int_eprior +from pycombat import check_mean_only, define_batchmod, check_ref_batch, treat_batches, treat_covariates, check_NAs +from pycombat import calculate_mean_var, calculate_stand_mean +from pycombat import standardise_data, fit_model, adjust_data +from pycombat import pycombat ########## print("\n#### Unit Testing for pyComBat ####\n") From 9cabb302662efeaa41fbdca6c9203120c0631e73 Mon Sep 17 00:00:00 2001 From: coldfire79 Date: Mon, 19 Sep 2022 16:00:41 -0700 Subject: [PATCH 2/2] revert the test file --- combat/test_unit.py | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/combat/test_unit.py b/combat/test_unit.py index 29c13db..b0a148c 100644 --- a/combat/test_unit.py +++ b/combat/test_unit.py @@ -34,13 +34,13 @@ ########## # import function used for unit testing -from pycombat import model_matrix, all_1 +from .pycombat import model_matrix, all_1 #covariate_model_matrix -from pycombat import compute_prior, postmean, postvar, it_sol, int_eprior -from pycombat import check_mean_only, define_batchmod, check_ref_batch, treat_batches, treat_covariates, check_NAs -from pycombat import calculate_mean_var, calculate_stand_mean -from pycombat import standardise_data, fit_model, adjust_data -from pycombat import pycombat +from .pycombat import compute_prior, postmean, postvar, it_sol, int_eprior +from .pycombat import check_mean_only, define_batchmod, check_ref_batch, treat_batches, treat_covariates, check_NAs +from .pycombat import calculate_mean_var, calculate_stand_mean +from .pycombat import standardise_data, fit_model, adjust_data +from .pycombat import pycombat ########## print("\n#### Unit Testing for pyComBat ####\n")