test of stratified permutation test using the sample mean

testPerm.ipynb

{
 "metadata": {
  "name": "",
  "signature": "sha256:b08bfa817da53374d2b46a5e268f47f5bbc118f487251147b77854d3fa7a3f38"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numpy as np\n",
      "def stratifiedPermutationTest(group, condition, response, iterations):\n",
      "    '''\n",
      "    Stratified permutation test using the sum of the differences in means between two conditions in\n",
      "    each group (stratum) as the test statistic.\n",
      "    The test statistic is\n",
      "        \\sum_{g in groups} [\n",
      "                            mean(response for cases in group g assigned to first condition) -\n",
      "                            mean(response for cases in group g assigned to second condition)\n",
      "                           ].\n",
      "    There should be at least one group and no more than two conditions.\n",
      "    Under the null hypothesis, all assignments to the two conditions that preserve the number of\n",
      "    cases assigned to the conditions are equally likely.\n",
      "    Groups in which all cases are assigned to the same condition are skipped; they do not contribute \n",
      "    to the p-value since all randomizations give the same contribution to the difference in means.\n",
      "    \n",
      "    Dependencies: numpy (as np)\n",
      "    '''   \n",
      "    groups = np.unique(group)\n",
      "    conditions = np.unique(condition) # to do: test that number of conditions <= 2\n",
      "    def testStatistic(group, condition, response, groups, conditions):\n",
      "        tst = 0.0\n",
      "        for g in groups:\n",
      "             x = (group == g) & (condition == conditions[0])\n",
      "             y = (group == g) & (condition == conditions[1])\n",
      "             if (any(x) & any(y)):\n",
      "                 tst = tst + \\\n",
      "                      np.mean(response[(group == g) & (condition == conditions[0])]) - \\\n",
      "                      np.mean(response[(group == g) & (condition == conditions[1])])\n",
      "        return tst\n",
      "#   \n",
      "    def permuteWithinGroups(group, condition, groups):\n",
      "        permuted = condition\n",
      "        for g in groups:\n",
      "            permuted[group == g] = np.random.permutation(condition[group == g])      \n",
      "        return permuted\n",
      "#\n",
      "#  (alternate approaches: assign randoms to all, then sort, or use pandas.groupby and transform(np.random.permutation))\n",
      "    tst = testStatistic(group, condition, response, groups, conditions)\n",
      "    dist = np.zeros(iterations)\n",
      "    for i in np.arange(iterations):\n",
      "         dist[i] = testStatistic( group, \n",
      "                                  permuteWithinGroups(group, condition, groups), \\\n",
      "                                  response, groups, conditions\\\n",
      "                                )\n",
      "    pLeft = float(np.sum(dist <= tst))/float(iterations)\n",
      "    pRight = float(np.sum(dist >= tst))/float(iterations)\n",
      "    pBoth = float(np.sum(abs(dist) >= abs(tst)))/float(iterations)\n",
      "    return pLeft, pRight, pBoth, tst, dist\n"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "group = np.array([0, 0, 0, 1, 1, 1])\n",
      "condition = np.array([0, 0, 1, 0, 1, 1])\n",
      "response = np.array([0, 0, 1, 0, 1, 1])\n",
      "# for these data, left-sided and two-sided p-values should be 1/9 = .111\n",
      "iterations = 10**4\n",
      "leftP, rightP, bothP, teststat, dist = stratifiedPermutationTest(group, condition, response, iterations)\n",
      "print \"left, right, both p:\", leftP, rightP, bothP, \"test statistic:\", teststat"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "left, right, both p: 0.1135 1.0 0.1135 test statistic: -2.0\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}