Adding results for logistic regression

python/logistic/data_utils.py

@@ -0,0 +1,58 @@
+import numpy as np
+
+"""
+  Extract data fom dictionaries for specific organ
+
+  Input: 
+    organ: string in ["lung", "bowel", "thyroid"]
+    perc: int, percentiles [1, ... , 100]
+    suv_dict: dict containing percentiles of suv
+    flags_dict: dict containing states
+"""
+
+def get_data(organ, perc, suv_dict, flags_dict, nr_patient = 58):
+    fname = "get_data"
+
+    # get concrete data set
+    suv = suv_dict[organ + '_SUVperc_COMBINED'][:nr_patient,:,:]  # suv percentiles
+    
+    # index of percentile
+    perc_idx = perc -1
+
+    # computing max SUV percentile per patient, ignoring nans
+    x = np.nanmax(suv[:,:,perc_idx], axis = 1)
+
+    # determining index in the flags based on organ
+    match organ:
+        case "lung":
+            flags_idx = 3
+        case "bowel":
+            flags_idx = 1
+        case "thyroid":
+            flags_idx = 5
+        case _:
+            assert False, f"{fname}::this organ {organ = } is not supported"
+
+    # getting state of patients: 0 == NC, 1  == AE
+    y = flags_dict['flags'][:nr_patient, flags_idx]  
+
+    return x, y
+
+
+"""
+    Check if a vector lies within the specified bounds for each dimension.
+
+    Parameters:
+    - vector (np.ndarray): 1D array representing the point to check. Shape: (n,)
+    - bounds (np.ndarray): 2D array of shape (n, 2), where each row is (min, max) for a dimension.
+
+    Returns:
+    - bool: True if the vector is within bounds in all dimensions, False otherwise.
+"""
+def within_bounds(vector: np.ndarray, bounds: np.ndarray) -> bool:
+
+    if vector.shape[-1] != bounds.shape[0]:
+        raise ValueError("Dimension mismatch: vector length and bounds rows must be equal.")
+    
+    return np.apply_along_axis(lambda x: np.all((x >= bounds[:, 0]) & (x <= bounds[:, 1])), -1, vector)
+

python/logistic/logit_utils.py

@@ -0,0 +1,245 @@
+import numpy as np
+import scipy
+import scipy.stats
+
+
+"""
+    Model function 
+
+        p(x|b) = 1/(1 + exp(-F(x|b)))
+    
+    with log odds of polynomial form:
+
+        log(p(x|b)/(1 - p(x|b))) = F(x|b)
+
+        F(x|b) = sum_{i=0}^degree b_i x^i
+                
+    with decision function (aka logit) F and parameters 
+        
+        b = [b_i]_{i=0}^degree
+
+    Input:
+        x: scalar value or a array of values
+        b = [b_i]_{i=0}^degree: array of r = degree + 1 floats
+           model parameters array of floats
+
+    Return:
+        model function values 
+
+"""
+def logit_poly_model(x, b):
+
+    X = np.column_stack([x**i for i in range(len(b))])
+    F = X @ b
+
+    return 1/(1 + np.exp(-F))
+
+
+"""
+    Performing logistic regression with log odds of polynomial form:
+
+        log(p(x|b)/(1 - p(x|b))) = F(x|b)
+    
+        F(x|b) = sum_{i=0}^degree b_i x^i
+                
+    with decision function (aka logit) F and b = [b_i]_{i=0}^degree. 
+    This gives 
+
+        p(x|b) = 1/(1 + exp(-F(x|b)))
+    
+    Input:
+        lm: instance linear_model.LogisticRegression
+        x: array of n floats
+        y: array of n floats
+        d: degree of decision function
+    
+    Return:
+        params = [b_0, ..., b_degree], array of r = degree + 1 floats
+"""
+
+def logit_poly_fit(lm, x, y, degree = 1):
+   
+   X_feature = np.column_stack([x**i for i in range(1, degree+1)])
+
+   lm.fit(X_feature, y)
+
+   return np.r_[lm.intercept_[0], lm.coef_[0,:]]
+
+"""
+    Producing goodness of fit measures:
+    
+        LLF = log_likelihood function
+        AIC = Akaike information criterion
+        BIC = Bayesian information criterion
+
+    Input:
+        x: array of n floats
+        y: array of n int in {0,1}
+        b: array of r = degree+1 floats, model parameters
+    
+    Return:
+        {"n": n, "k":k, "dof":n-k, "LLF": log_likelihood, "AIC": AIC, "BIC": BIC}
+    
+    Ref:
+        https://en.wikipedia.org/wiki/Logistic_regression
+        https://en.wikipedia.org/wiki/Akaike_information_criterion
+        https://www.medicine.mcgill.ca/epidemiology/joseph/courses/epib-621/logfit.pdf
+"""
+def logit_poly_goodness_of_fit(x, y, b):
+    
+    # model probabilities
+    p = logit_poly_model(x, b)
+
+    # log likelihood
+    eps = 1e-20  # prevent log(0)
+    llf = np.sum(y*np.log(p + eps) + (1 - y)*np.log(1 - p + eps))
+    
+    # information criteria
+    k, n = len(b), len(x)
+    AIC = 2*k - 2*llf
+    BIC = k*np.log(n) - 2*llf
+
+    # chi2
+    dof = n - k
+    r = (y - p)/np.sqrt(p*(1-p))
+    chi2 = np.sum(r**2)
+    p_val = scipy.stats.chi2.sf(chi2, dof)
+
+    return {"LLF": llf, "AIC": AIC, "BIC": BIC, 
+            "chi2": chi2, "p-value(chi2)": p_val,
+            "n": n, "k": k, "dof": dof}
+
+"""
+    Calculation of asymptotic variance-covariance matrix for the 
+    logistic regression of the polynomial model:
+    
+        log(p(x)/(1 - p(x))) ~ sum_{i=0}^degree b_i x^i
+
+    Input:
+        x: array of n floats
+        b: array of r = degree+1 floats, model parameters
+        
+    Return:
+        array of rxr floats; r = degree + 1
+    
+    Ref:
+        https://stats.stackexchange.com/questions/89484/how-to-compute-the-standard-errors-of-a-logistic-regressions-coefficients
+"""
+def logit_poly_cov(x, b):
+    
+    # Calculate matrix of predicted class probabilities.
+    probs = logit_poly_model(x, b)
+    
+    # Design matrix -- add column of 1's at the beginning of your X_train matrix
+    X = np.column_stack([x**i for i in range(len(b))])
+
+    # Initiate matrix of 0's, fill diagonal with each predicted observation's variance
+    V = np.diagflat(probs*(1 - probs)) # dig.matrix where each element is p*(1-p)
+
+    # Covariance matrix C_params = (X^T V X)^-1
+    return np.linalg.inv(X.T@V@X)
+
+
+"""
+    Calculating quantiles of the model parameters at given probabilities p 
+    for normal distribution of parameters:
+
+        b ~ N(mean_b, cov)
+    
+    Input:
+        probs: array of m floats, probabilities
+        mean_b: array of r = degree+1 floats, mean model parameters
+        cov: array of rxr floats, variance-covariance matrix of parameters
+    
+    Return:
+        array of mxn floats
+"""
+
+def logit_poly_pars_quantiles_normal(probs, mean_b, cov):
+    
+    # mean and standard variance parameters
+    locs = mean_b
+    scales = np.sqrt(np.diag(cov))
+
+    # computing quantiles of parameters
+    return locs + np.outer(scipy.stats.norm.ppf(probs), scales)
+
+"""
+    Calculating quantiles of the model values 
+    
+        p(x|b) = 1/(1 + exp(-F(x|b)))
+    
+    with
+
+        F(x|b) = sum_{i=0}^degree x^i b_i
+    
+    at given probabilities p and values x assuming 
+    normal distribution of parameters:
+
+        b ~ N(mean_b, cov)
+    
+    Input:
+        x: array of n float
+        probs: array of m floats, probabilities
+        mean_b: array of r = degree+1 floats, mean model parameters
+        cov: array of rxr floats, variance-covariance matrix of parameters
+    
+    Return:
+        array of mxn floats
+"""
+
+def logit_poly_model_quantiles_normal(x, probs, mean_b, cov):
+    
+    X = np.column_stack([x**i for i in range(len(mean_b))])
+
+    # mean and standard variance of logit (aka log of odds)
+    locs = X@mean_b
+    scales = np.sqrt(np.diag(X@cov@X.T))
+
+    # computing quantiles of logit
+    Q = locs + np.outer(scipy.stats.norm.ppf(probs),scales)
+    
+    # convert logit to expit
+    return 1/(1 + np.exp(-Q))
+
+"""
+    Calculating quantiles using delta method of the model values 
+    
+        p(x|b) = 1/(1 + exp(-F(x|b)))
+
+    with
+
+        F(x|b) = sum_{i=0}^degree x^i b_i
+    
+    at given probabilities p and values x assuming 
+    normal distribution of parameters:
+
+        b ~ N(mean_b, cov)
+
+    We approximate exact model with linear expansion
+
+        p(x|b) = p(x|b_mean) + dp/db (x| b_mean) (b - b_mean)
+    
+    and the last term is normally distributed.
+
+    Input:
+        x: array of n float
+        probs: array of m floats, probabilities
+        mean_b: array of r = degree+1 floats, mean model parameters
+        cov: array of rxr floats, variance-covariance matrix of parameters 
+    
+    Return:
+        array of mxn floats
+"""
+def logit_poly_model_quantiles_delta(x, probs, mean_b, cov):
+
+    X = np.column_stack([x**i for i in range(len(mean_b))])
+    F = X@mean_b
+
+    locs = 1/(1 + np.exp(-F))
+    scales = np.sqrt(np.diag(X@cov@X.T))/(4*np.cosh(F/2)**2)
+
+    # computing quantiles of logit
+    Q = locs + np.outer(scipy.stats.norm.ppf(probs), scales)
+
+    return np.clip(Q, a_min = 0, a_max = 1)

