diff --git a/python/src/manocca.py b/python/src/manocca.py
index 21c8f0e1a592b3120c14ac0266656561a16f92be..7ee039733fb31d48a3524b1522dae230053422de 100644
--- a/python/src/manocca.py
+++ b/python/src/manocca.py
@@ -110,11 +110,13 @@ class MANOCCA:
         self.predictors = predictors
         self.cols_predictors = cols_predictors
         self.predictors, self.cols_predictors = pt._extract_cols(self.predictors, self.cols_predictors)
+        self.predictors = scale(self.predictors)
 
         self.covariates = covariates
         self.cols_covariates = cols_covariates
         if covariates is not None :
             self.covariates, self.cols_covariates = pt._extract_cols(self.covariates, self.cols_covariates)  
+            self.covariates = scale(self.covariates)
 
         # Parameters
         if n_comp :
@@ -141,6 +143,9 @@ class MANOCCA:
         # Compute prodV_red if not given
         if prodV_red is not None:
             self.prodV_red = prodV_red
+            if use_resid and isinstance(covariates, np.ndarray):
+                print("adjusting predictor")
+                self.predictors = np.apply_along_axis(lambda x : pt.adjust_covariates(x,covariates), axis = 0, arr = self.predictors)
 
         else : # If not we preprocess the data and compute prodV and prodV_red
             self.prodV = self.get_prodV_para_wrap(self.outputs)
@@ -158,7 +163,13 @@ class MANOCCA:
 
             if use_resid and isinstance(covariates, np.ndarray):
                 print("computing residuals")
+                print("adjusting outputs")
                 self.prodV_red = np.apply_along_axis(lambda x : pt.adjust_covariates(x,covariates), axis = 0, arr = self.prodV_red)
+                print("adjusting predictor")
+                self.predictors = np.apply_along_axis(lambda x : pt.adjust_covariates(x,covariates), axis = 0, arr = self.predictors)
+                # print("adjusting predictor")
+                # self.predictors = np.apply_along_axis(lambda x : pt.adjust_covariates(x,covariates), axis = 0, arr = self.predictors)
+                # print(self.predictors.mean())
             
     ### Functions to generate prodV and prodV_red ###
 
@@ -287,6 +298,7 @@ class MANOCCA:
         if isinstance(var, type(None)) :
             var = self.predictors
 
+
         n_preds = var.shape[1]
         p = np.empty((n_preds,1))
         for i_pred in range(n_preds):
@@ -297,7 +309,7 @@ class MANOCCA:
                 p[i_pred] = custom_mancova(scale(self.prodV_red[~nan_mask,min_comp:n_comp]), scale(var[~nan_mask,i_pred].reshape(-1,1)), scale(self.covariates[~nan_mask,:]))
             else :
                 # print("Without covariates")
-                p[i_pred] = custom_manova(self.prodV_red[~nan_mask,min_comp:n_comp], scale(var[~nan_mask,i_pred].reshape(-1,1)))
+                p[i_pred] = custom_manova(self.prodV_red[~nan_mask,min_comp:n_comp], var[~nan_mask,i_pred].reshape(-1,1))
         self.p = p
 
 
diff --git a/python/src/manova.py b/python/src/manova.py
index 345a191d7758a7a034de55a2ddcbb8faa0f8f2b8..a8d56eebdf73951b5ac0c05a0d9fa4ad724425d5 100644
--- a/python/src/manova.py
+++ b/python/src/manova.py
@@ -35,6 +35,8 @@ class MANOVA :
         if use_resid and isinstance(covariates, np.ndarray):
             print("computing residuals")
             self.outputs = np.apply_along_axis(lambda x : pt.adjust_covariates(x,covariates), axis = 0, arr = self.outputs) 
+            print("adjusting predictor")
+            self.predictors = np.apply_along_axis(lambda x : pt.adjust_covariates(x,covariates), axis = 0, arr = self.predictors)
 
 
 
diff --git a/python/src/simulation.py b/python/src/simulation.py
index 2530f17e4e0fdb8a7c7c208eaf1b8af1570b7afa..6206e9ca50732466f94c0f5c0d71c388960edff1 100644
--- a/python/src/simulation.py
+++ b/python/src/simulation.py
@@ -77,7 +77,6 @@ class Simu :
         return self.L_models[num]
 
     def run_simu(self, n_comp, var = None):
-
         self.p = None
         for m in self.L_models :
             m.test(var = var, n_comp = n_comp)
@@ -113,6 +112,10 @@ class Simu :
     def _null_loop(self, m, var_null, N_iters, n_comp = None ) :
         p_null = np.empty((N_iters, 1))
         for i in range(N_iters):
