2_calculate_baselines.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compute age-specific baselines for the US data using the Karlinsky & Kobak, 2021 model\n",
    "Adapted from https://github.com/dkobak/excess-mortality/blob/main/baselines-stmf.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-21T22:19:37.974585Z",
     "start_time": "2022-04-21T22:19:35.610718Z"
    }
   },
   "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>Jurisdiction</th>\n",
       "      <th>Year</th>\n",
       "      <th>Age</th>\n",
       "      <th>Week</th>\n",
       "      <th>Number of Deaths</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Alabama</td>\n",
       "      <td>2015</td>\n",
       "      <td>&lt;65</td>\n",
       "      <td>1</td>\n",
       "      <td>345</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Alabama</td>\n",
       "      <td>2015</td>\n",
       "      <td>&lt;65</td>\n",
       "      <td>2</td>\n",
       "      <td>333</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Alabama</td>\n",
       "      <td>2015</td>\n",
       "      <td>&lt;65</td>\n",
       "      <td>3</td>\n",
       "      <td>303</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Alabama</td>\n",
       "      <td>2015</td>\n",
       "      <td>&lt;65</td>\n",
       "      <td>4</td>\n",
       "      <td>294</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Alabama</td>\n",
       "      <td>2015</td>\n",
       "      <td>&lt;65</td>\n",
       "      <td>5</td>\n",
       "      <td>316</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Jurisdiction  Year  Age  Week  Number of Deaths\n",
       "0      Alabama  2015  <65     1               345\n",
       "1      Alabama  2015  <65     2               333\n",
       "2      Alabama  2015  <65     3               303\n",
       "3      Alabama  2015  <65     4               294\n",
       "4      Alabama  2015  <65     5               316"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np, pandas as pd\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "df_stmf = pd.read_csv(\"data/tidy_deaths_age_state_CDC.csv\"\n",
    "                   ).rename(columns = {\"Age Group\":\"Age\"})\n",
    "\n",
    "df_stmf.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-21T22:19:37.996869Z",
     "start_time": "2022-04-21T22:19:37.978944Z"
    }
   },
   "outputs": [],
   "source": [
    "years_with_53_weeks = [2009, 2015, 2020, 2026]\n",
    "\n",
    "def predict(X):    \n",
    "    # Fit regression model on pre-2020 data from 2015 on \n",
    "    ind = (X[:,0] >= 2015) & (X[:,0] < 2020) & (X[:,1]<53)\n",
    "\n",
    "    m = int(np.max(X[ind,1]))\n",
    "    onehot = np.zeros((np.sum(ind), m))\n",
    "    for i,k in enumerate(X[ind,1]):\n",
    "        onehot[i,int(k)-1] = 1\n",
    "    predictors = np.concatenate((X[ind,:1], onehot), axis=1)\n",
    "    reg = LinearRegression(fit_intercept=False).fit(predictors, X[ind,2])\n",
    "        \n",
    "    # Compute 2020 baseline\n",
    "    ind2 = X[:,0] == 2020\n",
    "    predictors2020 = np.concatenate((np.ones((m,1))*2020, np.eye(m)), axis=1)\n",
    "    baseline = reg.predict(predictors2020)\n",
    "    \n",
    "    # Week 53 usually does not have enough data, so we'll use \n",
    "    # the same baseline value as for week 52\n",
    "    baseline = np.concatenate((baseline, [baseline[-1]]))\n",
    "    \n",
    "    df = pd.DataFrame()\n",
    "    \n",
    "    for Year in range(2015, 2023):\n",
    "        \n",
    "        Week = [i for i in range(1, 53)]\n",
    "        \n",
    "        predictors_year = np.concatenate((np.ones((m,1))*Year, np.eye(m)), axis=1)\n",
    "        baseline_year = reg.predict(predictors_year)\n",
    "        \n",
    "        df_tmp = pd.DataFrame(data = baseline_year, \n",
    "                              index = Week, \n",
    "                              columns = [\"Number of Deaths\"]) \n",
    "        df_tmp[\"Year\"] = Year\n",
    "        df_tmp = df_tmp.reset_index().rename(columns = {\"index\":\"Week\"})\n",
    "        \n",
    "        if Year in years_with_53_weeks:\n",
    "            # Week 53 again\n",
    "            last_row = df_tmp.iloc[-1].copy()\n",
    "            last_row[0] = 53\n",
    "            df_tmp = df_tmp.append(last_row)\n",
    "            \n",
    "        df = df.append(df_tmp)\n",
    "            \n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2022-04-21T22:20:02.210381Z",
     "start_time": "2022-04-21T22:19:38.000785Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "....................................................."
     ]
    },
    {
     "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>Week</th>\n",
       "      <th>Number of Deaths</th>\n",
       "      <th>Year</th>\n",
       "      <th>Jurisdiction</th>\n",
       "      <th>Age Group</th>\n",
       "      <th>Type</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1.0</td>\n",
       "      <td>329.123077</td>\n",
       "      <td>2015.0</td>\n",
       "      <td>Alabama</td>\n",
       "      <td>&lt;65</td>\n",
       "      <td>Predicted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2.0</td>\n",
       "      <td>333.923077</td>\n",
       "      <td>2015.0</td>\n",
       "      <td>Alabama</td>\n",
       "      <td>&lt;65</td>\n",
       "      <td>Predicted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3.0</td>\n",
       "      <td>310.923077</td>\n",
       "      <td>2015.0</td>\n",
       "      <td>Alabama</td>\n",
       "      <td>&lt;65</td>\n",
       "      <td>Predicted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4.0</td>\n",
       "      <td>312.123077</td>\n",
       "      <td>2015.0</td>\n",
       "      <td>Alabama</td>\n",
       "      <td>&lt;65</td>\n",
       "      <td>Predicted</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5.0</td>\n",
       "      <td>321.523077</td>\n",
       "      <td>2015.0</td>\n",
       "      <td>Alabama</td>\n",
       "      <td>&lt;65</td>\n",
       "      <td>Predicted</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Week  Number of Deaths    Year Jurisdiction Age Group       Type\n",
       "0   1.0        329.123077  2015.0      Alabama       <65  Predicted\n",
       "1   2.0        333.923077  2015.0      Alabama       <65  Predicted\n",
       "2   3.0        310.923077  2015.0      Alabama       <65  Predicted\n",
       "3   4.0        312.123077  2015.0      Alabama       <65  Predicted\n",
       "4   5.0        321.523077  2015.0      Alabama       <65  Predicted"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "countries = np.unique(df_stmf['Jurisdiction'])\n",
    "agebands = [\"<65\", \"65-74\", \"75-84\", \"85+\"]\n",
    "            \n",
    "df = pd.DataFrame()\n",
    "\n",
    "for country in countries:\n",
    "    print('.', end='')\n",
    "    for ageband in agebands:\n",
    "        X = df_stmf[(df_stmf['Jurisdiction']==country) & (df_stmf['Age']==ageband)]\n",
    "        X = X[['Year','Week', 'Number of Deaths']].values\n",
    "        baselines  = predict(X)\n",
    "        baselines[\"Jurisdiction\"] = country\n",
    "        baselines[\"Age Group\"] = ageband\n",
    "        df = df.append(baselines)\n",
    "df[\"Type\"] = \"Predicted\"\n",
    "order = [\"Jurisdiction\", \"Year\", \"Age Group\", \"Week\", \"Number of Deaths\", \"Type\"]\n",
    "df[order].set_index(\"Jurisdiction\").to_csv(\"data/tidy_predicted_deaths_age_state_CDC.csv\")\n",
    "\n",
    "df.head()"
   ]
  }
 ],
 "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"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}