python/logistic/mvn.py

@@ -0,0 +1,173 @@
+import numpy as np
+import scipy
+
+from typing import Tuple
+
+
+# * * * STATISTICAL TESTS * * *
+
+def mardia_test(data: np.ndarray, cov: bool = True) -> Tuple[float, float, float, float]:
+    """
+    https://rdrr.io/cran/MVN/src/R/mvn.R
+    https://stats.stackexchange.com/questions/317147/how-to-get-a-single-p-value-from-the-two-p-values-of-a-mardias-multinormality-t
+    Mardia's multivariate skewness and kurtosis.
+    Calculates the Mardia's multivariate skewness and kurtosis coefficients
+    as well as their corresponding statistical test. For large sample size
+    the multivariate skewness is asymptotically distributed as a Chi-square
+    random variable; here it is corrected for small sample size. However,
+    both uncorrected and corrected skewness statistic are presented. Likewise,
+    the multivariate kurtosis it is distributed as a unit-normal.
+
+     Syntax: function [Mskekur] = Mskekur(X,c,alpha)
+
+     Inputs:
+          X - multivariate data matrix [Size of matrix must be n(data)-by-p(variables)].
+          cov - boolean to whether to normalize the covariance matrix by n (c=1[default]) or by n-1 (c~=1)
+
+     Outputs:
+          - skewness test statistic
+          - kurtosis test statistic
+          - significance value for skewness
+          - significance value for kurtosis
+    """
+    n, p = data.shape
+
+    # correct for small sample size
+    small: bool = True if n < 20 else False
+
+    if cov:
+        S = ((n - 1)/n) * np.cov(data.T)
+    else:
+        S = np.cov(data.T)
+
+    # calculate mean
+    data_mean = data.mean(axis=0)
+    # inverse - check if singular matrix
+    try:
+        iS = np.linalg.inv(S)
+    except Exception as e:
+        # print for now
+        print(e)
+        return 0.0, 0.0, 0.0, 0.0
+    # squared-Mahalanobis' distances matrix
+    D: np.ndarray = (data - data_mean) @ iS @ (data - data_mean).T
+    # multivariate skewness coefficient
+    g1p: float = np.sum(D**3)/n**2
+    # multivariate kurtosis coefficient
+    g2p: float = np.trace(D**2)/n
+    # small sample correction
+    k: float = ((p + 1)*(n + 1)*(n + 3))/(n*(((n + 1)*(p + 1)) - 6))
+    # degrees of freedom
+    df: float = (p * (p + 1) * (p + 2))/6
+
+    if small:
+        # skewness test statistic corrected for small sample: it approximates to a chi-square distribution
+        g_skew = (n * g1p * k)/6
+    else:
+        # skewness test statistic:it approximates to a chi-square distribution
+        g_skew = (n * g1p)/6
+
+    # significance value associated to the skewness corrected for small sample
+    p_skew: float = 1.0 - scipy.stats.chi2.cdf(g_skew, df)
+
+    # kurtosis test statistic: it approximates to a unit-normal distribution
+    g_kurt = (g2p - (p*(p + 2)))/(np.sqrt((8 * p * (p + 2))/n))
+    # significance value associated to the kurtosis
+    p_kurt: float = 2 * (1.0 - scipy.stats.norm.cdf(np.abs(g_kurt)))
+
+    return g_skew, g_kurt, p_skew, p_kurt
+
+
+def hz_test(data: np.ndarray, cov: bool = True) -> Tuple[float, float]:
+    """
+    Henze-Zirkler method for goodness of fit of data to a multivariate normal distribution.
+    Researchers tend to use this MVN test for larger samples (N > 100).
+    https://www.tandfonline.com/doi/abs/10.1080/03610929008830400
+
+    :param data: multivariate data matrix [Size of matrix must be n(data)-by-p(variables)].
+    :param cov: boolean to whether to normalize the covariance matrix by n (c=1[default]) or by n-1 (c~=1)
+    :return:
+        HZ - Henze-Zirkler test statistic
+        p_value - significance value
+    """
+    n, p = data.shape
+
+    if cov:
+        S = ((n - 1)/n) * np.cov(data.T)
+    else:
+        S = np.cov(data.T)
+
+    # calculate mean
+    data_mean = data.mean(axis=0)
+
+    try:
+        iS = np.linalg.inv(S)
+    except Exception as e:
+        print(e)
+        return 0.0, 0.0
+
+    Y = data @ iS @ data.T
+    Dj = np.diag((data - data_mean) @ iS @ (data - data_mean).T)
+
+    Djk = - 2 * Y.T + np.tensordot(np.diag(Y.T), np.ones(n), axes=0) + np.tensordot(np.ones(n), np.diag(Y.T), axes=0)
+    b: float = 1 / (np.sqrt(2)) * ((2 * p + 1) / 4) ** (1 / (p + 4)) * (n ** (1 / (p + 4)))
+
+    # calculate rank of matrix
+    S_rank = np.linalg.matrix_rank(S)
+
+    if S_rank == p:
+        HZ = n * (1 / (n ** 2) * np.sum(np.sum(np.exp(- (b ** 2) / 2 * Djk))) - 2 * ((1 + (b ** 2)) ** (- p / 2)) * (1 / n) * (np.sum(np.exp(- ((b ** 2) / (2 * (1 + (b ** 2)))) * Dj))) + ((1 + (2 * (b ** 2))) ** (- p / 2)))
+    else:
+        HZ = n * 4
+
+    wb = (1 + b ** 2) * (1 + 3 * b ** 2)
+    a = 1 + 2 * b ** 2
+
+    # HZ mean
+    mu = 1 - a ** (- p / 2) * (1 + p * b ** 2 / a + (p * (p + 2) * (b ** 4)) / (2 * a ** 2))  # HZ mean
+
+    # HZ variance
+    si2 = 2 * (1 + 4 * b ** 2) ** (- p / 2) + 2 * a ** (- p) * (1 + (2 * p * b ** 4) / a ** 2 + (3 * p * (p + 2) * b ** 8) / (4 * a ** 4)) - 4 * wb ** (- p / 2) * (1 + (3 * p * b ** 4) / (2 * wb) + (p * (p + 2) * b ** 8) / (2 * wb ** 2))
+
+    pmu = np.log(np.sqrt(mu ** 4 / (si2 + mu ** 2)))  # lognormal HZ mean
+    psi = np.sqrt(np.log((si2 + mu ** 2) / mu ** 2))  # lognormal HZ standard deviation
+
+    # calculate p-value
+    p_value = 1.0 - scipy.stats.lognorm.cdf(HZ, psi, scale=np.exp(pmu))
+
+    return HZ, p_value
+
+
+import numpy as np
+from scipy.stats import shapiro, chi2
+
+def royston_test(X):
+    """
+    Royston's Multivariate Normality Test using Fisher's method on Shapiro-Wilk p-values.
+    
+    Parameters:
+        X (ndarray): 2D array (n_samples x n_variables)
+
+    Returns:
+        stat (float): Fisher's combined test statistic
+        p_value (float): p-value for overall multivariate normality
+    """
+    X = np.asarray(X)
+    n, p = X.shape
+
+    if n < 3:
+        raise ValueError("At least 3 observations are required.")
+    if p < 2:
+        raise ValueError("At least 2 variables required.")
+
+    p_values = []
+    for i in range(p):
+        _, pval = shapiro(X[:, i])
+        p_values.append(pval)
+
+    p_values = np.clip(p_values, 1e-16, 1.0)  # avoid log(0)
+    stat = -2 * np.sum(np.log(p_values))
+    df = 2 * p
+    p_combined = 1 - chi2.cdf(stat, df)
+    
+    return stat, p_combined