# -*- coding: utf-8 -*-
"""
Created on Thu Jun 6 13:28:15 2024
@author: Polina Arsenteva
"""
import itertools
import warnings
import re
from tqdm import tqdm
import numpy as np
import pandas as pd
from statsmodels.iolib import summary2
from patsy import dmatrices, DesignInfo, ContrastMatrix
from scipy.stats import chi2, f, gaussian_kde, beta
import matplotlib.pyplot as plt
from matplotlib import rc, colormaps
from torch import cdist, exp, matmul, diag, trace, sqrt, pow, cat
from torch import eye, ones, tensor, float64, from_numpy, Tensor, finfo
from torch.linalg import multi_dot, eigh, inv
import sys
sys.setrecursionlimit(10 ** 6)
[docs]
def ordered_eigsy(matrix, eps=None, clip=True):
"""
Calculates the eigendecomposition of the matrix, using a solver
implemented in C++.
Eigen values and eigen vectors are stored by decreasing eigen value order.
Note 1: only positive non-null eigen values and corresponding eigen
vectors are returned. Eigen values lower than `eps` threshold a
clipped to zeros.
Note 2: input matrix is assumed to be symmetrical and positive
semi-definite, such that its eigen values are positive or null. It
may happen due to numerical issues that the eigen decomposition finds
negative eigen values for such matrix anyway. These are also clipped
to 0.
Parameters
----------
matrix : 2-d array_like
Matrix to decompose.
eps : float, optional
minimum threshold value to clip lower eigen values to zeros.
If `None` (default), then machine precision (given by
`torch.finfo()`) for matrix dtype is used as threshold.
clip : boolean,
flag to enable/disable eigen value clipping.
Returns
-------
sp : torch.Tensor
Eigenvalues.
ev : torch.Tensor
Eigenvectors.
"""
# eigen values, eigen vectors
# following ascending eigen value order (no need for sorting)
sp, ev = eigh(matrix)
# revert to get decreasing order for eigen values
sp = sp.flip(dims=(0,))
ev = ev.flip(dims=(1,))
# clipping?
if clip:
# eigen value clipping threshold
if eps is None:
eps = finfo(matrix.dtype).eps
# select non null eigen val
sp_mask = sp >= eps
# reduc dimension
reduc_dim = sp_mask.sum()
# warning
if reduc_dim < len(sp):
warnings.warn(
f"Clipping last {len(sp) - reduc_dim} eigen values " +
f"out of {len(sp)} dimensions, " +
f"that are lower than threshold {eps}."
)
# output
return sp[sp_mask], ev[:, sp_mask]
else:
# no clipping
# output
return sp, ev
[docs]
def convert_to_torch(A):
"""
Converts A to torch.Tensor.
Parameters
----------
A : array_like
Container to convert.
Returns
-------
B : torch.Tensor
A converted to torch.
"""
if isinstance(A, pd.Series):
B = from_numpy(A.to_numpy().reshape(-1, 1)).double()
if isinstance(A, pd.DataFrame):
B = from_numpy(A.to_numpy()).double()
if isinstance(A, Tensor):
B = A.double()
else:
try:
X = A.to_numpy() if not isinstance(A, np.ndarray) else A.copy()
B = from_numpy(X).double()
except AttributeError:
print(f'Unknown data type {type(A)}')
return B
[docs]
def distances(x, y=None):
"""
Computes the distances between each pair of the two collections of row
vectors of x and y, or x and itself if y is not provided.
Parameters
----------
x : torch.Tensor
Input 2-d tensor.
y : None or torch.Tensor
Input 2-d tensor. Replaced by `x` if None.
Returns
-------
sq_dists : torch.Tensor
Distance matrix.
"""
if y is None:
y = x.clone()
if x.ndim != 2 or y.ndim != 2:
raise ValueError("2-dimensional input is expected.")
if x.shape[1] != y.shape[1]:
raise ValueError("`x` and `y` must have the same second dimension.")
sq_dists = cdist(x, y,
compute_mode='use_mm_for_euclid_dist_if_necessary').pow(2)
return sq_dists
def _median(x, y=None):
if y is None:
dtot = distances(x)
else:
dtot = distances(cat((x, y)))
median = dtot.median()
if median == 0:
warnings.warn('The median is null. To avoid a kernel with zero bandwidth, the median is replaced by the mean.')
mean = dtot.mean()
if mean == 0:
warnings.warn('The whole dataset is null.')
return mean
else:
return median
[docs]
def linear_kernel(x, y=None):
"""
Computes the standard linear kernel k(x,y)= <x,y>.
Parameters
----------
x : torch.Tensor
2-d tensor containing the data to kernalize.
y : None or torch.Tensor
2-d tensor containing the data to kernalize. Replaced by `x` if None.
Returns
-------
K : torch.Tensor
Kernel matrix (gram).
"""
if y is None:
y = x.clone()
K = matmul(x, y.T)
return K
[docs]
def gauss_kernel(x, y=None, sigma=1):
"""
Computes the standard Gaussian kernel k(x,y)=exp(- ||x-y||**2 / (2 * sigma**2)).
Parameters
----------
x : torch.Tensor
2-d tensor containing the data to kernalize.
y : None or torch.Tensor
2-d tensor containing the data to kernalize. Replaced by `x` if None.
sigma : int
Standard deviation for Gaussian kernel, 1 by default.
Returns
-------
K : torch.Tensor
Kernel matrix (gram).
"""
d = distances(x, y) # [sq_dists]_ij=||X_j - Y_i \\^2
K = exp(-d / (2 * sigma**2)) # Gram matrix
return K
def _calculate_XXinv_and_ProjImX(X):
XX = matmul(X.T, X)
sp, ev = ordered_eigsy(XX)
# Cut off the spectrum since the matrix is not full rank:
cutoff = np.linalg.matrix_rank(X)
sp = sp[: cutoff]
ev = ev[:, : cutoff]
_XXinv = multi_dot([ev, diag(sp ** -1), ev.T])
_ProjImX = multi_dot([X, _XXinv, X.T])
return _XXinv, _ProjImX
[docs]
class OneHot(object):
"""
Class defining One Hot encoding: a coding scheme for linear models, along
with other well known ones such as Treatment or Difference coding. It is to
be integrated in a formula defining the linear model, provided by patsy's
formula interface. It is recommended to use one hot encoding with the
testing framework implemented in AOV, especially with more than one factor.
Example of a formula with OneHot::
'y1 + y2 ~ C(x1, OneHot) + C(x2, OneHot)'
See Readme and tutorials for kAOV for more details and examples.
"""
def __init__(self, reference=-1):
self.reference = reference
[docs]
def code_with_intercept(self, levels):
return ContrastMatrix(np.eye(len(levels)),
["[%s]" % (level,) for level in levels])
[docs]
def code_without_intercept(self, levels):
return self.code_with_intercept(levels)
[docs]
class Data:
"""
Class containing data-related structures for the kernel analysis of variance.
Parameters
----------
endog : 2-d array_like
An array_like with dimensions nobs x nvar containing nobs values of
nvar dependent variables.
exog : 2-d array_like
An array_like with dimensions nobs x nlvl containing nobs values of
nlvl independent variables.
meta : None or 2-d array_like, optional
An array_like with the metadata for the dataset, i.e. containing
information on factors. Used for visualizations.
endog_names : None or 1-d array_like, optional
A 1-dimensional array_like containing names of the dependent variables.
If not specified (default), will be retrieved from `endog` or assigned
to numbers with respect to the order.
exog_names : None or 1-d array_like, optional
A 1-dimensional array_like containing names of the independent variables.
If not specified (default), will be retrieved from `exog` or assigned
to numbers with respect to the order.
nystrom : bool, optional
If True, computes the Nystrom landmarks, in which case the observations
in all attributes correspond to the landmarks and not the original data.
The default is False.
n_landmarks: int, optional
Number of landmarks used in the Nystrom method. If unspecified, one
fifth of the observations are selected as landmarks.
random_gen : int, Generator, RandomState instance or None
Determines random number generation for the landmarks selection.
If None (default), the generator is the RandomState instance used
by `np.random`. To ensure the results are reproducible, pass an int
to instanciate the seed, or a Generator/RandomState instance (recommended).
Attributes
----------
endog : 2-d torch.tensor
A tensor with dimensions nobs x nvar containing nobs values of
nvar dependent variables.
exog : 2-d torch.tensor
A tensor with dimensions nobs x nlvl containing nobs values of
nlvl independent variables.
meta : None or 2-d array_like
An array_like with the metadata for the dataset, i.e. containing
information on factors. Used for visualizations. The default is None.
nobs : int
Number of observations. If `nystrom=True`, corresponds to the number
of landmarks `n_landmarks`.
nlvl : int
Number of independent variables (i.e. levels of all factors).
nvar : int
Number of dependent variables.
index : 1-d array_like
Observation labels.
endog_names : 1-d array_like
A 1-dimensional array_like containing names of the dependent variables.
exog_names : 1-d array_like
A 1-dimensional array_like containing names of the independent variables.
nystrom : bool
False by default, True if the Nystrom approximation is performed, in
which case the observations in all attributes correspond to the
landmarks and not the original data.
"""
def __init__(self, endog, exog, meta=None, endog_names=None, exog_names=None,
nystrom=False, n_landmarks=None, random_gen=None):
self.exog = convert_to_torch(exog)
self.endog = convert_to_torch(endog)
self.meta = meta
self.index = exog.index if hasattr(exog, "index") else range(self.nobs)
self.nobs = self.exog.shape[0]
self.nlvl = self.exog.shape[1]
self.nvar = self.endog.shape[1]
if endog_names is not None:
self.endog_names = endog_names
elif hasattr(endog, "columns"):
self.endog_names = endog.columns
else:
self.endog_names = ['y' + str(x) for x in range(self.nvar)]
if exog_names is not None:
self.exog_names = exog_names
elif hasattr(exog, "columns"):
self.exog_names = exog.columns
else:
if hasattr(self, '_factor_info') and 'Intercept' in self._factor_info:
self.exog_names = ['x' + str(x + 1) for x in range(self.nlvl - 1)]
self.exog_names.insert(0, 'Intercept')
else:
self.exog_names = ['x' + str(x + 1) for x in range(self.nlvl)]
if self.exog.shape[1] != len(self.exog_names):
raise ValueError("The length of `exog_names` should be equal to "\
"the number of columns in `exog`.")
if self.endog.shape[1] != len(self.endog_names):
raise ValueError("The length of `endog_names` should be equal to "\
"the number of columns in `endog`.")
# Calculate useful matrices:
self._XXinv, self._ProjImX = _calculate_XXinv_and_ProjImX(self.exog)
self.nystrom = nystrom
if self.nystrom:
self.nobs = (n_landmarks if n_landmarks is not None
else min(self.nobs, self.nlvl * 30))
generators = (np.random.RandomState, np.random.Generator)
if isinstance(random_gen, generators):
rnd_gen = random_gen
elif isinstance(random_gen, int):
rnd_gen = np.random.default_rng(random_gen)
else:
rnd_gen = np.random
h_ii = diag(self._ProjImX)
ny_ind = rnd_gen.choice(np.arange(self.exog.shape[0]), size=self.nobs,
p=h_ii/h_ii.sum(), replace=False)
ny_ind.sort()
self.exog = self.exog[ny_ind]
self.endog = self.endog[ny_ind]
if isinstance(self.meta, (pd.Series, pd.DataFrame)):
self.meta = self.meta.iloc[ny_ind]
else:
self.meta = self.meta[ny_ind]
self.index = self.index[ny_ind]
self._XXinv, self._ProjImX = _calculate_XXinv_and_ProjImX(self.exog)
self._ProjImXorthogonal = eye(self.nobs) - self._ProjImX
def _diagonalize_residual_covariance(self, K):
"""
Computes the spectral decomposition of a matrix that shares the
spectrum with the residual covariance operator. A normalized version of
the eigenvectors is stored since it is used in all calculations.
"""
Kresidual = 1 / self.nobs * multi_dot([self._ProjImXorthogonal,
K, self._ProjImXorthogonal])
sp, ev = ordered_eigsy(Kresidual)
sp_power = -1/2 if self.nystrom else -1
sp12 = sp ** sp_power * self.nobs ** (-1/2)
to_ignore = sp12.isnan()
sp12 = sp12[~to_ignore]
U = sp12 * ev[:, ~to_ignore]
U_norm = matmul(self._ProjImXorthogonal, U)
return sp, U_norm
[docs]
class AOV:
"""
Class implementing Kernel Analysis Of Variance.
Parameters
----------
endog : 2-d array_like
An array_like with dimensions nobs x nvar containing nobs values of
nvar dependent variables.
exog : 2-d array_like
An array_like with dimensions nobs x nlvl containing nobs values of
nlvl independent variables.
meta : None or 2-d array_like, optional
An array_like with the metadata for the dataset, i.e. containing
information on factors. Used for visualizations.
endog_names : None or 1-d array_like, optional
A 1-dimensional array_like containing names of the dependent variables.
If not specified (default), will be retrieved from `endog` or assigned
to numbers with respect to the order.
exog_names : None or 1-d array_like, optional
A 1-dimensional array_like containing names of the independent variables.
If not specified (default), will be retrieved from `exog` or assigned
to numbers with respect to the order.
kernel_function : callable or str, optional
Specifies the kernel function. Acceptable values in the form of a
string are 'gauss' (default) and 'linear'. Pass a callable for a
user-defined kernel function.
kernel_bandwidth : 'median' or float, optional
Value of the bandwidth for kernels using a bandwidth. If 'median'
(default), the bandwidth will be set as the median or its multiple,
depending on the value of the parameter `median_coef`. Pass a float
for a user-defined value of the bandwidth.
kernel_median_coef : float, optional
Multiple of the median to compute bandwidth if `kernel_bandwidth='median'`.
The default is 1.
nystrom : bool, optional
If True, computes the Nystrom landmarks, in which case the observations
in all attributes correspond to the landmarks and not the original data.
The default is False.
n_landmarks: int, optional
Number of landmarks used in the Nystrom method. If unspecified, one
fifth of the observations are selected as landmarks.
random_gen : int, Generator, RandomState instance or None
Determines random number generation for the landmarks selection.
If None (default), the generator is the RandomState instance used
by `np.random`. To ensure the results are reproducible, pass an int
to instanciate the seed, or a Generator/RandomState instance (recommended).
Attributes
----------
data : instance of class Data
Contains various information on the original dataset, see the
documentation of the class Data for more details.
data_nystrom : None or instance of class Data
Contains various information on the Nystrom dataset, see the
documentation of the class Data for more details. If None, Nystrom is
not taken into account in the computations.
kernel_function : callable or str, optional
Specifies the kernel function.
kernel_bandwidth : 'median' or float, optional
Value of the bandwidth for kernels using a bandwidth.
kernel_median_coef : float, optional
Multiple of the median to compute bandwidth if `kernel_bandwidth='median'`.
The default is 1.
kernel: callable
Kernel function used for calculations.
computed_bandwidth : float
The value of the kernel bandwidth.
verbose : int, optional
The higher the verbosity, the more messages keeping track of
computations. The default is 0.
- < 1: no messages,
- 1: warnings are printed once,
- 2: warnings are printed every time they appear.
Notes
-----
The `from_formula` interface is the recommended method to specify
a model.
"""
def __init__(self, endog, exog, meta=None, endog_names=None, exog_names=None,
nystrom=False, n_landmarks=None, random_gen=None,
kernel_function='gauss', kernel_bandwidth='median',
kernel_median_coef=1, verbose=1):
if verbose < 1:
warnings.simplefilter("ignore")
elif verbose == 1:
warnings.simplefilter("once")
else:
warnings.simplefilter("always")
self.data = Data(endog, exog, meta=meta, endog_names=endog_names,
exog_names=exog_names)
### Nystrom:
self.data_nystrom = None
if nystrom:
self.data_nystrom = Data(endog, exog, meta=meta, endog_names=endog_names,
exog_names=exog_names, nystrom=True,
n_landmarks=n_landmarks, random_gen=random_gen)
### Kernel:
self.kernel_function = kernel_function
self.kernel_bandwidth = kernel_bandwidth
self.kernel_median_coef = kernel_median_coef
if self.kernel_function == 'gauss':
self.kernel, self.computed_bandwidth = gauss_kernel_median(
x=self.data.endog if not nystrom else self.data_nystrom.endog,
bandwidth=kernel_bandwidth,
median_coef=kernel_median_coef,
return_bandwidth=True
)
elif self.kernel_function == 'linear':
self.kernel = linear_kernel
else:
self.kernel = self.kernel_function
# Diagnostics:
self.diagnostics = None
def __str__(self):
s = "Kernel linear model"
s + f" with {self.data.nobs} observations and {self.data.nvar} features."
if self.data_nystrom is not None:
n_n, g = self.data.endog.shape
s += "\nNystrom approximation with"
s += f" {self.data_nystrom.nobs} landmarks."
return s
[docs]
def compute_diagnostics(self, n_trunc=100, n_anchors=None):
"""
Calculates diagnostics associated with the model and saves them in the
`diagnostics` attribute. The latter is a dictionary containing the
following quantities:
- Embeddings: projections of the embeddings on the first
eigenfunctions of the residual covariance operator.
- Predictions: projections of the predictions of the embeddings on
the first eigenfunctions of the residual covariance operator.
- Residuals: projections of the residuals on the first
eigenfunctions of the residual covariance operator.
Parameters
----------
n_trunc : int, optional
Maximal truncation for projections calculation, the default is 100.
n_anchors : int, optional
Number of anchors used in the Nystrom method. If None, the value is
set at `n_trunc`.
"""
data = self.data_nystrom if self.data_nystrom is not None else self.data
K = self.kernel(data.endog)
sp, U_norm = data._diagonalize_residual_covariance(K)
if self.data_nystrom is not None:
if n_anchors is None:
n_anchors = n_trunc
norm_anchors = U_norm[:, : n_anchors]
sp, U_a = self._diagonalize_covariance_of_anchor_projected_residuals(norm_anchors)
U_norm = multi_dot([norm_anchors, U_a])
K = self.kernel(self.data.endog, self.data_nystrom.endog)
else:
U_norm = sp ** (1/2) * U_norm
n_trunc = min(n_trunc, (sp > 0).sum())
U_norm_T = U_norm[:, : n_trunc]
columns = list(range(1, n_trunc + 1))
embeddings = pd.DataFrame(multi_dot([K, U_norm_T]),
index=self.data.index, columns=columns)
predictions = pd.DataFrame(multi_dot([self.data._ProjImX, K, U_norm_T]),
index=self.data.index, columns=columns)
residuals = pd.DataFrame(multi_dot([self.data._ProjImXorthogonal, K, U_norm_T]),
index=self.data.index, columns=columns)
self.diagnostics = {'Embeddings': embeddings,
'Predictions': predictions,
'Residuals': residuals}
[docs]
def plot_diagnostics(self, trunc=100, diagnostic='residuals', n_trunc=100,
n_anchors=None, colormap='viridis', alpha=.75,
legend_fontsize=12, font_family='serif', figsize=None):
"""
Plots diagnostics associated with the model. The graph contains sublots
for each factor, with the projections of either the residuals
(`diagnostic='residuals'`) or the embeddings (`diagnostic='embeddings'`)
on the first eigenfunctions of the residual covariance operator,
plotted against the projections of the predictions on these eigenfunctions.
Parameters
----------
trunc : int, optional
Truncation to plot, by default equal to the maximal truncation.
diagnostic : str, optional
The type of diagnostic to plot against the predictions.
The default is 'residuals', the alternative is 'embeddings'.
n_trunc : int, optional
Maximal truncation for projections calculation, the default is 100.
n_anchors : int, optional
Number of anchors used in the Nystrom method. If None, the value is
set at `n_trunc`.
colormap : str, optional
The name of a matplotlib colormap to be used for different factor
levels. The default is 'viridis'.
alpha : float, optional
The alpha blending value, between 0 (transparent) and 1 (opaque).
The default is 0.5.
legend_fontsize : int, optional
Legend font size. The default is 15.
font_family : str, optional
Legend and labels' font family name accepted by matplotlib
(e.g., 'serif', 'sans-serif', 'monospace', 'fantasy' or 'cursive'),
the default is 'serif'.
figsize : tuple, optional
The size of the figure. If not specified, is set to
(8 * nb_factors, 6).
Returns
-------
fig : matplotlib.figure.Figure
A Figure object of the plot.
axs : numpy.ndarray of matplotlib.axes._axes.Axes
An Axes object of the plot.
"""
if not self.diagnostics or trunc not in self.diagnostics['Predictions']:
self.compute_diagnostics(n_trunc=max(trunc, n_trunc), n_anchors=n_anchors)
factors = self.data.meta.columns
nb_factors = len(factors)
n_trunc = len(self.diagnostics['Predictions'].columns)
t = min(trunc, n_trunc)
pred = self.diagnostics['Predictions']
a, b = pred[t].min(), pred[t].max()
x = np.arange(a - (b - a) / 10, b + (b - a) / 10, (b - a) / 12)
combs = list(itertools.product(*[self.data.meta[c].unique() for c in self.data.meta.columns]))
if diagnostic == 'residuals':
diagn = self.diagnostics['Residuals']
figtitle = 'Residual plot'
y = [0] * len(x)
elif diagnostic == 'embeddings':
diagn = self.diagnostics['Embeddings']
figtitle = 'Response plot'
y = x.copy()
figtitle += f' (trunc. {t})'
rc('font', **{'family': font_family})
figsize = figsize if figsize is not None else (8 * nb_factors, 6)
fig, axs = plt.subplots(figsize=figsize, ncols=nb_factors)
medianprops = dict(linewidth=0.75, color='black', alpha=alpha)
whiskerprops = dict(linewidth=0.75, alpha=alpha)
for fct, factor in enumerate(factors):
factor_lvls = self.data.meta[factor].unique()
cmap = colormaps[colormap]
colors = cmap(np.linspace(0.1, 0.9, len(factor_lvls)))
ax = axs if nb_factors == 1 else axs[fct]
ax.plot(x, y, color='black', lw=.8, linestyle='--')
ax.set_xlim(a - (b - a) / 10, b + (b - a) / 10)
for c in combs:
c_obs = (self.data.meta == c).all(axis=1)
if c_obs.sum() > 0:
color_ind = np.intersect1d(factor_lvls, c, return_indices=True)[1]
bp = ax.boxplot(diagn[t][c_obs], positions=pred.round(5)[t][c_obs].unique(),
widths=(b - a) / len(combs), patch_artist=True,
manage_ticks=False, medianprops=medianprops,
whiskerprops=whiskerprops)
bp['boxes'][0].set_alpha(alpha)
bp['boxes'][0].set_facecolor(colors[color_ind])
markers = [plt.Line2D([0, 0], [0, 0], color=color, marker='o', linestyle='')
for color in colors]
ax.legend(markers, factor_lvls, numpoints=1, fontsize=legend_fontsize,
bbox_to_anchor=(1.01, 0.5), loc='center left')
ax.set_title(factor, fontsize=18)
ax.set_xlabel('Predictions', fontsize=14)
ax.set_ylabel('Residuals', fontsize=14)
plt.tight_layout()
fig.suptitle(figtitle, fontsize=25, y=1.05)
plt.show()
return fig, axs
[docs]
def set_hypotheses(self, hypotheses=None, by_level=False,
test_intercept=False, true_proportions=False):
"""
Set hypotheses to be tested.
Parameters
----------
hypotheses : str or None or list[tuple]
Hypotheses to be tested.
- if str: either 'pairwise' (default for OneHot) or 'one-vs-all'.
Recommended options in combination with OneHot encoding.
If 'pairwise', the levels of a factor are compared one to another
in the pairwise way. If 'one-vs-all', each level is compared to
the factor mean.
- if None: produces an identity contrast matrix for each factor.
Intended for other coding schemes (e.g. Treatment, Sum, etc).
- if list[tuple]: custom hypothesis option. Each element of the
list should be a tuple of size 2: `(name, contrast_L)`, where `name`
is a string and contrast_L is a contrast matrix in the form of a
torch.tensor (`dtype=torch.float64`).
by_level : bool, optional
If False (default), computes the global test. If True, computes the
test by level or by a pair of levels.
test_intercept : bool, optional
If True, adds a test for the intercept, which is set as the
grand mean (or the actual mean if `true_proportions=True`) of all the
level effects. The default is False. It is unnecessary to add an
intercept test manually eith this option if an intercept is present
in the design matrix.
true_proportions : bool, optional
Relevant for the calculation of the factor mean, i.e. if
`hypotheses='one-vs-all'` or `test_intercept=True`. If False (default),
the factor mean is the grand mean (mean of means). If True, the
true level proportions are taken into account, so the factor mean
is the actual global mean of the factor.
Returns
-------
hypotheses : list[tuple]
List of hypotheses to be tested. Each element of the list is a
tuple of size 2: `(name, contrast_L)`, where `name` is a string and
contrast_L is a contrast matrix in the form of a torch.tensor.
"""
if hypotheses is None and hasattr(self, 'formula') and 'OneHot' in self.formula:
hypotheses = 'pairwise'
if isinstance(hypotheses, str):
if not hasattr(self, 'formula') or 'OneHot' not in self.formula:
warnings.warn("'pairwise' and 'one-vs-all' options for `hypotheses` "
"are intended for OneHot encoding. Otherwise, "
"the test might be difficult to interpret.")
hypotheses_type = hypotheses
hypotheses = []
if not hasattr(self, '_factor_info'):
warnings.warn("The absence of `_factor_info` attribute was "
"caused by not using the `from_formula` interface. "
"With multiple factors, using `from_formula` "
"is highly recommended.")
self._factor_info = {'Factor': slice(0, None, None)}
for _name, _slice in self._factor_info.items():
lvls = self.data.exog_names[_slice]
DI = DesignInfo(self.data.exog_names)
if hypotheses_type == 'one-vs-all' or (test_intercept and _name != 'Intercept'):
if true_proportions:
proportions = (np.array(self.data.exog.sum(axis=0, dtype=int)[_slice])
.astype(str))
lvls_p = [proportions[i] + ' * ' + lvls[i] for i in range(len(lvls))]
factor_mean = '(' + ' + '.join(lvls_p) + f') / {self.data.nobs}'
else:
factor_mean = '(' + ' + '.join(lvls) + f') / {len(lvls)}'
if not by_level:
if hypotheses_type == 'pairwise':
formula = ' = '.join(lvls)
elif hypotheses_type == 'one-vs-all':
formula = [f + ' = ' + factor_mean for f in lvls[:-1]]
L = DI.linear_constraint(formula).coefs
hypotheses.append([_name, convert_to_torch(L)])
else:
if hypotheses_type == 'pairwise':
combs = list(itertools.combinations(lvls, 2))
formula = [' = '.join(comb) for comb in combs]
hyp_name = formula
elif hypotheses_type == 'one-vs-all':
formula = [l + ' = ' + factor_mean for l in lvls]
hyp_name = [l + ' = ' + _name + ' Grand Mean' for l in lvls]
L = DI.linear_constraint(formula).coefs
hypotheses.extend(list(zip(hyp_name, convert_to_torch(L))))
if test_intercept and _name != 'Intercept':
if len(self._factor_info.items()) == 1:
L = DI.linear_constraint(factor_mean).coefs
hypotheses.insert(0, ['Intercept', convert_to_torch(L)])
else:
raise ValueError("`test_intercept=True` is only "
"accepted in the 1-factor setting.")
elif hypotheses is None:
if not by_level:
if not hasattr(self, '_factor_info'):
warnings.warn("If `hypotheses` is not specified, "
"and `from_formula` interface not used to define the model, "
"all `exog_names` are considered as levels of one factor.")
self._factor_info = {'Factor': slice(0, None, None)}
hypotheses = []
for _name, _slice in self._factor_info.items():
L = eye(self.data.nlvl, dtype=float64)[_slice, :]
if 'Intercept' not in self._factor_info:
L = L[1:]
hypotheses.append([_name, L])
if test_intercept and 'Intercept' not in self._factor_info:
raise ValueError("Incompatible combination: "
"`test_intercept=True` and `hypotheses=None`.")
else:
hypotheses = list(zip(self.data.exog_names, eye(self.data.nlvl, dtype=float64)))
if test_intercept:
raise ValueError("Incompatible combination: "
"`test_intercept=True` and `hypotheses=None`.")
return hypotheses
def _compute_LXXLinv(self, L):
"""
Intermediate quantity used in statistics calculations, based on a
transformation of the design matrix X and the contrast matrix L.
"""
LXXL = multi_dot([L, self.data._XXinv, L.T])
if len(LXXL.shape) == 0:
LXXLinv = tensor([[LXXL**-1]])
else:
cutoff = np.linalg.matrix_rank(LXXL)
if cutoff < len(LXXL):
# Generalized inverse:
sp, ev = ordered_eigsy(LXXL)
sp = sp[: cutoff]
ev = ev[:, : cutoff]
LXXLinv = multi_dot([ev, diag(sp ** -1), ev.T])
else:
LXXLinv = LXXL.inverse()
return LXXLinv
def _compute_D(self, L):
"""
Intermediate quantity used in statistics calculations, based on a
transformation of the design matrix X and the contrast matrix L.
"""
LXXLinv = self._compute_LXXLinv(L)
A = multi_dot([L.T, LXXLinv, L])
D = multi_dot([self.data.exog, self.data._XXinv, A, self.data._XXinv, self.data.exog.T])
return D
def _diagonalize_covariance_of_anchor_projected_residuals(self, anchors):
"""
Computes the kernel trick version of the covariance operator associated
with the residuals projected onto the Nystrom anchors. Returns its
eigendecomposition.
"""
K_YZ = self.kernel(self.data.endog, self.data_nystrom.endog)
K_anchor_projected = multi_dot([anchors.T, K_YZ.T,
self.data._ProjImXorthogonal,
K_YZ, anchors])
K_anchor_projected *= (1/self.data.nobs)
return ordered_eigsy(K_anchor_projected)
def _compute_K_T(self, n_trunc=100, n_anchors=None):
"""
Intermediate quantity used in statistics calculations, based on a
transformation of the gram matrix K.
"""
data = self.data_nystrom if self.data_nystrom is not None else self.data
K = self.kernel(data.endog)
sp, U_norm = data._diagonalize_residual_covariance(K)
if self.data_nystrom is not None:
norm_anchors = U_norm[:, : n_anchors]
sp, U_a = self._diagonalize_covariance_of_anchor_projected_residuals(norm_anchors)
U_norm = multi_dot([norm_anchors, (sp ** (-1/2) * U_a)])
K = self.kernel(self.data_nystrom.endog, self.data.endog)
n_trunc = min(n_trunc, (sp > 0).sum())
U_norm_T = U_norm[:, : n_trunc]
K_T = multi_dot([U_norm_T.T, K])
return K_T
def _compute_kernel_HL_test_statistic(self, K_T, D, d, n_trunc=100, f_norm=True):
"""
Computes the truncated kernel Hotelling-Lawley test statistic based on
some intermediate quantities.
"""
K_T_D_K_T = multi_dot([K_T, D, K_T.T])
stats = []
pvals = []
for t in range(1, n_trunc):
norm_factor = ((self.data.nobs - self.data.nlvl) / (self.data.nobs * d * t)
if f_norm else 1)
stat = trace(K_T_D_K_T[:t, :t]).item() * norm_factor
stats += [stat]
if f_norm:
pvals += [f.sf(stat, t * d, self.data.nobs - t * d)]
else:
pvals += [chi2.sf(stat, t * d)]
return stats, pvals
[docs]
def project_on_discriminant(self, K, K_T, D, center=True):
"""
Computes of embeddings of each observation on the discriminant axes
associated with the test. The latter are correspond to the
eigenfunctions of the test statistic operator.
Parameters
----------
K : torch.tensor
The gram matrix.
K_T : torch.tensor
A transformed gram matrix, obtained with `_compute_K_T`.
D : torch.tensor
A quantity containing the information on the contrast matrix,
obtained with `_compute_D`.
center : bool, optional
If True (default), the projections are centered with respect to
the factor mean.
Returns
-------
proj : pandas.DataFrame
Contains projection values, with rows corresponding to
observations, and columns to truncations of the residual covariance
operator.
"""
K_T_D_K_T = multi_dot([K_T, D, K_T.T])
sp, ev = ordered_eigsy(K_T_D_K_T)
# Cut off the spectrum since the matrix is not full rank:
cutoff = np.linalg.matrix_rank(K_T_D_K_T)
sp = sp[: cutoff]
ev = ev[:, : cutoff]
norm = diag(multi_dot([ev.T, K_T, D, K, D, K_T.T, ev]))
norm = pow(norm, -1/2)
centering_mat = ones((self.data.nobs, self.data.nobs), dtype=float64) / self.data.nobs
K_ = K - matmul(K, centering_mat) if center else K
proj = norm * multi_dot([ev.T, K_T, D, K_]).T
proj = pd.DataFrame(proj, index=self.data.index)
proj.columns += 1
return proj
[docs]
def compute_cook_distances(self, L, n_trunc=100, normalize=True):
"""
Computes influence measures in the form of Cook's distances for each
observation as well as the corresponding p-values. The p-values are
obtained based on the Beta distribution associated with normalized
Cook's distances.
Parameters
----------
L : torch.tensor
Contrast matrix, defining a statistical test of interest.
n_trunc : int, optional
Maximal truncation of the residual covariance operator,
the default is 100.
normalize : bool, optional
If True (default), the distances are normalized with respect to
the design and contrasts.
Returns
-------
cook_distances : pandas.DataFrame
Contains influence measure values, with rows corresponding to
observations, and columns to truncations of the residual covariance
operator.
cook_pvalues : pandas.DataFrame
Contains p-values associated with the influence measure, with rows
corresponding to observations, and columns to truncations of the
residual covariance operator.
"""
W = multi_dot([self.data.exog, self.data._XXinv, L.T])
LXXLinv = self._compute_LXXLinv(L)
coef_D = diag(multi_dot([W, LXXLinv, W.T]))
one_minus_hii = (1 - diag(self.data._ProjImX))
coef_D /= (one_minus_hii ** 2 * L.shape[0])
LL = matmul(L, L.T)
L_inv = matmul(L.T, inv(LL))
XX = matmul(self.data.exog.T, self.data.exog)
hii_star = diag(multi_dot([self.data.exog, self.data._XXinv, L_inv, L, XX, L_inv, L,
self.data._XXinv, self.data.exog.T]))
C_i = hii_star * (self.data.nobs - L.shape[0]) / (L.shape[0] * one_minus_hii)
K_T = self._compute_K_T(n_trunc=n_trunc)
cook_distances = {}
cook_pvalues = {}
for t in range(1, n_trunc + 1):
# Cook distances:
cook_traces = diag(multi_dot([self.data._ProjImXorthogonal, K_T[:t].T,
K_T[:t], self.data._ProjImXorthogonal]))
cook_distances_norm = (cook_traces * coef_D).numpy() / np.array(C_i)
cook_distances[t] = (cook_distances_norm if normalize else
(cook_traces * coef_D).numpy())
# P-values:
a = t / 2
b = (self.data.nobs - L.shape[0] - t) / 2
cook_pvalues[t] = beta(a, b).sf(cook_distances_norm)
cook_distances = pd.DataFrame(cook_distances, index=self.data.index)
cook_pvalues = pd.DataFrame(cook_pvalues, index=self.data.index)
return cook_distances, cook_pvalues
[docs]
def correct_pvalues(self, pvalues, correction='bonferroni', by_level=False):
"""
Corrects the p-values according to the chosen correction strategy.
Parameters
----------
pvalues : pandas.DataFrame
Data frame with p-values to correct. Columns indicate hypotheses
that are tested and lines indicate truncation levels.
correction : str, optional
Relevent for multiple test comparisons, in particular when
`by_level=True`. If 'bonferroni' (default), permorms the Bonferroni
correction of the p-values. If 'BH', perfoms the Benjamini–Hochberg
correction.
by_level : by_level : bool, optional
If False (default), computes the global test. If True, computes the
test by level or by a pair of levels.
"""
if hasattr(self, '_factor_info') and by_level:
factor_cols = [pvalues.columns.str.contains(k)
for k in self._factor_info.keys()]
else:
factor_cols = [[True,] * len(pvalues.columns),]
for fct in factor_cols:
pvalues.loc[:, fct] *= np.sum(fct)
if correction == 'BH':
pval_ranks = pvalues.loc[:, fct].rank(axis=1, method='dense')
pvalues.loc[:, fct] /= pval_ranks
pvalues.clip(upper=1, inplace=True)
[docs]
def test(self, hypotheses=None, hypotheses_subset=None,
by_level=False, n_trunc=100, correction=None, test_intercept=False,
true_proportions=False, center_projections=True, verbose=1,
n_anchors=None, f_norm=True, skip_projections_and_cook=False,
norm_cook=True):
"""
Performs kernel hypothesis tests for the given model. Simultaneously
calculates projections on the associated discriminant axes as well as
influences of observations with respect to the test with their p-values.
Parameters
----------
hypotheses : str or None or list[tuple]
Hypotheses to be tested.
- if str: either 'pairwise' (default for OneHot) or 'one-vs-all'.
Recommended options in combination with OneHot encoding.
- if None: produces an identity contrast matrix for each factor.
Intended for other coding schemes (e.g. Treatment, Sum, etc).
- if list[tuple]: custom hypothesis option. Each element of the
list should be a tuple of size 2: `(name, contrast_L)`, where `name`
is a string and contrast_L is a contrast matrix in the form of a
torch.tensor (`dtype=torch.float64`).
hypotheses_subset : list of strings
Names of tests to perform, subset of all the tests in the
`hypotheses` variable (particularly useful with the by_level
testing option, when the total number of hypotheses is high and
the interest lies in the subset). The default is None, i.e. all
hypotheses are tested.
by_level : bool, optional
If False (default), computes the global test. If True, computes the
test by level or by a pair of levels.
n_trunc : int, optional
Maximal truncation for statistics calculation, the default is 100.
correction : str or None, optional
Relevent for multiple test comparisons, in particular when
`by_level=True`. If 'bonferroni', permorms the Bonferroni correction
of the p-values. If 'BH', perfoms the Benjamini–Hochberg
correction. If None (default), the p-values remain uncorrected.
test_intercept : bool, optional
If True, adds a test for the intercept, which is set as the
grand mean (or the actual mean if `true_proportions=True`) of all the
level effects. The default is False. It is unnecessary to add an
intercept test manually eith this option if an intercept is present
in the design matrix.
true_proportions : bool, optional
Relevant for the calculation of the factor mean, i.e. if
`hypotheses='one-vs-all'` or `test_intercept=True`. If False (default),
the factor mean is the grand mean (mean of means). If True, the
true level proportions are taken into account, so the factor mean
is the actual global mean of the factor.
center_projections : bool, optional
If True (default), the projections are centered with respect to
the factor mean.
n_anchors : int, optional
Number of anchors used in the Nystrom method. If None, the value is
set at `n_trunc`.
verbose : int, optional
The higher the verbosity, the more messages keeping track of
computations. The default is 0.
- < 1: no messages,
- 1: progress bar with computation time,
- 2: print tested hypothesis' name, warnings are printed once,
- 3: warnings are printed every time they appear.
f_norm : bool, optional
If True (default), the test statistic is normalized and
asymptotically follows an f-distribution. Otherwise, the original
chi-2 version is returned.
skip_projections_and_cook : bool, optional
If False (default), projections on the discriminant axes as well as
Cook's distances will be computed. Set to True to avoid computing
them if they are not needed, in order to reduce computation time.
norm_cook : bool, optional
If True (default), the Cook's distances are normalized with
respect to the design and contrasts.
Returns
-------
KernelAOVResults object
See the documentation for KernelAOVResults.
Examples
--------
Importing data and create an intsance of AOV:
>>> import pandas as pd
>>> from kaov import AOV
>>> url = "https://raw.githubusercontent.com/LMJL-Alea/kAOV/refs/heads/main/Data/reversion_kAOV.csv"
>>> data = pd.read_csv(url, index_col=0)
>>> kfit = AOV.from_formula('AACS + ACSL6 + ACSS1 ~ C(Medium, OneHot)', data=data)
Test for factor effects:
>>> res = kfit.test()
>>> print(res)
Kernel Analysis of Variance (trunc. 1):
====================================
------------------------------------
Factor test | factor stat pval
------------------------------------
| Medium 34.0842 0.0000
====================================
"""
if verbose > 0:
print('-Computing the Gram matrix...')
if verbose <= 1:
warnings.simplefilter("ignore")
elif verbose == 2:
warnings.simplefilter("once")
else:
warnings.simplefilter("always")
K = self.kernel(self.data.endog)
if n_anchors is None:
n_anchors = n_trunc
K_T = self._compute_K_T(n_trunc=n_trunc, n_anchors=n_anchors)
if verbose > 0:
print('-Testing hypotheses:')
if hypotheses is None and hasattr(self, 'formula') and 'OneHot' in self.formula:
hypotheses = 'pairwise'
hyps = self.set_hypotheses(hypotheses=hypotheses, by_level=by_level,
test_intercept=test_intercept,
true_proportions=true_proportions)
if hypotheses_subset is not None:
hyps = [h for h in hyps if h[0] in hypotheses_subset]
results = {}
projections = {}
cook_distances = {}
cook_pvalues = {}
if correction is not None:
corrected_pvals_dict = {}
# Get factor dummies to add the factor information to the data frames:
factor_dummies = pd.DataFrame(self.data.exog.numpy().astype(int),
columns=self.data.exog_names, index=self.data.index)
it = tqdm(range(len(hyps))) if verbose > 0 else range(len(hyps))
for i in it:
name, L = hyps[i]
if verbose > 1:
print('\n')
print(f'-Testing {name}...')
if any(isinstance(l, str) for l in L):
L = DesignInfo(self.data.exog_names).linear_constraint(L).coefs
L = convert_to_torch(L)
L = L.unsqueeze(0) if L.dim() == 1 else L
D = self._compute_D(L)
stats, pvals = self._compute_kernel_HL_test_statistic(K_T, D, len(L),
n_trunc=n_trunc,
f_norm=f_norm)
results_dict = {'stat': stats,
'pval': pvals}
results[name] = pd.DataFrame(results_dict)
results[name].index.name = 'Truncation'
results[name].index += 1
if correction is not None:
corrected_pvals_dict[name] = pvals
# Projections on discriminant axes and Cook's distances:
if not skip_projections_and_cook:
projections[name] = self.project_on_discriminant(K, K_T, D,
center=center_projections)
(cook_distances[name],
cook_pvalues[name]) = self.compute_cook_distances(L, n_trunc=n_trunc,
normalize=norm_cook)
if not hasattr(self, 'formula') or hypotheses not in [None, 'pairwise', 'one-vs-all']:
projections[name][name] = np.nan
cook_distances[name][name] = np.nan
elif name == 'Intercept':
pass
else:
# Adding factor information:
pd.options.mode.chained_assignment = None # default='warn'
if not by_level: # specify all levels for all factors
_slice = self._factor_info[name]
dummies_factor_i = factor_dummies.iloc[:, _slice]
else: # specify only those levels that are relevant for a given test
if hypotheses == 'one-vs-all':
_slice = self._factor_info[name.split(' = ')[1][:-11]]
else:
_slice = [i for i, en in enumerate(self.data.exog_names) if en in name]
dummies_factor_i = factor_dummies.iloc[:, _slice]
# Case of interaction effects:
if ':' in name:
interaction_cols = dummies_factor_i.columns.str.contains(':')
dummies_factor_i = dummies_factor_i.loc[:, interaction_cols]
# Put NA for the irrelevant levels (needed for visualizations)
nan_obs = (dummies_factor_i.max(axis=1) == 0)
dummies_factor_i['NA'] = 0
dummies_factor_i.loc[nan_obs, 'NA'] = 2
projections[name][name] = dummies_factor_i.idxmax(axis=1)
cook_distances[name][name] = dummies_factor_i.idxmax(axis=1)
cook_pvalues[name][name] = dummies_factor_i.idxmax(axis=1)
# Correct p-values:
if correction is not None:
if verbose > 0:
print('-Correcting p-values for multiple tests...')
corrected_pvals_df = pd.DataFrame(corrected_pvals_dict)
corrected_pvals_df.index += 1
self.correct_pvalues(corrected_pvals_df, correction=correction, by_level=by_level)
for hyp in hyps:
name, L = hyp
results[name]['pval'] = corrected_pvals_df[name]
hypothesis_type = (hypotheses if hypotheses is None or type(hypotheses) == str
else 'custom')
return KernelAOVResults(hyps, results, projections,
cook_distances, cook_pvalues,
hypothesis_type=hypothesis_type,
by_level=by_level, factor_info=self._factor_info)
[docs]
class KernelAOVResults():
"""
Class implementing Kernel Analysis Of Variance.
Parameters
----------
hypotheses : list
All the consideres hypotheses. Each element of the list represents a
hypothesis in the form of a list with two elements: a name and a
congtrast matrix associated with the test.
stats : dict
A dictionary with keys corresponding to hypothesis names, and values
to instances of pandas.DataFrame with the results of the corresponding
tests. Each data frame contains two columns, the first containing the
truncated kernel Hotelling-Lawley test statistic values, and the
second containing the associated p-values, indexed by truncations of the
residual covariance operator used in the calculations.
projections : dict
A dictionary with keys corresponding to hypothesis names, and values
to instances of pandas.DataFrame with the projections obtained with
`AOV.project_on_discriminant`. In each data frame, rows correspond to
observations, and columns to truncations of the residual covariance
operator.
cook_distances : dict
A dictionary with keys corresponding to hypothesis names, and values
to instances of pandas.DataFrame with the influence measure values
obtained with `AOV.compute_cook_distances`. In each data frame, rows
correspond to observations, and columns to truncations of the
residual covariance operator.
cook_pvalues : dict
A dictionary with keys corresponding to hypothesis names, and values
to instances of pandas.DataFrame with the p-values associated with the
influence measures obtained with `AOV.compute_cook_distances`.
In each data frame, rows correspond to observations, and columns to
truncations of the residual covariance operator.
hypothesis_type : str or None
Types of hypotheses that were tested. If not None, expected options
are 'pairwise', 'one-vs-all' and 'custom'.
by_level : bool, optional
False if global (factor-wise) tests were perfromed, True if tests by
level or by a pair of levels were perfromed.
factor_info : dict
_factor_info attribute of the AOV class.
Attributes
----------
hypotheses : list
See Parameters.
stats : dict
See Parameters.
projections : dict
See Parameters.
cook_distances : dict
See Parameters.
cook_pvalues : dict
See Parameters.
hypothesis_type : str or None
See Parameters.
by_level : bool
See Parameters.
"""
def __init__(self, hypotheses, stats, projections, cook_distances,
cook_pvalues, hypothesis_type, by_level, factor_info):
self.hypotheses = hypotheses
self.stats = stats
self.projections = projections
self.cook_distances = cook_distances
self.cook_pvalues = cook_pvalues
self.hypothesis_type = hypothesis_type
self.by_level = by_level
self._factor_info = factor_info
[docs]
def summary(self, trunc, factor=None):
"""
Creates a pandas.DataFrame or a distionary with a summary of the test
for a given truncation and factor (in the by-level case).
Parameters
----------
trunc : int
Truncation for which to return the test results.
factor : str or None
None by default, in which case the results of tests for each factor
are returned. If the factor is specified, returns the results of
tests on comparisons related with the chosen factor.
Returns
-------
sum_df : pandas.DataFrame or pandas.Series or dict
A data frame with the summary of test results. If not by_level and
factor is specified, returns a row of this data frame corresponding
to the chosen factor. If by_level and factor is not specified,
returns a dictionary with data frames for each factor.
"""
if not self.by_level or self.hypothesis_type == 'custom':
if factor is not None:
summ = pd.Series(index=['factor', 'stat', 'pval'], dtype=object)
summ['factor'] = factor
summ['stat'] = self.stats[factor].loc[trunc, 'stat']
summ['pval'] = self.stats[factor].loc[trunc, 'pval']
else:
summ = pd.DataFrame(columns=['factor', 'stat', 'pval'],
index=np.arange(1, len(self.stats) + 1),
dtype=object)
for i, (hyp, stat) in enumerate(self.stats.items()):
summ.loc[i + 1, 'factor'] = hyp
summ.loc[i + 1, 'stat'] = self.stats[hyp].loc[trunc, 'stat']
summ.loc[i + 1, 'pval'] = self.stats[hyp].loc[trunc, 'pval']
else:
summ = {}
if factor is None:
factors = list(self._factor_info.keys())
else:
factors = [factor, ]
for fct in factors:
nb_factors = fct.count(':') + 1
if self.hypothesis_type is None:
factor_cols = [f'factor_{f + 1}' for f in range(nb_factors)]
else:
factor_cols = [f'factor_{f + 1}_1' for f in range(nb_factors)]
factor_cols.extend([f'factor_{f + 1}_2' for f in range(nb_factors)])
sum_df = pd.DataFrame(columns=factor_cols + ['stat', 'pval'],
index=np.arange(1, len(self.stats) + 1),
dtype=object)
for i, (hyp, stat) in enumerate(self.stats.items()):
is_factor_hyp = ((nb_factors == 1 and fct in hyp and ':' not in hyp)
or (nb_factors > 1 and np.all([f in hyp for f in fct.split(':')])))
if is_factor_hyp:
levels = re.findall(r"\[(.*?)\]", hyp)
for j, lvl in enumerate(levels):
sum_df.iloc[i, j] = lvl
sum_df.loc[i + 1, 'stat'] = self.stats[hyp].loc[trunc, 'stat']
sum_df.loc[i + 1, 'pval'] = self.stats[hyp].loc[trunc, 'pval']
if self.hypothesis_type == 'one-vs-all':
sum_df.iloc[i, nb_factors : 2 * nb_factors] = 'Grand Mean'
sum_df.dropna(axis=0, inplace=True)
sum_df.index = np.arange(1, len(sum_df) + 1)
summ[fct] = sum_df
if factor is not None:
summ = summ[factor]
else:
# Pop emty factors from the dictionary:
empty_fct = []
for key, val in summ.items():
if val.empty:
empty_fct.append(key)
for key in empty_fct:
summ.pop(key)
return summ
def _summary_obj(self):
"""
Creates a summary object to display a summary of the test based on the
`summary` method.
Returns
-------
An instance of statsmodels.iolib.summary2
"""
t = 1
float_format = '%.3f'
summ = self.summary(trunc=t)
if type(summ) != dict:
summ = {'Factor test' : summ}
summ_print = summary2.Summary()
summ_print.add_title(f'Kernel Analysis of Variance (trunc. {t}):')
for i, (key, df) in enumerate(summ.items()):
summ_print.add_dict({'' : ''})
df.index = [' |',] * len(df)
df = df.reset_index()
c = list(df.columns)
c[0] = key + ' |'
df.columns = c
df.index = ['',] * len(df)
summ_print.add_df(df, float_format=float_format)
return summ_print
def __str__(self):
return self._summary_obj().__str__()
[docs]
def plot_density(self, comp=1, tests=None, colormap='viridis', alpha=.5,
legend_fontsize=12, font_family='serif', figsize=None):
"""
Plots kernel-densities of projections of the embeddings on the chosen
discriminant axis, associated with the tests underlying the
KernelAOVResults object. Produces separate subplots for each test.
Parameters
----------
comp : int, optional
Component to plot, i.e. the embeddings are projected on the comp-th
eigenfunction.
tests : list of strings or None
List containing names of tests to plot, out of all the tests in
the KernelAOVResults object (particularly useful with the by_level
testing option). The default is None, i.e. all tests are plotted.
colormap : str, optional
The name of a matplotlib colormap to be used for different factor
levels. The default is 'viridis'.
alpha : float, optional
The alpha blending value, between 0 (transparent) and 1 (opaque).
The default is 0.5.
legend_fontsize : int, optional
Legend font size. The default is 15.
font_family : str, optional
Legend and labels' font family name accepted by matplotlib
(e.g., 'serif', 'sans-serif', 'monospace', 'fantasy' or 'cursive'),
the default is 'serif'.
figsize : tuple, optional
The size of the figure. If not specified, is set to
(8 * nb_factors, 6).
Returns
-------
fig : matplotlib.figure.Figure
A Figure object of the plot.
axs : numpy.ndarray of matplotlib.axes._axes.Axes
An Axes object of the plot.
"""
warnings.simplefilter("always")
tests = self.projections.keys() if tests is None else tests
nb_tests = len(tests)
rc('font', **{'family': font_family})
figsize = figsize if figsize is not None else (8 * nb_tests, 6)
fig, axs = plt.subplots(ncols=nb_tests, figsize=figsize)
for j, test in enumerate(tests):
ax = axs if nb_tests == 1 else axs[j]
T_max = len(self.projections[test].columns)
if test in self.projections[test]:
T_max -= 1
t = min(comp, T_max)
if t != comp:
warnings.warn(f'In {test}: comp={t} will be plotted, since {comp} ' +
'is larger than the maximal number of components.')
proj_j = self.projections[test]
test_lvls = proj_j[test].unique()
test_lvls = test_lvls[test_lvls != 'NA'] # extract relevant observations
cmap = colormaps[colormap]
colors = cmap(np.linspace(0.1, 0.9, len(test_lvls)))
no_lvl_info = proj_j[test].isnull().all()
for i, test_lvl in enumerate(test_lvls):
if no_lvl_info:
lvl_proj = proj_j[t]
else:
lvl_proj = proj_j[proj_j[test] == test_lvl][t]
min_proj, max_proj = lvl_proj.min(), lvl_proj.max()
min_scaled = min_proj - 0.1 * (max_proj - min_proj)
max_scaled = max_proj + 0.1 * (max_proj - min_proj)
x = np.linspace(min_scaled, max_scaled, 200)
try:
density = gaussian_kde(lvl_proj, bw_method=.2)
y = density(x)
ax.plot(x, y, color=colors[i], lw=2)
ax.fill_between(x, y, y2=0, color=colors[i], label=test_lvl, alpha=alpha)
except ValueError:
pass
if not no_lvl_info:
ax.legend(bbox_to_anchor=(1.01, 0.5), loc='center left',
fontsize=legend_fontsize)
ax.set_title(test, fontsize=18)
ax.set_xlabel('Discriminant axis', fontsize=14)
ax.set_ylabel('Density', fontsize=14)
plt.tight_layout()
fig.suptitle(f'Discriminant axis projection density (trunc. {t})',
fontsize=25, y=1.05)
plt.draw()
return fig, axs
[docs]
def plot_mean_embedding_projections(self, comp1=1, comp2=2, tests=None,
figsize=None, ylim=None, xlim=None,
alpha=1, s=50, marker='o',
colors=None, colormap='viridis',
font_family='serif', legend=True,
legend_fontsize=15, figtitle=None):
"""
Plots projections of mean embeddings of relevat groups on the chosen
discriminant axis, associated with the tests underlying the
KernelAOVResults object. Produces separate subplots for each test.
Parameters
----------
comp1 : int, optional
Component x of the plot, i.e. the mean embeddings are projected on the
comp1-th eigenfunction.
comp2 : int, optional
Component y of the plot, i.e. the mean embeddings are projected on the
comp2-th eigenfunction.
tests : list of strings or None
List containing names of tests to plot, out of all the tests in
the KernelAOVResults object (particularly useful with the by_level
testing option). The default is None, i.e. all tests are plotted.
figsize : tuple, optional
The size of the figure. If not specified, is set to
(8 * nb_factors, 6).
alpha : float, optional
The alpha blending value, between 0 (transparent) and 1 (opaque).
The default is 1.
s : int or dict
Marker sizes of the mean embedding projections. The default is 50.
If int, same sizes for all mean embeddings. To pass different
values for different tests and levels, pass a dictionary with keys
corresponding to test names, and values that are dictionaries with
keys corresponding to the levels of the test and values that are
marker size integers.
marker : str or dict
Marker styles of the mean embedding projections. The default is 'o'.
If string, same marker for all mean embeddings. To pass different
values for different tests and levels, pass a dictionary with keys
corresponding to test names, and values that are dictionaries with
keys corresponding to the levels of the test and values that are
marker style strings.
colors : None or dict
Colors of the mean embedding projections. If None (default), colors
are chosen from the specified colormap. To customize, pass a
dictionary with keys corresponding to test names, and values that
are dictionaries with keys corresponding to the levels of the test
and values that are colors.
colormap : str, optional
The name of a matplotlib colormap to be used for different factor
levels (if colors are not provided). The default is 'viridis'.
font_family : str, optional
Legend and labels' font family name accepted by matplotlib
(e.g., 'serif', 'sans-serif', 'monospace', 'fantasy' or 'cursive'),
the default is 'serif'.
legend : bool, optional
If True (default), legend is plotted automatically.
legend_fontsize : int, optional
Legend font size. The default is 15.
figtitle : str, optional
The title of the figure.
Returns
-------
fig : matplotlib.figure.Figure
A Figure object of the plot.
axs : numpy.ndarray of matplotlib.axes._axes.Axes
An Axes object of the plot.
"""
warnings.simplefilter("always")
tests = self.projections.keys() if tests is None else tests
nb_tests = len(tests)
rc('font', **{'family': font_family})
figsize = figsize if figsize is not None else (8 * nb_tests, 6)
fig, axs = plt.subplots(ncols=nb_tests, figsize=figsize,
sharex=True, sharey=True)
for j, test in enumerate(tests):
ax = axs if nb_tests == 1 else axs[j]
ax.yaxis.set_tick_params(labelleft=True)
ax.axvline(0, color='grey', linestyle='--', alpha=0.5)
ax.axhline(0, color='grey', linestyle='--', alpha=0.5)
proj_j = self.projections[test]
t1_j = comp1
t2_j = t1_j if comp2 not in proj_j.columns else comp2 # sometimes only one axis available
if t2_j != comp2:
warnings.warn(f'In {test}: comp2={t2_j} will be plotted, since ' +
f'{comp2} is larger than the maximal number of components.')
test_lvls = proj_j[test].unique()
test_lvls = test_lvls[test_lvls != 'NA'] # extract relevant observations
nb_lvls = len(test_lvls)
# Create dictionaries for colors, markers and marker sizes:
if colors is None:
cmap = colormaps[colormap]
cmap_colors = cmap(np.linspace(0.1, 0.9, nb_lvls))
color = {test_lvls[i] : cmap_colors[i] for i in range(nb_lvls)}
elif type(colors) is dict:
color = colors[test]
if type(marker) is str:
markers = {test_lvls[i] : marker for i in range(nb_lvls)}
elif type(marker) is dict:
markers = marker[test]
if type(s) is int:
size = {test_lvls[i] : s for i in range(nb_lvls)}
elif type(s) is dict:
size = s[test]
no_lvl_info = proj_j[test].isnull().all()
# Scatter plots for each test:
for i, test_lvl in enumerate(test_lvls):
if no_lvl_info:
lvl_proj_1 = proj_j[t1_j]
lvl_proj_2 = proj_j[t2_j]
else:
lvl_proj_1 = proj_j[proj_j[test] == test_lvl][t1_j]
lvl_proj_2 = proj_j[proj_j[test] == test_lvl][t2_j]
try:
ax.scatter(lvl_proj_1.mean(), lvl_proj_2.mean(), alpha=alpha,
linewidths=3, color=color[test_lvl], label=test_lvl,
s=size[test_lvl], marker=markers[test_lvl])
except ValueError:
pass
if legend and not no_lvl_info:
ax.legend(bbox_to_anchor=(1.01, 0.5), loc='center left',
fontsize=legend_fontsize)
ax.set_title(test, fontsize=18)
ax.set_xlabel(f'Discriminant axis {t1_j}', fontsize=14)
ax.set_ylabel(f'Discriminant axis {t2_j}', fontsize=14)
plt.tight_layout()
if figtitle is not None:
fig.suptitle(figtitle)
else:
fig.suptitle('Discriminant axis projection of mean embeddings',
fontsize=25, y=1.05)
plt.draw()
return fig, axs
[docs]
def plot_influence(self, trunc=1, comp=1, tests=None,
marker='o', colors=None, colormap='viridis',
font_family='serif', legend=True, alpha=.5,
legend_fontsize=12, figsize=None):
"""
Plots influences (Cook's distances) of the embeddings, associated with
the tests underlying the KernelAOVResults object, against their
projections on the chosen discriminant axis. Produces separate
subplots for each test.
Parameters
----------
trunc : int, optional
Truncation of the resdual covariance operator used for the Cook's
distance calculation.
comp : int, optional
Component of the projections, i.e. the embeddings are projected on
the comp-th eigenfunction.
tests : list of strings
List containing a list of tests to plot, out of all the tests in
the KernelAOVResults object (particularly useful with the by_level
testing option).
marker : str or dict
Marker styles of the embedding projections. The default is 'o'.
If string, same marker for all mean embeddings. To pass different
values for different tests and levels, pass a dictionary with keys
corresponding to test names, and values that are dictionaries with
keys corresponding to the levels of the test and values that are
marker style strings.
colors : None or dict
Colors of the mean embedding projections. If None (default), colors
are chosen from the specified colormap. To customize, pass a
dictionary with keys corresponding to test names, and values that
are dictionaries with keys corresponding to the levels of the test
and values that are colors.
colormap : str, optional
The name of a matplotlib colormap to be used for different factor
levels (if colors are not provided). The default is 'viridis'.
alpha : float, optional
The alpha blending value, between 0 (transparent) and 1 (opaque).
The default is 0.5.
legend_fontsize : int, optional
Legend font size. The default is 15.
font_family : str, optional
Legend and labels' font family name accepted by matplotlib
(e.g., 'serif', 'sans-serif', 'monospace', 'fantasy' or 'cursive'),
the default is 'serif'.
figsize : tuple, optional
The size of the figure. If not specified, is set to
(8 * nb_factors, 6).
Returns
-------
fig : matplotlib.figure.Figure
A Figure object of the plot.
axs : numpy.ndarray of matplotlib.axes._axes.Axes
An Axes object of the plot.
"""
warnings.simplefilter("always")
tests = self.cook_distances.keys() if tests is None else tests
nb_tests = len(tests)
rc('font', **{'family': font_family})
figsize = figsize if figsize is not None else (8 * nb_tests, 6)
fig, axs = plt.subplots(ncols=nb_tests, figsize=figsize)
for j, test in enumerate(tests):
ax = axs if nb_tests == 1 else axs[j]
t1 = min(trunc, len(self.stats[test].index))
if t1 != trunc:
warnings.warn(f'In {test}: trunc={t1} will be plotted, since {trunc} ' +
'is larger than the maximal number of truncations.')
t2 = min(comp, len(self.projections[test].columns) - 1)
if t2 != comp:
warnings.warn(f'In {test}: comp={t2} will be plotted, since {comp} ' +
'is larger than the maximal number of components.')
cook_j = self.cook_distances[test]
proj_j = self.projections[test]
test_lvls = cook_j[test].unique()
test_lvls = test_lvls[test_lvls != 'NA'] # extract relevant observations
nb_lvls = len(test_lvls)
# Create dictionaries for colors and markers:
if colors is None:
cmap = colormaps[colormap]
cmap_colors = cmap(np.linspace(0.1, 0.9, nb_lvls))
color = {test_lvls[i] : cmap_colors[i] for i in range(nb_lvls)}
elif type(colors) is dict:
color = colors[test]
if type(marker) is str:
markers = {test_lvls[i] : marker for i in range(nb_lvls)}
elif type(marker) is dict:
markers = marker[test]
no_lvl_info = cook_j[test].isnull().all()
for i, test_lvl in enumerate(test_lvls):
if no_lvl_info:
lvl_cook = cook_j[t1]
lvl_proj = proj_j[t2]
else:
lvl_cook = cook_j[cook_j[test] == test_lvl][t1]
lvl_proj = proj_j[proj_j[test] == test_lvl][t2]
ax.scatter(lvl_proj, lvl_cook, color=color[test_lvl],
marker=markers[test_lvl], alpha=alpha, label=test_lvl)
if legend and not no_lvl_info:
ax.legend(bbox_to_anchor=(1.01, 0.5), loc='center left',
fontsize=legend_fontsize)
ax.set_title(test, fontsize=18)
plt.tight_layout()
fig.suptitle(f"Cook's distances (trunc. {t1}) against projections (comp. {t2})",
fontsize=25, y=1.05)
plt.draw()
return fig, axs
[docs]
def get_projections(self, n_comp=None, factor=None, hypothesis=None):
"""
Creates a pandas.DataFrame with the discriminant axis projections
associated with the test for a given number of components, factor
(for all tests of type 'pairwise' and 'one-vs-all') and
hypothesis (in the by-level cases and for user-specified tests).
Parameters
----------
n_comp : int or None
Number of discriminant axis components for which to return the
results. If None (default), projections on all available components
are returned.
factor : str or None
None by default (acceptable for user-specified tests) only. A factor
has to be specified in all other cases, then returns the results of
tests on comparisons related with the chosen factor. Factor names
are keys of the `_factor_info` attribute of the AOV class.
hypothesis : str or None
None by default, which is acceptable if the test is global or of type
'one-vs-all'. In the by-level pairwise or custom test cases a
hypothesis has to be specified. The list of possible hypothesis
names is accessible through the `hypotheses` attribute of the
KernelAOVResults class (first element of each tuple).
Returns
-------
sum_df : pandas.DataFrame
A data frame with the summary of projections.
"""
predef = (self.hypothesis_type in ['pairwise', 'one-vs-all'] or
(self.hypothesis_type is None and not self.by_level))
### Pre-defined test cases (compatible with the OneHot encoding):
if predef:
if factor is None:
raise ValueError("Set the factor parameter.")
nb_factors = factor.count(':') + 1
factor_cols = [f'factor_{f + 1}' for f in range(nb_factors)]
if not self.by_level:
t_cols = self.projections[factor].columns[:-1]
if n_comp is not None:
t_cols = t_cols[:n_comp]
proj = self.projections[factor]
sum_df = pd.DataFrame(columns=factor_cols + [f'proj_{i}' for i in t_cols],
index=proj.index, dtype=object)
for i in t_cols:
sum_df[f'proj_{i}'] = proj[i]
factors_split = proj[factor].str.split(':')
for i, fct in enumerate(factor_cols):
sum_df[fct] = factors_split.str.get(i).str.split('[').str.get(1).str.split(']').str.get(0)
elif self.hypothesis_type == 'one-vs-all':
# Only one direction in the by-level case:
t_cols = [1, ]
sum_df = pd.DataFrame(columns=factor_cols + [f'proj_{i}' for i in t_cols],
index=list(self.projections.values())[0].index,
dtype=object)
for i, (hyp, _) in enumerate(self.hypotheses):
is_factor_hyp = ((nb_factors == 1 and factor in hyp and ':' not in hyp)
or (nb_factors > 1 and np.all([f in hyp for f in factor.split(':')])))
if is_factor_hyp:
proj = self.projections[hyp]
where_hyp = proj[hyp] == hyp.split(' = ')[0]
for i in t_cols:
sum_df.loc[where_hyp, f'proj_{i}'] = proj.loc[where_hyp, i]
levels = re.findall(r"\[(.*?)\]", hyp)
for i, fct in enumerate(factor_cols):
sum_df.loc[where_hyp, fct] = levels[i]
else:
if hypothesis is None:
error_message = "Set the hypothesis parameter (necessary in "
error_message += "the by-level pairwise case)."
raise ValueError(error_message)
else:
t_cols = [1, ]
proj = self.projections[hypothesis]
where_hyp = proj[hypothesis] != 'NA'
sum_df = pd.DataFrame(columns=factor_cols + [f'proj_{i}' for i in t_cols],
index=proj.index,
dtype=object)
for i in t_cols:
sum_df.loc[where_hyp, f'proj_{i}'] = proj.loc[where_hyp, i]
pw_groups = hypothesis.split(' = ')
where_left = proj[hypothesis] == pw_groups[0]
where_right = proj[hypothesis] == pw_groups[1]
levels = re.findall(r"\[(.*?)\]", hypothesis)
for i, fct in enumerate(factor_cols):
sum_df.loc[where_left, fct] = levels[i]
sum_df.loc[where_right, fct] = levels[nb_factors + i]
sum_df.dropna(axis=1, how='all', inplace=True)
sum_df.dropna(axis=0, how='all', inplace=True)
### Other cases (custom contrasts, other coding schemes...):
else:
if hypothesis is None:
error_message = "Set the hypothesis parameter."
raise ValueError(error_message)
if self.by_level:
t_cols = [1, ]
else:
t_cols = self.projections[hypothesis].columns
if n_comp is not None:
t_cols = t_cols[:n_comp]
proj = self.projections[hypothesis]
sum_df = pd.DataFrame(columns=[f'proj_{i}' for i in t_cols],
index=proj.index,
dtype=object)
for i in t_cols:
sum_df.loc[:, f'proj_{i}'] = proj.loc[:, i]
sum_df.dropna(axis=1, how='all', inplace=True)
return sum_df
[docs]
def get_cook(self, n_trunc=None, factor=None, hypothesis=None):
"""
Creates a pandas.DataFrame with Cook's distances and their
p-values associated with the test for given truncations,
factor (for all tests of type 'pairwise' and 'one-vs-all') and
hypothesis (in the by-level cases and for user-specified tests).
Parameters
----------
n_trunc : int or None
Number of truncations for which to return the results. If None
(default), Cook's distances are returned for all available truncations.
factor : str or None
None by default (acceptable for user-specified tests) only. A factor
has to be specified in all other cases, then returns the results of
tests on comparisons related with the chosen factor. Factor names are
keys of the `_factor_info` attribute of the AOV class.
hypothesis : str or None
None by default, which is acceptable if the test is global or of type
'one-vs-all'. In the by-level pairwise or custom test cases a hypothesis
has to be specified. The list of possible hypothesis names is accessible
through the `hypotheses` attribute of the KernelAOVResults class (first
element of each tuple).
Returns
-------
sum_df : pandas.DataFrame
A data frame with the summary of projections.
"""
predef = (self.hypothesis_type in ['pairwise', 'one-vs-all'] or
(self.hypothesis_type is None and not self.by_level))
### Pre-defined test cases (compatible with the OneHot encoding):
if predef:
if factor is None:
raise ValueError("Set the factor parameter.")
nb_factors = factor.count(':') + 1
factor_cols = [f'factor_{f + 1}' for f in range(nb_factors)]
if not self.by_level:
t_cols = self.cook_distances[factor].columns[:-1]
if n_trunc is not None:
t_cols = t_cols[:n_trunc]
cook = self.cook_distances[factor]
cook_pval = self.cook_pvalues[factor]
sum_df = pd.DataFrame(columns=(factor_cols+ [f'cook_{i}' for i in t_cols]
+ [f'cook_pval_{i}' for i in t_cols]),
index=cook.index, dtype=object)
for i in t_cols:
sum_df[f'cook_{i}'] = cook[i]
sum_df[f'cook_pval_{i}'] = cook_pval[i]
factors_split = cook[factor].str.split(':')
for i, fct in enumerate(factor_cols):
sum_df[fct] = factors_split.str.get(i).str.split('[').str.get(1).str.split(']').str.get(0)
elif self.hypothesis_type == 'one-vs-all':
t_cols = list(self.cook_distances.values())[0].columns[:-1]
if n_trunc is not None:
t_cols = t_cols[:n_trunc]
sum_df = pd.DataFrame(columns=(factor_cols + [f'cook_{i}' for i in t_cols]
+ [f'cook_pval_{i}' for i in t_cols]),
index=list(self.cook_distances.values())[0].index,
dtype=object)
for i, (hyp, _) in enumerate(self.hypotheses):
is_factor_hyp = ((nb_factors == 1 and factor in hyp and ':' not in hyp)
or (nb_factors > 1 and np.all([f in hyp for f in factor.split(':')])))
if is_factor_hyp:
cook = self.cook_distances[hyp]
cook_pval = self.cook_pvalues[hyp]
where_hyp = cook[hyp] == hyp.split(' = ')[0]
for i in t_cols:
sum_df.loc[where_hyp, f'cook_{i}'] = cook.loc[where_hyp, i]
sum_df.loc[where_hyp, f'cook_pval_{i}'] = cook_pval.loc[where_hyp, i]
levels = re.findall(r"\[(.*?)\]", hyp)
for i, fct in enumerate(factor_cols):
sum_df.loc[where_hyp, fct] = levels[i]
else:
if hypothesis is None:
error_message = "Set the hypothesis parameter (necessary in "
error_message += "the by-level pairwise case)."
raise ValueError(error_message)
else:
t_cols = list(self.cook_distances.values())[0].columns[:-1]
if n_trunc is not None:
t_cols = t_cols[:n_trunc]
cook = self.cook_distances[hypothesis]
cook_pval = self.cook_pvalues[hypothesis]
where_hyp = cook[hypothesis] != 'NA'
sum_df = pd.DataFrame(columns=(factor_cols + [f'cook_{i}' for i in t_cols]
+ [f'cook_pval_{i}' for i in t_cols]),
index=cook.index,
dtype=object)
for i in t_cols:
sum_df.loc[where_hyp, f'cook_{i}'] = cook.loc[where_hyp, i]
sum_df.loc[where_hyp, f'cook_pval_{i}'] = cook_pval.loc[where_hyp, i]
pw_groups = hypothesis.split(' = ')
where_left = cook[hypothesis] == pw_groups[0]
where_right = cook[hypothesis] == pw_groups[1]
levels = re.findall(r"\[(.*?)\]", hypothesis)
for i, fct in enumerate(factor_cols):
sum_df.loc[where_left, fct] = levels[i]
sum_df.loc[where_right, fct] = levels[nb_factors + i]
sum_df.dropna(axis=1, how='all', inplace=True)
sum_df.dropna(axis=0, how='all', inplace=True)
### Other cases (custom contrasts, other coding schemes...):
else:
if hypothesis is None:
error_message = "Set the hypothesis parameter."
raise ValueError(error_message)
t_cols = self.stats[hypothesis].index
if n_trunc is not None:
t_cols = t_cols[:n_trunc]
cook = self.cook_distances[hypothesis]
cook_pval = self.cook_pvalues[hypothesis]
sum_df = pd.DataFrame(columns=([f'cook_{i}' for i in t_cols]
+ [f'cook_pval_{i}' for i in t_cols]),
index=cook.index,
dtype=object)
for i in t_cols:
sum_df[f'cook_{i}'] = cook[i]
sum_df[f'cook_pval_{i}'] = cook_pval[i]
sum_df.dropna(axis=1, how='all', inplace=True)
return sum_df