+            # print(var_null)
+            # print(var_null.flags)
+            # var_null.setflags(write=1)
+            # print(var_null.flags)
             np.random.shuffle(var_null)
             # print(var_null.shape)
             m.test(var_null, n_comp = n_comp)
@@ -120,10 +123,10 @@ class Simu :
         return p_null
 
     def simu_null(self, var_name, model, Nsimu = 1000, n_comp = None):
-        var_pos = self.cols_predictors.index(var_name)
-        var_null = self.predictors[:,var_pos].copy()
-        var_null = var_null.reshape(-1,1)
         m = self.L_models[self.methods.index(model)]
+        var_pos = m.cols_predictors.index(var_name)
+        var_null = m.predictors[:,var_pos].copy()
+        var_null = var_null.reshape(-1,1)
 
         if self.n_jobs == 1 : 
             return self._null_loop(m, var_null, Nsimu, n_comp)
diff --git a/python/src/tools/compute_manova.py b/python/src/tools/compute_manova.py
index 34d1cbc471af1f1a5e1c08f4eb75dd53734ff781..589199aeac6c8c9a2468580d1c07b171aaa441f9 100644
--- a/python/src/tools/compute_manova.py
+++ b/python/src/tools/compute_manova.py
@@ -18,9 +18,11 @@ def custom_manova(Y,X,C=None, return_beta = False):
     beta = linear_regression(Y, X, C)
     dot_Y = np.dot(Y.T,Y)
 
-    (sign_num, num) = np.linalg.slogdet(dot_Y-X.shape[0]*np.dot(beta,beta.T))
+    # (sign_num, num) = np.linalg.slogdet(dot_Y-X.shape[0]*np.dot(beta,beta.T))
+    (sign_num, num) = np.linalg.slogdet(dot_Y - np.dot( np.dot(beta, X.T) , np.dot(X, beta.T)))
     (sign_denom, denom) = np.linalg.slogdet(dot_Y)
     lamb = np.exp(sign_num*num-(sign_denom*denom))
+    # print(lamb)
     # lamb = np.exp(num-denom)
 
     # print(sign_num)
@@ -53,6 +55,7 @@ def custom_mancova(Y,X,C=None):
     #     beta = np.dot(Y.T,X)/X.shape[0]
     # print(Y.mean())
     # print(X.mean())
+    print("!!! DEPRECIATED !!! ")
     Y, beta = linear_regression_test(Y, X, C)
     dot_Y = np.dot(Y.T,Y)
 
@@ -74,7 +77,10 @@ def linear_regression(Y, X, C=None):
         beta = np.dot(np.linalg.inv(np.dot(X_C.T,X_C)), np.dot(X_C.T, Y))
         beta = beta[-1,:].reshape(-1,1) #-1
     else :
-        beta = np.dot(Y.T,X)/X.shape[0]
+        # beta = np.dot(Y.T,X)/X.shape[0]
+        beta = np.dot(np.linalg.inv(np.dot(X.T,X)) , np.dot(X.T,Y))
+        beta = beta.T
+        # beta = np.dot(X.T,Y)/X.shape[0]
     return beta
 
 
diff --git a/python/src/tools/preprocessing_tools.py b/python/src/tools/preprocessing_tools.py
index 9e77695821e3d179f661ad2ae6f5924e016c9cc5..06a539702f9f99f32fd718610106151c477b7f5e 100644
--- a/python/src/tools/preprocessing_tools.py
+++ b/python/src/tools/preprocessing_tools.py
@@ -103,8 +103,10 @@ def adjust_covariates(data_, covariates_):
     covariates = covariates_.copy()
 
     nan_mask = ~np.isnan(data) #[True]* data.shape[0]#
-    cov = scale(covariates[nan_mask])
+    cov = scale(covariates[nan_mask,:])
     dat = scale(data[nan_mask])
+    # cov = covariates[nan_mask,:]
+    # dat = data[nan_mask]
 
     lr = LinearRegression()
     lr.fit(cov, dat)
diff --git a/python/src/univariate.py b/python/src/univariate.py
index 84f9904e7b73b62aaedb5f95be379058c1fc106c..ef63592b98967dfad483868e7bbb00ad5e7c4ffb 100644
--- a/python/src/univariate.py
+++ b/python/src/univariate.py
@@ -24,7 +24,7 @@ class UNIVARIATE :
         self.outputs, self.cols_outputs = pt._extract_cols(self.outputs, self.cols_outputs)
         if len(L_preproc)>0:
             self.outputs = pt.pipeline(self.outputs, L_pipe = L_preproc)
-        self.outputs = scale(self.outputs)
+        # self.outputs = scale(self.outputs)
 
         self.predictors = predictors
         self.cols_predictors = cols_predictors