run_anova_tests.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"toc\"></a>\n",
    "# RUN ANOVA TESTS\n",
    "\n",
    "This notebook works through process of performing ANOVA and various statistical tests on each simulation set.\n",
    "\n",
    "---\n",
    "- [WORKSPACE VARIABLES](#workspace-variables)\n",
    "- [PREPARE DATA](#prepare-data)\n",
    "- **[TEST ASSUMPTIONS](#test-assumptions)**\n",
    "- **[RUN ANOVA](#run-anova)**\n",
    "- **[PLOT SUMMARY](#plot-summary)**\n",
    "---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. State the null hypotheses\n",
    "\n",
    "__Factor A Main Effect__\n",
    "\n",
    "$H_0: \\mu_{1\\cdot} = \\mu_{2\\cdot} =~...~= \\mu_{n\\cdot}$\n",
    "\n",
    "There is no difference between or among the factor A row means.\n",
    "\n",
    "__Factor B Main Effect__\n",
    "\n",
    "$H_0: \\mu_{\\cdot 1} = \\mu_{\\cdot 2} =~...~= \\mu_{\\cdot m}$\n",
    "\n",
    "There is no difference between or among the factor B column means.\n",
    "\n",
    "__Interaction Effect__\n",
    "\n",
    "$H_0: \\mu_{jk} - \\mu_{j\\cdot} - \\mu_{\\cdot k} + \\mu = 0$ for all $j$, $k$\n",
    "\n",
    "There is no difference in the factor A by factor B (JK) cell means that cannot be explained by the differences among the factor A (row) means, the factor B (column) means, or both.\n",
    "In other words, there is no interaction between the two independent variables (factor A and factor B)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2. State the alternative hypotheses\n",
    "\n",
    "__Factor A Main Effect__\n",
    "\n",
    "$H_A: \\mu_{j\\cdot} \\neq \\mu_{k\\cdot}$ for some $j$, $k$\n",
    "\n",
    "The factor A row means differ (are not equal) OR at least one pair or combination of means differ.\n",
    "\n",
    "__Factor B Main Effect__\n",
    "\n",
    "$H_A: \\mu_{\\cdot j} \\neq \\mu_{\\cdot k}$ for some $j$, $k$\n",
    "\n",
    "The factor B column means differ (are not equal) OR at least one pair or combination of means differ.\n",
    "\n",
    "__Interaction Effect__\n",
    "\n",
    "$H_A: \\mu_{jk} - \\mu_{j\\cdot} - \\mu_{\\cdot k} + \\mu \\neq 0$ for all $j$, $k$\n",
    "\n",
    "There are differences along the cell population means that cannot be attributed to main effects.\n",
    "In other words, there is an interaction between the two independent variables (factor A and factor B)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3. Check underlying assumptions\n",
    "\n",
    "**Check for assumption of independence**\n",
    "\n",
    "Assess how samples are taken to ensure groups are independent.\n",
    "\n",
    "**Check for assumption of normality**\n",
    "\n",
    "_Check for skewness_\n",
    "\n",
    "$$skewness = \\frac{\\sqrt{n(n - 1)}}{n - 2}\n",
    "\\cdot \\frac{\\frac{1}{n}\\sum_i^n (y_i - \\bar{y})^3}{\\left(\\frac{1}{n} \\sum_i^n (y_i - \\bar{y})^2\\right)^{3/2}}$$\n",
    "\n",
    "$$SE(skewness) = \\sqrt{\\frac{6n(n - 1)}{(n - 2)(n + 1)(n + 3)}}$$\n",
    "\n",
    "$$standardized~skewness = \\frac{skewness}{SE(skewness)}$$\n",
    "\n",
    "Standardized skewness is approzimately normal for $n > 20$ so we perform a two-tailed test and reject the null hypothesis if the test statistic is beyond the critical z-value for our selected significance level.\n",
    "\n",
    "> D. N. Joanes and C. A. Gill. Comparing measures of sample skewness and kurtosis. *The Statistician*, 47:183-189, 1998.\n",
    "\n",
    "_Check Shapiro-Wilks test_\n",
    "\n",
    "The Shapiro–Wilk test tests the null hypothesis that a sample $x_1, ..., x_n$ came from a normally distributed population.\n",
    "We can reject this null hypothesis if p value is less than or equal to significance level.\n",
    "\n",
    "**Check for assumption of homogeneity of variance**\n",
    "\n",
    "_Check Levene's test_\n",
    "\n",
    "Levene's test is used to assess the equality of variances for a variable calculated for two or more groups.\n",
    "It tests the null hypothesis that population variances are equal.\n",
    "If the resulting p-value is less than some significance level, the obtained differences in sample variances are unlikely to have occurred based on random sampling from a population with equal variances.\n",
    "Thus, the null hypothesis of equal variances is rejected and it is concluded that there is a difference between the variances in the population.\n",
    "\n",
    "In the case that variances are different, note that ANOVA is robust to different variances if the sample size is equal.\n",
    "\n",
    "> MJ Blanca, A Alarcon, J Arnau, R Bono, & R Bendayan. Effect of variance ratio on ANOVA robustness: Might 1.5 be the limit? _Behavior Research Methods_, 50, 937–962 (2018).\n",
    ">\n",
    "> _With a variance ratio as large as 9, F-test can, at least for the number of groups and sample sizes considered here, still be used without the Type I error rate being affected by heterogeneity when the design is balanced._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4. Calculate main and interaction effects\n",
    "\n",
    "For a two-way ANOVA, there are two _main effects_ (one for each factor).\n",
    "An _interaction_ between the two factors occurs if the effect of the level of one factor is not the same across all levels of the other factor (and vice versa)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5. Perform additional tests based on effects\n",
    "\n",
    "__If the interaction IS significant__, we can conclude that both factors affect the outcome and that the effect of changes in one factor depends on the level of the other factor, _i.e._, there is an interaction between the explanatory variables.\n",
    "The significance results of the main effects are ignored.\n",
    "\n",
    "We perform Simple Main Effects testing on each factor at each level of the other factor.\n",
    "The error term is still calculated from the overall data.\n",
    "To control for Type I error, we use a Bonferroni correction.\n",
    "\n",
    "__If the interaction IS NOT significant__, then we conclude that an additive model is sufficient, _i.e._, effects of changes in one factor are the same at every fixed level of the other factor.\n",
    "In other words, the effects on the outcome of a particular level change for one explanatory variable does not depend on the level of the other explanatory variable. \n",
    "\n",
    "We re-run the ANOVA without the interaction term and interpret the p-values as in one-way ANOVA. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6. Follow up with ad hoc tests\n",
    "\n",
    "For significant effects with more than two levels, we can use Tukey multiple comparison to identify which pairs of level are significant.\n",
    "\n",
    "__If the interaction IS significant__, do comparisons for a given factor at constant levels of the other factor.\n",
    "\n",
    "__If the interaction IS NOT significant__, do comparisons on full data set. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"workspace-variables\"></a>\n",
    "\n",
    "### WORKSPACE VARIABLES\n",
    "<span style=\"float:right;\">[back to top](#toc)</span>\n",
    "\n",
    "Set up workspace variables for running ANOVA.\n",
    "\n",
    "- **`PATH`** is the path for analysis files (`.json` and `.csv` files, `.tar.xz` compressed archives)\n",
    "- **`SIMULATIONS`** is the dictionary of simulations and simulation factors\n",
    "- **`METRICS`** is the dictionary of metrics and metric scaling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "PATH = \"/path/to/analysis/files/\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "SIMULATIONS = {\n",
    "    \"SR\": [\"geometry\", \"dimension\"],\n",
    "    \"CV\": [\"volume\", \"age\"],\n",
    "    \"ND\": [\"profile\", \"level\"],\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "METRICS = {\n",
    "    \"GROWTH\": 30,     # hexagon/day => um/day,\n",
    "    \"SYMMETRY\": 1,    # unitless\n",
    "    \"CYCLES\": 1/60.0, # minutes => hours\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"prepare-data\"></a>\n",
    "\n",
    "### PREPARE DATA\n",
    "<span style=\"float:right;\">[back to top](#toc)</span>\n",
    "\n",
    "Load and process data for ANOVA."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scripts.anova import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scripts.SYSTEM_REPRESENTATION import SYSTEM_REPRESENTATION as SR\n",
    "from scripts.CELL_VARIABILITY import CELL_VARIABILITY as CV\n",
    "from scripts.NUTRIENT_DYNAMICS import NUTRIENT_DYNAMICS as ND"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**LOAD DATA**\n",
    "\n",
    "Load metrics data from the `.SEEDS` (and `.LAYERS.SEEDS` for `SYSTEM_REPRESENTATION`) output files  produced in the basic analysis step."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = {}\n",
    "\n",
    "for metric in METRICS:\n",
    "    extension = f\".SEEDS.{metric}\"\n",
    "    data[f\"SR_{metric}\"] = SR.get(PATH, load_anova_data, (lambda out: pd.concat(out)), extension, timepoint=15.0)\n",
    "    data[f\"CV_{metric}\"] = CV.get(PATH, load_anova_data, (lambda out: pd.concat(out)), extension, timepoint=15.0)\n",
    "    data[f\"ND_{metric}\"] = ND.get(PATH, load_anova_data, (lambda out: pd.concat(out)), extension, timepoint=15.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**FORMAT DATA**\n",
    "\n",
    "Format loaded data into single dataframe with consistent columns names."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "columns = [\"Y\", \"X1\", \"X2\", \"SIMULATION\", \"METRIC\", \"CONTEXT\"]\n",
    "all_dfs = []\n",
    "\n",
    "for key, values in data.items():\n",
    "    # Get specific factors and scale for given simulation and metric.\n",
    "    simulation, metric = key.split(\"_\")\n",
    "    factors = SIMULATIONS[simulation]\n",
    "    scale = METRICS[metric]\n",
    "    \n",
    "    # Population data frame from loaded data.\n",
    "    df = pd.DataFrame(columns=columns)\n",
    "\n",
    "    df[\"Y\"] = values[\"metric\"] * scale\n",
    "    df[\"X1\"] = values[factors[0]] \n",
    "    df[\"X2\"] = values[factors[1]]\n",
    "    df[\"CONTEXT\"] = values[\"context\"]\n",
    "    df[\"SIMULATION\"] = simulation\n",
    "    df[\"METRIC\"] = metric\n",
    "    \n",
    "    all_dfs.append(df)\n",
    "    \n",
    "all_data = pd.concat(all_dfs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_data[all_data.SIMULATION == \"SR\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"test-assumptions\"></a>\n",
    "\n",
    "### TEST ASSUMPTIONS\n",
    "<span style=\"float:right;\">[back to top](#toc)</span>\n",
    "\n",
    "Test ANOVA assumptions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**GET DATA MEANS**\n",
    "\n",
    "Calculate the cell and effect means for each simulation set, metric, and context."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_cell_means = []\n",
    "all_effect_means = []\n",
    "\n",
    "for name, group in all_data.groupby([\"SIMULATION\", \"METRIC\", \"CONTEXT\"]):\n",
    "    cell_means = calculate_cell_means(group)\n",
    "    cell_means[[\"simulation\", \"metric\", \"context\"]] = list(name)\n",
    "    all_cell_means.append(cell_means)\n",
    "    \n",
    "    effect_means = calculate_effect_means(group)\n",
    "    effect_means[[\"simulation\", \"metric\", \"context\"]] = list(name)\n",
    "    all_effect_means.append(effect_means)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell_means_df = pd.concat(all_cell_means)\n",
    "cell_means_df.to_csv(f\"{PATH}_/ANOVA_cell_means.csv\", index=False)\n",
    "\n",
    "effect_means_df = pd.concat(all_effect_means)\n",
    "effect_means_df.to_csv(f\"{PATH}_/ANOVA_effect_means.csv\", index=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TEST NORMALITY**\n",
    "\n",
    "Check normality using skewness and Shapiro-Wilks for each simulation set, metric, and context."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, group in all_data.groupby([\"SIMULATION\", \"METRIC\", \"CONTEXT\"]):\n",
    "    print(name)\n",
    "    skewness = test_skewness(group)\n",
    "    shapiro = test_shapiro(group)\n",
    "    print(skewness.to_markdown(index=False))\n",
    "    print(shapiro.to_markdown(index=False))\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**TEST HOMOGENEITY OF VARIANCE**\n",
    "\n",
    "Check homogeneity of variance using Levene's for each simulation set, metric, and context."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for name, group in all_data.groupby([\"SIMULATION\", \"METRIC\", \"CONTEXT\"]):\n",
    "    print(name)\n",
    "    levene = test_levene(group)\n",
    "    print(levene.to_markdown(index=False))\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"run-anova\"></a>\n",
    "\n",
    "### RUN ANOVA\n",
    "<span style=\"float:right;\">[back to top](#toc)</span>\n",
    "\n",
    "Run ANOVA and ad hoc tests."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_out = {}\n",
    "\n",
    "for (simulation, metric, context), group in all_data.groupby([\"SIMULATION\", \"METRIC\", \"CONTEXT\"]):\n",
    "    out = {}\n",
    "    \n",
    "    # Run two-way ANOVA with interaction.\n",
    "    anova = run_two_way_anova_with_interaction(group, alpha=0.05)\n",
    "    out[\"anova_int\"] = anova.to_dict()\n",
    "\n",
    "    # Depending on interaction significance, perform additional tests.\n",
    "    if anova[\"is_significant\"][\"X1:X2\"]:\n",
    "        simple_main_1 = run_simple_main_effects(group, \"X1\", alpha=0.05)\n",
    "        simple_main_2 = run_simple_main_effects(group, \"X2\", alpha=0.05)\n",
    "        tukey_1 = run_tukey_tests(group, \"X1\", True, alpha=0.01)\n",
    "        tukey_2 = run_tukey_tests(group, \"X2\", True, alpha=0.01)\n",
    "        \n",
    "        out[\"simple_main\"] = {\n",
    "            \"X1\": simple_main_1.to_dict(),\n",
    "            \"X2\": simple_main_2.to_dict(),\n",
    "        }\n",
    "        \n",
    "        out[\"tukey\"] = {\n",
    "            \"X1\": tukey_1.to_dict(\"records\") if tukey_1 is not None else [],\n",
    "            \"X2\": tukey_2.to_dict(\"records\") if tukey_2 is not None else [],\n",
    "        }\n",
    "    else:\n",
    "        anova_noint = run_two_way_anova_without_interaction(group, alpha=0.05)\n",
    "        tukey_1 = run_tukey_tests(group, \"X1\", False, alpha=0.01)\n",
    "        tukey_2 = run_tukey_tests(group, \"X2\", False, alpha=0.01)\n",
    "        \n",
    "        out[\"anova_noint\"] =  anova_noint.to_dict()\n",
    "        \n",
    "        out[\"tukey\"] = {\n",
    "            \"X1\": tukey_1.to_dict(\"records\") if tukey_1 is not None else [],\n",
    "            \"X2\": tukey_2.to_dict(\"records\") if tukey_2 is not None else [],\n",
    "        }\n",
    "    \n",
    "    all_out[f\"{simulation}_{metric}_{context}\"] = out\n",
    "\n",
    "with open(f\"{PATH}_/ANOVA_testing.json\", \"w\") as f:\n",
    "    jn = json.dumps(all_out, indent = 4, separators = (',', ':'), sort_keys=True)\n",
    "    f.write(jn.replace(\"NaN\", '\"nan\"'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a class=\"anchor\" id=\"plot-summary\"></a>\n",
    "\n",
    "### PLOT SUMMARY\n",
    "<span style=\"float:right;\">[back to top](#toc)</span>\n",
    "\n",
    "Plot ANOVA results in summary interaction plots\n",
    "\n",
    "Summary figure for ANOVA results shows interaction plots of cell means with markers indicating:\n",
    "\n",
    "- significant interaction effect\n",
    "- if significant interaction effect, significant simple main effect(s)\n",
    "- if not significant interaction effect, significant main effect(s)\n",
    "\n",
    "[go to figure](http://0.0.0.0:8000/figures/anova_summary.html)"
   ]
  }
 ],
 "metadata": {
  "language_info": {
   "name": "python",
   "pygments_lexer": "ipython3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}