Games Howell post hoc test.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Games-Howell test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_The Games-Howell post-hoc test is another nonparametric approach to compare combinations of groups or treatments. Although rather similar to Tukey’s test in its formulation, the Games-Howell test does not assume equal variances and sample sizes. The test was designed based on Welch’s degrees of freedom correction and uses Tukey’s studentized range distribution, denoted q. The Games-Howell test is performed on the ranked variables similar to other nonparametric tests. Since the Games-Howell test does not rely on equal variances and sample sizes, it is often recommended over other approaches such as Tukey’s test._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import math\n",
    "from statsmodels.stats.libqsturng import psturng, qsturng\n",
    "from itertools import combinations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_dicts(names, *columns):\n",
    "    \"\"\"\n",
    "    names: np.array, names of the groups\n",
    "    columns: np.array, corresponding data    \n",
    "    \"\"\"\n",
    "    \n",
    "    data = dict(zip(names, columns))\n",
    "    n = []\n",
    "    Mean = []\n",
    "    std = []\n",
    "    for col in columns:\n",
    "        n.append(len(col))\n",
    "        Mean.append(col.mean())\n",
    "        std.append((np.var(col, ddof=1))**2)\n",
    "    nd = dict(zip(names, n))\n",
    "    Md = dict(zip(names, Mean))\n",
    "    stdd = dict(zip(names, std))\n",
    "    return data,nd,Md,stdd\n",
    "\n",
    "def mean_dif(dctry_m, grp):\n",
    "    key1 = grp[0]\n",
    "    key2 = grp[1]\n",
    "    return dctry_m[key1] - dctry_m[key2]\n",
    "\n",
    "def df(dctry_s,dctry_n, grp):\n",
    "    key1=grp[0]\n",
    "    key2=grp[1]\n",
    "    num = (dctry_s[key1]/dctry_n[key1] + dctry_s[key2]/dctry_n[key2])**2\n",
    "    den = ((dctry_s[key1]/dctry_n[key1])**2)/(dctry_n[key1]-1) + ((dctry_s[key2]/dctry_n[key2])**2)/(dctry_n[key2]-1)\n",
    "    return num/den\n",
    "\n",
    "def t(dctry_m, dctry_s, dctry_n, grp):\n",
    "    key1 = grp[0]\n",
    "    key2 = grp[1]\n",
    "    num = abs(mean_dif(dctry_m, grp))\n",
    "    den = math.sqrt(dctry_s[key1]/dctry_n[key1] + dctry_s[key2]/dctry_n[key2])\n",
    "    return num/den\n",
    "\n",
    "def sigma(dctry_s, dctry_n, grp):\n",
    "    key1 = grp[0]\n",
    "    key2 = grp[1]\n",
    "    return math.sqrt(0.5**(dctry_s[key1]/dctry_n[key1] + dctry_s[key2]/dctry_n[key2]))\n",
    "\n",
    "def p(dctry_m, dctry_s, dctry_n, grp):\n",
    "    K = len(dctry_m)\n",
    "    T = t(dctry_m, dctry_s, dctry_n, grp) * math.sqrt(2)\n",
    "    DF = df(dctry_s,dctry_n, grp)\n",
    "    return psturng(T, K, DF)\n",
    "\n",
    "def upper(dctry_m, dctry_s, dctry_n, grp):\n",
    "#     P = p(dctry_m, dctry_s, dctry_n, grp)\n",
    "    P = 0.95\n",
    "    K = len(dctry_m)\n",
    "    DF = df(dctry_s,dctry_n, grp)\n",
    "    return mean_dif(dctry_m, grp) + qsturng(P, K, DF)*sigma(dctry_s, dctry_n, grp)\n",
    "\n",
    "def lower(dctry_m, dctry_s, dctry_n, grp):\n",
    "#     P = p(dctry_m, dctry_s, dctry_n, grp)\n",
    "    P = 0.95\n",
    "    K = len(dctry_m)\n",
    "    DF = df(dctry_s,dctry_n, grp)\n",
    "    return mean_dif(dctry_m, grp) - qsturng(P, K, DF)*sigma(dctry_s, dctry_n, grp)\n",
    "\n",
    "\n",
    "#data = dctry[0]\n",
    "# nd = dctry[1]\n",
    "# Md = dctry[2]\n",
    "# stdd = dctry[3]\n",
    "\n",
    "def Games_Howell(names, *columns):\n",
    "    \"\"\"\n",
    "    names: np.array, names of the groups\n",
    "    columns: np.array, corresponding data   \n",
    "    \"\"\"\n",
    "    dctry = make_dicts(names, *columns)\n",
    "    groups = [list(x) for x in combinations(dctry[0].keys(), 2)]\n",
    "    N = len(groups)\n",
    "    MEAN_DIF = np.zeros(N)\n",
    "    SIGMA = np.zeros(N)\n",
    "    T = np.zeros(N)\n",
    "    DF = np.zeros(N)\n",
    "    P = np.zeros(N)\n",
    "    UPPER = np.zeros(N)\n",
    "    LOWER = np.zeros(N)\n",
    "    \n",
    "    for i in range(N):\n",
    "        grp = groups[i]\n",
    "        MEAN_DIF[i] = mean_dif(dctry[2], grp)\n",
    "        SIGMA[i] = sigma(dctry[3], dctry[1], grp)\n",
    "        T[i] = t(dctry[2], dctry[3], dctry[1], grp)\n",
    "        DF[i] = df(dctry[3], dctry[1], grp)\n",
    "        P[i] = p(dctry[2], dctry[3], dctry[1], grp)\n",
    "        UPPER[i] = upper(dctry[2], dctry[3], dctry[1], grp)\n",
    "        LOWER[i] = lower(dctry[2], dctry[3], dctry[1], grp)\n",
    "    \n",
    "    \n",
    "    data ={'group':groups,'Mean Difference':MEAN_DIF, 'Standard Error':SIGMA, 't':T, 'df':DF, \n",
    "           'p':P, 'Lower limit':LOWER, 'Upper limit':UPPER,}\n",
    "    \n",
    "    return pd.DataFrame(data=data)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Applications"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### *Example 1*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>group</th>\n",
       "      <th>Mean Difference</th>\n",
       "      <th>Standard Error</th>\n",
       "      <th>t</th>\n",
       "      <th>df</th>\n",
       "      <th>p</th>\n",
       "      <th>Lower limit</th>\n",
       "      <th>Upper limit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[a, b]</td>\n",
       "      <td>-0.5</td>\n",
       "      <td>0.990419</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>3.000000</td>\n",
       "      <td>0.112472</td>\n",
       "      <td>-6.347195</td>\n",
       "      <td>5.347195</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[a, c]</td>\n",
       "      <td>-4.8</td>\n",
       "      <td>0.966606</td>\n",
       "      <td>15.333038</td>\n",
       "      <td>4.000000</td>\n",
       "      <td>0.001000</td>\n",
       "      <td>-9.665080</td>\n",
       "      <td>0.065080</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[b, c]</td>\n",
       "      <td>-4.3</td>\n",
       "      <td>0.957345</td>\n",
       "      <td>12.124574</td>\n",
       "      <td>5.951411</td>\n",
       "      <td>0.001000</td>\n",
       "      <td>-8.459073</td>\n",
       "      <td>-0.140927</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    group  Mean Difference  Standard Error          t        df         p  \\\n",
       "0  [a, b]             -0.5        0.990419   3.000000  3.000000  0.112472   \n",
       "1  [a, c]             -4.8        0.966606  15.333038  4.000000  0.001000   \n",
       "2  [b, c]             -4.3        0.957345  12.124574  5.951411  0.001000   \n",
       "\n",
       "   Lower limit  Upper limit  \n",
       "0    -6.347195     5.347195  \n",
       "1    -9.665080     0.065080  \n",
       "2    -8.459073    -0.140927  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "names = np.array(['a','b','c'])\n",
    "a = np.array([1,1,1])\n",
    "b = np.array([1,2,2,1])\n",
    "c = np.array([5,6,7,6,5])\n",
    "d = make_dicts(names,a,b,c)\n",
    "\n",
    "Games_Howell(names, a,b,c)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### *Example 2*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>group</th>\n",
       "      <th>Mean Difference</th>\n",
       "      <th>Standard Error</th>\n",
       "      <th>t</th>\n",
       "      <th>df</th>\n",
       "      <th>p</th>\n",
       "      <th>Lower limit</th>\n",
       "      <th>Upper limit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[a1, a2]</td>\n",
       "      <td>-99.707387</td>\n",
       "      <td>0.910375</td>\n",
       "      <td>191.556101</td>\n",
       "      <td>85.191721</td>\n",
       "      <td>0.001000</td>\n",
       "      <td>-102.778551</td>\n",
       "      <td>-96.636223</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[a1, a3]</td>\n",
       "      <td>-0.591627</td>\n",
       "      <td>0.910583</td>\n",
       "      <td>1.138005</td>\n",
       "      <td>84.794873</td>\n",
       "      <td>0.494916</td>\n",
       "      <td>-3.663743</td>\n",
       "      <td>2.480490</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[a2, a3]</td>\n",
       "      <td>99.115760</td>\n",
       "      <td>0.993159</td>\n",
       "      <td>704.267358</td>\n",
       "      <td>157.826005</td>\n",
       "      <td>0.001000</td>\n",
       "      <td>95.792325</td>\n",
       "      <td>102.439196</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      group  Mean Difference  Standard Error           t          df  \\\n",
       "0  [a1, a2]       -99.707387        0.910375  191.556101   85.191721   \n",
       "1  [a1, a3]        -0.591627        0.910583    1.138005   84.794873   \n",
       "2  [a2, a3]        99.115760        0.993159  704.267358  157.826005   \n",
       "\n",
       "          p  Lower limit  Upper limit  \n",
       "0  0.001000  -102.778551   -96.636223  \n",
       "1  0.494916    -3.663743     2.480490  \n",
       "2  0.001000    95.792325   102.439196  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "names1 = np.array(['a1', 'a2', 'a3'])\n",
    "a1 = np.random.normal(0,2,80)\n",
    "a2 = np.random.normal(100,1,80)\n",
    "a3 = np.random.normal(1,1,80)\n",
    "w = {'a1':a1,'a2':a2,'a3':a3}\n",
    "df1 = pd.DataFrame(w)\n",
    "Games_Howell(names1, df1['a1'], df1['a2'], df1['a3'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### *Example 3*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>group</th>\n",
       "      <th>Mean Difference</th>\n",
       "      <th>Standard Error</th>\n",
       "      <th>t</th>\n",
       "      <th>df</th>\n",
       "      <th>p</th>\n",
       "      <th>Lower limit</th>\n",
       "      <th>Upper limit</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>[r1, r2]</td>\n",
       "      <td>-1.878403</td>\n",
       "      <td>0.524573</td>\n",
       "      <td>1.376731</td>\n",
       "      <td>71.590454</td>\n",
       "      <td>0.515989</td>\n",
       "      <td>-3.829859</td>\n",
       "      <td>0.073053</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>[r1, r3]</td>\n",
       "      <td>-14.828503</td>\n",
       "      <td>0.009278</td>\n",
       "      <td>4.035240</td>\n",
       "      <td>77.286931</td>\n",
       "      <td>0.001000</td>\n",
       "      <td>-14.862960</td>\n",
       "      <td>-14.794047</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>[r1, r4]</td>\n",
       "      <td>-14.851720</td>\n",
       "      <td>0.527355</td>\n",
       "      <td>10.930121</td>\n",
       "      <td>70.473807</td>\n",
       "      <td>0.001000</td>\n",
       "      <td>-16.814264</td>\n",
       "      <td>-12.889175</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>[r2, r3]</td>\n",
       "      <td>-12.950101</td>\n",
       "      <td>0.017266</td>\n",
       "      <td>3.784092</td>\n",
       "      <td>59.350892</td>\n",
       "      <td>0.002001</td>\n",
       "      <td>-13.014649</td>\n",
       "      <td>-12.885553</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>[r2, r4]</td>\n",
       "      <td>-12.973317</td>\n",
       "      <td>0.981361</td>\n",
       "      <td>55.679162</td>\n",
       "      <td>62.412315</td>\n",
       "      <td>0.001000</td>\n",
       "      <td>-16.636826</td>\n",
       "      <td>-9.309807</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>[r3, r4]</td>\n",
       "      <td>-0.023216</td>\n",
       "      <td>0.017358</td>\n",
       "      <td>0.006788</td>\n",
       "      <td>59.197198</td>\n",
       "      <td>0.900000</td>\n",
       "      <td>-0.088111</td>\n",
       "      <td>0.041679</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "      group  Mean Difference  Standard Error          t         df         p  \\\n",
       "0  [r1, r2]        -1.878403        0.524573   1.376731  71.590454  0.515989   \n",
       "1  [r1, r3]       -14.828503        0.009278   4.035240  77.286931  0.001000   \n",
       "2  [r1, r4]       -14.851720        0.527355  10.930121  70.473807  0.001000   \n",
       "3  [r2, r3]       -12.950101        0.017266   3.784092  59.350892  0.002001   \n",
       "4  [r2, r4]       -12.973317        0.981361  55.679162  62.412315  0.001000   \n",
       "5  [r3, r4]        -0.023216        0.017358   0.006788  59.197198  0.900000   \n",
       "\n",
       "   Lower limit  Upper limit  \n",
       "0    -3.829859     0.073053  \n",
       "1   -14.862960   -14.794047  \n",
       "2   -16.814264   -12.889175  \n",
       "3   -13.014649   -12.885553  \n",
       "4   -16.636826    -9.309807  \n",
       "5    -0.088111     0.041679  "
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "names2 = np.array(['r1','r2','r3', 'r4'])\n",
    "r1 = np.random.normal(0,3,70)\n",
    "r2 = np.random.normal(2, 1, 30)\n",
    "r3 = np.random.normal(15, 5, 60)\n",
    "r4 = np.random.normal(15, 1, 70)\n",
    "Games_Howell(names2, r1,r2,r3,r4)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}