compute_manova.py

import numpy as np

from scipy.stats import f

# to remove, only used for tests
# from sklearn.preprocessing import scale
from .preprocessing_tools import scale

def custom_manova(Y,X,C=None, return_beta = False):
    ### /!\ data needs to be scaled, X needs to be (k,1)
    # if C is not None:
    #     X_C = np.hstack([X,C]) #[C,X]
    #     # beta = np.dot(np.dot(Y.T,X_C),np.linalg.inv(np.dot(X_C.T,X_C)))
    #     beta = np.dot(np.linalg.inv(np.dot(X_C.T,X_C)), np.dot(X_C.T, Y))
    #     beta = beta[0,:].reshape(-1,1) #-1
    # else :
    #     beta = np.dot(Y.T,X)/X.shape[0]
    beta = linear_regression(Y, X, C)
    dot_Y = np.dot(Y.T,Y)

    # print(Y)

    # (sign_num, num) = np.linalg.slogdet(dot_Y-X.shape[0]*np.dot(beta,beta.T))
    # print('ou')
    (sign_num, num) = np.linalg.slogdet(dot_Y - np.dot( np.dot(beta, X.T) , np.dot(X, beta.T)))
    # (sign_num, num) = np.linalg.slogdet((Y - X @ beta.T).T@(Y - 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)
    # print(sign_denom)
    # print(num)
    # print(denom)
    # print(lamb)
    # print(num-denom)

    fisher = (Y.shape[0]-Y.shape[1]-1)/(Y.shape[1])*(1-lamb)/lamb

    # print(fisher)
    # print(Y.shape[1])
    # print(Y.shape[0])
    # print(f.sf(fisher, Y.shape[1], Y.shape[0]-Y.shape[1]-1))
    if return_beta == False :
        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


def custom_mancova(Y,X,C=None):
    ### /!\ data needs to be scaled, X needs to be (k,1)
    # if C is not None:
    #     X_C = np.hstack([X,C]) #[C,X]
    #     # beta = np.dot(np.dot(Y.T,X_C),np.linalg.inv(np.dot(X_C.T,X_C)))
    #     beta = np.dot(np.linalg.inv(np.dot(X_C.T,X_C)), np.dot(X_C.T, Y))
    #     beta = beta[0,:].reshape(-1,1) #-1
    # else :
    #     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)

    (sign_num, num) = np.linalg.slogdet(dot_Y-X.shape[0]*np.dot(beta,beta.T))
    (sign_denom, denom) = np.linalg.slogdet(dot_Y)
    lamb = np.exp(sign_num*num-(sign_denom*denom))

    fisher = (Y.shape[0]-Y.shape[1]-1)/(Y.shape[1])*(1-lamb)/lamb
    return f.sf(fisher, Y.shape[1], Y.shape[0]-Y.shape[1]-1)



def linear_regression(Y, X, C=None):
    ones = np.ones((X.shape[0],1))
    if C is not None:
        # X_C = np.hstack([X,C]) #,ones]) #[C,X]
        X_C = np.hstack([ones,C,X]) #,ones]) #[C,X]
        # beta = np.dot(np.dot(Y.T,X_C),np.linalg.inv(np.dot(X_C.T,X_C)))
        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]
        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 = beta.T
        # print(beta)
        beta = beta[-1,:].reshape(-1,1) 

        # beta = np.dot(X.T,Y)/X.shape[0]
        # print(beta.shape)
        # beta = beta.T
    return beta


def linear_regression_test(Y, X, C=None):
    ones = np.ones((X.shape[0],1))
    if C is not None:
        X_C = np.hstack([ones,C,X]) #,ones]) #[C,X]
        # beta = np.dot(np.dot(Y.T,X_C),np.linalg.inv(np.dot(X_C.T,X_C)))
        beta = np.dot(np.linalg.inv(np.dot(X_C.T,X_C)), np.dot(X_C.T, Y))
        
        Y = Y- np.dot(X_C[:,:-1],beta[:-1,:]) #Y - X_C[:,:-1] @ beta[:-1,:] # the @ option is slightly slower 13.5 us vs 23 us

        beta = beta[-1,:].reshape(-1,1) #-1
    else :
        beta = np.dot(Y.T,X)/X.shape[0]

    return Y, beta 



from scipy.stats import t

def ttest(beta, X, Y):
    L_p = []

    for i in range(beta.shape[0]) :
        se_b = np.sqrt(1/(X.shape[0]-2)*np.sum((Y[:,i].reshape(-1,1)-beta[i]*X)**2)/np.sum((X-X.mean())**2))
        L_p += [2*(t.sf(np.abs(beta[i])/se_b,X.shape[0]-2))[0]]
        # L_p += [1 - 2 * np.abs(0.5 - np.vectorize(t.cdf)(beta[i]/se_b, X.shape[0] - 2))[0]]

    return L_p

### TESTS ###
# Y = scale(np.random.randint(0,3,(100,50)))
# X = scale(np.random.randint(0,3,(100,1)))
# C = scale(np.random.randint(0,3,(100,3)))
# #Y = np.random.randint(0,3,(10,5))
# # print(Y)
# # print(X)
# # print(np.dot(X.T,X))


# X_C = np.hstack([X,C])
# print(np.dot(X_C.T, Y).shape)


# print(custom_manova(Y,X,C))


# import sys
# import os
# # os.listdir('/Documents/python/Micmat/tools')
# sys.path.insert(0, '\\Users\\cboetto\\Documents\\python\\Micmat\\tools')
# import preprocessing_tools as pt
# from manova import custom_manova as cm

# help(cm)
# print(cm(Y,X))

# Y_C = np.apply_along_axis(lambda x : pt.adjust_covariates(x,C), axis = 0, arr = Y)
# print(cm(Y_C,X))







### TRASH ###
# def custom_manova_old(Y,X_,C=None):
#     ### /!\ data needs to be scaled
#     n_var = X_.shape[1]
#     if C :
#         X = np.hstack([X_,C])

#     beta = np.dot(Y.T,X)/X.shape[0]
#     beta = []
#     dot_Y = np.dot(Y.T,Y)

#     (sign_num, num) = np.linalg.slogdet(dot_Y-X.shape[0]*np.dot(beta,beta.T))
#     (sign_denom, denom) = np.linalg.slogdet(dot_Y)
#     lamb = np.exp(sign_num*num-(sign_denom*denom))

#     fisher = (Y.shape[0]-Y.shape[1]-1)/(Y.shape[1])*(1-lamb)/lamb
#     return f.sf(fisher, Y.shape[1], Y.shape[0]-Y.shape[1]-1)