|
| 1 | +import sys |
| 2 | +import pytest |
| 3 | +import lqrt |
| 4 | +import numpy as np |
| 5 | +import scipy as sp |
| 6 | +import pandas as pd |
| 7 | +from .._stats_tools import effsize |
| 8 | +from .._classes import TwoGroupsEffectSize, PermutationTest, Dabest |
| 9 | + |
| 10 | + |
| 11 | +# Data for tests. |
| 12 | +# See: Asheber Abebe. Introduction to Design and Analysis of Experiments |
| 13 | +# with the SAS, from Example: Two-way RM Design Pg 137. |
| 14 | +hr = [72, 78, 71, 72, 66, 74, 62, 69, 69, 66, 84, 80, 72, 65, 75, 71, |
| 15 | + 86, 83, 82, 83, 79, 83, 73, 75, 73, 62, 90, 81, 72, 62, 69, 70] |
| 16 | + |
| 17 | +# Add experiment column |
| 18 | +e1 = np.repeat('Treatment1', 8).tolist() |
| 19 | +e2 = np.repeat('Control', 8).tolist() |
| 20 | +experiment = e1 + e2 + e1 + e2 |
| 21 | + |
| 22 | +# Add a `Drug` column as the first variable |
| 23 | +d1 = np.repeat('AX23', 8).tolist() |
| 24 | +d2 = np.repeat('CONTROL', 8).tolist() |
| 25 | +drug = d1 + d2 + d1 + d2 |
| 26 | + |
| 27 | +# Add a `Time` column as the second variable |
| 28 | +t1 = np.repeat('T1', 16).tolist() |
| 29 | +t2 = np.repeat('T2', 16).tolist() |
| 30 | +time = t1 + t2 |
| 31 | + |
| 32 | +# Add an `id` column for paired data plotting. |
| 33 | +id_col = [] |
| 34 | +for i in range(1, 9): |
| 35 | + id_col.append(str(i)+"a") |
| 36 | +for i in range(1, 9): |
| 37 | + id_col.append(str(i)+"c") |
| 38 | +id_col.extend(id_col) |
| 39 | + |
| 40 | +# Combine samples and gender into a DataFrame. |
| 41 | +df_test = pd.DataFrame({'ID' : id_col, |
| 42 | + 'Drug' : drug, |
| 43 | + 'Time' : time, |
| 44 | + 'Experiment': experiment, |
| 45 | + 'Heart Rate': hr |
| 46 | + }) |
| 47 | + |
| 48 | + |
| 49 | +df_test_control = df_test[df_test["Experiment"]=="Control"] |
| 50 | +df_test_control = df_test_control.pivot(index="ID", columns="Time", values="Heart Rate") |
| 51 | + |
| 52 | + |
| 53 | +df_test_treatment1 = df_test[df_test["Experiment"]=="Treatment1"] |
| 54 | +df_test_treatment1 = df_test_treatment1.pivot(index="ID", columns="Time", values="Heart Rate") |
| 55 | + |
| 56 | + |
| 57 | +# kwargs for Dabest class init. |
| 58 | +dabest_default_kwargs = dict(ci=95, |
| 59 | + resamples=5000, random_seed=12345, |
| 60 | + idx=None, proportional=False |
| 61 | + ) |
| 62 | + |
| 63 | + |
| 64 | +# example of unpaired delta-delta calculation |
| 65 | +unpaired = Dabest(data = df_test, x = ["Time", "Drug"], y = "Heart Rate", |
| 66 | + delta2 = True, experiment = "Experiment", |
| 67 | + experiment_label=None, x1_level=None, paired=None, id_col=None, |
| 68 | + **dabest_default_kwargs) |
| 69 | + |
| 70 | + |
| 71 | +# example of paired delta-delta calculation |
| 72 | +paired = Dabest(data = df_test, x = ["Time", "Drug"], y = "Heart Rate", |
| 73 | + delta2 = True, experiment = "Experiment", paired="sequential", id_col="ID", |
| 74 | + experiment_label=None, x1_level=None, |
| 75 | + **dabest_default_kwargs) |
| 76 | + |
| 77 | + |
| 78 | +# example of paired data with specified experiment/x1 level |
| 79 | +paired_specified_level = Dabest(data = df_test, x = ["Time", "Drug"], y = "Heart Rate", |
| 80 | + delta2 = True, experiment = "Experiment", paired="sequential", id_col="ID", |
| 81 | + experiment_label=["Control", "Treatment1"], x1_level=["T2", "T1"], |
| 82 | + **dabest_default_kwargs) |
| 83 | + |
| 84 | + |
| 85 | +def test_mean_diff_delta_unpaired(): |
| 86 | + mean_diff_results = unpaired.mean_diff.results |
| 87 | + all_mean_diff = mean_diff_results['difference'].to_list() |
| 88 | + diff1 = np.mean(df_test_treatment1["T2"])-np.mean(df_test_treatment1["T1"]) |
| 89 | + diff2 = np.mean(df_test_control["T2"])-np.mean(df_test_control["T1"]) |
| 90 | + diff3 = np.mean(mean_diff_results['bootstraps'][1]) - np.mean(mean_diff_results['bootstraps'][0]) |
| 91 | + np_result = [diff1, diff2, diff3] |
| 92 | + assert all_mean_diff == pytest.approx(np_result) |
| 93 | + |
| 94 | + |
| 95 | +def test_mean_diff_delta_paired(): |
| 96 | + mean_diff_results = paired.mean_diff.results |
| 97 | + all_mean_diff = mean_diff_results['difference'].to_list() |
| 98 | + diff1 = np.mean(df_test_treatment1["T2"]-df_test_treatment1["T1"]) |
| 99 | + diff2 = np.mean(df_test_control["T2"]-df_test_control["T1"]) |
| 100 | + diff3 = np.mean(mean_diff_results['bootstraps'][1] - mean_diff_results['bootstraps'][0]) |
| 101 | + np_result = [diff1, diff2, diff3] |
| 102 | + assert all_mean_diff == pytest.approx(np_result) |
| 103 | + |
| 104 | + |
| 105 | +def test_mean_diff_delta_paired_specified_level(): |
| 106 | + mean_diff_results = paired_specified_level.mean_diff.results |
| 107 | + all_mean_diff = mean_diff_results['difference'].to_list() |
| 108 | + diff1 = np.mean(df_test_control["T1"]-df_test_control["T2"]) |
| 109 | + diff2 = np.mean(df_test_treatment1["T1"]-df_test_treatment1["T2"]) |
| 110 | + diff3 = np.mean(mean_diff_results['bootstraps'][1] - mean_diff_results['bootstraps'][0]) |
| 111 | + np_result = [diff1, diff2, diff3] |
| 112 | + assert all_mean_diff == pytest.approx(np_result) |
| 113 | + |
| 114 | + |
| 115 | +def test_median_diff_unpaired(): |
| 116 | + all_median_diff = unpaired.median_diff.results |
| 117 | + median_diff = all_median_diff['difference'].to_list() |
| 118 | + diff1 = np.median(df_test_treatment1["T2"])-np.median(df_test_treatment1["T1"]) |
| 119 | + diff2 = np.median(df_test_control["T2"])-np.median(df_test_control["T1"]) |
| 120 | + np_result = [diff1, diff2] |
| 121 | + assert median_diff == pytest.approx(np_result) |
| 122 | + |
| 123 | + |
| 124 | +def test_median_diff_paired(): |
| 125 | + all_median_diff = paired.median_diff.results |
| 126 | + median_diff = all_median_diff['difference'].to_list() |
| 127 | + diff1 = np.median(df_test_treatment1["T2"]-df_test_treatment1["T1"]) |
| 128 | + diff2 = np.median(df_test_control["T2"]-df_test_control["T1"]) |
| 129 | + np_result = [diff1, diff2] |
| 130 | + assert median_diff == pytest.approx(np_result) |
| 131 | + |
| 132 | + |
| 133 | +def test_median_diff_paired_specified_level(): |
| 134 | + all_median_diff = paired_specified_level.median_diff.results |
| 135 | + median_diff = all_median_diff['difference'].to_list() |
| 136 | + diff1 = np.median(df_test_control["T1"]-df_test_control["T2"]) |
| 137 | + diff2 = np.median(df_test_treatment1["T1"]-df_test_treatment1["T2"]) |
| 138 | + np_result = [diff1, diff2] |
| 139 | + assert median_diff == pytest.approx(np_result) |
| 140 | + |
| 141 | + |
| 142 | +def test_cohens_d_unpaired(): |
| 143 | + all_cohens_d = unpaired.cohens_d.results |
| 144 | + cohens_d = all_cohens_d['difference'].to_list() |
| 145 | + diff1 = np.mean(df_test_treatment1["T2"])-np.mean(df_test_treatment1["T1"]) |
| 146 | + diff1 = diff1/np.sqrt((np.var(df_test_treatment1["T2"], ddof=1)+np.var(df_test_treatment1["T1"], ddof=1))/2) |
| 147 | + diff2 = np.mean(df_test_control["T2"])-np.mean(df_test_control["T1"]) |
| 148 | + diff2 = diff2/np.sqrt((np.var(df_test_control["T2"], ddof=1)+np.var(df_test_control["T1"], ddof=1))/2) |
| 149 | + np_result = [diff1, diff2] |
| 150 | + assert cohens_d == pytest.approx(np_result) |
| 151 | + |
| 152 | + |
| 153 | +def test_cohens_d_paired(): |
| 154 | + all_cohens_d = paired.cohens_d.results |
| 155 | + cohens_d = all_cohens_d['difference'].to_list() |
| 156 | + diff1 = np.mean(df_test_treatment1["T2"]-df_test_treatment1["T1"]) |
| 157 | + diff1 = diff1/np.sqrt((np.var(df_test_treatment1["T2"], ddof=1)+np.var(df_test_treatment1["T1"], ddof=1))/2) |
| 158 | + diff2 = np.mean(df_test_control["T2"]-df_test_control["T1"]) |
| 159 | + diff2 = diff2/np.sqrt((np.var(df_test_control["T2"], ddof=1)+np.var(df_test_control["T1"], ddof=1))/2) |
| 160 | + np_result = [diff1, diff2] |
| 161 | + assert cohens_d == pytest.approx(np_result) |
| 162 | + |
| 163 | + |
| 164 | +def test_cohens_d_paired_specified_level(): |
| 165 | + all_cohens_d = paired_specified_level.cohens_d.results |
| 166 | + cohens_d = all_cohens_d['difference'].to_list() |
| 167 | + diff1 = np.mean(df_test_control["T1"])-np.mean(df_test_control["T2"]) |
| 168 | + diff1 = diff1/np.sqrt((np.var(df_test_control["T2"], ddof=1)+np.var(df_test_control["T1"], ddof=1))/2) |
| 169 | + diff2 = np.mean(df_test_treatment1["T1"])-np.mean(df_test_treatment1["T2"]) |
| 170 | + diff2 = diff2/np.sqrt((np.var(df_test_treatment1["T2"], ddof=1)+np.var(df_test_treatment1["T1"], ddof=1))/2) |
| 171 | + np_result = [diff1, diff2] |
| 172 | + assert cohens_d == pytest.approx(np_result) |
| 173 | + |
| 174 | + |
| 175 | +def test_hedges_g_unpaired(): |
| 176 | + from math import gamma |
| 177 | + hedges_g = unpaired.hedges_g.results['difference'].to_list() |
| 178 | + a = 8*2-2 |
| 179 | + fac = gamma(a/2)/(np.sqrt(a/2)*gamma((a-1)/2)) |
| 180 | + diff1 = (np.mean(df_test_treatment1["T2"])-np.mean(df_test_treatment1["T1"]))*fac |
| 181 | + diff1 = diff1/np.sqrt((np.var(df_test_treatment1["T2"], ddof=1)+np.var(df_test_treatment1["T1"], ddof=1))/2) |
| 182 | + diff2 = (np.mean(df_test_control["T2"])-np.mean(df_test_control["T1"]))*fac |
| 183 | + diff2 = diff2/np.sqrt((np.var(df_test_control["T2"], ddof=1)+np.var(df_test_control["T1"], ddof=1))/2) |
| 184 | + np_result=[diff1, diff2] |
| 185 | + assert hedges_g == pytest.approx(np_result) |
| 186 | + |
| 187 | + |
| 188 | +def test_hedges_g_paired(): |
| 189 | + from math import gamma |
| 190 | + hedges_g = paired.hedges_g.results['difference'].to_list() |
| 191 | + a = 8*2-2 |
| 192 | + fac = gamma(a/2)/(np.sqrt(a/2)*gamma((a-1)/2)) |
| 193 | + diff1 = (np.mean(df_test_treatment1["T2"]-df_test_treatment1["T1"]))*fac |
| 194 | + diff1 = diff1/np.sqrt((np.var(df_test_treatment1["T2"], ddof=1)+np.var(df_test_treatment1["T1"], ddof=1))/2) |
| 195 | + diff2 = (np.mean(df_test_control["T2"]-df_test_control["T1"]))*fac |
| 196 | + diff2 = diff2/np.sqrt((np.var(df_test_control["T2"], ddof=1)+np.var(df_test_control["T1"], ddof=1))/2) |
| 197 | + np_result=[diff1, diff2] |
| 198 | + assert hedges_g == pytest.approx(np_result) |
| 199 | + |
| 200 | + |
| 201 | +def test_hedges_g_paired_specified_level(): |
| 202 | + from math import gamma |
| 203 | + hedges_g = paired_specified_level.hedges_g.results['difference'].to_list() |
| 204 | + a = 8*2-2 |
| 205 | + fac = gamma(a/2)/(np.sqrt(a/2)*gamma((a-1)/2)) |
| 206 | + diff1 = (np.mean(df_test_control["T1"]-df_test_control["T2"]))*fac |
| 207 | + diff1 = diff1/np.sqrt((np.var(df_test_control["T2"], ddof=1)+np.var(df_test_control["T1"], ddof=1))/2) |
| 208 | + diff2 = (np.mean(df_test_treatment1["T1"]-df_test_treatment1["T2"]))*fac |
| 209 | + diff2 = diff2/np.sqrt((np.var(df_test_treatment1["T2"], ddof=1)+np.var(df_test_treatment1["T1"], ddof=1))/2) |
| 210 | + np_result=[diff1, diff2] |
| 211 | + assert hedges_g == pytest.approx(np_result) |
| 212 | + |
| 213 | + |
| 214 | +def test_paired_stats_unpaired(): |
| 215 | + np_result = unpaired.mean_diff.results |
| 216 | + |
| 217 | + p1 = sp.stats.ttest_rel(np_result["bootstraps"][0], \ |
| 218 | + np_result["bootstraps"][1], nan_policy='omit').pvalue |
| 219 | + assert np_result["pvalue_paired_students_t"].to_list()[2] == pytest.approx(p1) |
| 220 | + |
| 221 | + p2 = sp.stats.wilcoxon(np_result["bootstraps"][0], \ |
| 222 | + np_result["bootstraps"][1],).pvalue |
| 223 | + assert np_result["pvalue_wilcoxon"].to_list()[2] == pytest.approx(p2) |
| 224 | + |
| 225 | + |
| 226 | +def test_paired_stats_paired(): |
| 227 | + np_result = paired.mean_diff.results |
| 228 | + p1 = sp.stats.ttest_rel(df_test_treatment1["T1"], \ |
| 229 | + df_test_treatment1["T2"], nan_policy='omit').pvalue |
| 230 | + p2 = sp.stats.ttest_rel(df_test_control["T1"], \ |
| 231 | + df_test_control["T2"], nan_policy='omit').pvalue |
| 232 | + p3 = sp.stats.ttest_rel(np_result["bootstraps"][0], \ |
| 233 | + np_result["bootstraps"][1], nan_policy='omit').pvalue |
| 234 | + assert np_result["pvalue_paired_students_t"].to_list() == pytest.approx([p1, p2, p3]) |
| 235 | + |
| 236 | + p1 = sp.stats.wilcoxon(df_test_treatment1["T1"], \ |
| 237 | + df_test_treatment1["T2"]).pvalue |
| 238 | + p2 = sp.stats.wilcoxon(df_test_control["T1"], \ |
| 239 | + df_test_control["T2"]).pvalue |
| 240 | + p3 = sp.stats.wilcoxon(np_result["bootstraps"][0], \ |
| 241 | + np_result["bootstraps"][1]).pvalue |
| 242 | + assert np_result["pvalue_wilcoxon"].to_list() == pytest.approx([p1, p2, p3]) |
| 243 | + |
| 244 | + |
| 245 | +def test_paired_stats_paired_speficied_level(): |
| 246 | + np_result = paired_specified_level.mean_diff.results |
| 247 | + p1 = sp.stats.ttest_rel(df_test_control["T2"], \ |
| 248 | + df_test_control["T1"], nan_policy='omit').pvalue |
| 249 | + p2 = sp.stats.ttest_rel(df_test_treatment1["T2"], \ |
| 250 | + df_test_treatment1["T1"], nan_policy='omit').pvalue |
| 251 | + p3 = sp.stats.ttest_rel(np_result["bootstraps"][0], \ |
| 252 | + np_result["bootstraps"][1], nan_policy='omit').pvalue |
| 253 | + assert np_result["pvalue_paired_students_t"].to_list() == pytest.approx([p1, p2, p3]) |
| 254 | + |
| 255 | + p1 = sp.stats.wilcoxon(df_test_control["T2"], \ |
| 256 | + df_test_control["T1"]).pvalue |
| 257 | + p2 = sp.stats.wilcoxon(df_test_treatment1["T2"], \ |
| 258 | + df_test_treatment1["T1"]).pvalue |
| 259 | + p3 = sp.stats.wilcoxon(np_result["bootstraps"][0], \ |
| 260 | + np_result["bootstraps"][1]).pvalue |
| 261 | + assert np_result["pvalue_wilcoxon"].to_list() == pytest.approx([p1, p2, p3]) |
| 262 | + |
| 263 | + |
| 264 | +def test_lqrt_delta_unpaired(): |
| 265 | + all_mean_diff = unpaired.mean_diff |
| 266 | + lqrt_result = all_mean_diff.lqrt["pvalue_paired_lqrt"].to_list()[2] |
| 267 | + |
| 268 | + p1 = lqrt.lqrtest_rel(all_mean_diff.results["bootstraps"][0], |
| 269 | + all_mean_diff.results["bootstraps"][1], |
| 270 | + random_state=12345).pvalue |
| 271 | + |
| 272 | + assert lqrt_result == pytest.approx(p1) |
| 273 | + |
| 274 | + |
| 275 | +def test_lqrt_delta_paired(): |
| 276 | + all_mean_diff = paired.mean_diff |
| 277 | + lqrt_result = all_mean_diff.lqrt["pvalue_paired_lqrt"].to_list() |
| 278 | + |
| 279 | + p1 = lqrt.lqrtest_rel(df_test_treatment1["T2"], |
| 280 | + df_test_treatment1["T1"], |
| 281 | + random_state=12345).pvalue |
| 282 | + p2 = lqrt.lqrtest_rel(df_test_control["T2"], |
| 283 | + df_test_control["T1"], |
| 284 | + random_state=12345).pvalue |
| 285 | + p3 = lqrt.lqrtest_rel(all_mean_diff.results["bootstraps"][0], |
| 286 | + all_mean_diff.results["bootstraps"][1], |
| 287 | + random_state=12345).pvalue |
| 288 | + |
| 289 | + assert lqrt_result == pytest.approx([p1, p2, p3]) |
| 290 | + |
| 291 | +def test_lqrt_delta_paired_specified_level(): |
| 292 | + all_mean_diff = paired_specified_level.mean_diff |
| 293 | + lqrt_result = all_mean_diff.lqrt["pvalue_paired_lqrt"].to_list() |
| 294 | + |
| 295 | + p1 = lqrt.lqrtest_rel(df_test_control["T2"], |
| 296 | + df_test_control["T1"], |
| 297 | + random_state=12345).pvalue |
| 298 | + p2 = lqrt.lqrtest_rel(df_test_treatment1["T2"], |
| 299 | + df_test_treatment1["T1"], |
| 300 | + random_state=12345).pvalue |
| 301 | + p3 = lqrt.lqrtest_rel(all_mean_diff.results["bootstraps"][0], |
| 302 | + all_mean_diff.results["bootstraps"][1], |
| 303 | + random_state=12345).pvalue |
| 304 | + |
| 305 | + assert lqrt_result == pytest.approx([p1, p2, p3]) |
| 306 | + |
0 commit comments