diff --git a/.gitignore b/.gitignore index ae8d3eb1..fbe06e35 100644 --- a/.gitignore +++ b/.gitignore @@ -125,3 +125,6 @@ tex.cache/ testtt.py real.py 0to2_beforeduringafter.csv + +TEST.ipynb +TrhCsCh.csv diff --git a/dabest/_api.py b/dabest/_api.py index fa2433f0..efb2ee98 100644 --- a/dabest/_api.py +++ b/dabest/_api.py @@ -4,8 +4,10 @@ # Email : joseshowh@gmail.com -def load(data, idx, x=None, y=None, paired=None, id_col=None, - ci=95, resamples=5000, random_seed=12345, proportional=False): +def load(data, idx=None, x=None, y=None, paired=None, id_col=None, + ci=95, resamples=5000, random_seed=12345, proportional=False, + delta2 = False, experiment = None, experiment_label = None, + x1_level = None): ''' Loads data in preparation for estimation statistics. @@ -18,7 +20,10 @@ def load(data, idx, x=None, y=None, paired=None, id_col=None, List of column names (if 'x' is not supplied) or of category names (if 'x' is supplied). This can be expressed as a tuple of tuples, with each individual tuple producing its own contrast plot - x : string, default None + x : string or list, default None + Column name(s) of the independent variable. This can be expressed as + a list of 2 elements if and only if 'delta2' is True; otherwise it + can only be a string. y : string, default None Column names for data to be plotted on the x-axis and y-axis. paired : string, default None @@ -37,6 +42,19 @@ def load(data, idx, x=None, y=None, paired=None, id_col=None, reported are replicable. proportional : boolean, default False. TO INCLUDE MORE DESCRIPTION ABOUT DATA FORMAT + delta2 : boolean, default False + Indicator of delta-delta experiment + experiment : String, default None + The name of the column of the dataframe which contains the label of + experiments + experiment_lab : list, default None + A list of String to specify the order of subplots for delta-delta plots. + This can be expressed as a list of 2 elements if and only if 'delta2' + is True; otherwise it can only be a string. + x1_level : list, default None + A list of String to specify the order of subplots for delta-delta plots. + This can be expressed as a list of 2 elements if and only if 'delta2' + is True; otherwise it can only be a string. Returns ------- @@ -65,4 +83,4 @@ def load(data, idx, x=None, y=None, paired=None, id_col=None, ''' from ._classes import Dabest - return Dabest(data, idx, x, y, paired, id_col, ci, resamples, random_seed, proportional) + return Dabest(data, idx, x, y, paired, id_col, ci, resamples, random_seed, proportional, delta2, experiment, experiment_label, x1_level) diff --git a/dabest/_classes.py b/dabest/_classes.py index 8395f5fa..15bc1530 100644 --- a/dabest/_classes.py +++ b/dabest/_classes.py @@ -9,8 +9,9 @@ class Dabest(object): Class for estimation statistics and plots. """ - def __init__(self, data, idx, x, y, paired, id_col, ci, resamples, - random_seed, proportional): + def __init__(self, data, idx, x, y, paired, id_col, ci, + resamples, random_seed, proportional, delta2, + experiment, experiment_label, x1_level): """ Parses and stores pandas DataFrames in preparation for estimation @@ -23,9 +24,10 @@ def __init__(self, data, idx, x, y, paired, id_col, ci, resamples, import pandas as pd import seaborn as sns + self.__delta2 = delta2 + self.__experiment = experiment self.__ci = ci self.__data = data - self.__idx = idx self.__id_col = id_col self.__is_paired = paired self.__resamples = resamples @@ -39,6 +41,98 @@ def __init__(self, data, idx, x, y, paired, id_col, ci, resamples, + # check if this is a 2x2 ANOVA case and x & y are valid columns + # create experiment_label and x1_level + if delta2: + # idx should not be specified + if idx: + err0 = '`idx` should not be specified when `delta2` is True.'.format(len(x)) + raise ValueError(err0) + + # check if x is valid + if len(x) != 2: + err0 = '`delta2` is True but the number of variables indicated by `x` is {}.'.format(len(x)) + raise ValueError(err0) + else: + for i in x: + if i not in data_in.columns: + err = '{0} is a column in `data`. Please check.'.format(i) + raise IndexError(err) + + # check if y is valid + if not y: + err0 = '`delta2` is True but `y` is not indicated.' + raise ValueError(err0) + elif y not in data_in.columns: + err = '{0} is not a column in `data`. Please check.'.format(y) + raise IndexError(err) + + # check if experiment is valid + if experiment not in data_in.columns: + err = '{0} is not a column in `data`. Please check.'.format(experiment) + raise IndexError(err) + + # check if experiment_label is valid and create experiment when needed + if experiment_label: + if len(experiment_label) != 2: + err0 = '`experiment_label` does not have a length of 2.' + raise ValueError(err0) + else: + for i in experiment_label: + if i not in data_in[experiment].unique(): + err = '{0} is an element in the column `{1}` of `data`. Please check.'.format(i, experiment) + raise IndexError(err) + else: + experiment_label = data_in[experiment].unique() + + # check if x1_level is valid + if x1_level: + if len(x1_level) != 2: + err0 = '`x1_level` does not have a length of 2.' + raise ValueError(err0) + else: + for i in x1_level: + if i not in data_in[x[0]].unique(): + err = '{0} is an element in the column `{1}` of `data`. Please check.'.format(i, experiment) + raise IndexError(err) + + else: + x1_level = data_in[x[0]].unique() + self.__experiment_label = experiment_label + self.__x1_level = x1_level + + + # check if idx is specified + if not delta2 and not idx: + err = '`idx` is not a column in `data`. Please check.' + raise IndexError(err) + + + # create new x & idx and record the second variable if this is a valid 2x2 ANOVA case + if delta2: + # add a new column which is a combination of experiment and the first variable + new_col_name = experiment+x[0] + while new_col_name in data_in.columns: + new_col_name += "_" + data_in[new_col_name] = data_in[x[0]].apply(lambda x: str(x)) + " " + data_in[experiment].apply(lambda x: str(x)) + + #create idx + idx = [] + for i in experiment_label: + temp = [] + for j in x1_level: + temp.append(j + " " + i) + idx.append(temp) + self.__idx = idx + self.__x1 = x[0] + self.__x2 = x[1] + # record the second variable and create idx + x = new_col_name + else: + self.__idx = idx + self.__x1 = None + self.__x2 = None + # Determine the kind of estimation plot we need to produce. if all([isinstance(i, str) for i in idx]): # flatten out idx. @@ -200,7 +294,11 @@ def __init__(self, data, idx, x, y, paired, id_col, ci, resamples, EffectSizeDataFrame_kwargs = dict(ci=ci, is_paired=paired, random_seed=random_seed, resamples=resamples, - proportional=proportional) + proportional=proportional, + delta2=delta2, + experiment_label=self.__experiment_label, + x1_level=self.__x1_level, + x2=self.__x2) self.__mean_diff = EffectSizeDataFrame(self, "mean_diff", **EffectSizeDataFrame_kwargs) @@ -269,6 +367,9 @@ def __repr__(self): for ix, test_name in enumerate(current_tuple[1:]): comparisons.append("{} minus {}".format(test_name, control_name)) + if self.__delta2: + comparisons.append("{} minus {} (only for mean difference)".format(self.__experiment_label[1], self.__experiment_label[0])) + for j, g in enumerate(comparisons): out.append("{}. {}".format(j+1, g)) @@ -495,7 +596,6 @@ def cliffs_delta(self): return self.__cliffs_delta - @property def data(self): """ @@ -503,6 +603,7 @@ def data(self): """ return self.__data + @property def idx(self): """ @@ -510,13 +611,36 @@ def idx(self): """ return self.__idx - # Removed due to the deprecation fo `is_paired` - #@property - #def is_paired(self): - # """ - # Returns True if the dataset was declared as paired to `dabest.load()`. - # """ - # return self.__is_paired + + @property + def x1(self): + return self.__x1 + + + @property + def x1_level(self): + return self.__x1_level + + + @property + def x2(self): + return self.__x2 + + + @property + def experiment(self): + return self.__experiment + + + @property + def experiment_label(self): + return self.__experiment_label + + + @property + def delta2(self): + return self.__delta2 + @property def is_paired(self): @@ -525,6 +649,7 @@ def is_paired(self): """ return self.__is_paired + @property def id_col(self): """ @@ -532,6 +657,7 @@ def id_col(self): """ return self.__id_col + @property def ci(self): """ @@ -539,6 +665,7 @@ def ci(self): """ return self.__ci + @property def resamples(self): """ @@ -546,6 +673,7 @@ def resamples(self): """ return self.__resamples + @property def random_seed(self): """ @@ -562,6 +690,7 @@ def x(self): """ return self.__x + @property def y(self): """ @@ -569,6 +698,7 @@ def y(self): """ return self.__y + @property def _xvar(self): """ @@ -576,6 +706,7 @@ def _xvar(self): """ return self.__xvar + @property def _yvar(self): """ @@ -583,6 +714,7 @@ def _yvar(self): """ return self.__yvar + @property def _plot_data(self): """ @@ -599,6 +731,7 @@ def proportional(self): """ return self.__proportional + @property def _all_plot_groups(self): """ @@ -621,7 +754,8 @@ def __init__(self, control, test, effect_size, is_paired=None, ci=95, resamples=5000, permutation_count=5000, - random_seed=12345): + random_seed=12345, + delta2=False): """ Compute the effect size between two groups. @@ -648,6 +782,9 @@ def __init__(self, control, test, effect_size, `random_seed` is used to seed the random number generator during bootstrap resampling. This ensures that the confidence intervals reported are replicable. + delta2 : boolean, default False + Indicate if the control and test data are boostrap deltas that can be + used to calculate delta-delta. Returns @@ -785,21 +922,25 @@ def __init__(self, control, test, effect_size, self.__random_seed = random_seed self.__ci = ci self.__alpha = ci2g._compute_alpha_from_ci(ci) - + self.__delta2 = delta2 self.__difference = es.two_group_difference( control, test, is_paired, effect_size) - + self.__jackknives = ci2g.compute_meandiff_jackknife( control, test, is_paired, effect_size) self.__acceleration_value = ci2g._calc_accel(self.__jackknives) - bootstraps = ci2g.compute_bootstrapped_diff( + if not delta2: + bootstraps = ci2g.compute_bootstrapped_diff( control, test, is_paired, effect_size, resamples, random_seed) - self.__bootstraps = npsort(bootstraps) + self.__bootstraps = bootstraps + else: + self.__bootstraps = self.__test-self.__control + sorted_bootstraps = npsort(self.__bootstraps) # Added in v0.2.6. # Raises a UserWarning if there are any infiinities in the bootstraps. num_infinities = len(self.__bootstraps[isinf(self.__bootstraps)]) @@ -825,8 +966,8 @@ def __init__(self, control, test, effect_size, self.__bca_interval_idx = (bca_idx_low, bca_idx_high) if ~isnan(bca_idx_low) and ~isnan(bca_idx_high): - self.__bca_low = self.__bootstraps[bca_idx_low] - self.__bca_high = self.__bootstraps[bca_idx_high] + self.__bca_low = sorted_bootstraps[bca_idx_low] + self.__bca_high = sorted_bootstraps[bca_idx_high] err1 = "The $lim_type limit of the interval" err2 = "was in the $loc 10 values." @@ -858,14 +999,19 @@ def __init__(self, control, test, effect_size, self.__bca_high = self.__difference warnings.warn(err_temp.substitute(lim_type="upper"), stacklevel=0) + if not self.__delta2: + # Compute percentile intervals. + pct_idx_low = int((self.__alpha/2) * resamples) + pct_idx_high = int((1-(self.__alpha/2)) * resamples) - # Compute percentile intervals. - pct_idx_low = int((self.__alpha/2) * resamples) - pct_idx_high = int((1-(self.__alpha/2)) * resamples) - - self.__pct_interval_idx = (pct_idx_low, pct_idx_high) - self.__pct_low = self.__bootstraps[pct_idx_low] - self.__pct_high = self.__bootstraps[pct_idx_high] + self.__pct_interval_idx = (pct_idx_low, pct_idx_high) + self.__pct_low = sorted_bootstraps[pct_idx_low] + self.__pct_high = sorted_bootstraps[pct_idx_high] + + else: + self.__pct_interval_idx = None + self.__pct_low = None + self.__pct_high = None # Perform statistical tests. @@ -1001,7 +1147,10 @@ def __repr__(self, show_resample_count=True, define_pval=True, sigfig=3): "es" : self.__EFFECT_SIZE_DICT[self.__effect_size], "paired_status": PAIRED_STATUS[str(self.__is_paired)]} - out1 = "The {paired_status} {es} {rm_status}".format(**first_line) + if self.__delta2: + out1 = "The delta-delta " + else: + out1 = "The {paired_status} {es} {rm_status}".format(**first_line) base_string_fmt = "{:." + str(sigfig) + "}" if "." in str(self.__ci): @@ -1038,7 +1187,9 @@ def __repr__(self, show_resample_count=True, define_pval=True, sigfig=3): # pval_rounded) - pvalue = "The p-value of the two-sided permutation t-test is {}. ".format(pval_rounded) + p1 = "The p-value of the two-sided permutation t-test is {}, ".format(pval_rounded) + p2 = "calculated for legacy purposes only. " + pvalue = p1 + p2 bs1 = "{} bootstrap samples were taken; ".format(self.__resamples) bs2 = "the confidence interval is bias-corrected and accelerated." @@ -1075,6 +1226,9 @@ def to_dict(self): return out + @property + def delta2(self): + return self.__delta2 @property @@ -1388,7 +1542,9 @@ def __init__(self, dabest, effect_size, is_paired, ci=95, proportional=False, resamples=5000, permutation_count=5000, - random_seed=12345): + random_seed=12345, + x1_level=None, x2=None, + delta2=False, experiment_label=None): """ Parses the data from a Dabest object, enabling plotting and printing capability for the effect size of interest. @@ -1402,6 +1558,10 @@ def __init__(self, dabest, effect_size, self.__permutation_count = permutation_count self.__random_seed = random_seed self.__proportional = proportional + self.__x1_level = x1_level + self.__experiment_label = experiment_label + self.__x2 = x2 + self.__delta2 = delta2 def __pre_calc(self): @@ -1435,21 +1595,22 @@ def __pre_calc(self): self.__permutation_count, self.__random_seed) r_dict = result.to_dict() - r_dict["control"] = cname r_dict["test"] = tname r_dict["control_N"] = int(len(control)) r_dict["test_N"] = int(len(test)) - out.append(r_dict) - if j == len(idx)-1 and ix == len(current_tuple)-2: - resamp_count = True - def_pval = True + if not self.__delta2: + resamp_count = True + def_pval = True + else: + resamp_count = False + def_pval = False else: resamp_count = False def_pval = False - + text_repr = result.__repr__(show_resample_count=resamp_count, define_pval=def_pval) @@ -1458,6 +1619,34 @@ def __pre_calc(self): reprs.append(text_repr) + if self.__delta2 and self.__effect_size == "mean_diff": + delta = TwoGroupsEffectSize(out[0]["bootstraps"], + out[1]["bootstraps"], + self.__effect_size, + "baseline", + self.__ci, + self.__resamples, + self.__permutation_count, + self.__random_seed, + self.__delta2 + ) + + r_dict = delta.to_dict() + r_dict["control"] = self.__experiment_label[0] + r_dict["test"] = self.__experiment_label[1] + r_dict["control_N"] = self.__resamples + r_dict["test_N"] = self.__resamples + out.append(r_dict) + resamp_count = True + def_pval = True + text_repr = delta.__repr__(show_resample_count=resamp_count, + define_pval=def_pval) + to_replace = "between {} and {} is".format(self.__experiment_label[0], self.__experiment_label[1]) + text_repr = text_repr.replace("is", to_replace, 1) + reprs.append(text_repr) + else: + self.__delta2 = False + varname = get_varname(self.__dabest_obj) lastline = "To get the results of all valid statistical tests, " +\ "use `{}.{}.statistical_tests`".format(varname, self.__effect_size) @@ -1504,7 +1693,6 @@ def __pre_calc(self): self.__results = out_.reindex(columns=columns_in_order) self.__results.dropna(axis="columns", how="all", inplace=True) - @@ -1526,6 +1714,7 @@ def __calc_lqrt(self): dat = db_obj._plot_data xvar = db_obj._xvar yvar = db_obj._yvar + delta2 = self.__delta2 out = [] @@ -1573,7 +1762,18 @@ def __calc_lqrt(self): "pvalue_lqrt_unequal_var" : lqrt_unequal_var_result.pvalue, "statistic_lqrt_unequal_var" : lqrt_unequal_var_result.statistic, }) - + if delta2: + lqrt_result = lqrt.lqrtest_rel(self.results["bootstraps"][0], + self.results["bootstraps"][1], + random_state=rnd_seed) + + out.append({"control": self.__experiment_label[0], + "test": self.__experiment_label[1], + "control_N": self.__resamples, + "test_N": self.__resamples, + "pvalue_paired_lqrt": lqrt_result.pvalue, + "statistic_paired_lqrt": lqrt_result.statistic + }) self.__lqrt_results = pd.DataFrame(out) @@ -1581,14 +1781,15 @@ def plot(self, color_col=None, raw_marker_size=6, es_marker_size=9, - swarm_label=None, barchart_label=None, contrast_label=None, - swarm_ylim=None, barchart_ylim=None, contrast_ylim=None, + swarm_label=None, barchart_label=None, contrast_label=None, delta_label=None, + swarm_ylim=None, barchart_ylim=None, contrast_ylim=None, custom_palette=None, swarm_desat=0.5, barchart_desat=0.5, halfviolin_desat=1, halfviolin_alpha=0.8, float_contrast=True, show_pairs=True, + show_delta2=True, group_summaries=None, group_summaries_offset=0.1, @@ -1618,15 +1819,16 @@ def plot(self, color_col=None, es_marker_size : float, default 9 The size (in points) of the effect size points on the difference axes. - swarm_label, contrast_label : strings, default None + swarm_label, contrast_label, delta_label : strings, default None Set labels for the y-axis of the swarmplot and the contrast plot, respectively. If `swarm_label` is not specified, it defaults to "value", unless a column name was passed to `y`. If `contrast_label` is not specified, it defaults to the effect size - being plotted. + being plotted. If `delta_label` is not specifed, it defaults to + "delta - delta" swarm_ylim, contrast_ylim : tuples, default None - The desired y-limits of the raw data (swarmplot) axes and the - difference axes respectively, as a tuple. These will be autoscaled + The desired y-limits of the raw data (swarmplot) axes, the + difference axes and the delta-delta axes respectively, as a tuple. These will be autoscaled to sensible values if they are not specified. custom_palette : dict, list, or matplotlib color palette, default None This keyword accepts a dictionary with {'group':'color'} pairings, @@ -1780,10 +1982,12 @@ def plot(self, color_col=None, from .plotter import EffectSizeDataFramePlotter, ProportionalDataFramePlotter - if hasattr(self, "results") is False: self.__pre_calc() + if self.__delta2: + color_col = self.__x2 + if self.__proportional: raw_marker_size = 0.01 @@ -1866,6 +2070,26 @@ def ci(self): """ return self.__ci + @property + def x1_level(self): + return self.__x1_level + + + @property + def x2(self): + return self.__x2 + + + @property + def experiment_label(self): + return self.__experiment_label + + + @property + def delta2(self): + return self.__delta2 + + @property def resamples(self): """ diff --git a/dabest/plotter.py b/dabest/plotter.py index 746161bc..62a5561d 100644 --- a/dabest/plotter.py +++ b/dabest/plotter.py @@ -14,12 +14,13 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): **plot_kwargs: color_col=None raw_marker_size=6, es_marker_size=9, - swarm_label=None, contrast_label=None, + swarm_label=None, contrast_label=None, delta_label=None, swarm_ylim=None, contrast_ylim=None, custom_palette=None, swarm_desat=0.5, halfviolin_desat=1, halfviolin_alpha=0.8, float_contrast=True, show_pairs=True, + show_delta2=True, group_summaries=None, group_summaries_offset=0.1, fig_size=None, @@ -59,23 +60,30 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): ytick_color = plt.rcParams["ytick.color"] axes_facecolor = plt.rcParams['axes.facecolor'] - dabest_obj = EffectSizeDataFrame.dabest_obj - plot_data = EffectSizeDataFrame._plot_data - xvar = EffectSizeDataFrame.xvar - yvar = EffectSizeDataFrame.yvar - is_paired = EffectSizeDataFrame.is_paired - + dabest_obj = EffectSizeDataFrame.dabest_obj + plot_data = EffectSizeDataFrame._plot_data + xvar = EffectSizeDataFrame.xvar + yvar = EffectSizeDataFrame.yvar + is_paired = EffectSizeDataFrame.is_paired + delta2 = EffectSizeDataFrame.delta2 + effect_size = EffectSizeDataFrame.effect_size all_plot_groups = dabest_obj._all_plot_groups idx = dabest_obj.idx + if effect_size != "mean_diff" or not delta2: + show_delta2 = False + else: + show_delta2 = plot_kwargs["show_delta2"] + + # Disable Gardner-Altman plotting if any of the idxs comprise of more than - # two groups. + # two groups or if it is a delta-delta plot. float_contrast = plot_kwargs["float_contrast"] effect_size_type = EffectSizeDataFrame.effect_size if len(idx) > 1 or len(idx[0]) > 2: @@ -84,11 +92,16 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): if effect_size_type in ['cliffs_delta']: float_contrast = False + if show_delta2: + float_contrast = False + if not is_paired: show_pairs = False else: show_pairs = plot_kwargs["show_pairs"] + + # Set default kwargs first, then merge with user-dictated ones. default_swarmplot_kwargs = {'size': plot_kwargs["raw_marker_size"]} if plot_kwargs["swarmplot_kwargs"] is None: @@ -109,6 +122,7 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): plot_kwargs["violinplot_kwargs"]) + # slopegraph kwargs. default_slopegraph_kwargs = {'lw':1, 'alpha':0.5} if plot_kwargs["slopegraph_kwargs"] is None: @@ -119,7 +133,6 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): - # Zero reference-line kwargs. default_reflines_kwargs = {'linestyle':'solid', 'linewidth':0.75, 'zorder': 2, @@ -217,10 +230,14 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): plot_palette_contrast = dict(zip(names, contrast_colors)) + # Infer the figsize. fig_size = plot_kwargs["fig_size"] if fig_size is None: all_groups_count = np.sum([len(i) for i in dabest_obj.idx]) + # Increase the width for delta-delta graph + if show_delta2: + all_groups_count += 2 if is_paired and show_pairs is True: frac = 0.75 else: @@ -236,7 +253,6 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): fig_size = (width_inches, height_inches) - # Initialise the figure. # sns.set(context="talk", style='ticks') init_fig_kwargs = dict(figsize=fig_size, dpi=plot_kwargs["dpi"], @@ -295,7 +311,6 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): gridspec_kw={"width_ratios": width_ratios_ga, "wspace": 0}, **init_fig_kwargs) - else: fig, axx = plt.subplots(nrows=2, gridspec_kw={"hspace": 0.3}, @@ -306,12 +321,12 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): contrast_ax_ylim_low = list() contrast_ax_ylim_high = list() contrast_ax_ylim_tickintervals = list() - rawdata_axes = axx[0] - contrast_axes = axx[1] - + rawdata_axes = axx[0] + contrast_axes = axx[1] rawdata_axes.set_frame_on(False) contrast_axes.set_frame_on(False) + redraw_axes_kwargs = {'colors' : ytick_color, 'facecolors' : ytick_color, 'lw' : 1, @@ -359,20 +374,20 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): current_pair = pivoted_plot_data else: current_pair = pivoted_plot_data[yvar] - grp_count = len(current_tuple) # Iterate through the data for the current tuple. for ID, observation in current_pair.iterrows(): x_points = [t for t in range(x_start,x_start+grp_count)] y_points = observation.tolist() - + if color_col is None: slopegraph_kwargs['color'] = ytick_color else: color_key = pivoted_plot_data[color_col, current_tuple[0]].loc[ID] - slopegraph_kwargs['color'] = plot_palette_raw[color_key] - slopegraph_kwargs['label'] = color_key + if not pd.isna(color_key): + slopegraph_kwargs['color'] = plot_palette_raw[color_key] + slopegraph_kwargs['label'] = color_key rawdata_axes.plot(x_points, y_points, **slopegraph_kwargs) x_start = x_start + grp_count @@ -450,7 +465,6 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): ticks_to_plot = np.arange(1, len(temp_all_plot_groups), 2).tolist() ticks_to_skip_contrast = np.cumsum([(len(t)-1)*2 for t in idx])[:-1].tolist() ticks_to_skip_contrast.insert(0, 0) - else: ticks_to_skip = np.cumsum([len(t) for t in idx])[:-1].tolist() ticks_to_skip.insert(0, 0) @@ -458,6 +472,12 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): # Then obtain the ticks where we have to plot the effect sizes. ticks_to_plot = [t for t in range(0, len(all_plot_groups)) if t not in ticks_to_skip] + + ticks_to_plot_violin = ticks_to_plot.copy() + ticks_to_skip_violin = ticks_to_skip.copy() + if show_delta2: + ticks_to_plot_violin.append(max(ticks_to_plot)+2) + ticks_to_skip_violin.append(max(ticks_to_skip)+2) # Plot the bootstraps, then the effect sizes and CIs. @@ -467,8 +487,8 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): results = EffectSizeDataFrame.results contrast_xtick_labels = [] - - for j, tick in enumerate(ticks_to_plot): + + for j, tick in enumerate(ticks_to_plot_violin): current_group = results.test[j] current_control = results.control[j] current_bootstrap = results.bootstraps[j] @@ -507,19 +527,35 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): current_control)) - # Make sure the contrast_axes x-lims match the rawdata_axes xlims. - contrast_axes.set_xticks(rawdata_axes.get_xticks()) + # Make sure the contrast_axes x-lims match the rawdata_axes xlims, + # and add an extra violinplot tick for delta-delta plot. + if not show_delta2: + contrast_axes.set_xticks(rawdata_axes.get_xticks()) + else: + temp = rawdata_axes.get_xticks() + temp = np.append(temp, [max(temp)+1, max(temp)+2]) + contrast_axes.set_xticks(temp) + + if show_pairs is True: max_x = contrast_axes.get_xlim()[1] rawdata_axes.set_xlim(-0.375, max_x) if float_contrast is True: contrast_axes.set_xlim(0.5, 1.5) + elif show_delta2: + # Increase the xlim of raw data by 2 + temp = rawdata_axes.get_xlim() + if show_pairs: + rawdata_axes.set_xlim(temp[0], temp[1]+0.25) + else: + rawdata_axes.set_xlim(temp[0], temp[1]+2) + contrast_axes.set_xlim(rawdata_axes.get_xlim()) else: contrast_axes.set_xlim(rawdata_axes.get_xlim()) # Properly label the contrast ticks. - for t in ticks_to_skip: + for t in ticks_to_skip_violin: contrast_xtick_labels.insert(t, "") contrast_axes.set_xticklabels(contrast_xtick_labels) @@ -693,7 +729,7 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): contrast_axes.set_ylim(low, high) else: contrast_axes.set_ylim(custom_contrast_ylim) - + # If 0 lies within the ylim of the contrast axes, # draw a zero reference line. contrast_axes_ylim = contrast_axes.get_ylim() @@ -720,8 +756,9 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): ax.set_ylim(ylim) del redraw_axes_kwargs['y'] - - rightend_ticks_contrast = np.array([(len(i)-1)*2-1 for i in idx]) + np.array(ticks_to_skip_contrast) + + temp_length = [(len(i)-1)*2-1 for i in idx] + rightend_ticks_contrast = np.array(temp_length) + np.array(ticks_to_skip_contrast) for ax in [contrast_axes]: sns.despine(ax=ax, bottom=True) @@ -741,11 +778,11 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): rightend_ticks = np.array([len(i)-1 for i in idx]) + np.array(ticks_to_skip) for ax in [rawdata_axes, contrast_axes]: sns.despine(ax=ax, bottom=True) - + ylim = ax.get_ylim() xlim = ax.get_xlim() redraw_axes_kwargs['y'] = ylim[0] - + for k, start_tick in enumerate(ticks_to_skip): end_tick = rightend_ticks[k] ax.hlines(xmin=start_tick, xmax=end_tick, @@ -754,8 +791,13 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): ax.set_ylim(ylim) del redraw_axes_kwargs['y'] - - + if show_delta2 is True: + ylim = contrast_axes.get_ylim() + redraw_axes_kwargs['y'] = ylim[0] + x_ticks = contrast_axes.get_xticks() + contrast_axes.hlines(xmin=x_ticks[-2], xmax=x_ticks[-1], + **redraw_axes_kwargs) + del redraw_axes_kwargs['y'] # Set raw axes y-label. swarm_label = plot_kwargs['swarm_label'] @@ -771,7 +813,7 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): 'hedges_g' : "Hedges' g", 'cliffs_delta' : "Cliff's delta"} default_contrast_label = contrast_label_dict[EffectSizeDataFrame.effect_size] - + if plot_kwargs['contrast_label'] is None: if is_paired: contrast_label = "paired\n{}".format(default_contrast_label) @@ -785,7 +827,6 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): if float_contrast is True: contrast_axes.yaxis.set_label_position("right") - # Set the rawdata axes labels appropriately rawdata_axes.set_ylabel(swarm_label) rawdata_axes.set_xlabel("") @@ -812,11 +853,26 @@ def EffectSizeDataFramePlotter(EffectSizeDataFrame, **plot_kwargs): **redraw_axes_kwargs) + if show_delta2 is True: + if plot_kwargs['delta_label'] is None: + delta_label = "delta - delta" + else: + delta_label = plot_kwargs['delta_label'] + delta_axes = contrast_axes.twinx() + delta_axes.set_frame_on(False) + delta_axes.set_ylabel(delta_label) + og_xlim_delta = contrast_axes.get_xlim() + og_ylim_delta = contrast_axes.get_ylim() + delta_axes.set_ylim(og_ylim_delta) + delta_axes.vlines(og_xlim_delta[1], + og_ylim_delta[0], og_ylim_delta[1], + **redraw_axes_kwargs) # Make sure no stray ticks appear! rawdata_axes.xaxis.set_ticks_position('bottom') rawdata_axes.yaxis.set_ticks_position('left') contrast_axes.xaxis.set_ticks_position('bottom') + if float_contrast is False: contrast_axes.yaxis.set_ticks_position('left') diff --git a/dabest/tests/baseline_images/test_47_cummings_unpaired_delta_delta_meandiff.png b/dabest/tests/baseline_images/test_47_cummings_unpaired_delta_delta_meandiff.png new file mode 100644 index 00000000..82d22b54 Binary files /dev/null and b/dabest/tests/baseline_images/test_47_cummings_unpaired_delta_delta_meandiff.png differ diff --git a/dabest/tests/baseline_images/test_48_cummings_sequential_delta_delta_meandiff.png b/dabest/tests/baseline_images/test_48_cummings_sequential_delta_delta_meandiff.png new file mode 100644 index 00000000..d4266e1d Binary files /dev/null and b/dabest/tests/baseline_images/test_48_cummings_sequential_delta_delta_meandiff.png differ diff --git a/dabest/tests/baseline_images/test_49_cummings_baseline_delta_delta_meandiff.png b/dabest/tests/baseline_images/test_49_cummings_baseline_delta_delta_meandiff.png new file mode 100644 index 00000000..d4266e1d Binary files /dev/null and b/dabest/tests/baseline_images/test_49_cummings_baseline_delta_delta_meandiff.png differ diff --git a/dabest/tests/baseline_images/test_50_delta_plot_ylabel.png b/dabest/tests/baseline_images/test_50_delta_plot_ylabel.png new file mode 100644 index 00000000..4d7f62ab Binary files /dev/null and b/dabest/tests/baseline_images/test_50_delta_plot_ylabel.png differ diff --git a/dabest/tests/baseline_images/test_51_delta_plot_change_palette_a.png b/dabest/tests/baseline_images/test_51_delta_plot_change_palette_a.png new file mode 100644 index 00000000..04ef61e2 Binary files /dev/null and b/dabest/tests/baseline_images/test_51_delta_plot_change_palette_a.png differ diff --git a/dabest/tests/baseline_images/test_52_delta_dot_sizes.png b/dabest/tests/baseline_images/test_52_delta_dot_sizes.png new file mode 100644 index 00000000..5c98b31c Binary files /dev/null and b/dabest/tests/baseline_images/test_52_delta_dot_sizes.png differ diff --git a/dabest/tests/baseline_images/test_53_delta_change_ylims.png b/dabest/tests/baseline_images/test_53_delta_change_ylims.png new file mode 100644 index 00000000..6b9a1c91 Binary files /dev/null and b/dabest/tests/baseline_images/test_53_delta_change_ylims.png differ diff --git a/dabest/tests/baseline_images/test_54_delta_invert_ylim.png b/dabest/tests/baseline_images/test_54_delta_invert_ylim.png new file mode 100644 index 00000000..a79a5e3a Binary files /dev/null and b/dabest/tests/baseline_images/test_54_delta_invert_ylim.png differ diff --git a/dabest/tests/baseline_images/test_55_delta_median_diff.png b/dabest/tests/baseline_images/test_55_delta_median_diff.png new file mode 100644 index 00000000..e20f95e9 Binary files /dev/null and b/dabest/tests/baseline_images/test_55_delta_median_diff.png differ diff --git a/dabest/tests/baseline_images/test_56_delta_cohens_d.png b/dabest/tests/baseline_images/test_56_delta_cohens_d.png new file mode 100644 index 00000000..3a6f2798 Binary files /dev/null and b/dabest/tests/baseline_images/test_56_delta_cohens_d.png differ diff --git a/dabest/tests/test_01_effsizes_pvals.py b/dabest/tests/test_01_effsizes_pvals.py index 3b2fba7d..6d2561a3 100644 --- a/dabest/tests/test_01_effsizes_pvals.py +++ b/dabest/tests/test_01_effsizes_pvals.py @@ -55,7 +55,10 @@ # kwargs for Dabest class init. dabest_default_kwargs = dict(x=None, y=None, ci=95, - resamples=5000, random_seed=12345) + resamples=5000, random_seed=12345, + proportional=False, delta2=False, experiment=None, + experiment_label=None, x1_level=None + ) @@ -79,7 +82,7 @@ def test_mean_diff_paired(): from numpy import mean as npmean mean_diff = effsize.func_difference(paired_wellbeing.pre, paired_wellbeing.post, - npmean, is_paired=True) + npmean, is_paired="baseline") assert mean_diff == pytest.approx(4.10) @@ -88,7 +91,7 @@ def test_median_diff_paired(): from numpy import median as npmedian median_diff = effsize.func_difference(paired_wellbeing.pre, paired_wellbeing.post, - npmedian, is_paired=True) + npmedian, is_paired="baseline") assert median_diff == pytest.approx(4.5) @@ -112,7 +115,7 @@ def test_hedges_g_unpaired(): def test_cohens_d_paired(): import numpy as np cohens_d = effsize.cohens_d(paired_wellbeing.pre, paired_wellbeing.post, - is_paired=True) + is_paired="baseline") assert np.round(cohens_d, 2) == pytest.approx(0.34) @@ -120,7 +123,7 @@ def test_cohens_d_paired(): def test_hedges_g_paired(): import numpy as np hedges_g = effsize.hedges_g(paired_wellbeing.pre, paired_wellbeing.post, - is_paired=True) + is_paired="baseline") assert np.round(hedges_g, 2) == pytest.approx(0.33) @@ -155,7 +158,7 @@ def test_paired_stats(): before = paired_wellbeing.pre after = paired_wellbeing.post - paired_es = TwoGroupsEffectSize(before, after, "mean_diff", is_paired=True) + paired_es = TwoGroupsEffectSize(before, after, "mean_diff", is_paired="baseline") p1 = sp.stats.ttest_rel(before, after, nan_policy='omit').pvalue assert paired_es.pvalue_paired_students_t == pytest.approx(p1) @@ -197,7 +200,7 @@ def test_paired_permutation_test(): perm_test = PermutationTest(paired_wellbeing.pre, paired_wellbeing.post, effect_size="mean_diff", - is_paired=True) + is_paired="baseline") assert perm_test.pvalue == pytest.approx(0.0124) @@ -223,7 +226,7 @@ def test_lqrt_unpaired(): def test_lqrt_paired(): paired_dabest = Dabest(paired_wellbeing, idx=("pre", "post"), - paired=True, id_col="ID", + paired="baseline", id_col="ID", **dabest_default_kwargs) lqrt_result = paired_dabest.mean_diff.lqrt diff --git a/dabest/tests/test_03_plotting.py b/dabest/tests/test_03_plotting.py index 78511321..818d8f86 100644 --- a/dabest/tests/test_03_plotting.py +++ b/dabest/tests/test_03_plotting.py @@ -25,7 +25,7 @@ two_groups_unpaired = load(df, idx=("Control 1", "Test 1")) two_groups_paired = load(df, idx=("Control 1", "Test 1"), - paired=True, id_col="ID") + paired="baseline", id_col="ID") multi_2group = load(df, idx=(("Control 1", "Test 1",), ("Control 2", "Test 2")) @@ -34,7 +34,7 @@ multi_2group_paired = load(df, idx=(("Control 1", "Test 1"), ("Control 2", "Test 2")), - paired=True, id_col="ID") + paired="baseline", id_col="ID") shared_control = load(df, idx=("Control 1", "Test 1", "Test 2", "Test 3", @@ -137,7 +137,7 @@ def test_11_inset_plots(): iris_dabest3 = load(data=iris_melt[iris_melt.species=="setosa"], x="metric", y="value", idx=("sepal_length", "sepal_width"), - paired=True, id_col="index") + paired="baseline", id_col="index") diff --git a/dabest/tests/test_04_repeated_measures_effsizes_pvals copy.py b/dabest/tests/test_04_repeated_measures_effsizes_pvals.py similarity index 97% rename from dabest/tests/test_04_repeated_measures_effsizes_pvals copy.py rename to dabest/tests/test_04_repeated_measures_effsizes_pvals.py index 5d2cc4e7..e3c052ce 100644 --- a/dabest/tests/test_04_repeated_measures_effsizes_pvals copy.py +++ b/dabest/tests/test_04_repeated_measures_effsizes_pvals.py @@ -37,7 +37,9 @@ # kwargs for Dabest class init. dabest_default_kwargs = dict(x=None, y=None, ci=95, - resamples=5000, random_seed=12345, proportional=False) + resamples=5000, random_seed=12345, proportional=False, + delta2 = False, experiment=None, + experiment_label=None, x1_level=None) # example of sequential repeated measures sequential = Dabest(df, id_col = "ID", diff --git a/dabest/tests/test_06_delta-delta_effsize_pvals.py b/dabest/tests/test_06_delta-delta_effsize_pvals.py new file mode 100644 index 00000000..c6aefbf5 --- /dev/null +++ b/dabest/tests/test_06_delta-delta_effsize_pvals.py @@ -0,0 +1,306 @@ +import sys +import pytest +import lqrt +import numpy as np +import scipy as sp +import pandas as pd +from .._stats_tools import effsize +from .._classes import TwoGroupsEffectSize, PermutationTest, Dabest + + +# Data for tests. +# See: Asheber Abebe. Introduction to Design and Analysis of Experiments +# with the SAS, from Example: Two-way RM Design Pg 137. +hr = [72, 78, 71, 72, 66, 74, 62, 69, 69, 66, 84, 80, 72, 65, 75, 71, + 86, 83, 82, 83, 79, 83, 73, 75, 73, 62, 90, 81, 72, 62, 69, 70] + +# Add experiment column +e1 = np.repeat('Treatment1', 8).tolist() +e2 = np.repeat('Control', 8).tolist() +experiment = e1 + e2 + e1 + e2 + +# Add a `Drug` column as the first variable +d1 = np.repeat('AX23', 8).tolist() +d2 = np.repeat('CONTROL', 8).tolist() +drug = d1 + d2 + d1 + d2 + +# Add a `Time` column as the second variable +t1 = np.repeat('T1', 16).tolist() +t2 = np.repeat('T2', 16).tolist() +time = t1 + t2 + +# Add an `id` column for paired data plotting. +id_col = [] +for i in range(1, 9): + id_col.append(str(i)+"a") +for i in range(1, 9): + id_col.append(str(i)+"c") +id_col.extend(id_col) + +# Combine samples and gender into a DataFrame. +df_test = pd.DataFrame({'ID' : id_col, + 'Drug' : drug, + 'Time' : time, + 'Experiment': experiment, + 'Heart Rate': hr + }) + + +df_test_control = df_test[df_test["Experiment"]=="Control"] +df_test_control = df_test_control.pivot(index="ID", columns="Time", values="Heart Rate") + + +df_test_treatment1 = df_test[df_test["Experiment"]=="Treatment1"] +df_test_treatment1 = df_test_treatment1.pivot(index="ID", columns="Time", values="Heart Rate") + + +# kwargs for Dabest class init. +dabest_default_kwargs = dict(ci=95, + resamples=5000, random_seed=12345, + idx=None, proportional=False + ) + + +# example of unpaired delta-delta calculation +unpaired = Dabest(data = df_test, x = ["Time", "Drug"], y = "Heart Rate", + delta2 = True, experiment = "Experiment", + experiment_label=None, x1_level=None, paired=None, id_col=None, + **dabest_default_kwargs) + + +# example of paired delta-delta calculation +paired = Dabest(data = df_test, x = ["Time", "Drug"], y = "Heart Rate", + delta2 = True, experiment = "Experiment", paired="sequential", id_col="ID", + experiment_label=None, x1_level=None, + **dabest_default_kwargs) + + +# example of paired data with specified experiment/x1 level +paired_specified_level = Dabest(data = df_test, x = ["Time", "Drug"], y = "Heart Rate", + delta2 = True, experiment = "Experiment", paired="sequential", id_col="ID", + experiment_label=["Control", "Treatment1"], x1_level=["T2", "T1"], + **dabest_default_kwargs) + + +def test_mean_diff_delta_unpaired(): + mean_diff_results = unpaired.mean_diff.results + all_mean_diff = mean_diff_results['difference'].to_list() + diff1 = np.mean(df_test_treatment1["T2"])-np.mean(df_test_treatment1["T1"]) + diff2 = np.mean(df_test_control["T2"])-np.mean(df_test_control["T1"]) + diff3 = np.mean(mean_diff_results['bootstraps'][1]) - np.mean(mean_diff_results['bootstraps'][0]) + np_result = [diff1, diff2, diff3] + assert all_mean_diff == pytest.approx(np_result) + + +def test_mean_diff_delta_paired(): + mean_diff_results = paired.mean_diff.results + all_mean_diff = mean_diff_results['difference'].to_list() + diff1 = np.mean(df_test_treatment1["T2"]-df_test_treatment1["T1"]) + diff2 = np.mean(df_test_control["T2"]-df_test_control["T1"]) + diff3 = np.mean(mean_diff_results['bootstraps'][1] - mean_diff_results['bootstraps'][0]) + np_result = [diff1, diff2, diff3] + assert all_mean_diff == pytest.approx(np_result) + + +def test_mean_diff_delta_paired_specified_level(): + mean_diff_results = paired_specified_level.mean_diff.results + all_mean_diff = mean_diff_results['difference'].to_list() + diff1 = np.mean(df_test_control["T1"]-df_test_control["T2"]) + diff2 = np.mean(df_test_treatment1["T1"]-df_test_treatment1["T2"]) + diff3 = np.mean(mean_diff_results['bootstraps'][1] - mean_diff_results['bootstraps'][0]) + np_result = [diff1, diff2, diff3] + assert all_mean_diff == pytest.approx(np_result) + + +def test_median_diff_unpaired(): + all_median_diff = unpaired.median_diff.results + median_diff = all_median_diff['difference'].to_list() + diff1 = np.median(df_test_treatment1["T2"])-np.median(df_test_treatment1["T1"]) + diff2 = np.median(df_test_control["T2"])-np.median(df_test_control["T1"]) + np_result = [diff1, diff2] + assert median_diff == pytest.approx(np_result) + + +def test_median_diff_paired(): + all_median_diff = paired.median_diff.results + median_diff = all_median_diff['difference'].to_list() + diff1 = np.median(df_test_treatment1["T2"]-df_test_treatment1["T1"]) + diff2 = np.median(df_test_control["T2"]-df_test_control["T1"]) + np_result = [diff1, diff2] + assert median_diff == pytest.approx(np_result) + + +def test_median_diff_paired_specified_level(): + all_median_diff = paired_specified_level.median_diff.results + median_diff = all_median_diff['difference'].to_list() + diff1 = np.median(df_test_control["T1"]-df_test_control["T2"]) + diff2 = np.median(df_test_treatment1["T1"]-df_test_treatment1["T2"]) + np_result = [diff1, diff2] + assert median_diff == pytest.approx(np_result) + + +def test_cohens_d_unpaired(): + all_cohens_d = unpaired.cohens_d.results + cohens_d = all_cohens_d['difference'].to_list() + diff1 = np.mean(df_test_treatment1["T2"])-np.mean(df_test_treatment1["T1"]) + diff1 = diff1/np.sqrt((np.var(df_test_treatment1["T2"], ddof=1)+np.var(df_test_treatment1["T1"], ddof=1))/2) + diff2 = np.mean(df_test_control["T2"])-np.mean(df_test_control["T1"]) + diff2 = diff2/np.sqrt((np.var(df_test_control["T2"], ddof=1)+np.var(df_test_control["T1"], ddof=1))/2) + np_result = [diff1, diff2] + assert cohens_d == pytest.approx(np_result) + + +def test_cohens_d_paired(): + all_cohens_d = paired.cohens_d.results + cohens_d = all_cohens_d['difference'].to_list() + diff1 = np.mean(df_test_treatment1["T2"]-df_test_treatment1["T1"]) + diff1 = diff1/np.sqrt((np.var(df_test_treatment1["T2"], ddof=1)+np.var(df_test_treatment1["T1"], ddof=1))/2) + diff2 = np.mean(df_test_control["T2"]-df_test_control["T1"]) + diff2 = diff2/np.sqrt((np.var(df_test_control["T2"], ddof=1)+np.var(df_test_control["T1"], ddof=1))/2) + np_result = [diff1, diff2] + assert cohens_d == pytest.approx(np_result) + + +def test_cohens_d_paired_specified_level(): + all_cohens_d = paired_specified_level.cohens_d.results + cohens_d = all_cohens_d['difference'].to_list() + diff1 = np.mean(df_test_control["T1"])-np.mean(df_test_control["T2"]) + diff1 = diff1/np.sqrt((np.var(df_test_control["T2"], ddof=1)+np.var(df_test_control["T1"], ddof=1))/2) + diff2 = np.mean(df_test_treatment1["T1"])-np.mean(df_test_treatment1["T2"]) + diff2 = diff2/np.sqrt((np.var(df_test_treatment1["T2"], ddof=1)+np.var(df_test_treatment1["T1"], ddof=1))/2) + np_result = [diff1, diff2] + assert cohens_d == pytest.approx(np_result) + + +def test_hedges_g_unpaired(): + from math import gamma + hedges_g = unpaired.hedges_g.results['difference'].to_list() + a = 8*2-2 + fac = gamma(a/2)/(np.sqrt(a/2)*gamma((a-1)/2)) + diff1 = (np.mean(df_test_treatment1["T2"])-np.mean(df_test_treatment1["T1"]))*fac + diff1 = diff1/np.sqrt((np.var(df_test_treatment1["T2"], ddof=1)+np.var(df_test_treatment1["T1"], ddof=1))/2) + diff2 = (np.mean(df_test_control["T2"])-np.mean(df_test_control["T1"]))*fac + diff2 = diff2/np.sqrt((np.var(df_test_control["T2"], ddof=1)+np.var(df_test_control["T1"], ddof=1))/2) + np_result=[diff1, diff2] + assert hedges_g == pytest.approx(np_result) + + +def test_hedges_g_paired(): + from math import gamma + hedges_g = paired.hedges_g.results['difference'].to_list() + a = 8*2-2 + fac = gamma(a/2)/(np.sqrt(a/2)*gamma((a-1)/2)) + diff1 = (np.mean(df_test_treatment1["T2"]-df_test_treatment1["T1"]))*fac + diff1 = diff1/np.sqrt((np.var(df_test_treatment1["T2"], ddof=1)+np.var(df_test_treatment1["T1"], ddof=1))/2) + diff2 = (np.mean(df_test_control["T2"]-df_test_control["T1"]))*fac + diff2 = diff2/np.sqrt((np.var(df_test_control["T2"], ddof=1)+np.var(df_test_control["T1"], ddof=1))/2) + np_result=[diff1, diff2] + assert hedges_g == pytest.approx(np_result) + + +def test_hedges_g_paired_specified_level(): + from math import gamma + hedges_g = paired_specified_level.hedges_g.results['difference'].to_list() + a = 8*2-2 + fac = gamma(a/2)/(np.sqrt(a/2)*gamma((a-1)/2)) + diff1 = (np.mean(df_test_control["T1"]-df_test_control["T2"]))*fac + diff1 = diff1/np.sqrt((np.var(df_test_control["T2"], ddof=1)+np.var(df_test_control["T1"], ddof=1))/2) + diff2 = (np.mean(df_test_treatment1["T1"]-df_test_treatment1["T2"]))*fac + diff2 = diff2/np.sqrt((np.var(df_test_treatment1["T2"], ddof=1)+np.var(df_test_treatment1["T1"], ddof=1))/2) + np_result=[diff1, diff2] + assert hedges_g == pytest.approx(np_result) + + +def test_paired_stats_unpaired(): + np_result = unpaired.mean_diff.results + + p1 = sp.stats.ttest_rel(np_result["bootstraps"][0], \ + np_result["bootstraps"][1], nan_policy='omit').pvalue + assert np_result["pvalue_paired_students_t"].to_list()[2] == pytest.approx(p1) + + p2 = sp.stats.wilcoxon(np_result["bootstraps"][0], \ + np_result["bootstraps"][1],).pvalue + assert np_result["pvalue_wilcoxon"].to_list()[2] == pytest.approx(p2) + + +def test_paired_stats_paired(): + np_result = paired.mean_diff.results + p1 = sp.stats.ttest_rel(df_test_treatment1["T1"], \ + df_test_treatment1["T2"], nan_policy='omit').pvalue + p2 = sp.stats.ttest_rel(df_test_control["T1"], \ + df_test_control["T2"], nan_policy='omit').pvalue + p3 = sp.stats.ttest_rel(np_result["bootstraps"][0], \ + np_result["bootstraps"][1], nan_policy='omit').pvalue + assert np_result["pvalue_paired_students_t"].to_list() == pytest.approx([p1, p2, p3]) + + p1 = sp.stats.wilcoxon(df_test_treatment1["T1"], \ + df_test_treatment1["T2"]).pvalue + p2 = sp.stats.wilcoxon(df_test_control["T1"], \ + df_test_control["T2"]).pvalue + p3 = sp.stats.wilcoxon(np_result["bootstraps"][0], \ + np_result["bootstraps"][1]).pvalue + assert np_result["pvalue_wilcoxon"].to_list() == pytest.approx([p1, p2, p3]) + + +def test_paired_stats_paired_speficied_level(): + np_result = paired_specified_level.mean_diff.results + p1 = sp.stats.ttest_rel(df_test_control["T2"], \ + df_test_control["T1"], nan_policy='omit').pvalue + p2 = sp.stats.ttest_rel(df_test_treatment1["T2"], \ + df_test_treatment1["T1"], nan_policy='omit').pvalue + p3 = sp.stats.ttest_rel(np_result["bootstraps"][0], \ + np_result["bootstraps"][1], nan_policy='omit').pvalue + assert np_result["pvalue_paired_students_t"].to_list() == pytest.approx([p1, p2, p3]) + + p1 = sp.stats.wilcoxon(df_test_control["T2"], \ + df_test_control["T1"]).pvalue + p2 = sp.stats.wilcoxon(df_test_treatment1["T2"], \ + df_test_treatment1["T1"]).pvalue + p3 = sp.stats.wilcoxon(np_result["bootstraps"][0], \ + np_result["bootstraps"][1]).pvalue + assert np_result["pvalue_wilcoxon"].to_list() == pytest.approx([p1, p2, p3]) + + +def test_lqrt_delta_unpaired(): + all_mean_diff = unpaired.mean_diff + lqrt_result = all_mean_diff.lqrt["pvalue_paired_lqrt"].to_list()[2] + + p1 = lqrt.lqrtest_rel(all_mean_diff.results["bootstraps"][0], + all_mean_diff.results["bootstraps"][1], + random_state=12345).pvalue + + assert lqrt_result == pytest.approx(p1) + + +def test_lqrt_delta_paired(): + all_mean_diff = paired.mean_diff + lqrt_result = all_mean_diff.lqrt["pvalue_paired_lqrt"].to_list() + + p1 = lqrt.lqrtest_rel(df_test_treatment1["T2"], + df_test_treatment1["T1"], + random_state=12345).pvalue + p2 = lqrt.lqrtest_rel(df_test_control["T2"], + df_test_control["T1"], + random_state=12345).pvalue + p3 = lqrt.lqrtest_rel(all_mean_diff.results["bootstraps"][0], + all_mean_diff.results["bootstraps"][1], + random_state=12345).pvalue + + assert lqrt_result == pytest.approx([p1, p2, p3]) + +def test_lqrt_delta_paired_specified_level(): + all_mean_diff = paired_specified_level.mean_diff + lqrt_result = all_mean_diff.lqrt["pvalue_paired_lqrt"].to_list() + + p1 = lqrt.lqrtest_rel(df_test_control["T2"], + df_test_control["T1"], + random_state=12345).pvalue + p2 = lqrt.lqrtest_rel(df_test_treatment1["T2"], + df_test_treatment1["T1"], + random_state=12345).pvalue + p3 = lqrt.lqrtest_rel(all_mean_diff.results["bootstraps"][0], + all_mean_diff.results["bootstraps"][1], + random_state=12345).pvalue + + assert lqrt_result == pytest.approx([p1, p2, p3]) + diff --git a/dabest/tests/test_07_delta-delta_plots.py b/dabest/tests/test_07_delta-delta_plots.py new file mode 100644 index 00000000..27ef53d8 --- /dev/null +++ b/dabest/tests/test_07_delta-delta_plots.py @@ -0,0 +1,90 @@ +#!/usr/bin/env python3 + +# -*- coding: utf-8 -*- + + +import pytest +import numpy as np +import pandas as pd +from scipy.stats import norm + + +import matplotlib as mpl +mpl.use('Agg') +import matplotlib.pyplot as plt +import matplotlib.ticker as Ticker + + +from .._api import load +from .utils import create_demo_dataset, create_demo_dataset_rm, create_demo_dataset_delta + + + +df = create_demo_dataset_delta() + +unpaired = load(data = df, x = ["Light", "Genotype"], y = "Y", delta2 = True, + experiment = "Experiment") + +baseline = load(data = df, x = ["Light", "Genotype"], y = "Y", delta2 = True, + experiment = "Experiment", + paired="baseline", id_col="ID") + +sequential = load(data = df, x = ["Light", "Genotype"], y = "Y", delta2 = True, + experiment = "Experiment", + paired="sequential", id_col="ID") + + +@pytest.mark.mpl_image_compare(tolerance=10) +def test_47_cummings_unpaired_delta_delta_meandiff(): + return unpaired.mean_diff.plot(fig_size=(12, 8), raw_marker_size=4); + + +@pytest.mark.mpl_image_compare(tolerance=10) +def test_48_cummings_sequential_delta_delta_meandiff(): + return sequential.mean_diff.plot(); + + +@pytest.mark.mpl_image_compare(tolerance=10) +def test_49_cummings_baseline_delta_delta_meandiff(): + return baseline.mean_diff.plot(); + + +@pytest.mark.mpl_image_compare(tolerance=10) +def test_50_delta_plot_ylabel(): + return baseline.mean_diff.plot(swarm_label="This is my\nrawdata", + contrast_label="The bootstrap\ndistribtions!", + delta_label="This is delta!"); + + +@pytest.mark.mpl_image_compare(tolerance=10) +def test_51_delta_plot_change_palette_a(): + return sequential.mean_diff.plot(custom_palette="Dark2"); + + +@pytest.mark.mpl_image_compare(tolerance=10) +def test_52_delta_dot_sizes(): + return sequential.mean_diff.plot(show_pairs=False,raw_marker_size=3, + es_marker_size=12); + + +@pytest.mark.mpl_image_compare(tolerance=10) +def test_53_delta_change_ylims(): + return sequential.mean_diff.plot(swarm_ylim=(0, 5), + contrast_ylim=(-2, 2), + fig_size=(15,6)); + + +@pytest.mark.mpl_image_compare(tolerance=10) +def test_54_delta_invert_ylim(): + return sequential.mean_diff.plot(contrast_ylim=(2, -2), + contrast_label="More negative is better!"); + + +@pytest.mark.mpl_image_compare(tolerance=10) +def test_55_delta_median_diff(): + return sequential.median_diff.plot(); + + +@pytest.mark.mpl_image_compare(tolerance=10) +def test_56_delta_cohens_d(): + return unpaired.cohens_d.plot(); \ No newline at end of file diff --git a/dabest/tests/utils.py b/dabest/tests/utils.py index f28f122b..79dc62e8 100644 --- a/dabest/tests/utils.py +++ b/dabest/tests/utils.py @@ -86,6 +86,53 @@ def create_demo_dataset_rm(seed=9999, N=20): return df +def create_demo_dataset_delta(seed=9999, N=20): + + import numpy as np + import pandas as pd + from scipy.stats import norm # Used in generation of populations. + + np.random.seed(seed) # Fix the seed so the results are replicable. + # pop_size = 10000 # Size of each population. + + from scipy.stats import norm # Used in generation of populations. + + # Create samples + y = norm.rvs(loc=3, scale=0.4, size=N*2) + + # Add experiment column + e1 = np.repeat('Control', N).tolist() + e2 = np.repeat('Test', N).tolist() + experiment = e1 + e2 + + # Add a `Light` column as the first variable + light = [] + for i in range(N): + light.append('L1') + light.append('L2') + + # Add a `genotype` column as the second variable + g1 = np.repeat('G1', N/2).tolist() + g2 = np.repeat('G2', N/2).tolist() + g3 = np.repeat('G3', N).tolist() + genotype = g1 + g2 + g3 + + # Add an `id` column for paired data plotting. + id_col = [] + for i in range(N): + id_col.append(i) + id_col.append(i) + + # Combine samples and gender into a DataFrame. + df = pd.DataFrame({'ID' : id_col, + 'Light' : light, + 'Genotype' : genotype, + 'Experiment': experiment, + 'Y' : y + }) + + return df + def get_swarm_yspans(coll, round_result=False, decimals=12): """