BCRA_R_to_Python.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 7,
   "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>ID</th>\n",
       "      <th>T1</th>\n",
       "      <th>T2</th>\n",
       "      <th>N_Biop</th>\n",
       "      <th>HypPlas</th>\n",
       "      <th>AgeMen</th>\n",
       "      <th>Age1st</th>\n",
       "      <th>N_Rels</th>\n",
       "      <th>Race</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>99</td>\n",
       "      <td>99</td>\n",
       "      <td>10</td>\n",
       "      <td>20</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>99</td>\n",
       "      <td>1</td>\n",
       "      <td>10</td>\n",
       "      <td>20</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>99</td>\n",
       "      <td>0</td>\n",
       "      <td>10</td>\n",
       "      <td>20</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>0</td>\n",
       "      <td>99</td>\n",
       "      <td>10</td>\n",
       "      <td>20</td>\n",
       "      <td>1</td>\n",
       "      <td>3</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>1</td>\n",
       "      <td>99</td>\n",
       "      <td>10</td>\n",
       "      <td>20</td>\n",
       "      <td>1</td>\n",
       "      <td>4</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   ID    T1    T2  N_Biop  HypPlas  AgeMen  Age1st  N_Rels  Race\n",
       "0   1  45.2  53.3      99       99      10      20       1     0\n",
       "1   2  45.2  53.3      99        1      10      20       1     1\n",
       "2   3  45.2  53.3      99        0      10      20       1     2\n",
       "3   4  45.2  53.3       0       99      10      20       1     3\n",
       "4   5  45.2  53.3       1       99      10      20       1     4"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np \n",
    "import pandas as pd\n",
    "data = pd.read_csv('BCRA_data.csv')\n",
    "data.head()"
   ]
  },
  {
   "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>Error_Ind</th>\n",
       "      <th>set_T1_missing</th>\n",
       "      <th>set_T2_missing</th>\n",
       "      <th>NB_Cat</th>\n",
       "      <th>AM_Cat</th>\n",
       "      <th>AF_Cat</th>\n",
       "      <th>NR_Cat</th>\n",
       "      <th>R_Hyp</th>\n",
       "      <th>set_HyperP_missing</th>\n",
       "      <th>set_R_Hyp_missing</th>\n",
       "      <th>set_Race_missing</th>\n",
       "      <th>CharRace</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>-300</td>\n",
       "      <td>??</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>1.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>-100.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-100</td>\n",
       "      <td>-100.00</td>\n",
       "      <td>1</td>\n",
       "      <td>Wh</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>1.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>-100.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-100</td>\n",
       "      <td>-100.00</td>\n",
       "      <td>2</td>\n",
       "      <td>AA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>0.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>3</td>\n",
       "      <td>HU</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>0.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>4</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>0.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>5</td>\n",
       "      <td>HF</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>0.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>6</td>\n",
       "      <td>Ch</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>0.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.82</td>\n",
       "      <td>1</td>\n",
       "      <td>1.82</td>\n",
       "      <td>7</td>\n",
       "      <td>Ja</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>1.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>-100.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-100</td>\n",
       "      <td>-100.00</td>\n",
       "      <td>8</td>\n",
       "      <td>Fi</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>0.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.93</td>\n",
       "      <td>0</td>\n",
       "      <td>0.93</td>\n",
       "      <td>9</td>\n",
       "      <td>Hw</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>1.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>-100.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-100</td>\n",
       "      <td>-100.00</td>\n",
       "      <td>10</td>\n",
       "      <td>oP</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>1.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>-100.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-100</td>\n",
       "      <td>-100.00</td>\n",
       "      <td>11</td>\n",
       "      <td>oA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>1.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>0.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>-300</td>\n",
       "      <td>??</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>1.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>-100.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-100</td>\n",
       "      <td>-100.00</td>\n",
       "      <td>-300</td>\n",
       "      <td>??</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>1.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>-100.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>-100</td>\n",
       "      <td>-100.00</td>\n",
       "      <td>1</td>\n",
       "      <td>Wh</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>0.0</td>\n",
       "      <td>45.2</td>\n",
       "      <td>53.3</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.93</td>\n",
       "      <td>0</td>\n",
       "      <td>0.93</td>\n",
       "      <td>2</td>\n",
       "      <td>AA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>0.0</td>\n",
       "      <td>35.0</td>\n",
       "      <td>40.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>3</td>\n",
       "      <td>HU</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>0.0</td>\n",
       "      <td>35.0</td>\n",
       "      <td>40.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>4</td>\n",
       "      <td>NA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>1.0</td>\n",
       "      <td>35.0</td>\n",
       "      <td>40.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>5</td>\n",
       "      <td>HF</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>1.0</td>\n",
       "      <td>35.0</td>\n",
       "      <td>40.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>6</td>\n",
       "      <td>Ch</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>0.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>90.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>7</td>\n",
       "      <td>Ja</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>0.0</td>\n",
       "      <td>27.0</td>\n",
       "      <td>90.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>8</td>\n",
       "      <td>Fi</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>26.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>9</td>\n",
       "      <td>Hw</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>1.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>10</td>\n",
       "      <td>oP</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>1.0</td>\n",
       "      <td>85.0</td>\n",
       "      <td>NaN</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>11</td>\n",
       "      <td>oA</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>1.0</td>\n",
       "      <td>86.0</td>\n",
       "      <td>90.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>1.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>1.00</td>\n",
       "      <td>99</td>\n",
       "      <td>1.00</td>\n",
       "      <td>-300</td>\n",
       "      <td>??</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    Error_Ind  set_T1_missing  set_T2_missing  NB_Cat  AM_Cat  AF_Cat  NR_Cat  \\\n",
       "0         1.0            45.2            53.3     0.0     2.0     1.0     1.0   \n",
       "1         1.0            45.2            53.3  -100.0     2.0     1.0     1.0   \n",
       "2         1.0            45.2            53.3  -100.0     1.0     0.0     1.0   \n",
       "3         0.0            45.2            53.3     0.0     0.0     1.0     1.0   \n",
       "4         0.0            45.2            53.3     1.0     2.0     1.0     1.0   \n",
       "5         0.0            45.2            53.3     1.0     0.0     0.0     1.0   \n",
       "6         0.0            45.2            53.3     0.0     0.0     0.0     1.0   \n",
       "7         0.0            45.2            53.3     1.0     0.0     0.0     1.0   \n",
       "8         1.0            45.2            53.3  -100.0     0.0     0.0     1.0   \n",
       "9         0.0            45.2            53.3     1.0     0.0     0.0     1.0   \n",
       "10        1.0            45.2            53.3  -100.0     0.0     0.0     1.0   \n",
       "11        1.0            45.2            53.3  -100.0     0.0     0.0     1.0   \n",
       "12        1.0            45.2            53.3     0.0     2.0     1.0     1.0   \n",
       "13        1.0            45.2            53.3  -100.0     2.0     1.0     1.0   \n",
       "14        1.0            45.2            53.3  -100.0     2.0     1.0     1.0   \n",
       "15        0.0            45.2            53.3     1.0     1.0     0.0     1.0   \n",
       "16        0.0            35.0            40.0     1.0     0.0     1.0     0.0   \n",
       "17        0.0            35.0            40.0     2.0     2.0     2.0     0.0   \n",
       "18        1.0            35.0            40.0     1.0     2.0     NaN     0.0   \n",
       "19        1.0            35.0            40.0     2.0     NaN     NaN     0.0   \n",
       "20        0.0            27.0            90.0     0.0     1.0     1.0     0.0   \n",
       "21        0.0            27.0            90.0     0.0     1.0     1.0     0.0   \n",
       "22        1.0             NaN            26.0     0.0     1.0     NaN     0.0   \n",
       "23        1.0             NaN             NaN     0.0     1.0     1.0     0.0   \n",
       "24        1.0            85.0             NaN     0.0     1.0     1.0     0.0   \n",
       "25        1.0            86.0            90.0     0.0     1.0     1.0     0.0   \n",
       "\n",
       "    R_Hyp  set_HyperP_missing  set_R_Hyp_missing  set_Race_missing CharRace  \n",
       "0    1.00                  99               1.00              -300       ??  \n",
       "1     NaN                -100            -100.00                 1       Wh  \n",
       "2     NaN                -100            -100.00                 2       AA  \n",
       "3    1.00                  99               1.00                 3       HU  \n",
       "4    1.00                  99               1.00                 4       NA  \n",
       "5    1.00                  99               1.00                 5       HF  \n",
       "6    1.00                  99               1.00                 6       Ch  \n",
       "7    1.82                   1               1.82                 7       Ja  \n",
       "8     NaN                -100            -100.00                 8       Fi  \n",
       "9    0.93                   0               0.93                 9       Hw  \n",
       "10    NaN                -100            -100.00                10       oP  \n",
       "11    NaN                -100            -100.00                11       oA  \n",
       "12   1.00                  99               1.00              -300       ??  \n",
       "13    NaN                -100            -100.00              -300       ??  \n",
       "14    NaN                -100            -100.00                 1       Wh  \n",
       "15   0.93                   0               0.93                 2       AA  \n",
       "16   1.00                  99               1.00                 3       HU  \n",
       "17   1.00                  99               1.00                 4       NA  \n",
       "18   1.00                  99               1.00                 5       HF  \n",
       "19   1.00                  99               1.00                 6       Ch  \n",
       "20   1.00                  99               1.00                 7       Ja  \n",
       "21   1.00                  99               1.00                 8       Fi  \n",
       "22   1.00                  99               1.00                 9       Hw  \n",
       "23   1.00                  99               1.00                10       oP  \n",
       "24   1.00                  99               1.00                11       oA  \n",
       "25   1.00                  99               1.00              -300       ??  "
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Make sure my version of NB_Cat is correct\n",
    "def recode_check(data, Raw_Ind=1):\n",
    "    ### set error indicator to default value of 0 for each subject\n",
    "    ## if mean not 0, implies ERROR in file\n",
    "    Error_Ind = np.zeros(data.shape[0])\n",
    "    ### test for consistency of T1 (initial age) and T2 (projection age)\n",
    "    set_T1_missing = data.T1.values.copy()\n",
    "\n",
    "    set_T2_missing = data.T2.values.copy()\n",
    "    set_T1_missing[np.where((data.T1 <20) | (data.T1 >= 90) | (data.T1 >= data.T2))] = np.nan\n",
    "    set_T2_missing[(data.T2.values > 90) | (data.T1.values >= data.T2.values)] = np.nan\n",
    "    Error_Ind[np.isnan(set_T1_missing)] = 1\n",
    "    Error_Ind[np.isnan(set_T2_missing)] = 1\n",
    " \n",
    "    ### RR covariates are in raw/original format  \n",
    "    if Raw_Ind == 1:\n",
    "        ### test for consistency of NumBiop (#biopsies) and Hyperplasia   \n",
    "        ## set NB_Cat to default value of -1\n",
    "        NB_Cat = np.repeat(-1., data.shape[0])\n",
    " \n",
    "        ## REQUIREMENT (A)\n",
    "        NB_Cat[((data.N_Biop == 0) | (data.N_Biop == 99)) & (data.HypPlas != 99)] = -100\n",
    "        Error_Ind[np.where(NB_Cat == -100)] = 1\n",
    " \n",
    "        ## REQUIREMENT (B)\n",
    "        NB_Cat[((data.N_Biop > 0) & (data.N_Biop < 99)) & ((data.HypPlas != 0) & (data.HypPlas != 1) & (data.HypPlas != 99))] = -200\n",
    "        Error_Ind[NB_Cat == -200] = 1\n",
    "     \n",
    "        ### editing and recoding for N_Biop\n",
    "        NB_Cat[(NB_Cat == -1) & ((data.N_Biop == 0) | (data.N_Biop == 99))] = 0\n",
    "        \n",
    "        NB_Cat[(NB_Cat == -1) & (data.N_Biop == 1)] = 1\n",
    "        NB_Cat[(NB_Cat == -1) & ((data.N_Biop >= 2) | (data.N_Biop != 99))] = 2\n",
    "        NB_Cat[NB_Cat == -1] = np.nan\n",
    "        \n",
    "        ### editing and recoding for AgeMen\n",
    "        AM_Cat = np.repeat(np.nan, data.shape[0])\n",
    "        AM_Cat[((data.AgeMen >= 14) & (data.AgeMen <= data.T1)) | (data.AgeMen ==99)] = 0\n",
    "        AM_Cat[(data.AgeMen >= 12) & (data.AgeMen < 14)] = 1\n",
    "        AM_Cat[(data.AgeMen > 0) & (data.AgeMen < 12)] = 2\n",
    "        AM_Cat[(data.AgeMen > data.T1) & (data.AgeMen !=99)] = np.nan\n",
    "        ## for African-Americans AgeMen code 2 (age <= 11) grouped with code 1(age == 12 or 13)\n",
    "        AM_Cat[(data.Race == 2) & (AM_Cat ==2)] = 1 \n",
    " \n",
    "        ### editing and recoding for Age1st\n",
    "        AF_Cat = np.repeat(np.nan, data.shape[0])\n",
    "        AF_Cat[(data.Age1st < 20) | (data.Age1st == 99)] = 0\n",
    "        AF_Cat[(data.Age1st >= 20) & (data.Age1st < 25)] = 1\n",
    "        AF_Cat[(((data.Age1st >= 25) & (data.Age1st < 30))) | (data.Age1st == 98)] = 2\n",
    "        AF_Cat[(data.Age1st >= 30) & (data.Age1st < 98)] = 3\n",
    "        AF_Cat[(data.Age1st < data.AgeMen) & (data.AgeMen != 99)] = np.nan\n",
    "        AF_Cat[(data.Age1st > data.T1) & (data.Age1st <98)] = np.nan\n",
    "        ## for African-Americans Age1st is not a RR covariate and not in RR model, set to 0\n",
    "        AF_Cat[data.Race == 2] = 0 \n",
    "        \n",
    "        ### editing and recoding for N_Rels\n",
    "        NR_Cat = np.repeat(np.nan, data.shape[0])\n",
    "        NR_Cat[(data.N_Rels == 0) | (data.N_Rels == 99)] = 0\n",
    "        NR_Cat[data.N_Rels == 1] = 1\n",
    "        NR_Cat[(data.N_Rels >= 2) & (data.N_Rels < 99)] = 2\n",
    "        ## for Asian-American NR_Cat=2 is pooled with NR_Cat=2\n",
    "        NR_Cat[((data.Race >= 6) & (data.Race <= 11)) & (NR_Cat == 2)] = 1\n",
    "    \n",
    " \n",
    "    ### Raw_Ind=0 means RR covariates have already been re-coded to 0, 1, 2 or 3 (when necessary)\n",
    "    ### edit/consistency checks for all relative four risk covariates not performed when Raw_Ind=0. (use this option at your own risk)\n",
    "    if Raw_Ind == 0:\n",
    "        NB_Cat = data.N_Biop\n",
    "        AM_Cat = data.AgeMen\n",
    "        AF_Cat = data.Age1st\n",
    "        NR_Cat = data.N_Rels\n",
    "     \n",
    "    ### setting RR multiplicative factor for atypical hyperplasia\n",
    "    R_Hyp = np.repeat(np.nan, data.shape[0])\n",
    "    R_Hyp[NB_Cat == 0] = 1.00\n",
    "    R_Hyp[((NB_Cat != -100) & (NB_Cat > 0)) & (data.HypPlas == 0)] = 0.93\n",
    "    R_Hyp[((NB_Cat != -100) & (NB_Cat > 0)) & (data.HypPlas == 1)] = 1.82\n",
    "    R_Hyp[((NB_Cat != -100) & (NB_Cat > 0)) & (data.HypPlas == 99)] = 1.00\n",
    "    \n",
    "    set_HyperP_missing = data.HypPlas.values\n",
    "    set_R_Hyp_missing = R_Hyp.copy()\n",
    "    set_HyperP_missing[NB_Cat == -100] = -100\n",
    "    set_R_Hyp_missing[NB_Cat == -100] = -100\n",
    "    set_HyperP_missing[NB_Cat == -200] = -200\n",
    "    set_R_Hyp_missing[NB_Cat == -200] = -200\n",
    " \n",
    "    set_Race_missing = data.Race.values\n",
    "    Race_range=np.array(range(1,12))\n",
    "    set_Race_missing[-data.Race.isin(Race_range)]=-300\n",
    " \n",
    "    Error_Ind[(np.isnan(NB_Cat)) | (np.isnan(AM_Cat)) | (np.isnan(AF_Cat)) | (np.isnan(NR_Cat)) | (set_Race_missing == -300)] = 1\n",
    " \n",
    "    ### african-american RR model from CARE study:(1) eliminates Age1st from model;(2) groups AgeMen=2 with AgeMen=1;\n",
    "    ## setting AF_Cat=0 eliminates Age1st and its interaction from RR model;\n",
    "    AF_Cat[data.Race == 2] = 0 \n",
    "    ## group AgeMen RR level 2 with 1;\n",
    "    AM_Cat[(data.Race == 2) & (AM_Cat ==2)] = 1 \n",
    "\n",
    "    \n",
    "    ### hispanic RR model from San Francisco Bay Area Breast Cancer Study (SFBCS):\n",
    "    ###         (1) groups N_Biop ge 2 with N_Biop eq 1\n",
    "    ###         (2) eliminates  AgeMen from model for US Born hispanic women\n",
    "    ###         (3) group Age1st=25-29 with Age1st=20-24 and code as 1\n",
    "    ###             for   Age1st=30+, 98 (nulliparous)       code as 2\n",
    "    ###         (4) groups N_Rels=2 with N_Rels=1;\n",
    "    NB_Cat[(data.Race.isin([3,5])) & (data.N_Biop.isin([0,99]))] = 0\n",
    "    NB_Cat[(data.Race.isin([3,5])) & (NB_Cat==2)] = 1\n",
    "    AM_Cat[data.Race==3] = 0\n",
    "   \n",
    "    AF_Cat[(data.Race.isin([3,5])) & (data.Age1st!=98) & (AF_Cat==2)] = 1\n",
    "    AF_Cat[(data.Race.isin([3,5])) & (AF_Cat==3)] = 2\n",
    "    NR_Cat[(data.Race.isin([3,5])) & (NR_Cat == 2)] = 1\n",
    "    \n",
    " \n",
    "    ### for asian-americans NR_Cat=2 is pooled with NR_Cat=1; \n",
    "    NR_Cat[(data.Race >= 6) & (data.Race <= 11) & (NR_Cat == 2)] = 1\n",
    "\n",
    "    CharRace = np.repeat('??', data.shape[0])\n",
    "    CharRace[data.Race == 1] = \"Wh\"      #white SEER 1983:87 BrCa Rate\n",
    "    CharRace[data.Race == 2] = \"AA\"      #african-american\n",
    "    CharRace[data.Race == 3] = \"HU\"      #hispanic-american (US born)\n",
    "    CharRace[data.Race == 4] = \"NA\"      #other (native american and unknown race)\n",
    "    CharRace[data.Race == 5] = \"HF\"      #hispanic-american (foreign born)\n",
    "    CharRace[data.Race == 6] = \"Ch\"      #chinese\n",
    "    CharRace[data.Race == 7] = \"Ja\"      #japanese\n",
    "    CharRace[data.Race == 8] = \"Fi\"      #filipino\n",
    "    CharRace[data.Race == 9] = \"Hw\"      #hawaiian\n",
    "    CharRace[data.Race == 10] = \"oP\"     #other pacific islander\n",
    "    CharRace[data.Race == 11] = \"oA\"     #other asian\n",
    "\n",
    "\n",
    "\n",
    "#     recode_check= cbind(Error_Ind, set_T1_missing, set_T2_missing, NB_Cat, AM_Cat, AF_Cat, NR_Cat, R_Hyp, set_HyperP_missing, set_R_Hyp_missing, set_Race_missing, CharRace)\n",
    "    recode_check = pd.DataFrame({'Error_Ind': Error_Ind, 'set_T1_missing':set_T1_missing, \n",
    "                                 'set_T2_missing':set_T2_missing, 'NB_Cat':NB_Cat, \n",
    "                                 'AM_Cat':AM_Cat, 'AF_Cat':AF_Cat, 'NR_Cat':NR_Cat, \n",
    "                                 'R_Hyp':R_Hyp, 'set_HyperP_missing':set_HyperP_missing, \n",
    "                                 'set_R_Hyp_missing':set_R_Hyp_missing, \n",
    "                                 'set_Race_missing':set_Race_missing, 'CharRace':CharRace})\n",
    "    return(recode_check)\n",
    "\n",
    "r_ch = recode_check(data, Raw_Ind=1)\n",
    "r_ch"
   ]
  },
  {
   "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>RR_Star1</th>\n",
       "      <th>RR_Star2</th>\n",
       "      <th>PatternNumber</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>29.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>1.489622</td>\n",
       "      <td>1.489622</td>\n",
       "      <td>5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.492604</td>\n",
       "      <td>4.117968</td>\n",
       "      <td>65.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>4.004947</td>\n",
       "      <td>4.004947</td>\n",
       "      <td>38.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>2.207490</td>\n",
       "      <td>2.207490</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>6.981955</td>\n",
       "      <td>6.981955</td>\n",
       "      <td>38.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>3.567702</td>\n",
       "      <td>3.567702</td>\n",
       "      <td>38.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>10</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>11</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>12</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>29.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>13</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>14</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>15</th>\n",
       "      <td>2.345782</td>\n",
       "      <td>2.097356</td>\n",
       "      <td>50.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>16</th>\n",
       "      <td>1.389460</td>\n",
       "      <td>1.389460</td>\n",
       "      <td>40.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>17</th>\n",
       "      <td>5.385991</td>\n",
       "      <td>3.027437</td>\n",
       "      <td>103.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>18</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>19</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>20</th>\n",
       "      <td>1.421020</td>\n",
       "      <td>1.421020</td>\n",
       "      <td>16.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>21</th>\n",
       "      <td>1.421020</td>\n",
       "      <td>1.421020</td>\n",
       "      <td>16.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>22</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>23</th>\n",
       "      <td>1.421020</td>\n",
       "      <td>1.421020</td>\n",
       "      <td>16.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>24</th>\n",
       "      <td>1.421020</td>\n",
       "      <td>1.421020</td>\n",
       "      <td>16.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25</th>\n",
       "      <td>NaN</td>\n",
       "      <td>NaN</td>\n",
       "      <td>16.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "    RR_Star1  RR_Star2  PatternNumber\n",
       "0        NaN       NaN           29.0\n",
       "1        NaN       NaN            NaN\n",
       "2        NaN       NaN            NaN\n",
       "3   1.489622  1.489622            5.0\n",
       "4   5.492604  4.117968           65.0\n",
       "5   4.004947  4.004947           38.0\n",
       "6   2.207490  2.207490            2.0\n",
       "7   6.981955  6.981955           38.0\n",
       "8        NaN       NaN            NaN\n",
       "9   3.567702  3.567702           38.0\n",
       "10       NaN       NaN            NaN\n",
       "11       NaN       NaN            NaN\n",
       "12       NaN       NaN           29.0\n",
       "13       NaN       NaN            NaN\n",
       "14       NaN       NaN            NaN\n",
       "15  2.345782  2.097356           50.0\n",
       "16  1.389460  1.389460           40.0\n",
       "17  5.385991  3.027437          103.0\n",
       "18       NaN       NaN            NaN\n",
       "19       NaN       NaN            NaN\n",
       "20  1.421020  1.421020           16.0\n",
       "21  1.421020  1.421020           16.0\n",
       "22       NaN       NaN            NaN\n",
       "23  1.421020  1.421020           16.0\n",
       "24  1.421020  1.421020           16.0\n",
       "25       NaN       NaN           16.0"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def relative_risk(data,Raw_Ind=1):\n",
    "    ## LN_RR, beta=lnRR, beta for NB, AM, AF, NR, AC*NB and AF*NR, from Gail/CARE model\n",
    "    White_Beta = np.array([[0.5292641686, 0.0940103059, \n",
    "                            0.2186262218, 0.9583027845, \n",
    "                            -0.2880424830, -0.1908113865]])\n",
    "    Black_Beta = np.array([[0.1822121131, 0.2672530336, \n",
    "                            0.0, 0.4757242578, -0.1119411682, 0.0]])\n",
    "    Hspnc_Beta = np.array([[0.0970783641, 0.0000000000, \n",
    "                            0.2318368334, 0.166685441, \n",
    "                            0.0000000000, 0.0000000000]])\n",
    "    FHspnc_Beta = np.array([[0.4798624017, 0.2593922322, \n",
    "                             0.4669246218, 0.9076679727, \n",
    "                             0.0000000000, 0.0000000000]])\n",
    "    Other_Beta = np.array([[0.5292641686, 0.0940103059, \n",
    "                            0.2186262218, 0.9583027845, \n",
    "                            -0.2880424830, -0.1908113865]])\n",
    "    Asian_Beta = np.array([[0.55263612260619, 0.07499257592975, \n",
    "                            0.27638268294593, 0.79185633720481, \n",
    "                            0.0, 0.0]])\n",
    "    Wrk_Beta_all = np.concatenate((White_Beta, Black_Beta, Hspnc_Beta, \n",
    "                                   Other_Beta, FHspnc_Beta, Asian_Beta, \n",
    "                                   Asian_Beta, Asian_Beta, Asian_Beta, \n",
    "                                   Asian_Beta, Asian_Beta))\n",
    "\n",
    "    ### define LP1 = Linear Predictor for woman of interest at ages < 50; LP2 = Linear Predictor for woman of interest at ages >= 50\n",
    "    LP1 = np.repeat(np.nan, data.shape[0])\n",
    "    LP2 = np.repeat(np.nan, data.shape[0])\n",
    "   \n",
    "    ### obtain covariates\n",
    "    check_cov=recode_check(data, Raw_Ind)\n",
    "\n",
    "    NB_Cat=check_cov.NB_Cat.values\n",
    "    \n",
    "    NB_Cat[(NB_Cat==-100) | (NB_Cat==-200)] = np.nan\n",
    "\n",
    "    AM_Cat    = check_cov.AM_Cat.values\n",
    "    AF_Cat    = check_cov.AF_Cat.values\n",
    "    NR_Cat    = check_cov.NR_Cat.values\n",
    "    R_Hyp     = check_cov.R_Hyp.values\n",
    "    CharRace  = check_cov.CharRace.values\n",
    "\n",
    "    ### define pattern number when NB_Cat, AM_Cat, AF_Cat, NR_Cat are meaningful\n",
    "    ### NB_Cat(3 levels), AM_Cat(3 levels), AF_Cat(4 levels), NR_Cat(3 levels), 3*3*4*3 = 108 patterns in total\n",
    "    ## let PNID be the ID numbers when all \"_Cat\" variables are numerical\n",
    "    PatternNumber = np.repeat(np.nan, data.shape[0])\n",
    "    PNID = np.argwhere((~np.isnan(NB_Cat)) & (~np.isnan(AM_Cat)) & (~np.isnan(AF_Cat)) & (~np.isnan(NR_Cat))).T[0]\n",
    "    PatternNumber[PNID] = NB_Cat[PNID]*36+AM_Cat[PNID]*12+AF_Cat[PNID]*3+NR_Cat[PNID]*1+1\n",
    "    for i in PNID:\n",
    "        if CharRace[i]!=\"??\":\n",
    "            Beta = Wrk_Beta_all[data.Race[i-1]]\n",
    "\n",
    "            # for woman at ages < 50\n",
    "            LP1[i] = NB_Cat[i]*Beta[0]+AM_Cat[i]*Beta[1]+AF_Cat[i]*Beta[2]+NR_Cat[i]*Beta[3]+AF_Cat[i]*NR_Cat[i]*Beta[5]+np.log(R_Hyp[i])\n",
    "            LP2[i] = LP1[i]+NB_Cat[i]*Beta[4]\n",
    "\n",
    "    ### define RR_Star1 = relative risk for woman of interest at ages < 50; RR_Star2 = relative risk for woman of interest at ages >= 50\n",
    "    RR_Star1 = np.exp(LP1)\n",
    "    RR_Star2 = np.exp(LP2)  \n",
    "    RR_Star = pd.DataFrame({'RR_Star1': RR_Star1, 'RR_Star2': RR_Star2, 'PatternNumber': PatternNumber}) \n",
    "    return(RR_Star)\n",
    "                                    \n",
    "relative_risk(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([       nan,        nan,        nan, 1.86859688, 4.44129446,\n",
       "       1.76935628, 1.24955047, 5.7756759 ,        nan, 3.90606306,\n",
       "              nan,        nan,        nan,        nan,        nan,\n",
       "       2.68990404, 0.31628271, 1.02298019,        nan,        nan,\n",
       "       8.82766887, 6.76784804,        nan,        nan,        nan,\n",
       "              nan])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def absolute_risk(data, Raw_Ind=1, Avg_White=0):\n",
    "    ### set up lambda1*, lambda2, beta & F(t) with known constants used in the nci brca risk disk\n",
    "    ## lambda1_Star, BrCa composite incidences\n",
    "    # SEER BrCa incidence rates (current) non-hispanic white women, SEER white 1983:87\n",
    "    White_lambda1 = np.array([[0.00001000, 0.00007600, 0.00026600, 0.00066100, 0.00126500, 0.00186600, 0.00221100, \n",
    "                     0.00272100, 0.00334800, 0.00392300, 0.00417800, 0.00443900, 0.00442100, 0.00410900]])\n",
    "    # SEER BrCa incidence rates for \"avg\" non-hispanic white women and \"avg\" other (native american) women, SEER white 1992:96\n",
    "    White_lambda1Avg = np.array([[0.00001220, 0.00007410, 0.00022970, 0.00056490, 0.00116450, 0.00195250, 0.00261540, \n",
    "                        0.00302790, 0.00367570, 0.00420290, 0.00473080, 0.00494250, 0.00479760, 0.00401060]])\n",
    "    # SEER BrCa indicdence rates (under study) for non-hispanic white women, SEER white 1995:2003\n",
    "    White_nlambda1 = np.array([[0.0000120469, 0.0000746893, 0.0002437767, 0.0005878291, 0.0012069622, 0.0019762053, 0.0026200977, \n",
    "                      0.0033401788, 0.0039743676, 0.0044875763, 0.0048945499, 0.0051610641, 0.0048268456, 0.0040407389]])\n",
    "    # SEER black 1994-98\n",
    "    Black_lambda1 = np.array([[0.00002696, 0.00011295, 0.00031094, 0.00067639, 0.00119444, 0.00187394, 0.00241504, \n",
    "                       0.00291112, 0.00310127, 0.00366560, 0.00393132, 0.00408951, 0.00396793, 0.00363712]])\n",
    "    # SEER Ca Hisp 1995-2004\n",
    "    Hspnc_lambda1 = np.array([[0.0000166, 0.0000741, 0.0002740, 0.0006099, 0.0012225, 0.0019027, 0.0023142, \n",
    "                       0.0028357, 0.0031144, 0.0030794, 0.0033344, 0.0035082, 0.0025308, 0.0020414]])\n",
    "    # SEER white 1983:87\n",
    "    Other_lambda1 = np.array([[0.00001000, 0.00007600, 0.00026600, 0.00066100, 0.00126500, 0.00186600, 0.00221100, \n",
    "                       0.00272100, 0.00334800, 0.00392300, 0.00417800, 0.00443900, 0.00442100, 0.00410900]])\n",
    "    # SEER Ca Hisp 1995-2004\n",
    "    FHspnc_lambda1 = np.array([[0.0000102, 0.0000531, 0.0001578, 0.0003602, 0.0007617, 0.0011599, 0.0014111,\n",
    "                        0.0017245,  0.0020619, 0.0023603, 0.0025575, 0.0028227, 0.0028295, 0.0025868]])\n",
    "    # seer18 chinese  1998:02\n",
    "    Chnes_lambda1 = np.array([[0.000004059636, 0.000045944465, 0.000188279352, 0.000492930493, 0.000913603501,\n",
    "                       0.001471537353, 0.001421275482, 0.001970946494, 0.001674745804, 0.001821581075,\n",
    "                       0.001834477198, 0.001919911972, 0.002233371071, 0.002247315779]])\n",
    "    # seer18 japanese 1998:02\n",
    "    Japns_lambda1 = np.array([[0.000000000001, 0.000099483924, 0.000287041681, 0.000545285759, 0.001152211095,\n",
    "                       0.001859245108, 0.002606291272, 0.003221751682, 0.004006961859, 0.003521715275,\n",
    "                       0.003593038294, 0.003589303081, 0.003538507159, 0.002051572909]])\n",
    "    # seer18 filipino 1998:02\n",
    "    Filip_lambda1 = np.array([[0.000007500161, 0.000081073945, 0.000227492565, 0.000549786433, 0.001129400541,\n",
    "                       0.001813873795, 0.002223665639, 0.002680309266, 0.002891219230, 0.002534421279,\n",
    "                       0.002457159409, 0.002286616920, 0.001814802825, 0.001750879130]])\n",
    "    # seer18 hawaiian 1998:02\n",
    "    Hawai_lambda1 = np.array([[0.000045080582, 0.000098570724, 0.000339970860, 0.000852591429, 0.001668562761,\n",
    "                       0.002552703284, 0.003321774046, 0.005373001776, 0.005237808549, 0.005581732512,\n",
    "                       0.005677419355, 0.006513409962, 0.003889457523, 0.002949061662]])\n",
    "    # seer18 otr pac isl 1998:02\n",
    "    OtrPI_lambda1 = np.array([[0.000000000001, 0.000071525212, 0.000288799028, 0.000602250698, 0.000755579402,\n",
    "                       0.000766406354, 0.001893124938, 0.002365580107, 0.002843933070, 0.002920921732,\n",
    "                       0.002330395655, 0.002036291235, 0.001482683983, 0.001012248203]])\n",
    "    # seer18 otr asian 1998:02\n",
    "    OtrAs_lambda1 = np.array([[0.000012355409, 0.000059526456, 0.000184320831, 0.000454677273, 0.000791265338,\n",
    "                       0.001048462801, 0.001372467817, 0.001495473711, 0.001646746198, 0.001478363563,\n",
    "                       0.001216010125, 0.001067663700, 0.001376104012, 0.000661576644]])\n",
    "    ## lambda2, Competing hazards\n",
    "    #nchs competing mortality (current) for non-hispanic white women, NCHS white 1985:87\n",
    "    White_lambda2 = np.array([[0.00049300, 0.00053100, 0.00062500, 0.00082500, 0.00130700, 0.00218100, 0.00365500, \n",
    "                       0.00585200, 0.00943900, 0.01502800, 0.02383900, 0.03883200, 0.06682800, 0.14490800]])\n",
    "    # nchs competing mortality for \"avg\" non-hispanic white women and \"avg\" other (native american) women, NCHS white 1992:96\n",
    "    White_lambda2Avg = np.array([[0.00044120, 0.00052540, 0.00067460, 0.00090920, 0.00125340, 0.00195700, 0.00329840, \n",
    "                          0.00546220, 0.00910350, 0.01418540, 0.02259350, 0.03611460, 0.06136260, 0.14206630]])\n",
    "    # nchs competing mortality (under study) for non-hispanic white women, NCHS white 1995:2003\n",
    "    White_nlambda2 = np.array([[0.0004000377, 0.0004280396, 0.0005656742, 0.0008474486, 0.0012752947, 0.0018601059, 0.0028780622, \n",
    "                        0.0046903348, 0.0078835252, 0.0127434461, 0.0208586233, 0.0335901145, 0.0575791439, 0.1377327125]])\n",
    "    # NCHS black 1996-00\n",
    "    Black_lambda2 = np.array([[0.00074354, 0.00101698, 0.00145937, 0.00215933, 0.00315077, 0.00448779, 0.00632281, \n",
    "                       0.00963037, 0.01471818, 0.02116304, 0.03266035, 0.04564087, 0.06835185, 0.13271262]])\n",
    "    # SEER Ca Hisp 1995-2004\n",
    "    Hspnc_lambda2 = np.array([[0.0003561, 0.0004038, 0.0005281, 0.0008875, 0.0013987, 0.0020769, 0.0030912,\n",
    "                       0.0046960, 0.0076050, 0.0120555, 0.0193805, 0.0288386, 0.0429634, 0.0740349]])                \n",
    "    # NCHS white 1985:87\n",
    "    Other_lambda2 = np.array([[0.00049300, 0.00053100, 0.00062500, 0.00082500, 0.00130700, 0.00218100, 0.00365500, \n",
    "                       0.00585200, 0.00943900, 0.01502800, 0.02383900, 0.03883200, 0.06682800, 0.14490800]])\n",
    "    # SEER Ca Hisp 1995-2004\n",
    "    FHspnc_lambda2 = np.array([[0.0003129, 0.0002908, 0.0003515, 0.0004943, 0.0007807, 0.0012840, 0.0020325,\n",
    "                        0.0034533, 0.0058674, 0.0096888, 0.0154429, 0.0254675, 0.0448037, 0.1125678]])\n",
    "    # NCHS mortality chinese  1998:02\n",
    "    Chnes_lambda2 = np.array([[0.000210649076, 0.000192644865, 0.000244435215, 0.000317895949, 0.000473261994,\n",
    "                       0.000800271380, 0.001217480226, 0.002099836508, 0.003436889186, 0.006097405623,\n",
    "                       0.010664526765, 0.020148678452, 0.037990796590, 0.098333900733]])\n",
    "    # NCHS mortality japanese 1998:02\n",
    "    Japns_lambda2 = np.array([[0.000173593803, 0.000295805882, 0.000228322534, 0.000363242389, 0.000590633044,\n",
    "                       0.001086079485, 0.001859999966, 0.003216600974, 0.004719402141, 0.008535331402,\n",
    "                       0.012433511681, 0.020230197885, 0.037725498348, 0.106149118663]])\n",
    "    # NCHS mortality filipino 1998:02\n",
    "    Filip_lambda2 = np.array([[0.000229120979, 0.000262988494, 0.000314844090, 0.000394471908, 0.000647622610,\n",
    "                       0.001170202327, 0.001809380379, 0.002614170568, 0.004483330681, 0.007393665092,\n",
    "                       0.012233059675, 0.021127058106, 0.037936954809, 0.085138518334]])\n",
    "    # NCHS mortality hawaiian 1998:02\n",
    "    Hawai_lambda2 = np.array([[0.000563507269, 0.000369640217, 0.001019912579, 0.001234013911, 0.002098344078,\n",
    "                       0.002982934175, 0.005402445702, 0.009591474245, 0.016315472607, 0.020152229069,\n",
    "                       0.027354838710, 0.050446998723, 0.072262026612, 0.145844504021]])\n",
    "    # NCHS mortality otr pac isl 1998:02\n",
    "    OtrPI_lambda2 = np.array([[0.000465500812, 0.000600466920, 0.000851057138, 0.001478265376, 0.001931486788,\n",
    "                       0.003866623959, 0.004924932309, 0.008177071806, 0.008638202890, 0.018974658371,\n",
    "                       0.029257567105, 0.038408980974, 0.052869579345, 0.074745721133]])\n",
    "    # NCHS mortality otr asian 1998:02\n",
    "    OtrAs_lambda2 = np.array([[0.000212632332, 0.000242170741, 0.000301552711, 0.000369053354, 0.000543002943,\n",
    "                       0.000893862331, 0.001515172239, 0.002574669551, 0.004324370426, 0.007419621918,\n",
    "                       0.013251765130, 0.022291427490, 0.041746550635, 0.087485802065]])\n",
    "    # F(t), 1-Attributable Risk=F(t) \n",
    "    White_1_AR = np.array([[0.5788413, 0.5788413]])\n",
    "    Black_1_AR = np.array([[0.72949880, 0.74397137]])\n",
    "    Hspnc_1_AR = np.array([[0.749294788397, 0.778215491668]])\n",
    "    Other_1_AR = np.array([[0.5788413, 0.5788413]])\n",
    "    FHspnc_1_AR = np.array([[0.428864989813, 0.450352338746]])\n",
    "    Asian_1_AR = np.array([[0.47519806426735, 0.50316401683903]])\n",
    "\n",
    "    # intialize \"avg white women\" and \"avg\" other (native american women) rate for each year in the 5yr age cat\n",
    "    Avg_lambda1 = np.zeros((14, 5))\n",
    "    Avg_lambda2 = np.zeros((14, 5))\n",
    "    for i in range(Avg_lambda1.shape[1]):\n",
    "        Avg_lambda1[:,i] = White_lambda1Avg\n",
    "        Avg_lambda2[:,i] = White_lambda2Avg\n",
    "    # initialize rate vectors with the correct rates for each woman under study based on her race\n",
    "    # for i=1 to 11, when Race=i, Wrk_lambda1=Wrk_lambda1_all[i], Wrk_lambda2=Wrk_lambda2_all[i], Wrk_Beta=Wrk_Beta_all[i], Wrk_1_AR=Wrk_1_AR_all[i]\n",
    "    Wrk_lambda1_all = np.concatenate((White_lambda1, Black_lambda1, Hspnc_lambda1, Other_lambda1, FHspnc_lambda1, Chnes_lambda1, Japns_lambda1, Filip_lambda1, Hawai_lambda1, OtrPI_lambda1, OtrAs_lambda1))\n",
    "    Wrk_lambda2_all = np.concatenate((White_lambda2, Black_lambda2, Hspnc_lambda2, Other_lambda2, FHspnc_lambda2, Chnes_lambda2, Japns_lambda2, Filip_lambda2, Hawai_lambda2, OtrPI_lambda2, OtrAs_lambda2)) \n",
    "    Wrk_1_AR_all = np.concatenate((White_1_AR, Black_1_AR, Hspnc_1_AR, Other_1_AR, FHspnc_1_AR, Asian_1_AR, Asian_1_AR, Asian_1_AR, Asian_1_AR, Asian_1_AR, Asian_1_AR))\n",
    "    AbsRisk = np.repeat(np.nan,data.shape[0])\n",
    "    ## obtain IDs without any error\n",
    "    check_cov = recode_check(data, Raw_Ind)\n",
    "\n",
    "    Error_Ind = check_cov.Error_Ind.values.copy()\n",
    "    IDwoERR = np.argwhere(Error_Ind==0).T[0]\n",
    "    for i in IDwoERR:\n",
    "        obs = data[data.index==i]\n",
    "        RR_Star = relative_risk(data,Raw_Ind)\n",
    "        rrstar1 = RR_Star.RR_Star1[i]\n",
    "        rrstar2 = RR_Star.RR_Star2[i]\n",
    "        One_AR_RR = np.repeat(np.nan, 70) \n",
    "        Strt_Intvl = int(np.floor(obs.T1)-20+1)\n",
    "        End_Intvl = int(np.ceil(obs.T2)-20+0)\n",
    "        NumbrIntvl = int(np.ceil(obs.T2)-np.floor(obs.T1))\n",
    "        RskWrk = 0\n",
    "        Cum_lambda = 0\n",
    "        lambda1_temp = np.zeros((14, 5))\n",
    "        lambda2_temp = np.zeros((14, 5))\n",
    "\n",
    "         ## calculate abs risk  \n",
    "        if Avg_White == 0:\n",
    "            One_AR1 = Wrk_1_AR_all[int(obs.Race.values)-1,0]\n",
    "            One_AR2 = Wrk_1_AR_all[int(obs.Race.values)-1,1]\n",
    "            \n",
    "            # (1-AR)*RR at ages < 50\n",
    "            One_AR_RR1 = One_AR1*rrstar1\n",
    "            # (1-AR)*RR at ages >= 50\n",
    "            One_AR_RR2 = One_AR2*rrstar2\n",
    "            # define One_AR_RR\n",
    "            One_AR_RR[0:30] = One_AR_RR1\n",
    "            One_AR_RR[30:70] = One_AR_RR2\n",
    "            for v in range(lambda1_temp.shape[1]):\n",
    "                lambda1_temp[:,v] = Wrk_lambda1_all[int(obs.Race.values)-1,:]\n",
    "                lambda2_temp[:,v] = Wrk_lambda2_all[int(obs.Race.values)-1,:]\n",
    "            lambda1 = lambda1_temp.flatten()\n",
    "            lambda2 = lambda2_temp.flatten()\n",
    "         # calculate avg abs risk\n",
    "        if Avg_White == 1:\n",
    "             # define One_AR_RR\n",
    "            One_AR_RR = np.repeat(1, 70)\n",
    "            for val in range(lambda1_temp.shape[1]):\n",
    "                lambda1_temp[:,val] = Wrk_lambda1_all[int(obs.Race.values)-1,:]\n",
    "                lambda2_temp[:,val] = Wrk_lambda2_all[int(obs.Race.values)-1,:]\n",
    "             # if Race=1, 4 or 5, lambda1.temp[race,1:5]=Avg_lambda1[race,],lambda2.temp[race,1:5]=Avg_lambda2[race,]\n",
    "            if obs.Race==1 or obs.Race==4:\n",
    "                lambda1_temp = Avg_lambda1\n",
    "                lambda2_temp = Avg_lambda2\n",
    "\n",
    "            lambda1 = lambda1_temp.flatten()\n",
    "            lambda2 = lambda1_temp.flatten()\n",
    "        for j in range(NumbrIntvl):\n",
    "            j_intvl = Strt_Intvl+j-1\n",
    "            if NumbrIntvl>1 and j>0 and j<NumbrIntvl-1:\n",
    "                IntgrlLngth = 1\n",
    "\n",
    "            if NumbrIntvl>1 and j==0:\n",
    "                IntgrlLngth = 1-float((obs.T1-np.floor(obs.T1))) \n",
    "            if NumbrIntvl>1 and j+1==NumbrIntvl:\n",
    "                z1 = np.where(obs.T2>np.floor(obs.T2).values[0], 1, 0)\n",
    "                z2 = np.where(obs.T2==np.floor(obs.T2).values[0], 1, 0)\n",
    "                IntgrlLngth = (float((obs.T2-np.floor(obs.T2)))*z1+z2)[0] \n",
    "\n",
    "            if NumbrIntvl==0:\n",
    "                IntgrlLngth = (obs.T2-obs.T1).values[0]\n",
    "            lambdaj = lambda1[j_intvl]*One_AR_RR[j_intvl]+lambda2[j_intvl]\n",
    "            \n",
    "            PI_j = ((One_AR_RR[j_intvl]*lambda1[j_intvl]/lambdaj)*np.exp(-Cum_lambda))*(1-np.exp(-lambdaj*IntgrlLngth))\n",
    "            RskWrk = RskWrk+PI_j \n",
    "            Cum_lambda = Cum_lambda+lambdaj*IntgrlLngth      \n",
    "        AbsRisk[i] = 100*RskWrk    \n",
    "     \n",
    "    return(AbsRisk)\n",
    "absolute_risk(data)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "evamodels",
   "language": "python",
   "name": "evamodels"
  },
  "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
}