Skip to content
Snippets Groups Projects
Commit 46b4049c authored by Christophe BOETTO's avatar Christophe BOETTO
Browse files

updated scripts to return fisher value

parent 0ab6e020
Branches master
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
python/src/manocca.py
+ 25
− 11
View file @ 46b4049c
...@@ -110,7 +110,7 @@ class MANOCCA: ...@@ -110,7 +110,7 @@ class MANOCCA:
def __init__(self, predictors, outputs, covariates=None, cols_outputs = None, def __init__(self, predictors, outputs, covariates=None, cols_outputs = None,
cols_predictors = None, cols_covariates = None, prodV_red=None, cols_predictors = None, cols_covariates = None, prodV_red=None,
n_comp = None, prod_to_keep = None, use_resid = True, use_pca = True, n_comp = None, prod_to_keep = None, use_resid = True, use_pca = True,
adj_pred = False, use_extended = False, corr_bias = False, n_jobs = 1, verbose = 1): adj_pred = False, use_extended = False, corr_bias = False, n_jobs = 1, apply_qt=True, return_fisher = False, verbose = 1):
### Initializing ### Initializing
self.outputs = outputs self.outputs = outputs
...@@ -143,15 +143,21 @@ class MANOCCA: ...@@ -143,15 +143,21 @@ class MANOCCA:
self.use_pca = use_pca self.use_pca = use_pca
self.use_extended = use_extended self.use_extended = use_extended
self.adj_pred = adj_pred self.adj_pred = adj_pred
self.apply_qt = apply_qt
self.return_fisher = return_fisher
self.bias = None self.bias = None
self.corr_bias = corr_bias self.corr_bias = corr_bias
if corr_bias : # if we want to use the correction, we can't ust the qt transformation
self.apply_qt = False
# Filled later # Filled later
self.prodV = None self.prodV = None
self.prodV_red = None self.prodV_red = None
self.pca = None self.pca = None
self.p = None self.p = None
self.fisher = None
# Filtering product columns # Filtering product columns
self.prod_to_keep = prod_to_keep self.prod_to_keep = prod_to_keep
...@@ -184,9 +190,9 @@ class MANOCCA: ...@@ -184,9 +190,9 @@ class MANOCCA:
print(self.prodV.shape) print(self.prodV.shape)
else : else :
if self.corr_bias: if self.corr_bias:
self.prodV = self.get_prodV_para_wrap(self.outputs, apply_qt=False) self.prodV = self.get_prodV_para_wrap(self.outputs, apply_qt=False) #Always false for corr_bias as qt can bring back the bias
else: else:
self.prodV = self.get_prodV_para_wrap(self.outputs, apply_qt=True) self.prodV = self.get_prodV_para_wrap(self.outputs, apply_qt=self.apply_qt)
if self.verbose>0 : if self.verbose>0 :
print(self.prodV.shape) print(self.prodV.shape)
if not isinstance(self.prod_to_keep, type(None)): # Filtering out some columns if not isinstance(self.prod_to_keep, type(None)): # Filtering out some columns
...@@ -233,12 +239,12 @@ class MANOCCA: ...@@ -233,12 +239,12 @@ class MANOCCA:
pred = self.predictors[:,i_pred] pred = self.predictors[:,i_pred]
beta = np.apply_along_axis(lambda x : linear_regression(Y_tmp, x.reshape(-1,1)).flatten(), axis = 0, arr = pred) beta = np.apply_along_axis(lambda x : linear_regression(Y_tmp, x.reshape(-1,1)).flatten(), axis = 0, arr = pred)
beta = beta.reshape(-1,1) beta = beta.reshape(-1,1)
prod_beta = self.get_prodV_para_wrap(beta.T, apply_qt = False, job_overide = True) prod_beta = self.get_prodV_para_wrap(beta.T, apply_qt = self.apply_qt, job_overide = True)
self.prod_beta = prod_beta self.prod_beta = prod_beta
self.bias = prod_beta*(pred*pred).reshape(-1,1) self.bias = prod_beta*(pred*pred).reshape(-1,1)
def get_prodV_para_effi(self, DD0, job_id, nb_compute, apply_qt = True): def get_prodV_para_effi(self, DD0, job_id, nb_compute, apply_qt):
""" Computes the product matrix for a set of outputs and a given number of jobs """ Computes the product matrix for a set of outputs and a given number of jobs
...@@ -264,14 +270,13 @@ class MANOCCA: ...@@ -264,14 +270,13 @@ class MANOCCA:
tmp = (DD0[:,(i+1):]*DD0[:,i].reshape(-1,1)) tmp = (DD0[:,(i+1):]*DD0[:,i].reshape(-1,1))
# std = tmp.std(axis = 0) # std = tmp.std(axis = 0)
if apply_qt: if apply_qt:
# print('coucou')
tmp = pt.get_qt(tmp) tmp = pt.get_qt(tmp)
# tmp = (std) * tmp # tmp = (std) * tmp
L_prodV += [tmp] L_prodV += [tmp]
return L_prodV return L_prodV
def get_prodV_para_wrap(self, DD0, apply_qt = True, job_overide = False ) : def get_prodV_para_wrap(self, DD0, apply_qt, job_overide = False ) :
""" Computes the product matrix for a set of outputs """ Computes the product matrix for a set of outputs
This function wraps the parallel processing to compute the product matrix. It parameters the get_prodV_para_effi function depending on the parralel status. This function wraps the parallel processing to compute the product matrix. It parameters the get_prodV_para_effi function depending on the parralel status.
...@@ -302,7 +307,7 @@ class MANOCCA: ...@@ -302,7 +307,7 @@ class MANOCCA:
if self.verbose>0 : if self.verbose>0 :
print("Parallel computation of prodV") print("Parallel computation of prodV")
nb_compute = cpu_count() nb_compute = cpu_count()
res = Parallel(n_jobs=n_jobs, verbose = self.verbose)(delayed(self.get_prodV_para_effi)(DD0,j, nb_compute, apply_qt = apply_qt) for j in range(nb_compute)) res = Parallel(n_jobs=n_jobs, verbose = self.verbose)(delayed(self.get_prodV_para_effi)(DD0,j, nb_compute, apply_qt = self.apply_qt) for j in range(nb_compute))
# reordering # reordering
res_ordered = [] res_ordered = []
for _ in range(len(res[0])): for _ in range(len(res[0])):
...@@ -453,6 +458,10 @@ class MANOCCA: ...@@ -453,6 +458,10 @@ class MANOCCA:
n_preds = var.shape[1] n_preds = var.shape[1]
p = np.empty((n_preds,1)) p = np.empty((n_preds,1))
if self.return_fisher:
fisher = np.empty((n_preds,1))
for i_pred in range(n_preds): for i_pred in range(n_preds):
nan_mask = np.isnan(var[:,i_pred]) nan_mask = np.isnan(var[:,i_pred])
...@@ -467,14 +476,19 @@ class MANOCCA: ...@@ -467,14 +476,19 @@ class MANOCCA:
if ((self.covariates is not None) and (self.use_resid == False)) : if ((self.covariates is not None) and (self.use_resid == False)) :
# print("With covariates") # print("With covariates")
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,:])) # need to dissociate case with and without bias #[i_pred,:].reshape(-1,1)) if self.return_fisher :
p[i_pred], fisher[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,:]), return_fisher=True)
else :
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 : else :
# print("Without covariates") # print("Without covariates")
if not isinstance(self.bias, type(None)): if self.return_fisher :
p[i_pred] = custom_manova(self.prodV_red[~nan_mask,min_comp:n_comp], var[~nan_mask,i_pred].reshape(-1,1))#[i_pred,:].reshape(-1,1)) p[i_pred], fisher[i_pred] = custom_manova(self.prodV_red[~nan_mask,min_comp:n_comp], var[~nan_mask,i_pred].reshape(-1,1), return_fisher=True)
else : else :
p[i_pred] = custom_manova(self.prodV_red[~nan_mask,min_comp:n_comp], 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 self.p = p
if self.return_fisher:
self.fisher = fisher
......
python/src/manova.py
+ 18
− 4
View file @ 46b4049c
...@@ -11,7 +11,7 @@ from tools.compute_manova import custom_manova ...@@ -11,7 +11,7 @@ from tools.compute_manova import custom_manova
class MANOVA : class MANOVA :
def __init__(self, predictors, outputs, covariates=None, cols_outputs = None, def __init__(self, predictors, outputs, covariates=None, cols_outputs = None,
cols_predictors = None, cols_covariates = None, use_resid = True): cols_predictors = None, cols_covariates = None, use_resid = True, return_fisher=False):
### Initializing ### Initializing
self.outputs = outputs self.outputs = outputs
self.cols_outputs = cols_outputs self.cols_outputs = cols_outputs
...@@ -21,7 +21,6 @@ class MANOVA : ...@@ -21,7 +21,6 @@ class MANOVA :
self.predictors = predictors self.predictors = predictors
self.cols_predictors = cols_predictors self.cols_predictors = cols_predictors
self.predictors, self.cols_predictors = pt._extract_cols(self.predictors, self.cols_predictors) self.predictors, self.cols_predictors = pt._extract_cols(self.predictors, self.cols_predictors)
self.predictors = scale(self.predictors) self.predictors = scale(self.predictors)
...@@ -38,7 +37,11 @@ class MANOVA : ...@@ -38,7 +37,11 @@ class MANOVA :
# print("adjusting predictor") # print("adjusting predictor")
# self.predictors = np.apply_along_axis(lambda x : pt.adjust_covariates(x,covariates), axis = 0, arr = self.predictors) # self.predictors = np.apply_along_axis(lambda x : pt.adjust_covariates(x,covariates), axis = 0, arr = self.predictors)
self.return_fisher = return_fisher
# For results
self.p = None
self.fisher = None
...@@ -50,15 +53,26 @@ class MANOVA : ...@@ -50,15 +53,26 @@ class MANOVA :
nb_tests = var.shape[1] nb_tests = var.shape[1]
p = np.empty((nb_tests,1)) p = np.empty((nb_tests,1))
if self.return_fisher:
fisher = np.empty((nb_tests,1))
for i_pred in range(nb_tests): for i_pred in range(nb_tests):
nan_mask = np.isnan(var[:,i_pred]) nan_mask = np.isnan(var[:,i_pred])
if ((self.covariates is not None) and (self.use_resid == False)) : if ((self.covariates is not None) and (self.use_resid == False)) :
# print("With covariates") # print("With covariates")
p[i_pred] = custom_manova(self.outputs[~nan_mask,:], var[~nan_mask,i_pred].reshape(-1,1), self.covariates[~nan_mask,:]) if self.return_fisher :
p[i_pred], fisher[i_pred] = custom_manova(self.outputs[~nan_mask,:], var[~nan_mask,i_pred].reshape(-1,1), self.covariates[~nan_mask,:], return_fisher=self.return_fisher)
else :
p[i_pred] = custom_manova(self.outputs[~nan_mask,:], var[~nan_mask,i_pred].reshape(-1,1), self.covariates[~nan_mask,:])
else : else :
# print("Without covariates") # print("Without covariates")
p[i_pred] = custom_manova(self.outputs[~nan_mask,:], var[~nan_mask,i_pred].reshape(-1,1)) if self.return_fisher :
p[i_pred], fisher[i_pred] = custom_manova(self.outputs[~nan_mask,:], var[~nan_mask,i_pred].reshape(-1,1), return_fisher=self.return_fisher)
else :
p[i_pred] = custom_manova(self.outputs[~nan_mask,:], var[~nan_mask,i_pred].reshape(-1,1))
self.p = p self.p = p
if self.return_fisher:
self.fisher = fisher
......
python/src/tools/compute_manova.py
+ 12
− 7
View file @ 46b4049c
...@@ -6,7 +6,7 @@ from scipy.stats import f ...@@ -6,7 +6,7 @@ from scipy.stats import f
# from sklearn.preprocessing import scale # from sklearn.preprocessing import scale
from .preprocessing_tools import scale from .preprocessing_tools import scale
def custom_manova(Y,X,C=None, return_beta = False, bias = None): def custom_manova(Y,X,C=None, return_beta = False, bias = None, return_fisher=False):
### /!\ data needs to be scaled, X needs to be (k,1) ### /!\ data needs to be scaled, X needs to be (k,1)
beta = linear_regression(Y, X, C) beta = linear_regression(Y, X, C)
dot_Y = np.dot(Y.T,Y) dot_Y = np.dot(Y.T,Y)
...@@ -18,10 +18,12 @@ def custom_manova(Y,X,C=None, return_beta = False, bias = None): ...@@ -18,10 +18,12 @@ def custom_manova(Y,X,C=None, return_beta = False, bias = None):
fisher = (Y.shape[0]-Y.shape[1]-1)/(Y.shape[1])*(1-lamb)/lamb fisher = (Y.shape[0]-Y.shape[1]-1)/(Y.shape[1])*(1-lamb)/lamb
if return_beta == False : if return_beta == True :
return f.sf(fisher, Y.shape[1], Y.shape[0]-Y.shape[1]-1)
else :
return f.sf(fisher, Y.shape[1], Y.shape[0]-Y.shape[1]-1), beta return f.sf(fisher, Y.shape[1], Y.shape[0]-Y.shape[1]-1), beta
elif return_fisher == True :
return f.sf(fisher, Y.shape[1], Y.shape[0]-Y.shape[1]-1), fisher
else :
return f.sf(fisher, Y.shape[1], Y.shape[0]-Y.shape[1]-1)
def custom_mancova(Y,X,C=None): def custom_mancova(Y,X,C=None):
...@@ -39,16 +41,19 @@ def custom_mancova(Y,X,C=None): ...@@ -39,16 +41,19 @@ def custom_mancova(Y,X,C=None):
def linear_regression(Y, X, C=None): def linear_regression(Y, X, C=None, return_all = False):
ones = np.ones((X.shape[0],1)) ones = np.ones((X.shape[0],1))
if C is not None: if C is not None:
X_C = np.hstack([ones,C,X]) X_C = np.hstack([ones,C,X])
beta = np.dot(np.linalg.inv(np.dot(X_C.T,X_C)), np.dot(X_C.T, Y)) 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 if return_all == False:
beta = beta[-1,:].reshape(-1,1) #-1
else : else :
X = np.hstack([ones,X]) #,ones]) #[C,X] X = np.hstack([ones,X]) #,ones]) #[C,X]
beta = np.dot(np.linalg.inv(np.dot(X.T,X)) , np.dot(X.T,Y)) beta = np.dot(np.linalg.inv(np.dot(X.T,X)) , np.dot(X.T,Y))
beta = beta[-1,:].reshape(-1,1) if return_all == False:
beta = beta[-1,:].reshape(-1,1)
return beta return beta
......
python/src/univariate.py
+ 4
− 5
View file @ 46b4049c
...@@ -31,8 +31,6 @@ class UNIVARIATE : ...@@ -31,8 +31,6 @@ class UNIVARIATE :
self.predictors, self.cols_predictors = pt._extract_cols(self.predictors, self.cols_predictors) self.predictors, self.cols_predictors = pt._extract_cols(self.predictors, self.cols_predictors)
self.predictors = scale(self.predictors) self.predictors = scale(self.predictors)
self.covariates = covariates self.covariates = covariates
self.cols_covariates = cols_covariates self.cols_covariates = cols_covariates
if covariates is not None : if covariates is not None :
...@@ -52,14 +50,15 @@ class UNIVARIATE : ...@@ -52,14 +50,15 @@ class UNIVARIATE :
self.p = [] self.p = []
for i in range(self.predictors.shape[1]): for i in range(self.predictors.shape[1]):
pred = self.predictors[:,i].reshape(-1,1) pred = self.predictors[:,i].reshape(-1,1)
# print(pred.shape)
nan_mask = np.isnan(pred).flatten() nan_mask = np.isnan(pred).flatten()
if ((self.covariates is not None) and (self.use_resid == False)) : if ((self.covariates is not None) and (self.use_resid == False)) :
print(self.outputs.shape) print(self.outputs.shape)
beta = linear_regression(self.outputs[~nan_mask], pred[~nan_mask], self.covariates[~nan_mask]) beta = linear_regression(self.outputs[~nan_mask,:], pred[~nan_mask,:], self.covariates[~nan_mask])
else : else :
beta = linear_regression(self.outputs[~nan_mask], pred[~nan_mask]) beta = linear_regression(self.outputs[~nan_mask,:], pred[~nan_mask,:])
self.beta = beta self.beta = beta
res = ttest(beta, pred[~nan_mask], self.outputs[~nan_mask]) res = ttest(beta, pred[~nan_mask,:], self.outputs[~nan_mask,:])
self.p += [res] self.p += [res]
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment