3_predictive_models_2023_10_02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "8d266f77-a697-401f-a89b-5e4f94730ed0",
   "metadata": {},
   "source": [
    "# Predictive models for changes in clinical variables related to T2D\n",
    "\n",
    "Includes regression and classification models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37f75449-8f26-4b36-84ff-edf37cd925b2",
   "metadata": {},
   "source": [
    "## 1. Load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ad4dd48f-16f7-4010-9d62-21a3e2e6c502",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pickle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "e781afb3-ebed-41ad-bad4-e7a067d57f37",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/tmp/ipykernel_642/3589811407.py:1: DtypeWarning: Columns (992) have mixed types. Specify dtype option on import or set low_memory=False.\n",
      "  combined_table = pd.read_csv('results_2023_10_18/combined_data_table.csv', dtype={'public_client_id': str})\n"
     ]
    }
   ],
   "source": [
    "combined_table = pd.read_csv('results_2023_10_18/combined_data_table.csv', dtype={'public_client_id': str})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "b1600516-c899-4403-bd54-670bb5e4a218",
   "metadata": {},
   "outputs": [],
   "source": [
    "combined_table.index = pd.MultiIndex.from_frame(combined_table[['public_client_id', 'days_in_program']])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "40393229-1c96-470e-8de7-dbf17fa1e699",
   "metadata": {},
   "outputs": [],
   "source": [
    "selected_columns_full = np.loadtxt('results_2023_10_18/selected_columns_full.txt', dtype=str, delimiter='\\t').tolist()\n",
    "chem_subset_cols = np.loadtxt('results_2023_10_18/chem_subset_cols.txt', dtype=str, delimiter='\\t').tolist()\n",
    "selected_chem_bp_cols = np.loadtxt('results_2023_10_18/selected_chem_bp_cols.txt', dtype=str, delimiter='\\t').tolist()\n",
    "selected_prot_cols = np.loadtxt('results_2023_10_18/selected_prot_cols.txt', dtype=str, delimiter='\\t').tolist()\n",
    "selected_met_cols = np.loadtxt('results_2023_10_18/selected_met_cols.txt', dtype=str, delimiter='\\t').tolist()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "da5f3106-b4c2-41f6-99ff-cab062c968ed",
   "metadata": {},
   "outputs": [],
   "source": [
    "chems_selected = ['HbA1C', 'GFR', 'GLUCOSE', 'INSULIN', 'HOMA-IR']\n",
    "chems_to_column = {\n",
    "    'HbA1C': 'GLYCOHEMOGLOBIN A1C',\n",
    "    'GFR': 'GFR, MDRD',\n",
    "    'GLUCOSE': 'GLUCOSE',\n",
    "    'INSULIN': 'INSULIN',\n",
    "    'HOMA-IR': 'HOMA-IR'\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "822278c4-0a35-4339-ac44-19dc3f75e6aa",
   "metadata": {},
   "source": [
    "### Create X, y arrays from data matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4c5b54f9-7b9a-4328-b98e-2f04aef3f351",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_all = combined_table[selected_columns_full].to_numpy()\n",
    "X_selected = combined_table[chem_subset_cols].to_numpy()\n",
    "X_chems = combined_table[selected_chem_bp_cols].to_numpy()\n",
    "X_prots = combined_table[selected_prot_cols].to_numpy()\n",
    "X_mets = combined_table[selected_met_cols].to_numpy()\n",
    "#X_prs = combined_table[prs_cols].to_numpy()\n",
    "\n",
    "y_current = {c: combined_table[chems_to_column[c]].to_numpy() for c in chems_selected}\n",
    "y_next = {c: combined_table['next_' + c].to_numpy() for c in chems_selected}\n",
    "y_next_1y = {c: combined_table['next_1y_' + c].to_numpy() for c in chems_selected}\n",
    "y_next_2y = {c: combined_table['next_2y_' + c].to_numpy() for c in chems_selected}\n",
    "y_delta = {c: combined_table['d_' + c].to_numpy() for c in chems_selected}\n",
    "y_delta_1y = {c: combined_table['d_1y_' + c].to_numpy() for c in chems_selected}\n",
    "y_delta_2y = {c: combined_table['d_2y_' + c].to_numpy() for c in chems_selected}\n",
    "# TODO: also try to predict age/other clinical variables?\n",
    "y_misc = {'age': combined_table['age'].to_numpy(), 'bmi': combined_table['bmi'].to_numpy()}\n",
    "y_delta_class = {c: combined_table['d_' + c + '_class'].to_numpy() for c in chems_selected}\n",
    "y_delta_class_1y = {c: combined_table['d_1y_' + c + '_class'].to_numpy() for c in chems_selected}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "4f785bc3-69c3-4279-81b6-7b2d343eda94",
   "metadata": {},
   "outputs": [],
   "source": [
    "# TODO: demographic-only model\n",
    "X_basic = combined_table[['age', 'is_m', 'bmi']].to_numpy()\n",
    "X_prots_demographics = np.concatenate([X_prots, X_basic], 1)\n",
    "X_mets_demographics = np.concatenate([X_mets, X_basic], 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "d3eae3a6-6ca6-4cbd-a6d9-1ea9b5a93b84",
   "metadata": {},
   "outputs": [],
   "source": [
    "client_ids_subset = np.array([x[0] for x in combined_table.index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "2a8f72b7-a5d0-4a7b-b9cf-c593fe76eb9d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All: (1131, 1042)\n",
      "selected clinical: (1131, 14)\n",
      "full clinical: (1131, 70)\n",
      "prots: (1131, 262)\n",
      "mets: (1131, 710)\n"
     ]
    }
   ],
   "source": [
    "# print data sizes\n",
    "print('All:', X_all.shape)\n",
    "print('selected clinical:', X_selected.shape)\n",
    "print('full clinical:', X_chems.shape)\n",
    "print('prots:', X_prots.shape)\n",
    "print('mets:', X_mets.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e15fb04d-75bd-4e9c-898c-4d31c1632559",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HbA1C 6 months: 185.0\n",
      "GFR 6 months: 299.0\n",
      "GLUCOSE 6 months: 250.0\n",
      "INSULIN 6 months: 453.0\n",
      "HOMA-IR 6 months: 455.0\n",
      "HbA1C 1 year: 87.0\n",
      "GFR 1 year: 204.0\n",
      "GLUCOSE 1 year: 176.0\n",
      "INSULIN 1 year: 247.0\n",
      "HOMA-IR 1 year: 253.0\n"
     ]
    }
   ],
   "source": [
    "# print \n",
    "for c, data in y_delta_class.items():\n",
    "    print(c, '6 months:', data.sum())\n",
    "    \n",
    "for c, data in y_delta_class_1y.items():\n",
    "    print(c, '1 year:', data[~np.isnan(data)].sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4d146b8-1e5c-4503-a134-ae175549469d",
   "metadata": {
    "tags": []
   },
   "source": [
    "## 2. Running machine learning predictions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "dfe1bac5-332f-4028-bd12-5056279bad9b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import ElasticNet, ElasticNetCV, Lasso, LassoCV\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.model_selection import cross_val_score\n",
    "from sklearn.metrics import r2_score\n",
    "\n",
    "from sklearn.model_selection import GroupKFold\n",
    "\n",
    "cv = GroupKFold(10)\n",
    "\n",
    "# scale data\n",
    "scaler = StandardScaler()\n",
    "X_all_scaled = scaler.fit_transform(X_all)\n",
    "X_selected_scaled = scaler.fit_transform(X_selected)\n",
    "X_chems_scaled = scaler.fit_transform(X_chems)\n",
    "X_prots_scaled = scaler.fit_transform(X_prots)\n",
    "X_mets_scaled = scaler.fit_transform(X_mets)\n",
    "\n",
    "\n",
    "en = ElasticNet()\n",
    "\n",
    "analyses = []\n",
    "dependent_vars = []\n",
    "current_vars = []\n",
    "for k in y_current.keys():\n",
    "    analyses.append('d_' + k.lower())\n",
    "    dependent_vars.append(y_delta[k])\n",
    "    current_vars.append(y_current[k])\n",
    "    \n",
    "for k in y_current.keys():\n",
    "    analyses.append('next_' + k.lower())\n",
    "    dependent_vars.append(y_next[k])\n",
    "    current_vars.append(y_current[k])\n",
    "    \n",
    "# TODO: add current\n",
    "for k in y_current.keys():\n",
    "    analyses.append(k)\n",
    "    dependent_vars.append(y_current[k])\n",
    "    current_vars.append(np.zeros(y_current[k].shape))\n",
    "    \n",
    "for k in y_current.keys():\n",
    "    analyses.append('next_1y_' + k.lower())\n",
    "    dependent_vars.append(y_next_1y[k])\n",
    "    current_vars.append(y_current[k])\n",
    "    analyses.append('d_1y_' + k.lower())\n",
    "    dependent_vars.append(y_delta_1y[k])\n",
    "    current_vars.append(y_current[k])\n",
    "    analyses.append('next_2y_' + k.lower())\n",
    "    dependent_vars.append(y_next_2y[k])\n",
    "    current_vars.append(y_current[k])\n",
    "    analyses.append('d_2y_' + k.lower())\n",
    "    dependent_vars.append(y_delta_2y[k])\n",
    "    current_vars.append(y_current[k])\n",
    "    \n",
    "# vars with no \"current\" variable\n",
    "for k in y_misc.keys():\n",
    "    analyses.append(k.lower())\n",
    "    dependent_vars.append(y_misc[k])\n",
    "    current_vars.append(np.zeros(y_misc[k].shape))\n",
    "    \n",
    "class_analyses = []\n",
    "class_dependent_vars = []\n",
    "class_current_vars = []\n",
    "for k in y_current.keys():\n",
    "    class_analyses.append('d_class_' + k.lower())\n",
    "    class_dependent_vars.append(y_delta_class[k])\n",
    "    class_current_vars.append(y_current[k])\n",
    "    class_analyses.append('d_class_1y_' + k.lower())\n",
    "    class_dependent_vars.append(y_delta_class_1y[k])\n",
    "    class_current_vars.append(y_current[k])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "169e170c-8aa5-43cb-aee7-dcd526ab375a",
   "metadata": {},
   "source": [
    "### Analysis functions"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "451c9212-1603-4150-a300-7a780c74bb17",
   "metadata": {},
   "source": [
    "#### Random Forest grid search function for RFR (unused here)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "303eac40-34ff-497b-b7e4-6599ccca1fc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.model_selection import RandomizedSearchCV, GridSearchCV\n",
    "\n",
    "param_distributions = {\n",
    "    \"max_features\": [1, 2, 3, 5, None],\n",
    "    \"max_leaf_nodes\": [10, 100, 500, None],\n",
    "    \"min_samples_leaf\": [1, 2, 5, 10, 20, 50],\n",
    "    \"max_depth\": [15, None],\n",
    "}\n",
    "\n",
    "def rf_grid_search(model, X_train, y_train, scoring='r2'):\n",
    "    # see: https://inria.github.io/scikit-learn-mooc/python_scripts/ensemble_hyperparameters.html\n",
    "    search_cv = RandomizedSearchCV(model, param_distributions=param_distributions, scoring=scoring)\n",
    "    search_cv.fit(X_train, y_train)\n",
    "    print(search_cv.best_params_, search_cv.best_score_)\n",
    "    return search_cv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3ce9c36f-e79d-4e61-9c63-42304d3f2341",
   "metadata": {},
   "source": [
    "#### Continuous regression analysis function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "4256cb83-7cd3-41ac-a7c6-853543a9c7cb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'HbA1C': 'GLYCOHEMOGLOBIN A1C',\n",
       " 'GFR': 'GFR, MDRD',\n",
       " 'GLUCOSE': 'GLUCOSE',\n",
       " 'INSULIN': 'INSULIN',\n",
       " 'HOMA-IR': 'HOMA-IR'}"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chems_to_column"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "d45518c5-615b-4c9f-9bc1-fb2af6c5e0d0",
   "metadata": {},
   "outputs": [],
   "source": [
    "# set of chems to remove just for the comparison\n",
    "related_chems = {'GFR': ['CREATININE ENZ, SER', 'BUN/CREAT RATIO'],\n",
    "                 'HOMA-IR': ['INSULIN'],\n",
    "                 'age': ['is_m']}\n",
    "\n",
    "# if the value of duplicate_current_vals for a given key is in the labels, set \"current\" to 0. This is because\n",
    "# in these cases, \"current\" is duplicated by something already present in the data.\n",
    "duplicate_current_vals = {}\n",
    "for analysis in analyses + class_analyses:\n",
    "    if 'gfr' in analysis:\n",
    "        duplicate_current_vals[analysis] = chems_to_column['GFR']\n",
    "    if 'hba1c' in analysis:\n",
    "        duplicate_current_vals[analysis] = chems_to_column['HbA1C']\n",
    "    if 'glucose' in analysis:\n",
    "        duplicate_current_vals[analysis] = chems_to_column['GLUCOSE']\n",
    "    if 'insulin' in analysis:\n",
    "        duplicate_current_vals[analysis] = chems_to_column['INSULIN']\n",
    "    if 'homa' in analysis:\n",
    "        duplicate_current_vals[analysis] = chems_to_column['HOMA-IR']\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "d48799d7-9fc2-472c-9e57-4c9db50ab99e",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_analyses(X, analyses, dependent_vars, current_vars, \n",
    "                 client_ids_subset, model, return_params=False, labels=None,\n",
    "                 use_random_search=False,\n",
    "                 save_models=False,\n",
    "                 save_model_prefix=''):\n",
    "    \"\"\"\n",
    "    X - input data\n",
    "    analyses - list of names of the analyses\n",
    "    dependent_vars - list of arrays\n",
    "    current_vars - list of arrays\n",
    "    client_ids_subset - array of strings indicating client ids\n",
    "    model - a scikit-learn model\n",
    "    \n",
    "    return_params - return the model parameters for each CV run. only for linear models.\n",
    "    labels - list or array of labels for each of the covariates\n",
    "    use_random_search - whether or not to use random parameter search\n",
    "    \"\"\"\n",
    "    all_cv_scores = []\n",
    "    cv = GroupKFold(10)\n",
    "    weights_results = []\n",
    "    predicted_values = []\n",
    "    predicted_r2_values = []\n",
    "    scaler = StandardScaler()\n",
    "    for analysis, var, current in zip(analyses, dependent_vars, current_vars):\n",
    "        # TODO: in some circumstances, set current to 0 if current duplicates one of the clinical variables present in the data.\n",
    "        cv_score = []\n",
    "        current = current.reshape((X.shape[0], 1))\n",
    "        if analysis in duplicate_current_vals:\n",
    "            if duplicate_current_vals[analysis] in labels:\n",
    "                print('analysis contains duplicate value', duplicate_current_vals[analysis])\n",
    "                current = np.zeros(current.shape)\n",
    "        X_new = np.concatenate([X, current], 1)\n",
    "        included_subset = (~np.isnan(var))\n",
    "        X_new = X_new[included_subset]\n",
    "        client_ids_subset_t = client_ids_subset[included_subset]\n",
    "        var = var[included_subset]\n",
    "        # this is a hack to remove the current variable from the list of variables.\n",
    "        if analysis in labels:\n",
    "            age_index = labels.index(analysis)\n",
    "            X_new[:, age_index] = 0\n",
    "        if analysis in chems_to_column and chems_to_column[analysis] in labels:\n",
    "            cc_index = labels.index(chems_to_column[analysis])\n",
    "            X_new[:, cc_index] = 0\n",
    "        if analysis in related_chems:\n",
    "            for c in related_chems[analysis]:\n",
    "                if c in labels:\n",
    "                    cc_index = labels.index(c)\n",
    "                    X_new[:, cc_index] = 0\n",
    "        predicted_values_analysis = np.zeros(var.shape)\n",
    "        if return_params:\n",
    "            cv_score = []\n",
    "            cv_round = 0\n",
    "            for train_index, test_index in cv.split(X_new, var, groups=client_ids_subset_t):\n",
    "                # TODO: move scaling step here. don't scale the data beforehand.\n",
    "                X_train = X_new[train_index]\n",
    "                X_test = X_new[test_index]\n",
    "                X_train = scaler.fit_transform(X_train)\n",
    "                X_test = scaler.transform(X_test)\n",
    "                filename = save_model_prefix + '_' + analysis + '_fold' + str(cv_round) + '.pkl'\n",
    "                model_exists = os.path.exists(filename)\n",
    "                if model_exists:\n",
    "                    with open(filename, 'rb') as f:\n",
    "                        model = pickle.load(f)\n",
    "                else:\n",
    "                    if use_random_search:\n",
    "                        cv_model = rf_grid_search(model, X_train, var[train_index])\n",
    "                        model = cv_model.best_estimator_\n",
    "                    else:\n",
    "                        model.fit(X_train, var[train_index])\n",
    "                if hasattr(model, 'coef_'):\n",
    "                    weights = list(model.coef_)\n",
    "                elif hasattr(model, 'feature_importances_'):\n",
    "                    weights = list(model.feature_importances_)\n",
    "                else:\n",
    "                    weights = [0]*(X_new.shape[1])\n",
    "                if hasattr(model, 'intercept_'):\n",
    "                    weights.append(model.intercept_)\n",
    "                if hasattr(model, 'alpha'):\n",
    "                    weights.append(model.alpha)\n",
    "                elif hasattr(model, 'alpha_'):\n",
    "                    weights.append(model.alpha_)\n",
    "                if use_random_search:\n",
    "                    for k in sorted(param_distributions.keys()):\n",
    "                        weights.append(model.__dict__[k])\n",
    "                weights.append(analysis)\n",
    "                predictions = model.predict(X_test)\n",
    "                predicted_values_analysis[test_index] = predictions\n",
    "                score = r2_score(var[test_index], predictions)\n",
    "                weights.append(score)\n",
    "                cv_score.append(score)\n",
    "                weights_results.append(weights)\n",
    "                if save_models and not model_exists:\n",
    "                    with open(filename, 'wb') as f:\n",
    "                        pickle.dump(model, f)\n",
    "                cv_round += 1\n",
    "            all_cv_scores.append(cv_score)\n",
    "            predicted_r2_values.append(r2_score(var, predicted_values_analysis))\n",
    "            predicted_values.append(predicted_values_analysis)\n",
    "            print(analysis, cv_score, np.mean(cv_score))\n",
    "            print('total r^2: ', predicted_r2_values[-1])\n",
    "        else:\n",
    "            cv_score = cross_val_score(model, X_new, var, groups=client_ids_subset_t, cv=cv, scoring='r2')\n",
    "            all_cv_scores.append(cv_score)\n",
    "            print(analysis, cv_score, np.mean(cv_score))\n",
    "    if return_params:\n",
    "        columns = list(labels) + ['current']\n",
    "        if hasattr(model, 'intercept_'):\n",
    "            columns.append('intercept')\n",
    "        if hasattr(model, 'alpha') or hasattr(model, 'alpha_'):\n",
    "            columns.append('alpha')\n",
    "        if use_random_search:\n",
    "            for k in sorted(param_distributions.keys()):\n",
    "                columns.append(k)\n",
    "        df = pd.DataFrame(weights_results, columns=columns + ['target', 'r2_score'])\n",
    "        return all_cv_scores, df, predicted_values, predicted_r2_values\n",
    "    else:\n",
    "        return all_cv_scores"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "22e032fd-899c-4592-95e6-57fcb5653a82",
   "metadata": {},
   "source": [
    "#### Classification analysis function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "94d9c53f-e19b-4c1b-bb13-9002bb05384f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.metrics import balanced_accuracy_score, f1_score\n",
    "\n",
    "def run_classification_analyses(X, analyses, dependent_vars, current_vars, \n",
    "                                client_ids_subset, model, return_params=False, labels=None,\n",
    "                                use_random_search=False,\n",
    "                                save_models=False,\n",
    "                                save_model_prefix=''):\n",
    "    \"\"\"\n",
    "    X - input data\n",
    "    analyses - list of names of the analyses\n",
    "    dependent_vars - list of arrays\n",
    "    current_vars - list of arrays\n",
    "    client_ids_subset - array of strings indicating client ids\n",
    "    model - a scikit-learn model\n",
    "    \n",
    "    return_params - return the model parameters for each CV run. only for linear models.\n",
    "    labels - list or array of labels for each of the covariates\n",
    "    \"\"\"\n",
    "    all_ba_scores = []\n",
    "    all_f1_scores = []\n",
    "    cv = GroupKFold(10)\n",
    "    weights_results = []\n",
    "    predicted_values = []\n",
    "    predicted_ba_scores = []\n",
    "    predicted_f1_scores = []\n",
    "    scaler = StandardScaler()\n",
    "    for analysis, var, current in zip(analyses, dependent_vars, current_vars):\n",
    "        cv_score = []\n",
    "        f1_scores = []\n",
    "        current = current.reshape((X.shape[0], 1))\n",
    "        if analysis in duplicate_current_vals:\n",
    "            if duplicate_current_vals[analysis] in labels:\n",
    "                print('analysis contains duplicate value', duplicate_current_vals[analysis])\n",
    "                current = np.zeros(current.shape)\n",
    "        X_new = np.concatenate([X, current], 1)\n",
    "        included_subset = (~np.isnan(var))\n",
    "        X_new = X_new[included_subset]\n",
    "        client_ids_subset_t = client_ids_subset[included_subset]\n",
    "        var = var[included_subset]\n",
    "        predicted_values_analysis = np.zeros(var.shape)\n",
    "        if return_params:\n",
    "            cv_score = []\n",
    "            f1_scores = []\n",
    "            cv_round = 0\n",
    "            for train_index, test_index in cv.split(X_new, var, groups=client_ids_subset_t):\n",
    "                X_train = X_new[train_index]\n",
    "                X_test = X_new[test_index]\n",
    "                X_train = scaler.fit_transform(X_train)\n",
    "                X_test = scaler.transform(X_test)\n",
    "                filename = save_model_prefix + '_' + analysis + '_fold' + str(cv_round) + '.pkl'\n",
    "                model_exists = os.path.exists(filename)\n",
    "                if model_exists:\n",
    "                    with open(filename, 'rb') as f:\n",
    "                        model = pickle.load(f)\n",
    "                else:\n",
    "                    if use_random_search:\n",
    "                        cv_model = rf_grid_search(model, X_train, var[train_index])\n",
    "                        model = cv_model.best_estimator_\n",
    "                    else:\n",
    "                        model.fit(X_train, var[train_index])\n",
    "                if hasattr(model, 'coef_'):\n",
    "                    weights = list(model.coef_[0,:])\n",
    "                elif hasattr(model, 'feature_importances_'):\n",
    "                    weights = list(model.feature_importances_)\n",
    "                else:\n",
    "                    weights = [0]*(X_new.shape[1])\n",
    "                if hasattr(model, 'intercept_'):\n",
    "                    weights.append(model.intercept_[0])\n",
    "                if hasattr(model, 'C'):\n",
    "                    weights.append(model.C)\n",
    "                elif hasattr(model, 'C_'):\n",
    "                    weights.append(model.C_[0])\n",
    "                weights.append(analysis)\n",
    "                predictions = model.predict(X_test)\n",
    "                predicted_values_analysis[test_index] = predictions\n",
    "                # TODO: change scoring method\n",
    "                score = balanced_accuracy_score(var[test_index], predictions)\n",
    "                f1 = f1_score(var[test_index], predictions)\n",
    "                f1_scores.append(f1)\n",
    "                weights.append(score)\n",
    "                weights.append(f1)\n",
    "                cv_score.append(score)\n",
    "                weights_results.append(weights)\n",
    "                if save_models and not model_exists:\n",
    "                    with open(filename, 'wb') as f:\n",
    "                        pickle.dump(model, f)\n",
    "                cv_round += 1\n",
    "            all_ba_scores.append(cv_score)\n",
    "            all_f1_scores.append(f1_scores)\n",
    "            predicted_ba_scores.append(balanced_accuracy_score(var, predicted_values_analysis))\n",
    "            predicted_f1_scores.append(f1_score(var, predicted_values_analysis))\n",
    "            predicted_values.append(predicted_values_analysis)\n",
    "            print(analysis, cv_score, np.mean(cv_score))\n",
    "            print('total balanced accuracy score: ', predicted_ba_scores[-1])\n",
    "            print('total F1 score: ', predicted_f1_scores[-1])\n",
    "        else:\n",
    "            cv_score = cross_val_score(model, X_new, var, groups=client_ids_subset_t, cv=cv, scoring='balanced_accuracy')\n",
    "            all_ba_scores.append(cv_score)\n",
    "            print(analysis, cv_score, np.mean(cv_score))\n",
    "    if return_params:\n",
    "        columns = list(labels)\n",
    "        if current_vars is not None:\n",
    "            columns.append('current')\n",
    "        if hasattr(model, 'intercept_'):\n",
    "            columns.append('intercept')\n",
    "        if hasattr(model, 'C') or hasattr(model, 'C_'):\n",
    "            columns.append('C')\n",
    "        df = pd.DataFrame(weights_results, columns=columns + ['target', 'balanced_accuracy_score', 'f1_score'])\n",
    "        return all_ba_scores, all_f1_scores, df, predicted_values, predicted_ba_scores\n",
    "    else:\n",
    "        return all_ba_scores"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5e8aadb-fdff-4841-bc56-b06420dc115d",
   "metadata": {},
   "source": [
    "### Regression Baseline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "86ac5d81-8392-43ed-a8a8-9e66e077e381",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "d_hba1c [-0.05445994727608383, -0.008291923016248592, -0.062461465926812565, -0.05943257785749245, -0.07063992359121318, -0.0640839284585728, -0.008337132807920256, -0.02805086384734401, -0.00771379493024682, -0.02650168877900172] -0.038997324649093626 0.205063842812146\n",
      "d_gfr [-0.008003258938506397, -0.008942276903969537, -5.814048558194784e-05, -0.0555466511088174, -0.014459631040406151, -0.00885745115794867, -0.13308034545930392, -0.01986480104574495, -2.5999049590330614e-05, -0.0003323715121705817] -0.024917092670203988 0.10987133400396662\n",
      "d_glucose [-9.258528991207271e-05, -0.002184116732982355, -0.026471986159432692, -0.01145041199771235, -0.0015797269829183058, -0.023904448579273208, -9.542212361934865e-05, -0.0558965181078761, -0.001807452125111686, -0.005237981861864416] -0.012872064996070253 0.1515667418801425\n",
      "d_insulin [-0.024715603125468277, -0.05551536241182542, -0.05398258572225423, -0.04159558827715992, -0.011342549192969509, -0.0197321401509587, -0.02666352169326469, -0.02792100562305322, -0.06451815483035572, -0.027229644113604845] -0.03532161551409145 0.2537297149475929\n",
      "d_homa-ir [-0.023515746965139694, -0.04027744194100524, -0.06029160185247573, -0.028360187629598244, -0.009989623268417347, -0.025388987119657758, -0.03534253231815332, -0.03398163464437354, -0.06245133749622722, -0.01505082469681307] -0.033464991793186114 0.2805136830359546\n",
      "next_hba1c [0.5073047210300428, 0.29102770305870584, 0.4610732452146401, 0.5282893303971153, 0.35656885705101726, 0.5181167716084718, 0.4444733513353488, 0.2594099753866159, 0.37602882703777296, 0.4028189815706672] 0.4145111763690398 0.5504057153685651\n",
      "next_gfr [0.5646358378451821, 0.6042225859247137, 0.6041150463800409, 0.641794584308157, 0.5034353163683869, 0.5632701173348189, 0.6396720297257584, 0.5032002788859324, 0.5958473941578186, 0.6274947509603239] 0.5847687941891133 0.6402927140415403\n",
      "next_glucose [0.572346699525351, 0.4889918406859356, 0.689097086260509, 0.5689115820525639, 0.35238376838986873, 0.5336361529314549, 0.7102564102564102, 0.11780352301875785, 0.32917698567467013, 0.490003513960673] 0.4852607562756194 0.5811744627398714\n",
      "next_insulin [0.14134672336184861, 0.31034446511808156, 0.524717476340215, 0.47216935872669685, 0.48907587369949046, 0.25157719094510045, 0.5107134176917038, 0.17943987218682944, 0.42050490969018095, 0.2665624859174849] 0.35664517736776324 0.5481719386588193\n",
      "next_homa-ir [-0.1687704425617369, 0.28632912265254806, 0.5236539688764378, 0.6318983318181048, 0.41284480308920446, 0.05281429129054027, 0.41857990371336007, 0.12604440414150975, 0.4624747171152429, 0.48573905161153863] 0.323160815174675 0.5636822690713679\n",
      "HbA1C [-138.4974543977473, -222.9849403854152, -135.03851144137772, -141.1801851731108, -203.24013095363816, -129.58824912724648, -143.65076920782408, -201.2747268378518, -169.98338871915809, -214.45697687347118] -169.98953331168406 -0.005659616740942864\n",
      "GFR [-32.24966420758381, -40.76999825723249, -39.071768859867305, -34.87636984560834, -35.572901644121444, -29.549669152803975, -34.30530675273345, -38.49105307375407, -34.86753109654644, -39.269479818131224] -35.90237427083825 -0.008915103156122739\n",
      "GLUCOSE [-45.659142524303945, -91.02015953964946, -67.4130334611788, -58.88304616460352, -62.102059955590946, -34.92520344677448, -47.83296917131582, -108.729494943527, -51.20651281293434, -96.99891835900614] -66.47705403788845 -0.004024792803333188\n",
      "INSULIN [-1.444994351606827, -1.7761829664596513, -1.9871763919388323, -3.0421559492472996, -2.29007188958192, -2.5415068803693615, -2.296293415150965, -1.9053568038195436, -2.851982219450398, -2.4734250115753693] -2.2609145879200168 -0.0057201328634412405\n",
      "HOMA-IR [-0.8866226532434225, -1.5271679094907693, -1.5211540160600814, -2.113581249691648, -1.7917665631105488, -1.2178962127532929, -1.3551530820623485, -1.4433071507081627, -1.6091183098408126, -1.655311380146817] -1.5121078527107905 -0.004305647219808262\n",
      "next_1y_hba1c [-0.5220455330124487, 0.30103667647375976, 0.42295164513540295, 0.7647529196303628, 0.38318410531596325, 0.4986733022712557, 0.06720331747774122, 0.6505604659193787, -0.03702248846116096, -0.6226824689040129] 0.19066119418462424 0.45661832043574846\n",
      "d_1y_hba1c [-0.13691045070465213, -0.29000710732054036, -0.20508691794551592, -0.135599551737019, -0.21050119331742234, -0.25092660056387994, -0.04423812124522142, -0.009735969146956291, -0.06290890790135895, -0.19011145057876822] -0.15360262704613345 0.13316847932895853\n",
      "next_2y_hba1c [0.24699085219065997, 0.4073935772964892, 0.033588761174967896, 0.09018987341772189, 0.10881104033970324, 0.15000000000000024, 0.742891252441936, 0.4662124736900414, -0.1091954022988495, 0.25755279377589135] 0.23944352220285614 0.49985774069016353\n",
      "d_2y_hba1c [-0.7074235807860254, -0.04960317460317443, -0.2208777024846722, -0.2940735183795953, -0.2520507084265473, -0.3344709897610927, -0.11377527033380352, -0.1273982218062697, -0.05595622979358361, -0.1866925064599485] -0.23423219028347128 0.1590190578628722\n",
      "next_1y_gfr [0.3897452803872057, 0.7366086018928957, 0.4295938298344153, 0.5152918444504728, 0.38192136984307923, 0.2823555493128407, 0.41008870272439124, 0.32155358527741096, 0.2977162211090283, 0.6296841493112346] 0.4394559134142975 0.5463770071753353\n",
      "d_1y_gfr [-0.0003880050994955475, -0.018153795995791366, -0.0005754843660079612, -0.01792551929788888, -0.0007092807794726319, -0.01611156371303979, -0.01987686473122663, -0.00012616931922648789, -0.03352348127551452, -0.00395246098176516] -0.011134262555942898 0.1664497467041743\n",
      "next_2y_gfr [0.7269042872477929, -0.0041498045133412376, -0.31238507508427826, 0.6518179053907303, -1.4926631359020455, 0.43916280466872626, 0.37547971774027156, 0.05545168498335551, 0.7255841564330519, 0.5317059656467402] 0.16969085066110037 0.4230992748355555\n",
      "d_2y_gfr [-0.007923673997412672, -0.04459666773698667, -0.002472640018727601, -2.079347916494889e-05, -8.772442832905014e-05, -0.023994202217259897, -0.019467834287638874, -0.036479183346677324, -0.012775842044134622, -0.03712035995500562] -0.01849389215113373 0.19033254203937527\n",
      "next_1y_glucose [-0.030388832186481363, 0.34197269221415416, -0.14238244255054489, 0.7378310355649045, 0.6976673938580299, 0.5343672423457811, -0.058942268873637405, 0.6185027388922703, 0.07187904128055178, 0.0393951257707742] 0.2809901726315802 0.48600074887110606\n",
      "d_1y_glucose [-0.0898959881129271, -0.017899909643916523, -0.0064450864837111155, -0.0016616021034552286, -0.0006720766056242944, -0.007308384973238713, -0.002340577979460301, -0.02359761574263075, -0.0013010892590243905, -0.0017643499245407757] -0.015288668082852919 0.13989161920656507\n",
      "next_2y_glucose [0.6279043579445216, -0.07306338028169002, -0.41641357027463655, 0.40032531215614997, 0.4761400856638952, 0.5181318610106666, -0.7505854800936769, 0.5404303463881397, 0.34719015896049066, 0.26512535655307] 0.19351850480269303 0.3845737878082197\n",
      "d_2y_glucose [-0.09766360314472378, -0.0172163907066627, -0.00036511341335399017, -0.004004805766920194, -0.00019241870309816278, -0.19016293735394885, -0.06600794540083532, -0.006040836051709242, -0.01228633161588899, -0.035759627592043985] -0.04297000097491852 0.1202330469858093\n",
      "next_1y_insulin [0.18288764709861094, 0.2956031426968976, 0.1217375442759806, 0.5967733310393093, 0.540064740129874, 0.6620920537246215, 0.5191886912223671, 0.6055905193351839, 0.3692483580179535, -0.08444335330516872] 0.38087426742356306 0.5264167107391649\n",
      "d_1y_insulin [-0.10969516106111787, -0.004532647365727005, -0.00326534123338984, -0.020749356873396296, -0.06546154792111025, -0.05348565670976613, -0.1714052826177499, -0.017917810424253888, -0.04095322361308629, -0.04444621211985811] -0.05319122399394556 0.15973021711083785\n",
      "next_2y_insulin [-1.981777998520379, 0.30132209229779605, 0.41505861635839725, -0.2961293660599198, 0.5803727131202789, -0.3834860915193208, 0.060527976858621546, 0.6333110482184141, -0.21071918666493605, -1.3689945303391475] -0.2250514726250195 0.3565768319063616\n",
      "d_2y_insulin [-0.2699269734472234, -0.14622158709220368, -0.009187712324430608, -0.19551609724790664, -0.0012620372887661713, -0.2624450516811214, -0.022516087670154228, -0.0056807772718587835, -0.2250773155001058, -0.14278430755712734] -0.12806179470808982 0.20323825112821456\n",
      "next_1y_homa-ir [-0.1929572050053554, 0.3910066754619874, -0.44002685171590095, 0.615198014404668, 0.5986164534215308, 0.6646785373378035, 0.3405054316425915, 0.540885919102909, 0.15563684539862843, -1.1898129659130916] 0.14837308541357705 0.41859366846703006\n",
      "d_1y_homa-ir [-0.12345628903510586, -0.011377577939336803, -0.0044269616762002695, -0.0036641130948280587, -0.05240189329709555, -0.03222047935700245, -0.12310998362862358, -0.0007960811674476709, -0.040331294941582385, -0.030325800065848618] -0.042211047420307124 0.20076026710724984\n",
      "next_2y_homa-ir [-3.7385286291636204, 0.42975511435677216, 0.37566233356461687, -0.4826411990528592, 0.5609391845863498, -0.6418271004038896, 0.23214907821689013, 0.49781652970151524, -0.05715764312282401, -2.084551135693388] -0.49083834670104365 0.32378847310723613\n",
      "d_2y_homa-ir [-0.18425426969295677, -0.10680045656733617, -0.0017357060926870282, -0.1874042299185541, -0.00014582060725931711, -0.26521311212090715, -0.020306065214483526, -0.004699357695883588, -0.18175616903307512, -0.13094626016765432] -0.10832614471107971 0.2410339084227231\n",
      "age [-20.24094343775218, -16.296646852411698, -21.571177775513025, -19.82415152420413, -22.05586472150941, -20.46681274915191, -19.802188921294693, -17.873193182055566, -17.432207587905477, -20.046374832141588] -19.56095615839397 -0.008436012461973763\n",
      "bmi [-21.610702190786036, -18.555490933850205, -20.190765681906612, -21.043230687605803, -21.98080140900663, -21.84213616921507, -21.522507276049584, -23.266288416695147, -19.492637064136037, -22.601491983445136] -21.21060518126962 -0.007251086099330162\n"
     ]
    }
   ],
   "source": [
    "from sklearn.model_selection import GroupKFold\n",
    "from sklearn.linear_model import LinearRegression\n",
    "\n",
    "lr = LinearRegression()\n",
    "\n",
    "# TODO: do 10-fold cv\n",
    "baseline_cv_scores = []\n",
    "baseline_lr_cv_scores = []\n",
    "for analysis, var, baseline in zip(analyses, dependent_vars, current_vars):\n",
    "    kf = GroupKFold(10)\n",
    "    X = baseline.reshape((len(baseline), 1))\n",
    "    if 'd_' in analysis:\n",
    "        baseline = np.zeros(baseline.shape)\n",
    "    included_subset = (~np.isnan(var))\n",
    "    X = X[included_subset]\n",
    "    baseline = baseline[included_subset]\n",
    "    client_ids_subset_t = client_ids_subset[included_subset]\n",
    "    var = var[included_subset]\n",
    "    scores = []\n",
    "    lr_scores = []\n",
    "    for train_indices, test_indices in kf.split(X, var, groups=client_ids_subset_t):\n",
    "        lr.fit(X[train_indices], var[train_indices])\n",
    "        var_test = lr.predict(X[test_indices])\n",
    "        scores.append(r2_score( var[test_indices], baseline[test_indices]))\n",
    "        lr_scores.append(r2_score(var[test_indices], var_test))\n",
    "    print(analysis, scores, np.mean(scores), np.mean(lr_scores))\n",
    "    baseline_cv_scores.append(scores)\n",
    "    baseline_lr_cv_scores.append(lr_scores)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f61ca0f-4846-4fbb-896b-002b265d664e",
   "metadata": {},
   "source": [
    "#### Saving baseline results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "8f9303c1-e99e-49b5-a727-0704766992c1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "output_cv = {'baseline': baseline_cv_scores, 'baseline_linear_regression': baseline_lr_cv_scores}\n",
    "output_cv = {k: {a: list(x) for a, x in zip(analyses, v)} for k, v in output_cv.items()}\n",
    "with open('results_2023_10_18/baseline_regression.json', 'w') as f:\n",
    "    json.dump(output_cv, f, indent=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82baa4ac-87b4-47f6-a776-284f4f6516b5",
   "metadata": {},
   "source": [
    "### Classification Baseline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "a94591a4-982a-4aa4-866f-4d385e601fbf",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "d_class_hba1c [0.5210526315789473, 0.5, 0.5260953608247423, 0.5282312925170068, 0.5405525846702317, 0.5333333333333333, 0.5, 0.5225146198830409, 0.5, 0.5573453608247423] [0.09523809523809525, 0.0, 0.1111111111111111, 0.11764705882352941, 0.15384615384615385, 0.125, 0.0, 0.09999999999999999, 0.0, 0.21052631578947367] 0.5229125183632044 0.09133687348083633\n",
      "d_class_1y_hba1c [0.5, 0.5, 0.49122807017543857, 0.5625, 0.5833333333333334, 0.5, 0.5, 0.5, 0.5, 0.5] [0.0, 0.0, 0.0, 0.2222222222222222, 0.2857142857142857, 0.0, 0.0, 0.0, 0.0, 0.0] 0.5137061403508772 0.050793650793650794\n",
      "d_class_gfr [0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.525, 0.5, 0.5, 0.5] [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.09523809523809523, 0.0, 0.0, 0.0] 0.5025000000000001 0.009523809523809523\n",
      "d_class_1y_gfr [0.5216450216450216, 0.4737171464330413, 0.5011074197120708, 0.5693877551020409, 0.5272727272727273, 0.5434782608695652, 0.5471794871794872, 0.5263157894736842, 0.5005537098560354, 0.5476190476190476] [0.15384615384615385, 0.14814814814814817, 0.14814814814814814, 0.28571428571428575, 0.16666666666666669, 0.16, 0.20689655172413793, 0.1, 0.08333333333333333, 0.23076923076923073] 0.5258276365162721 0.16835225183501046\n",
      "d_class_glucose [0.5161290322580645, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5] [0.0625, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 0.5016129032258064 0.00625\n",
      "d_class_1y_glucose [0.5, 0.5, 0.5, 0.5277777777777778, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5] [0.0, 0.0, 0.0, 0.10526315789473684, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0] 0.5027777777777778 0.010526315789473684\n",
      "d_class_insulin [0.6584375, 0.6744071146245059, 0.6822807017543859, 0.6436643835616438, 0.6645035691109669, 0.7145390070921986, 0.6965079365079365, 0.6156941649899397, 0.6588753387533876, 0.5956112852664577] [0.5842696629213484, 0.5822784810126581, 0.5882352941176471, 0.5476190476190476, 0.5641025641025642, 0.6511627906976744, 0.6190476190476191, 0.4864864864864864, 0.5526315789473684, 0.4705882352941177] 0.6604521001661422 0.564642176024653\n",
      "d_class_1y_insulin [0.43115942028985504, 0.6480938416422287, 0.5403940886699508, 0.5375, 0.5495951417004049, 0.6766766766766766, 0.667220376522702, 0.6400966183574879, 0.6161161161161162, 0.6247401247401247] [0.17142857142857143, 0.5217391304347827, 0.2285714285714286, 0.32432432432432434, 0.3243243243243243, 0.5499999999999999, 0.5161290322580645, 0.4666666666666666, 0.39999999999999997, 0.4878048780487804] 0.5931592404715547 0.39909883560569426\n",
      "d_class_homa-ir [0.6657329598506069, 0.6346528228423102, 0.6710616438356164, 0.6517123287671234, 0.6582849774339137, 0.7147955872809864, 0.6765079365079365, 0.5701219512195121, 0.6536246612466124, 0.5687106918238993] [0.5813953488372093, 0.5365853658536586, 0.575, 0.55, 0.5569620253164557, 0.6511627906976744, 0.5853658536585366, 0.4266666666666667, 0.5405405405405406, 0.41975308641975306] 0.6465205560808519 0.5423431677990495\n",
      "d_class_1y_homa-ir [0.5113636363636364, 0.65, 0.5436507936507936, 0.5572640509013785, 0.5335335335335336, 0.6724137931034483, 0.5916666666666667, 0.6085271317829457, 0.5925925925925926, 0.5689102564102564] [0.28571428571428575, 0.4615384615384615, 0.23529411764705882, 0.3428571428571428, 0.23529411764705882, 0.5128205128205129, 0.33333333333333337, 0.4242424242424242, 0.3125, 0.3783783783783784] 0.5829922455005252 0.35219727741786566\n"
     ]
    }
   ],
   "source": [
    "from sklearn.linear_model import LogisticRegression\n",
    "lr = LogisticRegression()\n",
    "\n",
    "# TODO: do 10-fold cv\n",
    "baseline_lr_cv_scores = []\n",
    "baseline_lr_cv_ba_scores = []\n",
    "for analysis, var, baseline in zip(class_analyses, class_dependent_vars, class_current_vars):\n",
    "    kf = GroupKFold(10)\n",
    "    X = baseline.reshape((len(baseline), 1))\n",
    "    if 'd_' in analysis:\n",
    "        baseline = np.zeros(baseline.shape)\n",
    "    included_subset = (~np.isnan(var))\n",
    "    X = X[included_subset]\n",
    "    baseline = baseline[included_subset]\n",
    "    client_ids_subset_t = client_ids_subset[included_subset]\n",
    "    var = var[included_subset]\n",
    "    scores = []\n",
    "    lr_scores = []\n",
    "    for train_indices, test_indices in kf.split(X, var, groups=client_ids_subset_t):\n",
    "        lr.fit(X[train_indices], var[train_indices])\n",
    "        var_test = lr.predict(X[test_indices])\n",
    "        scores.append(balanced_accuracy_score( var[test_indices], var_test))\n",
    "        lr_scores.append(f1_score(var[test_indices], var_test))\n",
    "    print(analysis, scores, lr_scores, np.mean(scores), np.mean(lr_scores))\n",
    "    baseline_lr_cv_scores.append(lr_scores)\n",
    "    baseline_lr_cv_ba_scores.append(scores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "3d102784-54c4-42de-81bb-6e90b57bd9b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "import json\n",
    "\n",
    "output_cv = {'baseline_logistic_regression': baseline_lr_cv_scores}\n",
    "output_cv = {k: {a: list(x) for a, x in zip(class_analyses, v)} for k, v in output_cv.items()}\n",
    "with open('results_2023_10_18/baseline_classification_f1.json', 'w') as f:\n",
    "    json.dump(output_cv, f, indent=1)\n",
    "    \n",
    "output_cv = {'baseline_logistic_regression': baseline_lr_cv_ba_scores}\n",
    "output_cv = {k: {a: list(x) for a, x in zip(class_analyses, v)} for k, v in output_cv.items()}\n",
    "with open('results_2023_10_18/baseline_classification_balanced_acc.json', 'w') as f:\n",
    "    json.dump(output_cv, f, indent=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d067976-e1b0-47e1-b938-7f4cda584783",
   "metadata": {},
   "source": [
    "### Run Regression models for all feature sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "f2246d79-85b8-41da-979f-d3c35da249f2",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.linear_model import ElasticNet, ElasticNetCV, Lasso, LassoCV, Ridge, RidgeCV\n",
    "\n",
    "import json\n",
    "import os"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "4b87d09b-a749-481a-8572-dd6d445b97d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "demo_columns = ['age', 'is_m', 'bmi']\n",
    "all_X_vals = [X_selected, X_chems, X_prots, X_mets, X_all, X_basic]\n",
    "all_labels = [chem_subset_cols, selected_chem_bp_cols, selected_prot_cols, selected_met_cols,\n",
    "              selected_columns_full, demo_columns]\n",
    "all_feature_names = ['selected_clinical', 'clinical', 'prots', 'mets_imputed', 'all', 'demo']\n",
    "all_models = [ElasticNet, LassoCV, ElasticNetCV, Ridge, RidgeCV]\n",
    "model_names = ['ElasticNet', 'LassoCV', 'ElasticNetCV', 'Ridge', 'RidgeCV']\n",
    "all_model_params = [{},\n",
    "                    {'selection': 'random', 'precompute': True, 'verbose': False, 'tol': 0.1},\n",
    "                    {'selection': 'random', 'precompute': True, 'verbose': False, 'tol': 0.1},\n",
    "                    {},\n",
    "                    {}]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "45b8398c-d967-4d00-beca-ef58f78592ed",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Feature set: selected_clinical ######\n",
      "results already exist\n",
      "Feature set: clinical ######\n",
      "results already exist\n",
      "Feature set: prots ######\n",
      "results already exist\n",
      "Feature set: mets_imputed ######\n",
      "results already exist\n",
      "Feature set: all ######\n",
      "results already exist\n",
      "Feature set: demo ######\n",
      "results already exist\n"
     ]
    }
   ],
   "source": [
    "for X, labels, feature_name in zip(all_X_vals, all_labels, all_feature_names):\n",
    "    model_cv_scores = {}\n",
    "    model_weights = {}\n",
    "    print('Feature set:', feature_name, '######')\n",
    "    if os.path.exists('results_2023_10_18/{0}_lassocv_weights.csv'.format(feature_name)) and\\\n",
    "       os.path.exists('results_2023_10_18/{0}_regression.json'.format(feature_name)):\n",
    "        print('results already exist')\n",
    "        continue\n",
    "    for Model, model_name, model_params in zip(all_models, model_names, all_model_params):\n",
    "        print(model_name)\n",
    "        m = Model(**model_params)\n",
    "        cv_scores, weights, predictions, r2 = run_analyses(X, \n",
    "                                                           analyses, \n",
    "                                                           dependent_vars, \n",
    "                                                           current_vars, client_ids_subset, m, True, labels, False,\n",
    "                                                           True, 'models_2023_10_18/' + feature_name + '_' + model_name)\n",
    "        model_cv_scores[model_name] = cv_scores\n",
    "        model_weights[model_name] = weights\n",
    "        print()\n",
    "    output_cv = {k: {a: list(x) for a, x in zip(analyses, v)} for k, v in model_cv_scores.items()}\n",
    "    with open('results_2023_10_18/{0}_regression.json'.format(feature_name), 'w') as f:\n",
    "        json.dump(output_cv, f, indent=1)\n",
    "    model_weights['LassoCV'].to_csv('results_2023_10_18/{0}_lassocv_weights.csv'.format(feature_name), index=None)\n",
    "    with open('results_2023_10_18/{0}_regression_model_weights.pkl'.format(feature_name), 'wb') as f:\n",
    "        pickle.dump(model_weights, f)\n",
    "    print('results saved\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "121c6956-cdda-41de-8f24-1c625f26c6ae",
   "metadata": {},
   "source": [
    "### Run Classification models for all feature sets"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "caea83f7-aff0-4a4d-8731-8b15f265ceab",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import LogisticRegression, LogisticRegressionCV, RidgeClassifier, RidgeClassifierCV\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.svm import SVC\n",
    "\n",
    "import json\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "b3cd48df-ceaa-45dc-b2c6-22614e6e0f8b",
   "metadata": {},
   "outputs": [],
   "source": [
    "demo_columns = ['age', 'is_m', 'bmi']\n",
    "all_X_vals = [X_selected, X_chems, X_prots, X_mets, X_all, X_basic]\n",
    "all_labels = [chem_subset_cols, selected_chem_bp_cols, selected_prot_cols, selected_met_cols,\n",
    "              selected_columns_full, demo_columns]\n",
    "all_feature_names = ['selected_clinical', 'clinical', 'prots', 'mets_imputed', 'all', 'demo']\n",
    "all_class_models = [LogisticRegression, LogisticRegression, LogisticRegressionCV, LogisticRegressionCV, RidgeClassifier, RidgeClassifierCV]\n",
    "class_model_names = ['LogReg', 'LogRegLasso', 'LogRegLassoCV', 'LogRegElasticNetCV', 'RidgeClassifier', 'RidgeClassifierCV']\n",
    "class_model_params = [{'max_iter': 1000},\n",
    "                      {'penalty': 'l1', 'solver': 'liblinear', 'max_iter': 1000},\n",
    "                      {'penalty': 'l1', 'solver': 'saga', 'scoring': 'f1', 'verbose': False},\n",
    "                      {'penalty': 'elasticnet', 'solver': 'saga', 'scoring': 'f1', 'verbose': False, 'l1_ratios': [0.5]},\n",
    "                      {},\n",
    "                      {}]"
   ]
  },
  {
   "cell_type": "raw",
   "id": "3ffbaa98-ce24-43a9-9c26-3fc9b397cc17",
   "metadata": {},
   "source": [
    "time for selected features with saga solver for LogRegLassoCV: 85.89411306381226\n",
    "clinical: 442.5704982280731\n",
    "mets: 5462.737760066986s\n",
    "all: 8479.637984275818s\n",
    "\n",
    "time for selected features with liblinear solver for LogRegLassoCV: 71.5005624294281\n",
    "clinical: way too long, dnf"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "0a70771c-2d74-4a76-9c1f-b4c8b6c512b1",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Feature set: selected_clinical #########\n",
      "results already exist\n",
      "Feature set: clinical #########\n",
      "results already exist\n",
      "Feature set: prots #########\n",
      "results already exist\n",
      "Feature set: mets_imputed #########\n",
      "results already exist\n",
      "Feature set: all #########\n",
      "results already exist\n",
      "Feature set: demo #########\n",
      "results already exist\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "\n",
    "for X, labels, feature_name in zip(all_X_vals, all_labels, all_feature_names):\n",
    "    t0 = time.time()\n",
    "    model_cv_scores = {}\n",
    "    model_weights = {}\n",
    "    model_f1_scores = {}\n",
    "    print('Feature set:', feature_name, '#########')\n",
    "    if os.path.exists('results_2023_10_18/{0}_classification_lassocv_weights.csv'.format(feature_name)) and\\\n",
    "       os.path.exists('results_2023_10_18/{0}_classification_lasso_weights.csv'.format(feature_name)) and \\\n",
    "       os.path.exists('results_2023_10_18/{0}_classification_f1.json'.format(feature_name)):\n",
    "        print('results already exist')\n",
    "        continue\n",
    "    for Model, model_name, model_params in zip(all_class_models, class_model_names, class_model_params):\n",
    "        print(model_name)\n",
    "        m = Model(**model_params)\n",
    "        ba_scores, f1_scores, weights, predictions, ba = run_classification_analyses(X,\n",
    "                                                                                 class_analyses, \n",
    "                                                                                 class_dependent_vars, \n",
    "                                                                                 class_current_vars,\n",
    "                                                                                 client_ids_subset, m, True, labels, False,\n",
    "                                                                True, 'models_2023_10_18/' + feature_name + '_classification_' + model_name)\n",
    "        model_cv_scores[model_name] = ba_scores\n",
    "        model_f1_scores[model_name] = f1_scores\n",
    "        model_weights[model_name] = weights\n",
    "        print('\\n')\n",
    "    output_cv = {k: {a: list(x) for a, x in zip(class_analyses, v)} for k, v in model_cv_scores.items()}\n",
    "    with open('results_2023_10_18/{0}_classification_balanced_acc.json'.format(feature_name), 'w') as f:\n",
    "        json.dump(output_cv, f, indent=1)\n",
    "    output_f1 = {k: {a: list(x) for a, x in zip(class_analyses, v)} for k, v in model_f1_scores.items()}\n",
    "    with open('results_2023_10_18/{0}_classification_f1.json'.format(feature_name), 'w') as f:\n",
    "        json.dump(output_f1, f, indent=1)\n",
    "    model_weights['LogRegLassoCV'].to_csv('results_2023_10_18/{0}_classification_lassocv_weights.csv'.format(feature_name), index=None)\n",
    "    model_weights['LogRegLasso'].to_csv('results_2023_10_18/{0}_classification_lasso_weights.csv'.format(feature_name), index=None)\n",
    "    with open('results_2023_10_18/{0}_classification_model_weights.pkl'.format(feature_name), 'wb') as f:\n",
    "        pickle.dump(model_weights, f)\n",
    "    print('results saved')\n",
    "    print('time elapsed:', time.time() - t0)\n",
    "    print('\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d638a9ef-b361-4e52-a280-e3532b93f616",
   "metadata": {},
   "source": [
    "### Advanced models for regression analysis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "2911fad2-6c6e-4eb9-96f0-991992b8063f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.svm import SVR, LinearSVR\n",
    "from sklearn.neighbors import KNeighborsRegressor\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "\n",
    "advanced_models = [RandomForestRegressor, SVR, LinearSVR, KNeighborsRegressor]\n",
    "advanced_model_names = ['RandomForestRegressor', 'SVR', 'LinearSVR', 'KNeighborsRegressor']\n",
    "advanced_model_params = [{'verbose': False, 'n_jobs': -1, 'max_depth': 15},\n",
    "                         {'verbose': False, 'C': 1},\n",
    "                         {'verbose': False, 'C': 1,  'max_iter': 10000},\n",
    "                         {}]"
   ]
  },
  {
   "cell_type": "raw",
   "id": "564adad7-64a4-494a-aabd-b761abbb572a",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "cdf657da-e901-4245-ac07-a2e3512282af",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Feature set: selected_clinical ######\n",
      "results already exist\n",
      "Feature set: clinical ######\n",
      "results already exist\n",
      "Feature set: prots ######\n",
      "results already exist\n",
      "Feature set: mets_imputed ######\n",
      "results already exist\n",
      "Feature set: all ######\n",
      "results already exist\n",
      "Feature set: demo ######\n",
      "results already exist\n"
     ]
    }
   ],
   "source": [
    "for X, labels, feature_name in zip(all_X_vals, all_labels, all_feature_names):\n",
    "    model_cv_scores = {}\n",
    "    model_weights = {}\n",
    "    print('Feature set:', feature_name, '######')\n",
    "    if os.path.exists('results_2023_10_18/{0}_rfr_weights.csv'.format(feature_name)) and\\\n",
    "       os.path.exists('results_2023_10_18/{0}_advanced_models_regression.json'.format(feature_name)):\n",
    "        print('results already exist')\n",
    "        continue\n",
    "    for Model, model_name, model_params in zip(advanced_models, advanced_model_names, advanced_model_params):\n",
    "        print(model_name)\n",
    "        m = Model(**model_params)\n",
    "        cv_scores, weights, predictions, r2 = run_analyses(X, \n",
    "                                                           analyses, \n",
    "                                                           dependent_vars, \n",
    "                                                           current_vars, client_ids_subset, m, True, labels, False,\n",
    "                                                          False, 'models_2023_10_18/' + feature_name + '_' + model_name)\n",
    "        model_cv_scores[model_name] = cv_scores\n",
    "        model_weights[model_name] = weights\n",
    "        print()\n",
    "    output_cv = {k: {a: list(x) for a, x in zip(analyses, v)} for k, v in model_cv_scores.items()}\n",
    "    with open('results_2023_10_18/{0}_advanced_models_regression.json'.format(feature_name), 'w') as f:\n",
    "        json.dump(output_cv, f, indent=1)\n",
    "    model_weights['RandomForestRegressor'].to_csv('results_2023_10_18/{0}_rfr_weights.csv'.format(feature_name), index=None)\n",
    "    print('results saved\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d070854b-8b79-48b0-be7e-22264959648b",
   "metadata": {},
   "source": [
    "### Advanced models for classification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "ab319aaa-f353-4e09-ba2f-49385b8d453b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.svm import SVC, LinearSVC\n",
    "from sklearn.neighbors import KNeighborsClassifier\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "\n",
    "advanced_models = [RandomForestClassifier, SVC, LinearSVC, KNeighborsClassifier]\n",
    "advanced_model_names = ['RandomForestClassifier', 'SVC', 'LinearSVC', 'KNeighborsClassifier']\n",
    "advanced_model_params = [{'verbose': False, 'n_jobs': -1, 'max_depth': 15},\n",
    "                         {'verbose': False, 'C': 1},\n",
    "                         {'verbose': False, 'C': 1, 'max_iter': 10000},\n",
    "                         {}]"
   ]
  },
  {
   "cell_type": "raw",
   "id": "7c2c4c6f-8c55-49b9-8d2b-388d05f7222f",
   "metadata": {},
   "source": [
    "Timing results:\n",
    "selected: 24.89202189445495\n",
    "clinical: 29.304314136505127\n",
    "prots: 61.24226665496826\n",
    "mets_zeros: 846.0482602119446 \n",
    "mets_imputed: 848.3962519168854\n",
    "all: 3818.5492985248566"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "a94de4c9-99ac-42a7-8811-1d333c752c85",
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Feature set: selected_clinical #########\n",
      "results already exist\n",
      "Feature set: clinical #########\n",
      "results already exist\n",
      "Feature set: prots #########\n",
      "results already exist\n",
      "Feature set: mets_imputed #########\n",
      "results already exist\n",
      "Feature set: all #########\n",
      "results already exist\n",
      "Feature set: demo #########\n",
      "results already exist\n"
     ]
    }
   ],
   "source": [
    "for X, labels, feature_name in zip(all_X_vals, all_labels, all_feature_names):\n",
    "    t0 = time.time()\n",
    "    model_cv_scores = {}\n",
    "    model_weights = {}\n",
    "    model_f1_scores = {}\n",
    "    print('Feature set:', feature_name, '#########')\n",
    "    if os.path.exists('results_2023_10_18/{0}_classification_rfc_weights.csv'.format(feature_name)) and\\\n",
    "       os.path.exists('results_2023_10_18/{0}_advanced_models_classification_f1.json'.format(feature_name)):\n",
    "        print('results already exist')\n",
    "        continue\n",
    "    for Model, model_name, model_params in zip(advanced_models, advanced_model_names, advanced_model_params):\n",
    "        print(model_name)\n",
    "        m = Model(**model_params)\n",
    "        ba_scores, f1_scores, weights, predictions, ba = run_classification_analyses(X,\n",
    "                                                                                 class_analyses, \n",
    "                                                                                 class_dependent_vars, \n",
    "                                                                                 class_current_vars,\n",
    "                                                                                 client_ids_subset, m, True, labels, False,\n",
    "                                                                                False, 'models_2023_10_18/' + feature_name + '_classification_' + model_name)\n",
    "        model_cv_scores[model_name] = ba_scores\n",
    "        model_f1_scores[model_name] = f1_scores\n",
    "        model_weights[model_name] = weights\n",
    "        print('\\n')\n",
    "    output_cv = {k: {a: list(x) for a, x in zip(class_analyses, v)} for k, v in model_cv_scores.items()}\n",
    "    with open('results_2023_10_18/{0}_advanced_models_classification_balanced_acc.json'.format(feature_name), 'w') as f:\n",
    "        json.dump(output_cv, f, indent=1)\n",
    "    output_f1 = {k: {a: list(x) for a, x in zip(class_analyses, v)} for k, v in model_f1_scores.items()}\n",
    "    with open('results_2023_10_18/{0}_advanced_models_classification_f1.json'.format(feature_name), 'w') as f:\n",
    "        json.dump(output_f1, f, indent=1)\n",
    "    model_weights['RandomForestClassifier'].to_csv('results_2023_10_18/{0}_classification_rfc_weights.csv'.format(feature_name), index=None)\n",
    "    print('results saved')\n",
    "    print('time elapsed:', time.time() - t0)\n",
    "    print('\\n')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dafb2d42-e722-4c72-97b4-7712436a1fa3",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  },
  "widgets": {
   "application/vnd.jupyter.widget-state+json": {
    "state": {},
    "version_major": 2,
    "version_minor": 0
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}