manuscript_demo_graphs.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 190,
   "id": "62fe277d",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T19:10:16.580277Z",
     "start_time": "2023-11-04T19:10:16.569201Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "import csv\n",
    "import paramiko\n",
    "import gzip\n",
    "\n",
    "import pymysql\n",
    "from sshtunnel import SSHTunnelForwarder\n",
    "\n",
    "from tqdm import tqdm\n",
    "\n",
    "from datetime import timedelta\n",
    "from datetime import date\n",
    "from datetime import datetime\n",
    "\n",
    "import pandas as pd\n",
    "from collections import defaultdict\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "import matplotlib.dates as mdates\n",
    "\n",
    "import re"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "id": "77c7bc7b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:24:59.149935Z",
     "start_time": "2023-11-03T20:24:59.050244Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "plt.rcParams['font.family'] = 'Arial'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "id": "65168265",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:24:59.326858Z",
     "start_time": "2023-11-03T20:24:59.154516Z"
    }
   },
   "outputs": [],
   "source": [
    "database = \"clinical_merge_v5_2020q4\"\n",
    "ped_data = \"/home/kel2158/Projects/Heritability_Base_Data/pedigree_v3_updated_results/\"\n",
    "output_folder = \"/home/kel2158/Projects/COVIDHerit/data/12_29_21/\"\n",
    "hhid_folder = \"/home/kel2158/Projects/hhid_update/data/\"\n",
    "\n",
    "demo_file = os.path.join(output_folder, \"demo.txt.gz\")\n",
    "hosp_days_file = os.path.join(output_folder, \"days_hosp.txt.gz\")\n",
    "covid_nopro_file = os.path.join(output_folder, \"covid_noproband.txt.gz\")\n",
    "covid_file = os.path.join(output_folder, \"covid.txt.gz\")\n",
    "\n",
    "hosp_stat_file = os.path.join(output_folder, \"hosp_stat.txt.gz\")\n",
    "hosp_stat_nopro_file = os.path.join(output_folder, \"hosp_stat_noproband.txt.gz\")\n",
    "high_deg_file = os.path.join(output_folder, \"covid_high_degfam.txt.gz\")\n",
    "high_deg_nopro_file = os.path.join(output_folder, \"covid_high_degfam_noproband.txt.gz\")\n",
    "\n",
    "time2pcr_file = os.path.join(output_folder, \"time2firstpcr.txt.gz\")\n",
    "\n",
    "cummulative_directory = os.path.join(output_folder, \"cummulative_classswitch\")\n",
    "\n",
    "# Input files\n",
    "finename_fam_input = os.path.join(ped_data, \"family_v3.txt.gz\")\n",
    "filename_hhid = os.path.join(hhid_folder, \"hhid_trimmed_notcln.txt.gz\")\n",
    "buildid_file = os.path.join(hhid_folder, \"buildid_notcln.txt.gz\")\n",
    "filename_ped_input = os.path.join(ped_data, \"pedigree_v3.txt.gz\")\n",
    "\n",
    "\n",
    "#Figure output directory\n",
    "figure_directory = '/home/kel2158/Projects/COVIDHerit/figures'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 120,
   "id": "9f1b5217",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:30:14.927165Z",
     "start_time": "2023-11-03T20:30:14.917567Z"
    }
   },
   "outputs": [],
   "source": [
    "def return_sign_in(login_loc = '/home/kel2158/login.txt'):\n",
    "    login_dict = dict()\n",
    "    with open(login_loc, \"r\") as file:\n",
    "        for line in file.readlines():\n",
    "            login_dict[ line.split(\",\")[0]] = line.split(\",\")[1].strip()\n",
    "    return(login_dict)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "id": "55a2f3a1",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:30:18.627931Z",
     "start_time": "2023-11-03T20:30:18.468418Z"
    }
   },
   "outputs": [],
   "source": [
    "login_dict = return_sign_in()\n",
    "\n",
    "ssh = paramiko.SSHClient()\n",
    "ssh.set_missing_host_key_policy(paramiko.AutoAddPolicy())\n",
    "ssh.connect('10.144.220.21', username=login_dict['user'], password=login_dict['mimir_pw'])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "id": "290e652b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:30:19.918469Z",
     "start_time": "2023-11-03T20:30:19.437748Z"
    }
   },
   "outputs": [],
   "source": [
    "server = SSHTunnelForwarder(\n",
    "    '10.144.220.21',\n",
    "    ssh_username='kel2158',\n",
    "    ssh_password=login_dict['mimir_pw'],\n",
    "    remote_bind_address=('127.0.0.1', 3306)\n",
    ")\n",
    "server.start()\n",
    "\n",
    "db = pymysql.connect(\n",
    "    host='127.0.0.1',\n",
    "    port=server.local_bind_port,\n",
    "    user=login_dict['user'],\n",
    "    password=login_dict[\"db_pw\"],\n",
    "    db=database\n",
    ")\n",
    "\n",
    "cur = db.cursor()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d7dacc67",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "446c04f8",
   "metadata": {},
   "source": [
    "## Pull Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 123,
   "id": "424eabcd",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:30:39.512255Z",
     "start_time": "2023-11-03T20:30:26.864096Z"
    },
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "#Read in time to first pcr trait file \n",
    "file_name = filename_ped_input\n",
    "gfile = gzip.open(file_name)\n",
    "ped_df = pd.read_csv(gfile,sep='\\t')\n",
    "gfile.close()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "id": "05a9a457",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:30:39.731239Z",
     "start_time": "2023-11-03T20:30:39.518358Z"
    }
   },
   "outputs": [],
   "source": [
    "#Read in demo df\n",
    "file_name = demo_file\n",
    "gfile = gzip.open(file_name)\n",
    "demo_df = pd.read_csv(gfile,sep='\\t')\n",
    "gfile.close()\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "id": "f13cf9b8",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:30:40.487677Z",
     "start_time": "2023-11-03T20:30:39.734326Z"
    }
   },
   "outputs": [],
   "source": [
    "#Read in family and id dataframe\n",
    "file_name = finename_fam_input\n",
    "gfile = gzip.open(file_name)\n",
    "fam_df = pd.read_csv(gfile,sep='\\t')\n",
    "gfile.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "id": "020de5ca",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:30:45.342858Z",
     "start_time": "2023-11-03T20:30:40.494025Z"
    }
   },
   "outputs": [],
   "source": [
    "#Read in HHID file \n",
    "file_name = filename_hhid\n",
    "gfile = gzip.open(file_name)\n",
    "hhid_df = pd.read_csv(gfile,sep='\\t')\n",
    "gfile.close()\n",
    "\n",
    "\n",
    "mrn2hhid = dict(zip( hhid_df[\"mrn\"], hhid_df[\"trimmed_hhid\"] ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 127,
   "id": "77b4b888",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:31:08.197165Z",
     "start_time": "2023-11-03T20:30:45.347364Z"
    }
   },
   "outputs": [],
   "source": [
    "#Read in BuildingID file \n",
    "file_name = buildid_file\n",
    "gfile = gzip.open(file_name)\n",
    "bid_df = pd.read_csv(gfile,sep='\\t')\n",
    "gfile.close()\n",
    "\n",
    "\n",
    "mrn2bid = dict(zip( bid_df[\"mrn\"], bid_df[\"building_id\"] ))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "id": "b421f195",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:31:08.308495Z",
     "start_time": "2023-11-03T20:31:08.201946Z"
    }
   },
   "outputs": [],
   "source": [
    "#Read in time to first pcr trait file \n",
    "file_name = time2pcr_file\n",
    "gfile = gzip.open(file_name)\n",
    "time2pcr_df = pd.read_csv(gfile,sep='\\t')\n",
    "gfile.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "id": "ef8b1e4a",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:31:08.433247Z",
     "start_time": "2023-11-03T20:31:08.313066Z"
    }
   },
   "outputs": [],
   "source": [
    "#Read in covid susceptibility file\n",
    "file_name = covid_file\n",
    "gfile = gzip.open(file_name)\n",
    "cov_df = pd.read_csv(gfile,sep='\\t')\n",
    "gfile.close()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "id": "a82a84db",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:31:08.536580Z",
     "start_time": "2023-11-03T20:31:08.436352Z"
    }
   },
   "outputs": [],
   "source": [
    "#Read in cov noproband\n",
    "file_name = covid_nopro_file\n",
    "gfile = gzip.open(file_name)\n",
    "covnopro_df = pd.read_csv(gfile,sep='\\t')\n",
    "gfile.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "id": "f54dcbce",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:31:08.616398Z",
     "start_time": "2023-11-03T20:31:08.539636Z"
    }
   },
   "outputs": [],
   "source": [
    "#Read in days_hosp\n",
    "file_name = hosp_days_file\n",
    "gfile = gzip.open(file_name)\n",
    "hospday_df = pd.read_csv(gfile,sep='\\t')\n",
    "\n",
    "gfile.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 132,
   "id": "4f1208fa",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:31:08.745926Z",
     "start_time": "2023-11-03T20:31:08.621447Z"
    }
   },
   "outputs": [],
   "source": [
    "#Read in hospitalization_stat\n",
    "file_name = hosp_stat_file\n",
    "gfile = gzip.open(file_name)\n",
    "hosp_df = pd.read_csv(gfile,sep='\\t')\n",
    "\n",
    "gfile.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 133,
   "id": "9764af51",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:31:08.822732Z",
     "start_time": "2023-11-03T20:31:08.749267Z"
    }
   },
   "outputs": [],
   "source": [
    "#Read in hospitalization_stat_noproband\n",
    "file_name = hosp_stat_nopro_file\n",
    "gfile = gzip.open(file_name)\n",
    "hosp_nopro_df = pd.read_csv(gfile,sep='\\t')\n",
    "\n",
    "gfile.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 134,
   "id": "14b83792",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:31:08.916972Z",
     "start_time": "2023-11-03T20:31:08.826190Z"
    }
   },
   "outputs": [],
   "source": [
    "#Read in high_deg\n",
    "file_name = high_deg_file\n",
    "gfile = gzip.open(file_name)\n",
    "high_deg_df = pd.read_csv(gfile,sep='\\t')\n",
    "\n",
    "gfile.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 135,
   "id": "593143de",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:31:09.021304Z",
     "start_time": "2023-11-03T20:31:08.920198Z"
    }
   },
   "outputs": [],
   "source": [
    "#Read in high_deg_noproband\n",
    "file_name = high_deg_nopro_file\n",
    "gfile = gzip.open(file_name)\n",
    "high_deg_nopro_df = pd.read_csv(gfile,sep='\\t')\n",
    "\n",
    "gfile.close()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 136,
   "id": "01d843e2",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:23.796212Z",
     "start_time": "2023-11-03T20:31:09.024178Z"
    }
   },
   "outputs": [],
   "source": [
    "pull_data_query = \"\"\"SELECT *\n",
    "FROM clinical_relationships_v3.all_relationships_to_generate_pedigree_file\"\"\"\n",
    "\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata = cur.fetchall()\n",
    "\n",
    "rel_df = pd.DataFrame(list(rawdata),\n",
    "                    columns=[\"mrn\",\"relationship\",\"relation_mrn\",\"provided_relationship\",\n",
    "                             ])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 137,
   "id": "5a9c38fd",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:23.983771Z",
     "start_time": "2023-11-03T20:32:23.800619Z"
    }
   },
   "outputs": [],
   "source": [
    "pull_data_query = \"\"\"SELECT *\n",
    "FROM clinical_relationships.relationships_degree\"\"\"\n",
    "\n",
    "\n",
    "cur.execute(pull_data_query)\n",
    "rawdata = cur.fetchall()\n",
    "\n",
    "reldegree_df = pd.DataFrame(list(rawdata),\n",
    "                    columns=[\"relationship\",\"degree\"\n",
    "                             ])\n",
    "\n",
    "rel2degree = dict( zip( reldegree_df[\"relationship\"],reldegree_df[\"degree\"] ) ) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 138,
   "id": "b73d69f3",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:24.082168Z",
     "start_time": "2023-11-03T20:32:23.987484Z"
    }
   },
   "outputs": [],
   "source": [
    "db.close()\n",
    "server.stop() "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c191475",
   "metadata": {},
   "source": [
    "## Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "id": "6f4b58b2",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:24.095486Z",
     "start_time": "2023-11-03T20:32:24.086467Z"
    }
   },
   "outputs": [],
   "source": [
    "# Initialize a list to keep track of all figures\n",
    "all_figures = []\n",
    "\n",
    "# Create a function to save a figure and add it to the list\n",
    "def save_and_track_figure(fig, filename, save_directory = figure_directory):\n",
    "    # Combine the save_directory and filename to create the full path\n",
    "    save_path = os.path.join(save_directory, filename)\n",
    "    \n",
    "    fig.patch.set_facecolor('white')\n",
    "    # Save the figure with the full path\n",
    "    fig.savefig(save_path, dpi=300, bbox_inches=\"tight\")\n",
    "    \n",
    "    # Append the saved file path to the list of saved figures\n",
    "    all_figures.append(save_path)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "id": "e5d936ad",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:24.225208Z",
     "start_time": "2023-11-03T20:32:24.099235Z"
    }
   },
   "outputs": [],
   "source": [
    "def pull_fam_count( pheno, df, print_res = True, \\\n",
    "                   return_fams_df=False, type_pheno=\"binary_withproband\",\n",
    "                   return_patient_limited_df = False ):\n",
    "    \n",
    "    \"\"\"Returns the number of families in a phenotype for SOLAR based on pheno \n",
    "     Needs phenotype column, famid column and ptid column\"\"\"    \n",
    "    \n",
    "    #subset on phenotype of interest and those in our pedigree\n",
    "    df_phenosub = df[ (df[\"pheno\"] == pheno ) \n",
    "                       &  (~df[\"famid\"].isnull())] \n",
    "    \n",
    "    #Get the total number of  members in each fam and then number of cases\n",
    "    df_phenogrouped = df_phenosub.groupby( 'famid' )\\\n",
    "        .agg({ \"ptid\":\"count\",\"value\":sum })\\\n",
    "        .reset_index()\\\n",
    "        .rename(columns={ \"ptid\":\"num_fam_members\",\"value\":\"num_cases\" })\n",
    "    \n",
    "    #Subset to fams that will be included in SOLAR\n",
    "    #ie has more than one family member an at least one case in each family\n",
    "    if type_pheno == \"binary_withproband\":\n",
    "        solar_fams_df = df_phenogrouped[ (df_phenogrouped[\"num_fam_members\"] > 1) \\\n",
    "                                        & (df_phenogrouped[\"num_cases\"] > 0) ]\n",
    "    \n",
    "    else:\n",
    "         solar_fams_df = df_phenogrouped[ (df_phenogrouped[\"num_fam_members\"] > 1)]\n",
    "    \n",
    "    solar_limited_df = df_phenosub[ df_phenosub[\"famid\"].isin(solar_fams_df[\"famid\"].unique())]\n",
    "    solar_limited_df.drop_duplicates(inplace=True)\n",
    "    solar_patients = set( solar_limited_df[\"ptid\"] )\n",
    "    \n",
    "    \n",
    "    \n",
    "    if print_res:\n",
    "            print(f'There are {len(solar_fams_df)} SOLAR families')\n",
    "    \n",
    "    if return_fams_df:\n",
    "        return solar_fams_df, solar_patients\n",
    "    \n",
    "    elif return_patient_limited_df:\n",
    "        return solar_limited_df\n",
    "    \n",
    "    else:\n",
    "        return(len( solar_fams_df) )\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d07dc8cc",
   "metadata": {},
   "source": [
    "## Process Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "id": "4f178eee",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:24.352553Z",
     "start_time": "2023-11-03T20:32:24.227961Z"
    }
   },
   "outputs": [],
   "source": [
    "#Map for differences in how SOLAR includes people\n",
    "#Options: binary_withproband OR quant_or_binary_noproband\n",
    "pheno2type = {   \n",
    "\"covid\":\"binary_withproband\",\n",
    "\"covid_noproband\":\"binary_noproband\",\n",
    "\"days_hosp\":\"quant\",\n",
    "    \"hosp_stat\":\"binary_withproband\",\n",
    "    \"hosp_stat_nopro\":\"binary_noproband\",\n",
    "    \"covid_highdeg_fam\":\"binary_withproband\",\n",
    "    \"covid_highdeg_nopro\":\"binary_noproband\",\n",
    "}\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "id": "c0f2b694",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:26.861362Z",
     "start_time": "2023-11-03T20:32:24.356755Z"
    }
   },
   "outputs": [],
   "source": [
    "#Info for mapping\n",
    "mrn2sex = dict(zip( demo_df[\"ptid\"],demo_df[\"sex\"] ))\n",
    "mrn2famid = dict(zip( fam_df[\"ptid\"],fam_df[\"famid\"] ))\n",
    "\n",
    "#Map famid to demo df\n",
    "demo_df[\"famid\"]  = demo_df[\"ptid\"].map(mrn2famid)\n",
    "\n",
    "pheno_dfs = [cov_df, covnopro_df, \n",
    "             hospday_df,hosp_df,hosp_nopro_df,\n",
    "             high_deg_df,high_deg_nopro_df]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "id": "ef12ae92",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:26.873559Z",
     "start_time": "2023-11-03T20:32:26.865884Z"
    }
   },
   "outputs": [],
   "source": [
    "def merge_limit2solar(df, pheno):\n",
    "    #Merge with demo data\n",
    "    df = df.merge(demo_df, how=\"left\", left_on=\"ptid\", right_on=\"ptid\") \n",
    "    \n",
    "    #Get patients included in solar runs\n",
    "    df = pull_fam_count(pheno, df,\n",
    "                                type_pheno=pheno2type[pheno],\n",
    "                                return_patient_limited_df=True)\n",
    "    return df\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "id": "5d1e987b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:28.158104Z",
     "start_time": "2023-11-03T20:32:26.876775Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "There are 5676 SOLAR families\n",
      "There are 1324 SOLAR families\n",
      "There are 158 SOLAR families\n",
      "There are 115 SOLAR families\n",
      "There are 158 SOLAR families\n",
      "There are 296 SOLAR families\n",
      "There are 962 SOLAR families\n"
     ]
    }
   ],
   "source": [
    "#Map demo and limit to SOLAR patients\n",
    "covnopro_limited_df = merge_limit2solar(covnopro_df,\"covid_noproband\")\n",
    "cov_limited_df = merge_limit2solar(cov_df,\"covid\")\n",
    "hospday_limited_df = merge_limit2solar(hospday_df,\"days_hosp\")\n",
    "\n",
    "hosp_limited_df = merge_limit2solar(hosp_df,\"hosp_stat\")\n",
    "hosp_nopro_limited_df = merge_limit2solar(hosp_nopro_df,\"hosp_stat_nopro\")\n",
    "highdeg_limited_df = merge_limit2solar(high_deg_df,\"covid_highdeg_fam\")\n",
    "highdeg_nopro_limited_df = merge_limit2solar(high_deg_nopro_df,\"covid_highdeg_nopro\")\n",
    "\n",
    "\n",
    "df_limited = [covnopro_limited_df,cov_limited_df,hospday_limited_df ,\n",
    "             hosp_limited_df,hosp_nopro_limited_df,highdeg_limited_df,highdeg_nopro_limited_df]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "id": "f4d24ecc",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:28.171691Z",
     "start_time": "2023-11-03T20:32:28.162255Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "12764"
      ]
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "len( covnopro_limited_df[\"ptid\"].unique() )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 146,
   "id": "30398be3",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:29.344886Z",
     "start_time": "2023-11-03T20:32:28.174788Z"
    },
    "code_folding": []
   },
   "outputs": [],
   "source": [
    "covnopro_limited_rel_df = covnopro_limited_df.copy()\n",
    "covnopro_limited_rel_df[\"ptid\"] = covnopro_limited_rel_df[\"ptid\"].astype(str)\n",
    "\n",
    "#Merge relationship and pheno and limit to solar families AND both sides of relationship must\n",
    "#have covid tests so limit rel file where relation_mrn is in pheno_df also\n",
    "rel_limit_df = covnopro_limited_rel_df.merge( rel_df[ rel_df[\"relation_mrn\"].isin( covnopro_limited_rel_df[\"ptid\"] ) ]\n",
    "                                         , how='inner', left_on='ptid', right_on='mrn')\n",
    " \n",
    "#Map the degree of relationship    \n",
    "rel_limit_df[\"rel_degree\"] = rel_limit_df[\"relationship\"].map(rel2degree)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5883515f",
   "metadata": {},
   "source": [
    "## Check Ethnicity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "id": "10622e96",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:29.359580Z",
     "start_time": "2023-11-03T20:32:29.348700Z"
    }
   },
   "outputs": [],
   "source": [
    "#Dictionary to map ethnicity to SOLARStrap friendly coding\n",
    "map_ethnic_dict = { \"NOT HISPANIC OR LATINO OR SPANISH ORIGIN\":'U',\n",
    "\"DECLINED\":'U',\n",
    "\"HISPANIC OR LATINO OR SPANISH ORIGIN\":'H',\n",
    "\"AFRICAN AMERICAN\":'B',\n",
    "\"MULTI-RACIAL\":\"Other\",\n",
    "\"UNKNOWN\":\"U\",\n",
    "\"CAUCASIAN\":\"W\",\n",
    "None:\"U\",\n",
    "\"(null)\":'U',\n",
    "\"ASIAN / PACIFIC ISLANDER\":\"Other\",\n",
    "\"AMERICAN INDIAN / ESKIMO\":\"Other\"\n",
    "                  }\n",
    "\n",
    "map_race_dict = {'WHITE':'W'\n",
    "                    , 'BLACK OR AFRICAN AMERICAN':'B'\n",
    "                    , 'ASIAN':'Other'\n",
    "                    ,'OTHER COMBINATIONS NOT DESCRIBED':'Other'\n",
    "                   ,'AMERICAN INDIAN OR ALASKA NATION':'Other'\n",
    "                   ,'NAT.HAWAIIAN/OTH.PACIFIC ISLAND':'Other'\n",
    "                    ,'SEPHARDIC JEWISH':'Other'\n",
    "                       ,'ASHKENAZI JEWISH':'Other'\n",
    "                    ,'NULL':'NA'\n",
    "                  , '(null)':'NA'\n",
    "                  , 'DECLINED':'NA'\n",
    "              }\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "id": "5a0f30d7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:29.506283Z",
     "start_time": "2023-11-03T20:32:29.363582Z"
    }
   },
   "outputs": [],
   "source": [
    "#Solar implementation on combining race and ethnicity\n",
    "def fe_ethnicity_race(mrn, race, ethnicity):\n",
    "    if race in ('NA','U', 'D'):\n",
    "        if ethnicity == 'W':\n",
    "            race = 'White'\n",
    "        elif ethnicity == 'B':\n",
    "            race = 'Black'\n",
    "        elif ethnicity in ('H', 'S', 'O'):\n",
    "            race = 'Hispanic'\n",
    "        elif ethnicity in ('NULL', '', 'U', 'D', '2'):\n",
    "            race = 'Unknown'\n",
    "        else:\n",
    "            race = 'Other'\n",
    "    else:\n",
    "        if race == 'W':\n",
    "            if ethnicity == 'H':\n",
    "                race = 'Hispanic'\n",
    "            else:\n",
    "                race = 'White'\n",
    "        elif race == 'B':\n",
    "            if ethnicity == 'H':\n",
    "                race = 'Hispanic'\n",
    "            else:\n",
    "                race = 'Black'\n",
    "        elif race == 'O':\n",
    "            race = 'Hispanic'\n",
    "        else:\n",
    "            race = 'Other'\n",
    "    \n",
    "    empi2race_eth[mrn] = race\n",
    "#    mrn2race_fe[race][mrn] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "id": "5631a244",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:29.626194Z",
     "start_time": "2023-11-03T20:32:29.510467Z"
    }
   },
   "outputs": [],
   "source": [
    "covnopro_limited_rel_df[\"ethnicity_mapped\"] = covnopro_limited_rel_df[\"ethnic_code\"].map(map_ethnic_dict).fillna(covnopro_limited_rel_df[\"ethnic_code\"])\n",
    "covnopro_limited_rel_df[\"race_mapped\"] = covnopro_limited_rel_df[\"race_code\"].map(map_race_dict).fillna(covnopro_limited_rel_df[\"race_code\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 150,
   "id": "9ef39af6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:31.144131Z",
     "start_time": "2023-11-03T20:32:29.633727Z"
    }
   },
   "outputs": [],
   "source": [
    "empi2race_eth = defaultdict(dict)\n",
    "\n",
    "for row in covnopro_limited_rel_df.iterrows():\n",
    "    \n",
    "    mrn = str(row[1][\"ptid\"])\n",
    "    race = row[1][\"race_mapped\"]\n",
    "    ethnicity = row[1][\"ethnicity_mapped\"]\n",
    "    \n",
    "    fe_ethnicity_race(mrn, race, ethnicity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "id": "cd4c84f9",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:31.159844Z",
     "start_time": "2023-11-03T20:32:31.147123Z"
    }
   },
   "outputs": [],
   "source": [
    "def assign_family_ethnicities(fam2empi, empi2demog, print_breakdown=True):\n",
    "    \"\"\"\n",
    "    Assign each family a single ethnicity. This makes the (incorrect) assumption that everyone in the same\n",
    "    family will have the same ethnicity. However, because there are so much missing data, this is better\n",
    "    than nothing at all.\n",
    "\n",
    "    fam2empi        dict, map from family id to a set of patients identifiers that belong in that family\n",
    "    empi2demog      dict, map from patient identifiers to the demographic data (including self-reported race/ethnicity data)\n",
    "    \"\"\"\n",
    "\n",
    "    print(\"Identifying most commonly reported ethnicity for each family...\")\n",
    "\n",
    "    eth2fam = defaultdict(set)\n",
    "    fam2eth = dict()\n",
    "   \n",
    "    for fid, members in tqdm(fam2empi.items()):\n",
    "        \n",
    "\n",
    "        family_ethnicities = [empi2race_eth[e] for e in members if empi2race_eth.get(e,0) != 0 if empi2race_eth[e] != 'Unknown' ]\n",
    "        if len(family_ethnicities) == 0:\n",
    "            mode = 'Unknown'\n",
    "        else:\n",
    "            mode = max(set(family_ethnicities), key=family_ethnicities.count)\n",
    "\n",
    "        eth2fam[mode].add(fid)\n",
    "        fam2eth[fid] = mode\n",
    "\n",
    "    if print_breakdown:\n",
    "        for eth, families in eth2fam.items():\n",
    "            print( \"\\t%s: %d\" % (eth, len(families)))\n",
    "\n",
    "    return eth2fam, fam2eth\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "id": "e8726430",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:32.870793Z",
     "start_time": "2023-11-03T20:32:31.162517Z"
    },
    "code_folding": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "12764it [00:01, 9098.18it/s]\n"
     ]
    }
   ],
   "source": [
    "#map family to patient ID\n",
    "fam2empi = defaultdict(list)\n",
    "\n",
    "for row in tqdm( fam_df[ fam_df[\"ptid\"].isin(covnopro_limited_rel_df[\"ptid\"].astype(int)) ]\\\n",
    "                .iterrows() ):\n",
    "    fam2empi[row[1][\"famid\"]].append(row[1][\"ptid\"].astype(str))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "id": "ee6a49cc",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:33.074025Z",
     "start_time": "2023-11-03T20:32:32.874439Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Identifying most commonly reported ethnicity for each family...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|████████████████████████████████████| 5676/5676 [00:00<00:00, 47960.30it/s]\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\tWhite: 1197\n",
      "\tBlack: 718\n",
      "\tOther: 2050\n",
      "\tHispanic: 1711\n"
     ]
    }
   ],
   "source": [
    "eth2fam, fam2eth = assign_family_ethnicities(fam2empi,empi2race_eth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "id": "f61d39b2",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:33.128327Z",
     "start_time": "2023-11-03T20:32:33.078095Z"
    }
   },
   "outputs": [],
   "source": [
    "covnopro_limited_rel_df[\"fam_race_eth\"] = covnopro_limited_rel_df[\"famid\"].map(fam2eth)\n",
    "covnopro_limited_rel_df[\"fe_race_eth\"] = covnopro_limited_rel_df[\"ptid\"].astype(str).map(empi2race_eth)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "id": "530d3e78",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:33.362047Z",
     "start_time": "2023-11-03T20:32:33.139928Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "covid_noproband\n",
      "covid\n",
      "days_hosp\n",
      "hosp_stat\n",
      "hosp_stat_nopro\n",
      "covid_highdeg_fam\n",
      "covid_highdeg_nopro\n"
     ]
    }
   ],
   "source": [
    "holder_list = []\n",
    "holder_list_cols = [\"pheno\", \"nfemale\", \"nnonfemales\", \"nfemale_cases\", \"nnon_female_cases\",\n",
    "                    \"percent_cases_female\", \"percent_cases_nonfemale\",\n",
    "                    \"percent_controls_female\", \"percent_controls_nonfemale\",\n",
    "                    \"percent_female_cases\", \"percent_nonfemale_cases\",\n",
    "                    \"percent_cases_overall\", \"variance\"]\n",
    "sample_dict = dict()\n",
    "\n",
    "for df in df_limited:\n",
    "    pheno = df[\"pheno\"].unique()[0]\n",
    "    print(pheno)\n",
    "\n",
    "    type_pheno = pheno2type[pheno]\n",
    "    \n",
    "    sample_dict[pheno] = dict()\n",
    "\n",
    "    \n",
    "    if type_pheno == \"quant\":\n",
    "        count_df = df.groupby([\"sex\"], dropna=False).count().reset_index()\n",
    "        total_f = count_df[count_df[\"sex\"] == \"F\"][\"ptid\"].sum()\n",
    "        total_nf = count_df[count_df[\"sex\"] != \"F\"][\"ptid\"].sum()\n",
    "\n",
    "        f_cases = nf_cases = None\n",
    "        f_controls = nf_controls=None\n",
    "        \n",
    "        percent_controls_overall = None\n",
    "        percent_cases_overall = None\n",
    "        \n",
    "\n",
    "        \n",
    "    else:\n",
    "        count_df = df.groupby([\"sex\", \"value\"], dropna=False).count().reset_index()\n",
    "        total_f = count_df[count_df[\"sex\"] == \"F\"][\"ptid\"].sum()\n",
    "        total_nf = count_df[count_df[\"sex\"] != \"F\"][\"ptid\"].sum()\n",
    "\n",
    "        f_cases = count_df[(count_df[\"sex\"] == \"F\") & (count_df[\"value\"] == 1)][\"ptid\"].sum()\n",
    "        nf_cases = count_df[(count_df[\"sex\"] != \"F\") & (count_df[\"value\"] == 1)][\"ptid\"].sum()\n",
    "        percent_cases_overall = round((f_cases + nf_cases) / (total_f + total_nf) * 100)\n",
    "\n",
    "        f_controls = count_df[(count_df[\"sex\"] == \"F\") & (count_df[\"value\"] == 0)][\"ptid\"].sum()\n",
    "        nf_controls = count_df[(count_df[\"sex\"] != \"F\") & (count_df[\"value\"] == 0)][\"ptid\"].sum()\n",
    "        percent_controls_overall = round((f_controls + nf_controls) / (total_f + total_nf) * 100)\n",
    "\n",
    "\n",
    "\n",
    "    percent_f = total_f / (total_f + total_nf)\n",
    "    percent_nf = total_nf / (total_f + total_nf)\n",
    "    percent_cases_f = round(f_cases / (f_cases + nf_cases) * 100,2) if f_cases is not None else None\n",
    "    percent_cases_nf = round(nf_cases / (f_cases + nf_cases) * 100,2) if nf_cases is not None else None\n",
    "    percent_f_cases = round(f_cases / total_f * 100,2) if f_cases is not None else None\n",
    "    percent_nf_cases = round(nf_cases / total_nf * 100,2) if nf_cases is not None else None\n",
    "\n",
    "    percent_controls_f = round(f_controls / (f_controls + nf_controls) * 100,2) if f_cases is not None else None\n",
    "    percent_controls_nf = round(nf_controls / (f_controls + nf_controls) * 100,2) if nf_cases is not None else None\n",
    "    percent_f_controls = round(f_controls / total_f * 100,2) if f_controls is not None else None\n",
    "    percent_nf_controls = round(nf_controls / total_nf * 100,2) if nf_cases is not None else None\n",
    "    \n",
    "    count_cases = f_cases + nf_cases if f_cases is not None else None\n",
    "    count_controls = f_controls + nf_controls if f_controls is not None else None\n",
    "    \n",
    "    x_variance = np.var(df[\"value\"])\n",
    "\n",
    "\n",
    "    sample_dict[pheno][\"Total Patients\"] = count_df\n",
    "    sample_dict[pheno][\"n Cases\"] = count_cases\n",
    "    sample_dict[pheno][\"n Controls\"] = count_controls\n",
    "    sample_dict[pheno][\"n Female\"] = total_f\n",
    "    sample_dict[pheno][\"n Not Female\"] = total_nf\n",
    "    sample_dict[pheno][\"Percent Female\"] = percent_f\n",
    "    sample_dict[pheno][\"Percent Non-female\"] = percent_nf\n",
    "    \n",
    "    sample_dict[pheno][\"Female Cases\"] = f_cases\n",
    "    sample_dict[pheno][\"Non-female Cases\"] = nf_cases\n",
    "    sample_dict[pheno][\"Female Controls\"] = f_controls\n",
    "    sample_dict[pheno][\"Non-female Controls\"] = nf_controls\n",
    "    sample_dict[pheno][\"Percent Cases\"] = percent_cases_overall\n",
    "    sample_dict[pheno][\"Percent Controls\"] = percent_controls_overall\n",
    "    sample_dict[pheno][\"Percent Cases (Female)\"] = percent_cases_f\n",
    "    sample_dict[pheno][\"Percent Cases (Non-female)\"] = percent_cases_nf\n",
    "    sample_dict[pheno][\"Percent Controls (Female)\"] = percent_controls_f\n",
    "    sample_dict[pheno][\"Percent Controls (Non-female)\"] = percent_controls_nf\n",
    "    sample_dict[pheno][\"Percent Female That Are Cases\"] = percent_f_cases\n",
    "    sample_dict[pheno][\"Percent Female That Are Controls\"] = percent_f_controls\n",
    "    sample_dict[pheno][\"Percent Non-female That Are Cases\"] = percent_nf_cases\n",
    "    sample_dict[pheno][\"Percent Non-female That Are Controls\"] = percent_nf_controls\n",
    "\n",
    "samp_df = pd.DataFrame(sample_dict).transpose()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "id": "edfabe8e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:33.381393Z",
     "start_time": "2023-11-03T20:32:33.365111Z"
    }
   },
   "outputs": [],
   "source": [
    "# set the max_columns option to 40\n",
    "pd.set_option('display.max_columns', 40)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "id": "533d5666",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:33.553052Z",
     "start_time": "2023-11-03T20:32:33.385515Z"
    },
    "scrolled": true
   },
   "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>Total Patients</th>\n",
       "      <th>n Cases</th>\n",
       "      <th>n Controls</th>\n",
       "      <th>n Female</th>\n",
       "      <th>n Not Female</th>\n",
       "      <th>Percent Female</th>\n",
       "      <th>Percent Non-female</th>\n",
       "      <th>Female Cases</th>\n",
       "      <th>Non-female Cases</th>\n",
       "      <th>Female Controls</th>\n",
       "      <th>Non-female Controls</th>\n",
       "      <th>Percent Cases</th>\n",
       "      <th>Percent Controls</th>\n",
       "      <th>Percent Cases (Female)</th>\n",
       "      <th>Percent Cases (Non-female)</th>\n",
       "      <th>Percent Controls (Female)</th>\n",
       "      <th>Percent Controls (Non-female)</th>\n",
       "      <th>Percent Female That Are Cases</th>\n",
       "      <th>Percent Female That Are Controls</th>\n",
       "      <th>Percent Non-female That Are Cases</th>\n",
       "      <th>Percent Non-female That Are Controls</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>covid_noproband</th>\n",
       "      <td>sex  value  ptid  pheno  birth_decade  race...</td>\n",
       "      <td>1580</td>\n",
       "      <td>11184</td>\n",
       "      <td>8554</td>\n",
       "      <td>4210</td>\n",
       "      <td>0.670166</td>\n",
       "      <td>0.329834</td>\n",
       "      <td>1091</td>\n",
       "      <td>489</td>\n",
       "      <td>7463</td>\n",
       "      <td>3721</td>\n",
       "      <td>12</td>\n",
       "      <td>88</td>\n",
       "      <td>69.05</td>\n",
       "      <td>30.95</td>\n",
       "      <td>66.73</td>\n",
       "      <td>33.27</td>\n",
       "      <td>12.75</td>\n",
       "      <td>87.25</td>\n",
       "      <td>11.62</td>\n",
       "      <td>88.38</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>covid</th>\n",
       "      <td>sex  value  ptid  pheno  birth_decade  race...</td>\n",
       "      <td>1580</td>\n",
       "      <td>1594</td>\n",
       "      <td>2133</td>\n",
       "      <td>1041</td>\n",
       "      <td>0.672023</td>\n",
       "      <td>0.327977</td>\n",
       "      <td>1091</td>\n",
       "      <td>489</td>\n",
       "      <td>1042</td>\n",
       "      <td>552</td>\n",
       "      <td>50</td>\n",
       "      <td>50</td>\n",
       "      <td>69.05</td>\n",
       "      <td>30.95</td>\n",
       "      <td>65.37</td>\n",
       "      <td>34.63</td>\n",
       "      <td>51.15</td>\n",
       "      <td>48.85</td>\n",
       "      <td>46.97</td>\n",
       "      <td>53.03</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>days_hosp</th>\n",
       "      <td>sex  ptid  pheno  value  birth_decade  race_...</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>214</td>\n",
       "      <td>122</td>\n",
       "      <td>0.636905</td>\n",
       "      <td>0.363095</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "      <td>None</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hosp_stat</th>\n",
       "      <td>sex  value  ptid  pheno  birth_decade  race_...</td>\n",
       "      <td>196</td>\n",
       "      <td>49</td>\n",
       "      <td>158</td>\n",
       "      <td>87</td>\n",
       "      <td>0.644898</td>\n",
       "      <td>0.355102</td>\n",
       "      <td>126</td>\n",
       "      <td>70</td>\n",
       "      <td>32</td>\n",
       "      <td>17</td>\n",
       "      <td>80</td>\n",
       "      <td>20</td>\n",
       "      <td>64.29</td>\n",
       "      <td>35.71</td>\n",
       "      <td>65.31</td>\n",
       "      <td>34.69</td>\n",
       "      <td>79.75</td>\n",
       "      <td>20.25</td>\n",
       "      <td>80.46</td>\n",
       "      <td>19.54</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>hosp_stat_nopro</th>\n",
       "      <td>sex  value  ptid  pheno  birth_decade  race_...</td>\n",
       "      <td>196</td>\n",
       "      <td>140</td>\n",
       "      <td>214</td>\n",
       "      <td>122</td>\n",
       "      <td>0.636905</td>\n",
       "      <td>0.363095</td>\n",
       "      <td>126</td>\n",
       "      <td>70</td>\n",
       "      <td>88</td>\n",
       "      <td>52</td>\n",
       "      <td>58</td>\n",
       "      <td>42</td>\n",
       "      <td>64.29</td>\n",
       "      <td>35.71</td>\n",
       "      <td>62.86</td>\n",
       "      <td>37.14</td>\n",
       "      <td>58.88</td>\n",
       "      <td>41.12</td>\n",
       "      <td>57.38</td>\n",
       "      <td>42.62</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>covid_highdeg_fam</th>\n",
       "      <td>sex  value  ptid  pheno  birth_decade  race_...</td>\n",
       "      <td>369</td>\n",
       "      <td>558</td>\n",
       "      <td>624</td>\n",
       "      <td>303</td>\n",
       "      <td>0.673139</td>\n",
       "      <td>0.326861</td>\n",
       "      <td>257</td>\n",
       "      <td>112</td>\n",
       "      <td>367</td>\n",
       "      <td>191</td>\n",
       "      <td>40</td>\n",
       "      <td>60</td>\n",
       "      <td>69.65</td>\n",
       "      <td>30.35</td>\n",
       "      <td>65.77</td>\n",
       "      <td>34.23</td>\n",
       "      <td>41.19</td>\n",
       "      <td>58.81</td>\n",
       "      <td>36.96</td>\n",
       "      <td>63.04</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>covid_highdeg_nopro</th>\n",
       "      <td>sex  value  ptid  pheno  birth_decade  race_...</td>\n",
       "      <td>369</td>\n",
       "      <td>2332</td>\n",
       "      <td>1832</td>\n",
       "      <td>869</td>\n",
       "      <td>0.678267</td>\n",
       "      <td>0.321733</td>\n",
       "      <td>257</td>\n",
       "      <td>112</td>\n",
       "      <td>1575</td>\n",
       "      <td>757</td>\n",
       "      <td>14</td>\n",
       "      <td>86</td>\n",
       "      <td>69.65</td>\n",
       "      <td>30.35</td>\n",
       "      <td>67.54</td>\n",
       "      <td>32.46</td>\n",
       "      <td>14.03</td>\n",
       "      <td>85.97</td>\n",
       "      <td>12.89</td>\n",
       "      <td>87.11</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                                                        Total Patients  \\\n",
       "covid_noproband         sex  value  ptid  pheno  birth_decade  race...   \n",
       "covid                   sex  value  ptid  pheno  birth_decade  race...   \n",
       "days_hosp              sex  ptid  pheno  value  birth_decade  race_...   \n",
       "hosp_stat              sex  value  ptid  pheno  birth_decade  race_...   \n",
       "hosp_stat_nopro        sex  value  ptid  pheno  birth_decade  race_...   \n",
       "covid_highdeg_fam      sex  value  ptid  pheno  birth_decade  race_...   \n",
       "covid_highdeg_nopro    sex  value  ptid  pheno  birth_decade  race_...   \n",
       "\n",
       "                    n Cases n Controls n Female n Not Female Percent Female  \\\n",
       "covid_noproband        1580      11184     8554         4210       0.670166   \n",
       "covid                  1580       1594     2133         1041       0.672023   \n",
       "days_hosp              None       None      214          122       0.636905   \n",
       "hosp_stat               196         49      158           87       0.644898   \n",
       "hosp_stat_nopro         196        140      214          122       0.636905   \n",
       "covid_highdeg_fam       369        558      624          303       0.673139   \n",
       "covid_highdeg_nopro     369       2332     1832          869       0.678267   \n",
       "\n",
       "                    Percent Non-female Female Cases Non-female Cases  \\\n",
       "covid_noproband               0.329834         1091              489   \n",
       "covid                         0.327977         1091              489   \n",
       "days_hosp                     0.363095         None             None   \n",
       "hosp_stat                     0.355102          126               70   \n",
       "hosp_stat_nopro               0.363095          126               70   \n",
       "covid_highdeg_fam             0.326861          257              112   \n",
       "covid_highdeg_nopro           0.321733          257              112   \n",
       "\n",
       "                    Female Controls Non-female Controls Percent Cases  \\\n",
       "covid_noproband                7463                3721            12   \n",
       "covid                          1042                 552            50   \n",
       "days_hosp                      None                None          None   \n",
       "hosp_stat                        32                  17            80   \n",
       "hosp_stat_nopro                  88                  52            58   \n",
       "covid_highdeg_fam               367                 191            40   \n",
       "covid_highdeg_nopro            1575                 757            14   \n",
       "\n",
       "                    Percent Controls Percent Cases (Female)  \\\n",
       "covid_noproband                   88                  69.05   \n",
       "covid                             50                  69.05   \n",
       "days_hosp                       None                   None   \n",
       "hosp_stat                         20                  64.29   \n",
       "hosp_stat_nopro                   42                  64.29   \n",
       "covid_highdeg_fam                 60                  69.65   \n",
       "covid_highdeg_nopro               86                  69.65   \n",
       "\n",
       "                    Percent Cases (Non-female) Percent Controls (Female)  \\\n",
       "covid_noproband                          30.95                     66.73   \n",
       "covid                                    30.95                     65.37   \n",
       "days_hosp                                 None                      None   \n",
       "hosp_stat                                35.71                     65.31   \n",
       "hosp_stat_nopro                          35.71                     62.86   \n",
       "covid_highdeg_fam                        30.35                     65.77   \n",
       "covid_highdeg_nopro                      30.35                     67.54   \n",
       "\n",
       "                    Percent Controls (Non-female)  \\\n",
       "covid_noproband                             33.27   \n",
       "covid                                       34.63   \n",
       "days_hosp                                    None   \n",
       "hosp_stat                                   34.69   \n",
       "hosp_stat_nopro                             37.14   \n",
       "covid_highdeg_fam                           34.23   \n",
       "covid_highdeg_nopro                         32.46   \n",
       "\n",
       "                    Percent Female That Are Cases  \\\n",
       "covid_noproband                             12.75   \n",
       "covid                                       51.15   \n",
       "days_hosp                                    None   \n",
       "hosp_stat                                   79.75   \n",
       "hosp_stat_nopro                             58.88   \n",
       "covid_highdeg_fam                           41.19   \n",
       "covid_highdeg_nopro                         14.03   \n",
       "\n",
       "                    Percent Female That Are Controls  \\\n",
       "covid_noproband                                87.25   \n",
       "covid                                          48.85   \n",
       "days_hosp                                       None   \n",
       "hosp_stat                                      20.25   \n",
       "hosp_stat_nopro                                41.12   \n",
       "covid_highdeg_fam                              58.81   \n",
       "covid_highdeg_nopro                            85.97   \n",
       "\n",
       "                    Percent Non-female That Are Cases  \\\n",
       "covid_noproband                                 11.62   \n",
       "covid                                           46.97   \n",
       "days_hosp                                        None   \n",
       "hosp_stat                                       80.46   \n",
       "hosp_stat_nopro                                 57.38   \n",
       "covid_highdeg_fam                               36.96   \n",
       "covid_highdeg_nopro                             12.89   \n",
       "\n",
       "                    Percent Non-female That Are Controls  \n",
       "covid_noproband                                    88.38  \n",
       "covid                                              53.03  \n",
       "days_hosp                                           None  \n",
       "hosp_stat                                          19.54  \n",
       "hosp_stat_nopro                                    42.62  \n",
       "covid_highdeg_fam                                  63.04  \n",
       "covid_highdeg_nopro                                87.11  "
      ]
     },
     "execution_count": 157,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "samp_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "id": "f04dd773",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:33.662728Z",
     "start_time": "2023-11-03T20:32:33.556336Z"
    }
   },
   "outputs": [],
   "source": [
    "#Get race/ethnicity data broken out by case and control status\n",
    "total = len( covnopro_limited_rel_df )\n",
    "\n",
    "race_eth_df = covnopro_limited_rel_df.groupby([\"fe_race_eth\",\"value\"],dropna=False ).count()\\\n",
    "[\"ptid\"].reset_index()\n",
    "\n",
    "race_eth_df.rename(columns={\"ptid\":\"n\"},inplace=True)\n",
    "\n",
    "race_eth_df[\"percent_race_eth\"] = round( (race_eth_df[\"n\"]/total)*100, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "id": "83540719",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:33.771071Z",
     "start_time": "2023-11-03T20:32:33.665555Z"
    }
   },
   "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>fe_race_eth</th>\n",
       "      <th>value</th>\n",
       "      <th>n</th>\n",
       "      <th>percent_race_eth</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Black</td>\n",
       "      <td>0</td>\n",
       "      <td>1226</td>\n",
       "      <td>9.61</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Black</td>\n",
       "      <td>1</td>\n",
       "      <td>192</td>\n",
       "      <td>1.50</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Hispanic</td>\n",
       "      <td>0</td>\n",
       "      <td>2180</td>\n",
       "      <td>17.08</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Hispanic</td>\n",
       "      <td>1</td>\n",
       "      <td>348</td>\n",
       "      <td>2.73</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Other</td>\n",
       "      <td>0</td>\n",
       "      <td>5730</td>\n",
       "      <td>44.89</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Other</td>\n",
       "      <td>1</td>\n",
       "      <td>882</td>\n",
       "      <td>6.91</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>White</td>\n",
       "      <td>0</td>\n",
       "      <td>2048</td>\n",
       "      <td>16.05</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>White</td>\n",
       "      <td>1</td>\n",
       "      <td>158</td>\n",
       "      <td>1.24</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  fe_race_eth  value     n  percent_race_eth\n",
       "0       Black      0  1226              9.61\n",
       "1       Black      1   192              1.50\n",
       "2    Hispanic      0  2180             17.08\n",
       "3    Hispanic      1   348              2.73\n",
       "4       Other      0  5730             44.89\n",
       "5       Other      1   882              6.91\n",
       "6       White      0  2048             16.05\n",
       "7       White      1   158              1.24"
      ]
     },
     "execution_count": 159,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "race_eth_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 160,
   "id": "f5f3ae00",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:33.855804Z",
     "start_time": "2023-11-03T20:32:33.773813Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "value\n",
       "0    38.852647\n",
       "1    46.291772\n",
       "Name: age, dtype: float64"
      ]
     },
     "execution_count": 160,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Average age of case and control groups for dataset\n",
    "covnopro_limited_rel_df.groupby([\"value\"],dropna=False).mean(\"age\")[\"age\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 161,
   "id": "934a034f",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:34.004838Z",
     "start_time": "2023-11-03T20:32:33.858571Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "value\n",
       "0    25.257443\n",
       "1    25.094659\n",
       "Name: age, dtype: float64"
      ]
     },
     "execution_count": 161,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Standard deviation for age of case and control groups for dataset\n",
    "covnopro_limited_rel_df.groupby([\"value\"],dropna=False).agg(np.std)[\"age\"]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 201,
   "id": "e76b5184",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T19:15:07.760859Z",
     "start_time": "2023-11-04T19:15:07.718539Z"
    }
   },
   "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>fe_race_eth</th>\n",
       "      <th>ptid</th>\n",
       "      <th>percent_raceeth</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Black</td>\n",
       "      <td>1418</td>\n",
       "      <td>11.109370</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Hispanic</td>\n",
       "      <td>2528</td>\n",
       "      <td>19.805704</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Other</td>\n",
       "      <td>6612</td>\n",
       "      <td>51.801943</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>White</td>\n",
       "      <td>2206</td>\n",
       "      <td>17.282983</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  fe_race_eth  ptid  percent_raceeth\n",
       "0       Black  1418        11.109370\n",
       "1    Hispanic  2528        19.805704\n",
       "2       Other  6612        51.801943\n",
       "3       White  2206        17.282983"
      ]
     },
     "execution_count": 201,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Totals by race/ethnicity\n",
    "race_eth_totals_df = covnopro_limited_rel_df.groupby([\"fe_race_eth\"],dropna=False ).count()\\\n",
    "[\"ptid\"].transpose().reset_index()\n",
    "\n",
    "race_eth_totals_df[\"ptid\"] = race_eth_totals_df[\"ptid\"].astype(int)\n",
    "\n",
    "race_eth_totals_df[\"percent_raceeth\"] = race_eth_totals_df[\"ptid\"]/race_eth_totals_df[\"ptid\"].sum()*100\n",
    "\n",
    "race_eth_totals_df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 200,
   "id": "6d89d4c8",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-04T19:15:00.116926Z",
     "start_time": "2023-11-04T19:15:00.107318Z"
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "ad74885c",
   "metadata": {},
   "source": [
    "## Graphs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 163,
   "id": "3d97c388",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:34.183421Z",
     "start_time": "2023-11-03T20:32:34.105632Z"
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib inline  \n",
    "import seaborn as sns\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "id": "44df96f8",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:34.315511Z",
     "start_time": "2023-11-03T20:32:34.185944Z"
    }
   },
   "outputs": [],
   "source": [
    "#Manuscript figure width requirements\n",
    "%matplotlib inline\n",
    "\n",
    "plt.rcParams['font.family'] = 'Arial'\n",
    "figsize_2colw = 175/25.4\n",
    "figsize_1colw= 85/25.4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "id": "c0ff55c7",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:34.447903Z",
     "start_time": "2023-11-03T20:32:34.319638Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.2487667371388302"
      ]
     },
     "execution_count": 165,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Average number of family members - use noproband covid susceptibility because includes max patients included\n",
    "count_fam = covnopro_limited_df.groupby([\"famid\"]).count()[\"ptid\"].reset_index()\n",
    "np.mean( count_fam[\"ptid\"] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 166,
   "id": "6afdd0b2",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:35.001842Z",
     "start_time": "2023-11-03T20:32:34.451967Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPMAAADPCAYAAAAgTufpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAckElEQVR4nO3daVhTZ9oH8H8AFa2ALLbVjtYFcUHqxhZDILigzKDWCtWigCCbdavWqi1vqw4UqtURO0UdbAHX4bq6XVqrIqhBQZGCihZlEQFxBSoKRZFAnveDFxkR4iGQjXD/PpGTw8l9on/OkjzPzWOMMRBCOj09TRdACFEOCjMhOoLCTIiOoDAToiMozIToCAozITqiS4S5vLwcAwYMQF5eHsrLyzFr1iw4OztDIBCgqKgIALBixQpMmDABIpEIIpEIjx8/Rm1tLXx9fSEUCuHg4IDMzEwN7wkh8ulEmKdPny73OYlEgpCQEPTs2RMAsGbNGsyfPx9nzpxBREQE8vLyAADZ2dlISkqCWCyGWCyGiYkJvv76a4wePRpnz57F7t27kZ+fr5b9IaQ9dCLMlZWVcp9bvXo1QkND0b9/fwBAeno6bt++jSlTpuDAgQMQiUSQSqUoLCxEcHAwBAIB4uLiAABJSUno3r07pk2bhvDwcEybNk0t+0NIe+hEmOVJSEhA3759m4WwpKQEpqamSElJwcCBA7Fp0ybU1tZi2bJl2L9/P44fP44dO3bgypUrqKysRFVVFZKSkjBjxgysXr1ag3tDyKvpdJjj4uKQnJwMkUiEy5cvw9fXF/r6+pg5cyYAYMaMGcjKykKvXr2wYsUK9OrVC0ZGRpg0aRJycnJgbm7eYl1CtJVOh/nMmTNITU2FWCzG2LFjsXfvXnh4eODo0aOy562trVFQUACBQIDGxkZIJBKkpaVh/PjxcHJyarEuIdrKQNMFqNvWrVsRGBiInTt3wsTEBAcPHoSpqSl8fHzg6OiIbt26wdfXF9bW1vjss88QGBgIPp+Pbt26Ye/evZounxC5eLowasrW1pZOgUmX12WOzIfSizu8jVmCwUqohBDV0OlrZkK6EgozITqCwkyIjqAwE6IjKMyE6AgKMyE6gsJMiI6gMBOiIyjMhOgICjMhOoLCTIiOoDAToiMozIToCAozITqCwkyIjqAwE6IjKMyE6AgKMyE6gsJMiI5QWZhf7O9048YNODk5QSgUYvHixZBKpQCAjRs3wt7eHhMnTpT1cZK3LiHk1VQS5pf7O61atQoRERE4e/YsGGM4dOgQLl68iNTUVFy4cAGJiYlYsmSJ3HUJIdxUEuaX+ztlZ2fDxcUFAODu7o6UlBSkpaXBzc0NPB4PAwcORENDAyoqKlpdtzWxsbGwtbWFra0tKioqVLEbhHQqSg9za/2dGGPg8XgAACMjIzx+/BjV1dUwMTGRrdO0vLV1WxMcHIysrCxkZWWhb9++yt4NQjodpc+bHRcXBx6Ph5SUFFl/p/LyctnzNTU16NOnD4yNjVFTU9NiuZ6eXotlhBBuSj8yt9bfyd3dHWKxGABw7NgxCIVCCAQCJCUlQSqV4tatW5BKpbCwsMC4ceNarEsI4aaWjhZbt25FUFAQ6uvrMXLkSHh6ekJfXx9CoRB8Ph9SqRQxMTFy1yWEcOsyvaaoPQ3RdfSlEUJ0BIWZEB1BYSZER1CYCdERFGZCdARnmO/cuYNr166hoKAAixYtwuXLl9VQFiFEUZxh9vb2xoMHD/DZZ59h6tSpWLlypTrqIoQoiDPMenp6cHZ2xqNHjzBv3rxmX7ckhGgPzmRKJBKsWbMGzs7OOH36NOrr69VRFyFEQZxhjo+Px9ChQ7Fu3TpUVFRgz5496qiLEKIgzjAPGTIE3bt3R0REBF5//XUYGxuroy5CiII4wxwSEoJbt24hOTkZNTU18PX1VUddhBAFcYa5qKgI//znP9GzZ0/MmDFD7mQBhBDN4gxzQ0MDKisrATyfLIDuZhOinTjHM0dEREAgEODevXtwdHREdHS0GsoihCiKM8wuLi7Iz89HRUUFLCwsZPNzEUK0i9wwL126FN9++y34fH6LAJ87d07lhRFCFCM3zJ9//jkAIDExUW3FEELaT26Yo6Oj5Z5SR0ZGqqwgQkj7yA3ziBEj1FkHIaSD5IbZ2toatra2OHHihDrrIYS0k9wwnzp1Cra2tvjvf//bbDmPx4Obm5vKCyOEKEZumNesWQPg+UCL6upq1NXVqa0oQojiOD9n9vPzQ1paGkxMTGR9oC5evKiO2gghCuAMc15eHoqKitRRCyGkAzi/aG1vb4/8/HyFNtrY2IiAgAAIBAI4OTnhjz/+oIbrhKgYZ5hNTExgZ2eH/v37o1+/frKey6/y66+/AgDS09MRERGBsLAwarhOiIpxnmafOnUKDx8+hIFB23vMvfvuu/Dw8AAAlJaWok+fPkhJSWnWRP3EiRMYPnx4mxqunzhxArNnz27P/hHSZXAema2srPDgwQOFN2xgYAA/Pz8sW7YM8+fPV3rD9djYWNja2sLW1hYVFRUK10eIruEMc3p6OgYNGoR+/fq1+TS7yZ49e1BQUICgoCA8ffpUtlwZDdeDg4ORlZWFrKws9O3bt801EaKrOMNcWFgIiUSCe/fu4d69e7h79y7nRvft24eoqCgAQK9evaCnpwdbW1tquE6ICnFeCGdkZCA+Ph4SiQSMMdy9exdJSUmv/J333nsP/v7+cHZ2hkQiQXR0NEaOHEkN1wlRIc5m6+PGjcOaNWvw448/wsbGBoWFhThw4IC66msTarZOSBtOsy0sLPDBBx/A2NgYGzZswO3bt9VRFyFEQW1qT5Obm4snT54gPz8fDx8+VEddhBAFcYb5X//6F3Jzc7F8+XJ4e3tj0aJF6qiLEKIguTfAsrKyYGtrC2tra1hbWwMAsrOz1VYYIUQxco/MTUMgAWDFihVqKYYQ0n5yw/ziTe6rV6+qpRhCSPvJDTPNj01I5yL3mvnOnTuIjY0FY0z2c5Pg4GC1FEcIaTu5Yfb29sa9e/da/ExHbEK0k9wwr1+/Xp11EEI6iFo6EqIj5IaZ+jAT0rnIDfM//vEPAMDixYvVVgwhpP3kXjN369YNdnZ2KCwsRE5ODgDIZgChLpCEaB+5YU5JScGdO3ewePFi7Ny5ExwjJQkhGiY3zPr6+hg4cCAOHTqE2NhY5ObmwsrKik67CdFSnHezQ0JCcOPGDUydOhUlJSUIDAxUR12EEAVxThtUWFiIM2fOAHg+he7EiRNVXhQhRHGcR+a6ujo8efIEAPD06VM0NjaqvChCiOI4j8wrVqzAmDFjMHr0aFy7dg0bN25UR12EEAVxhnn+/Plwd3fHzZs3MXjwYJibm6ujLkKIgtrUc8bMzAxmZmaqroUQ0gH03WxCdARnmLds2aKOOgghHcQZ5qNHj9IdbEI6Ac4wV1ZWon///nB0dASfz+f8nFkikcDHxwdCoRD29vY4fPgwNVonRA04b4A1NU5vq/3798Pc3Bz79u3Dw4cPMXbsWIwdOxYREREQiUQIDQ3FoUOH8Pbbb8sarZeVlWHOnDn4/fffZY3WX1yXejMTwo3zyGxgYICwsDAEBQXhxIkTuH///ivX9/LyQnh4OIDno6wMDAxaNE9PSUlBWlpamxqtp6SkdHQfCekSOMMcHByMgIAASCQSODs7c86h3bt3bxgZGaGmpgaenp6IiIhQeqN1gJqtE/IyzjA/ffoUkyZNAo/Hw/Dhw2FoaMi50bKyMri6usLHxwfe3t6tNk/vSKN1gJqtE/IyzjAbGhoiKSkJjY2NyMjI4AzzgwcP4Obmhk2bNiEgIAAAWm2eTo3WCVEuzhtgsbGxWL16NSorK7Flyxbs3LnzletHRkaiqqoK4eHhsmvn7du3Y/ny5dRonRAV4my2DgBXrlxBQUEBRo8ejREjRqijLoVQs3VC2nCaHRERgQ8//BDp6elYtGgRoqOj1VAWIURRnKfZv/32G9LT06Gnp4eGhgY4OTnho48+UkNphBBFcB6Z33jjDdnkBPX19XTnmBAtJffIzOfzwePxUF5ejmHDhmHMmDG4du0ajWcmREvJDXNiYqI66yCEdJDcML/99tsAgMzMTCQmJqKurk723I4dO1RfGSFEIZw3wPz8/LB27VqYmpqqox5CSDtxhnnYsGFYuHChGkohhHQEZ5jnzJmDefPmYdSoUbJlX3zxhUqLIoQojjPMMTExmDNnjtwBD4QQ7cAZZnNzc6xdu1YdtRBCOoAzzBYWFggJCcH48eNl44yDg4NVXhghRDGcYba0tAQAzhlGCCGaxRlmf39/ddRBCOkgzjDPnTsXPB4PUqkUxcXFGDZsGNLS0tRRGyFEAZxhPn/+vOznR48e0fUyIVpKofY0JiYmuHnzpqpqIYR0AOeRuWn0FGMMFRUVmDJlijrqIoQoiDPML46eMjQ0xBtvvKHSgggh7SM3zHv37pX7S76+viophhDSfnLDfP369WaPGWOIj49Hr169KMyEaCG5YY6KipL9XFRUBD8/P3h4eNCEfoRoqTYNtIiOjsa2bdvg4eGhjpoIIe0gN8x37tyBv78/zMzMkJmZSZMTEKLl5IbZ2toaPXr0wKRJk7BkyZJmzx08eFDlhRFCFCM3zIcOHVJnHYSQDpIb5qYeye114cIFrF27FmKxGDdu3MDChQvB4/EwevRoxMTEQE9PDxs3bsRvv/0GAwMDREdHw97eXu66hJBXU0lKNm/ejMDAQNmMnqtWrUJERATOnj0LxhgOHTqEixcvIjU1FRcuXEBiYqLsVL61dQkh3FQS5qFDh+Lnn3+WPc7OzpYd6d3d3ZGSkoK0tDS4ubmBx+Nh4MCBaGhoQEVFRavrtkZbmq0nJCRAJBJBJBLB0dERhoaGyMrKwpQpU+Ds7IypU6fizz//BAB88skn4PP5sLOzw+7duzVWM9FNnB9NtcecOXNQUlIie8wYk81SYmRkhMePH6O6urpZd4ym5a2t25rg4GDZCC5bW1tV7EabLFy4UDZ76ZIlSxAQEIA1a9YgMjISjo6O+Omnn1BQUIC6ujrcuHED58+fx7Nnz2BtbQ1PT0/6lIAojVouRl+85q2pqUGfPn1gbGyMmpqaFstbW7czyMrKQm5uLnx8fFBeXo5ff/0VIpEI58+fh729Pfh8PuLi4gAAPB4PjY2N6Natm4arJrpELWEeN24cxGIxAODYsWMQCoUQCARISkqCVCrFrVu3IJVKYWFh0eq6nUFkZCTWr1+Phw8fIjc3F1OmTMHp06dRVVWFPXv2wNDQEKamppBIJPDz80NwcDB69+6t6bKJDlHJafbLtm7diqCgINTX12PkyJHw9PSEvr4+hEIh+Hw+pFIpYmJi5K6r7R49eoT8/Hy4urri6dOnMDIygqurKwDAw8MDycnJCAgIQFVVFTw9PSESifDpp59quGqia3iMMabpIjrK1tYWWVlZr1znUHpxh19nlmBwq8sPHz6MlJQUfPPNN7J6tm3bBqFQiJUrV8LS0hIBAQEQCAT4+OOPMX/+/A7XQsjL6ANcJcjPz8eQIUNkj7///nt8+umncHR0xP379xEUFIRdu3bh5s2b2L17t+zud3Fxx//AENKEjswKkHdkJkQbqOWaWZfRHwmiLeg0mxAdQWEmREdQmAnRERRmQnQEhZkQHUFhJkRHUJgJ0REUZkJ0BIWZEB1BYSZER1CYCdERFGYtd+HCBYhEIgDAvHnzZCOuBg0ahHnz5gEAwsLC4ODgAEdHR9nEDqTroYEWWmzz5s3Yt28fXnvtNQD/a69bVVUFV1dXbNu2DZcuXUJGRgYyMjJQWlqKWbNmIScnR5NlEw2hI7MWe3mW0ybr16/HsmXL0K9fP4wbNw5JSUng8XgoLS3tNHOmEeWjMGuxOXPmtJj0r7y8HCdPnpTNCAoABgYGCAsLg4eHB/z9/dVcJdEWFOZO5scff4S3tzf09fWbLf/yyy9x9+5dfP311ygqKmrz9saPHy+7Dvf398cvv/yCoUOHypalpqYqexeIilCYO5mUlBS4u7vLHp86dUrWDcTQ0BDdunVrczufuro6MMYgFoshFosRHx+P7OxsbN68WbasvW2K2nLjrr2ePXsGb29vODo6ws3NDYWFhR3aXhOJRAIfHx8IhULY29vj8OHDWrlNeegGWCfz8nxjLi4u+OGHHyAQCNDY2IglS5Zg8OC2zVySk5ODJ0+ewM3NDQ0NDYiMjER2djYuXbok6/21adMmGBgo9t+kLTfuOmL37t3o3bs3MjIykJ+fj6VLlyIpKalD2wSA/fv3w9zcHPv27cPDhw8xduxYzJw5U+u2KQ+FWcsNGjQIGRkZsse5ubnNntfX18fOnTvbte1evXph9erVCAwMRGFhIdzd3RESEgJPT08MHjwYoaGh2LVrF5YuXarQdptu3Pn4+DRb/uKNu464du2a7Oxk+PDhuH79eoe218TLy0s2tTNjTOE/YurapjwUZi2kjHnFAO65xaysrGBpaQkejwcrKyuYm5vjgw8+wIABA57//qxZ+OmnnxR+3ZfbEwH/u3HX0aMyAIwdOxZHjhzBu+++iwsXLuDOnTtobGxscR9BUU1NCWpqauDp6YmIiIgO16qKbcpD18xdWFxcHD7++GMAwN27d/H48WM4ODjg9u3bAICTJ09iwoQJSnkteTfu2iMgIADGxsYQCoX45ZdfMGHCBKVsFwDKysrg6uoKHx8feHt7a+02W0NH5i5s0aJFWLhwIZycnMDj8RAfH4+//voL7733Hnr27IlRo0YhKChIKa+VkpKC//u//1PKtn7//XdMnjwZ27ZtQ1ZWFkpLS5Wy3QcPHsDNzQ3ffvstJk+erLXblEcrwyyVSvHhhx8iJycHPXr0wHfffQdLS0tNl6VzunfvjoMHD7ZY7ubmpvTXevnGXUcMGzYMn3/+Ob788kv06dMH33//vVK2GxkZiaqqKoSHhyM8PBzA835nPXv21KptyqOVk+D//PPPOHz4MBISEpCRkYGoqKhXNl3X5CT4qtiuqq6Z1XUtTjRDK4/MaWlpmD59OgDA0dGRM6hEszrTH0pdppVhrq6uhomJieyxvr4+Ghoamt3Wj42NRWxsLAAgLy9PKQ3XKyoq0LdvX7nPh3f4FZS3Xa5a27vdtlBFvZ2pVk2ysLDA8ePHW31OK8P8ciN2qVTa4vO54OBgBAcHK/V123K6ri06U61A56q3M9X6Iq38aEogEODo0aMAgIyMDNjY2Gi4IkK0n1YemWfPno3k5GRMnDgRjDHEx8druiRCtJ5WhllPTw+7du1S++sq+7RdlTpTrUDnqrcz1foirfxoihCiOK28ZiaEKI7CTIiO6PJhVufgcWUpLy/HgAEDkJeXp+lSOEVFRYHP52PChAlK+9qlqkgkEnh7e2PixIkQCoWd4v19UZcPc9Pg8bNnz+L48eMKj91VN4lEgpCQEJV8t1fZxGIxzp07h/T0dKSmpqKsrEzTJb3S0aNH0dDQgHPnzuGLL75AWFiYpktSSJcPs5eXl+wL8KoePK4Mq1evRmhoKPr376/pUjglJSXBxsYGs2fPxowZM+Dh4aHpkl7JysoKDQ0NkEqlqK6ubjGZorbT7v+5aqDOweMdlZCQgL59+2LatGmIiorSdDmcKisrUVpaiiNHjqC4uBgzZ85EXl4eeDyepktrVe/evVFSUoIRI0agsrISR44c0XRJCunyR2ZAfYPHOyouLg7JyckQiUS4fPkyfH19cf/+fU2XJZe5uTmmTZuG7t27Y/jw4TA0NERFRYWmy5Jr27ZtmDZtGgoKCpCTkwM/Pz/U1dVpuqw26/JHZnUOHu+oM2fOyH4WiUTYtWsX3nzzTQ1W9GpOTk7Yvn07Vq1ahXv37qG2thbm5uaaLksuU1NT2am1mZkZJBIJGhsbNVxV23X5MKtz8HhX4+HhgTNnzsDe3h5SqRQxMTFKm95HFVauXImAgAAIhULU19cjMjJSNsNoZ0DfACNER9A1MyE6gsJMiI6gMBOiIyjMhOgICjMhOoLCrARisRgmJibNvnu8bt06JCQktHubJSUlcHR0VEJ1LTU0NMDV1RUTJ05EVVWVbPmGDRtgZWUl69goEomQmZnZ7tf56quvkJmZiYSEBKxbt45zfbFYDB6PJ2s01+Sdd95p1o9aUap8L7VJl/+cWVl69OgBf39/JCcna+3XFZvcvXsX1dXVyM7ObvHcqlWrEBoaqpTXaQrwtWvX2vw7I0aMQGJioqzt69WrV1FbW6uUenQdHZmVZNKkSTAzM0NMTEyz5S8fFRwdHVFSUoINGzZgwYIFmD59Ouzs7JCQkIAZM2Zg+PDhsq6PFRUVmDlzJhwcHGRfaCkrK4O7uztEIhHc3d1RVlaGkpIS2NjYQCQSYfPmzc1e/8CBA7Czs4OTkxP8/f0hkUgQGhqKwsJChISEtGnf/vjjD7i5uWHy5MkYM2YMzp07BwCwtLSEr68vHB0dERYWhqVLl8Le3l7W/XHhwoXNpoWNjY3FJ598AgBobGyEjY1Ni69LjhkzBqWlpXj8+DGA56Pa5s+fL3v+hx9+AJ/Ph5OTk+yPhSrfyx07dsDBwQF8Ph/Lly9v0/ulMYx02OnTp9ncuXNZZWUlGzJkCCssLGRr165l8fHxrLi4mDk4OMjWdXBwYMXFxWz9+vUsMDCQMcZYVFQU8/LyYowxFhcXx1asWMGKi4vZ66+/zh49esQaGhoYn89nly9fZu+//z47evQoY4yxlJQU5u3tzYqLi5mFhQV79uxZs7oqKyvZ0KFDWXV1NWOMsY8++oj9+9//blFTk/Xr17Nhw4YxFxcX5uLiwpYuXcoYYywxMZFduXKFMcbYgQMHZHXr6+uz0tJSVl9fz3r37s1yc3OZVCplgwYNYlVVVczPz48dO3aMxcfHs7Vr17Lq6mpmaWnJGhoa2JEjR9iyZctafR/Dw8NZXFwck0qlTCQSsWPHjjE/Pz/2559/spEjR7La2lrGGGMLFixgJ06cUOl7aWtryzIzMxljjO3YsYNJJBIF/3eoD51mK5G5uTmio6Ph5+cHgUDQ6jrshS/cjR8/HgDQp08fjBo1CsDz7wc3Ha3GjBkjawZgb2+PgoICXL16FZGRkdi0aRMYY7LvEg8ePBjdu3dv9lo3b96EtbU1jIyMAADOzs44ceLEK4citnaa/dZbbyE8PBw9e/ZETU0NjI2NZfs7cOBAAMBrr70m2wcTE5NWBygYGRnBxcUFSUlJiI+PxxdffNFqDd7e3li8eDGGDBkCoVAoW37jxg1UVFTg73//O4DnI92KiopU+l7Gx8djy5YtKC4uBp/Pb/bvp23oNFvJmk7vmm5+GRoaory8HI2NjXj06BGKi//XcoXr2vr69ev466+/0NDQgAsXLsDa2hojRozApk2bIBaL8Z///AdeXl4Ans9o+rLBgwfj2rVrsmvO1NRUWFlZKbxPy5cvx8aNG7Fnzx7Y2NjI/kO3595AUFAQvvvuO5SXl+Odd95pdZ0hQ4agtrYW33zzDRYsWNBsfwYMGIDk5GSIxWIsW7ZMdgmjqvdy9+7d2LVrF1JTU3Hp0iXZJYY2oiOzCkRHR+PkyZMAgDfffBNTp06FnZ0dhg4dqlA3SzMzM8ydOxcVFRWYO3cuRo0ahS1btmDx4sWoq6vD06dPsX37drm/b2FhgY0bN8LV1RV6enqwtLTEV199pfCwyQULFsDLywumpqb429/+hsrKSoV+/0UODg64ceMGlixZ8sr15s6di3379sHKygo3b94EAPTt2xerVq2Ci4sLGhsbMWjQILz//vttet32vpc2NjYQCoUwMjLCW2+9BQcHB8V3Wk1ooAVRK6lUCoFAgKSkJNnpOlEOOs0malNcXIzx48dj3rx5FGQVoCMzITqCjsyE6AgKMyE6gsJMiI6gMBOiIyjMhOiI/wetjHiFQRmlwwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 240.945x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "hist_df = count_fam.groupby(\"ptid\").count().reset_index()\\\n",
    ".rename(columns={\"ptid\":\"Number of Family Members\",\"famid\":\"Number of Families\"})\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(figsize_1colw,3)\n",
    "bars = ax.bar( hist_df[\"Number of Family Members\"], hist_df[\"Number of Families\"],color=\"lightsteelblue\", )\n",
    "\n",
    "ax.bar_label(bars)\n",
    "sns.despine()\n",
    "\n",
    "plt.xlabel(\"Number of Family Members\")\n",
    "plt.ylabel(\"Number of Families\")\n",
    "\n",
    "save_and_track_figure(fig,\"number_of_family_members.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "id": "37148d1c",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:35.463154Z",
     "start_time": "2023-11-03T20:32:35.004858Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "17336\n",
      "17336\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPMAAADPCAYAAAAgTufpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAdhElEQVR4nO3de1SU5fr/8feAIBGCBmjZF4wUFA8BipxJPISSZJlmlqloipjunUZlpqsgSU1Tce0smlLKU66dVtouRdEAQZHQUldESeIhhZRyi4LEYeb3hz9mg8wwDDADjtdrLdeCmXsersH5cD8zz31QqNVqNUKI255FWxcghGgdEmYhzISEWQgzIWEWwkxImIUwExJmIczEHR3mZcuWERgYyKBBg1i/fr3m9vnz55OUlKT5fs2aNfj7++Pv7098fHy9Y+Tn5+Pg4EBFRYXJ6hZCG7MI86hRowx+TFpaGocOHSIrK4v09HTOnz/P5cuXiYiIYNeuXZp2p0+fZsuWLRw6dIjs7Gz27t3LiRMnACgtLSU2NpaOHTu22nMRornMIswlJSUGPyYlJYUBAwYwduxYHnvsMSIjI7l+/TpxcXFMnjxZ087FxYU9e/ZgaWmJQqGgqqoKGxsb1Go10dHRLF26FFtb29Z8OkI0S4e2LqCtlJSUcPbsWf7zn/9QWFjImDFjyM/Px83Njd27d2vaWVlZ4eTkhFqt5pVXXsHHxwcPDw/i4uIYPXo0Xl5ebfgshPgfs+iZm8PR0ZGRI0dibW1N7969sbGx4fLly1rbVlRUMGnSJK5du8b7778PwObNm1m/fj1hYWEUFxcTHh5uyvKFaOCODXNISAh79uxBrVZz8eJFysrKcHR0bNBOrVbz+OOP4+XlxYcffoilpSUABQUFpKWlkZaWxr333svevXtN/RSEqOeOPc2OjIwkIyMDPz8/VCoV69at0wS1rq+++or09HT+/vtvzel37afgQrQnCnOYNeXr60tubm5blyFEmzLbnrnw0vW2LqEBt652bV2CMGNGe89864CMgoICQkJCCA0NZfbs2ahUKgDi4+Px8/MjKCiInJwcAJ1thRC6GSXM2gZkvPTSSyQkJHDw4EHUajU7d+7k2LFjpKenc+TIEbZt28acOXMAtLYVQjTOKGHWNiDj6NGjDBkyBICIiAhSU1PJzMwkPDwchUKBq6sr1dXVXL58WWtbIUTjjPKeWduADJVKhUKhAKBTp05cvXqV0tLSepeDam9Xq9UN2t5KqVSiVCoBdF4fFuJOYpSeWduAjLqBvHbtGp07d8be3p5r1641uN3CwqLBbbeKjo4mNzeX3NxcnJ2djfE0hLitGCXM2gZkDB8+nLS0NAB2795NaGgowcHBpKSkoFKpOHfuHCqVCicnJ3x8fBq0FUI0ziin2doGZLi5uTFz5kwqKyvx9PRk/PjxWFpaEhoaSmBgoKYdwKpVqxq0FUI0zmwHjch1ZnGnuWPHZgthbiTMQpgJCbMQZkLCLISZkDALYSYkzEKYCQmzEGZCwiyEmZAwC2EmJMxCmAkJsxBmQsIshJnQG+aysjJ+//13/vjjD5YsWcLZs2dNUZcQwkB6wzx+/HiOHj3KK6+8gpWVFdHR0aaoSwhhIL1hLi8vZ8yYMfz++++89tpr1NTUmKIuIYSB9Ia5srKStWvXMmjQIPLy8igrKzNFXUIIA+kN87vvvktRURGLFy/mwIEDrF271hR1CSEMpDfMwcHBBAQEkJSUxAMPPICfn1+TDjxw4EDCwsIICwtj2rRpZGdn4+/vT3BwMPHx8QCoVCpiYmIIDAwkLCyMgoICAK1thRCN07sG2IwZM7h+/TqBgYFs3LiR/fv3s2bNmkYfU1FRgVqt1izKB+Dt7c2OHTt48MEHGT16ND/88AOFhYVUVFRw+PBhsrOziY2NZefOncTExDRo6+Pj0+InK4Q50xvmkydPcuTIEQBefPFFAgIC9B70+PHjlJeXEx4eTnV1NXFxcfz999/07NkTgJEjR5KamkpRURGjRo0CICAggNzcXEpLS7W2lTAL0Ti9p9m9evWisLAQgEuXLuHq6qr3oLa2trz88sukpKSQlJTEtGnTsLW11dxfdxF8BwcHze2WlpaUlpZib2/foO2tlEolvr6++Pr6yiL4QtCEnvnw4cP06dMHV1dXLly4QMeOHbnvvvtQKBRcvHhR62M8PDzo1asXCoUCDw8PHBwc+OuvvzT31y5sX15eXm8RfJVKpXNh/FtFR0drrnn7+vo2+QkLYa709synT5/m77//5tSpU5SXl3PlyhWKiop0Bhlgw4YNxMbGAnDx4kXKy8u5++67+e2331Cr1aSkpGgWwf/222+Bmx96DRgwAHt7e6ytrRu0FUI0TmfPnJCQwOLFi3nmmWc0+z7V2rp1a6MHff7554mKiiIkJASFQsGGDRuwsLBg0qRJ1NTUEB4ejr+/P4MHD2bfvn0EBQWhVqtJTk4GICkpqUFbIUTjdC6Cf/z4cby8vEhPT29wX+0Oje2FLIIvRCM9s5eXFwA9evRg+/btlJeXa+5rb2EWQjThPfMzzzxDWVkZ3bp10/wTQrQ/ej/NtrW15c033zRFLUKIFtAZ5l9//RWAbt268dlnnzFw4EDNB2EeHh6mqU4I0WQ6wzxr1izN10qlUvO1QqHgwIEDxq1KCGEwnWH+7rvvGtx2/vx5XFxcjFqQEKJ59L5nXrlyJZ07d+a///0vycnJjBo1itWrV5uiNiGEAfSGeceOHWRkZDBq1Cjy8vIYOnSoKeoSQhhI76UpS0tLiouLNZekbty4YfSihBCG09sz1y4wsHnzZubPn8/o0aNNUZcQwkA6h3NqU1VVhZWVlTHraRYZzilEE3rmffv2sWbNGioqKjS3yaUpIdofvWGeP38+iYmJcklKiHZOb5hdXV0ZMWKEKWoRQrSA3jB37dqVmJgYfHx8NMM5ZVcLIdofvWF2c3MDoLi42OjFCCGaT+915jfffBNfX1/uuusuvL29ZQaVEO2U3jAvXLiQ5ORkrK2t+fTTT3n55ZdNUZcQwkB6w5yRkcH27duZN28eO3bs4ODBg0068KVLl3BxcSE/P5+CggJCQkIIDQ1l9uzZqFQqAOLj4/Hz8yMoKIicnBwAnW2FEI3TG+aqqipNoNRqdYPF/XQ9ZtasWdx1110AvPTSSyQkJHDw4EHUajU7d+7k2LFjpKenc+TIEbZt28acOXN0thVC6Kc3zE8//TTBwcHMnz+fkJAQnn76ab0Hffnll4mJiaF79+4AHD16VLNuWEREBKmpqWRmZhIeHo5CocDV1ZXq6mouX76sta02sgi+EPXpDXNsbCwfffQRwcHBKJVK5s+f32j7Tz75BGdnZ0aOHKm5rW6Prms3i9rbtbXVJjo6mtzcXHJzc3F2dtb/TIUwczovTX388cfMmDGDhQsXasJ17NgxAJYuXarzgBs2bEChUJCamsqPP/7IlClTuHTpkub+2h0qdO1cYWFh0eA2IYR+Onvm2uGbffr0oXfv3pp/ffr0afSAGRkZpKenk5aWhre3Nxs3biQiIkKzI+Tu3bs1u1mkpKSgUqk4d+4cKpUKJycnfHx8GrQVQuins2euPU3+/vvvee+99zS3T5kyhSlTphj0Q1atWsXMmTOprKzE09OT8ePHY2lpSWhoKIGBgahUKtatW6ezrRBCP51TINetW0dCQgJXrlzhnnvuAW6+9+3bty/79+83aZH6yBRIIZown3np0qW8/vrrpqqnWSTMQjRhbHZMTAyfffYZVVVVqNVqLl68yMKFC01RmxDCAHrDPHbsWDw9PTlx4gR33XVXvU3ThRDth97rzGq1mqSkJPr06cO+ffvqbZouhGg/9Ia5Q4cOVFRUUFZWhkKhoLq62hR1CSEMpDfMc+bMYc2aNYSHh+Pi4qKZ3yyEaF8MWp2ztLQUe3t7Y9bTLPJpthCNfAAWGBioc4bUoUOHjFaQEKJ5dIZ527ZtpqxDCNFCOsPco0cPAC5cuMCCBQu4dOkSTz31FA899JDmPiFE+6H3A7Do6GimT59OVVUVDz/8MC+++KIp6hJCGEhvmG/cuMGwYcNQKBT07t0bGxsbU9QlhDCQ3jDb2NiQkpJCTU0N2dnZEmYh2im9YVYqlSQnJ1NSUsK7775LUlKSKeoSQhhI79js//u//6v3yfbu3bt54IEHjFmTEKIZdPbMn3zyCd27d6dnz5788MMPXL16lQkTJvDqq6+asj4hRBPp7JlXr17NTz/9RFFREbGxsVy8eJHHH3+cLVu2mLI+IUQT6QzzPffcQ5cuXejSpQt5eXkkJSURERHRpIPW1NQwc+ZMfvnlFxQKBUlJSdjY2BAVFYVCoaB///6sW7cOCwsL4uPj+eabb+jQoQOJiYn4+flRUFCgta0QQjedCakbnh49ejQ5yABff/01AFlZWSQkJLBo0SJZCF8II9PZM//555/s27cPlUpFaWkpe/fu1dwXHh7e6EGfeOIJIiMjATh79iydO3cmNTW13uL2e/fupXfv3k1aCH/v3r2MHTu2xU9WCHOmM8wDBw5k69atAPj4+PDZZ58BoFAo9IYZbs6Dnjp1Kl9++SXbt29n3759WhfCd3R01DzGkIXwlUolSqUSQHa0EIJGwpycnNzig3/66ae88847+Pv7c+PGDc3trbEQfnR0tGbTd19f3xbXKsTtziifKm3atIlly5YBYGtri4WFBb6+vrIQvhBGpLNnvnr1ar29oAzx5JNPMm3aNB5++GGqqqpITEzE09NTFsIXwoh0rjQSEhJCZmYms2fP5oMPPjB1XQaRlUaEaKRntrKyYvDgwZw6dYrjx48D/9vNUVYaEaL90Rnm1NRULly4oOmZDVgqTAjRBnSG2dLSEldXV3bu3IlSqeSnn37Cw8OD2bNnm7I+IUQT6f00e9asWRQUFPDII49w5swZZsyYYYq6hBAG0jsF8tSpU2RkZAA3R3YFBQUZvSghhOH09swVFRWUl5cDN5cQqqmpMXpRQgjD6e2ZX3zxRby8vOjfvz95eXnEx8eboi4hhIH0hnnSpElERERw+vRp3Nzc6o2lFkK0H3rDDDfnNt9zzz3GrkUI0QIy418IM6E3zO+++64p6hBCtJDeMH/77bfyCbYQtwG975lLSkro3r07bm5uKBQKGZstRDulN8y163kJIdo3vWHu0KGD7AIpxG1AdoEUwkzILpBCmIlW3wWyqqqKyZMnExoaip+fH7t27aKgoICQkBBCQ0OZPXs2KpUKgPj4ePz8/AgKCiInJwdAZ1shROMM3gVS3xJCmzdvxtHRkYMHD7Jnzx7mzp0rC+ALYQJN2gXy9ddf59dff6V///64ubk12v6pp57SLMCnVqvp0KGD1kXtZQF8IVqX3p45ISGBF154gaysLJ5//nkSExMbbW9nZ0enTp24du0a48ePJyEhQeui9qWlpfVW/zRkAXy4ecbg6+uLr6+vLIIvBE0I8zfffENGRgZr1qwhPT293l7Nupw/f56hQ4cyefJknn32Wa2L2rdkAXy4+Sl7bm4uubm5ODs7661JCHOnN8zdunXTLE5QWVmpNzh//PEH4eHhvPPOO0yfPh1A66L2sgC+EK1L53vmwMBAFAoFly5dwt3dHS8vL/Ly8vTOZ166dClXrlxhyZIlLFmyBIC1a9fyz3/+UxbAF8KIdC6Cf/bsWZ0Pam8jwGQRfCEa6ZlrA5uTk8O2bduoqKjQ3Pf+++8bvzIhhEH0XpqaOnUqCxYsoEuXLqaoRwjRTHrD7O7uTlRUlAlKEUK0hN4wjxs3jokTJ9K3b1/NbW+88YZRixJCGE5vmNetW8e4ceN0Xu8VQrQPesPs6OjIggULTFGLEKIF9IbZycmJWbNmMXDgQM0wy+joaKMXJoQwjN4w9+rVC4Di4mKjFyOEaD69YZ42bZop6hDNcOTIERYsWKAZ/vrll1/y+eefs3XrVgDCwsI0bfPz84mKimL58uVaHytuf3rD/PTTT6NQKFCpVBQWFuLu7k5mZqYpahONWLFiBZs2beLuu+8Gbu4JlpKSgre3t6ZNbVBPnz7NhAkTWLx4sdbHCvOgd6LF4cOHOXToENnZ2fzyyy90797dFHUJPXr27MkXX3yh+T4oKEjnwhHz5s3jnXfewc7OTutjhXkwaHsaBwcHTp8+baxahAHGjRuHlZWV5vvaM6hbnThxgtLSUoYPH67zscI86D3Nrp09pVaruXz5MiNGjDBFXaKVbN68mZkzZ7Z1GcIE9Ia57mIENjY2dOvWzagFida1f/9+GSdwh9AZ5o0bN+p80JQpU4xSjGh9xcXFsqf2HULnfOaFCxfW+16tVpOcnIytrS2FhYUmKa6pZD6zEI30zMuWLdN8/dtvvzF16lQiIyP1LugnWi4972Jbl9DAkL5yFaO9a9JEi8TERNasWUNkZKQpahJCNIPOS1MXLlwgPDycgwcPkpOTY3CQjxw5ohmBJDtaCGF8OsPcr18/jh8/jkKhYM6cOTz77LOaf/qsWLGCGTNmaJYakh0thDA+nafZLQlR7QijyZMnAxhlRwulUolSqQSQRfCFoJEw1waqOcaNG8eZM2c03+va0aLuJRNDd7SIjo7WTMX09fVtdq1CmAuDhnM2+4cYYUcLIUR9Jgmz7GghhPHpvTTVGrTtUiE7WgjRunSOALudmNsIMBk0IprDJKfZQgjjkzALYSYkzEKYCQmzEGZCwiyEmZAwC2EmJMxCmAkJs2hTdafKHjt2DD8/P0JDQ/nHP/7R6FRZ0ZCEWbSZW6fKRkdHk5iYyMGDB3FwcGDr1q06p8q2tbp/hPLy8ggJCSE4OJioqCiqq6uBmwt7DB48GD8/P/79738bvSYJs2gzty7G//vvvxMUFARAcHAwmZmZZGZmap0q25Zu/SP0+uuvs3TpUrKysgD4+uuvKSkp4YMPPuDQoUPs37+f2NhYjD3YUsIs2syti/E/+OCDpKenAzcDUVZWRmlpKQ4ODpo2jU2LNZVb/wjt2LGDhx9+mMrKSoqLi3FwcMDJyYkff/wRKysriouLsbGx0bpJQWuSMIt2Izk5mWXLljF8+HC6du2Kk5OTzqmybenWP0KWlpacPXuWfv36UVJSgpeXFwAdOnTgvffeIyAggOeee87odUmYRbvxzTffsGXLFvbv38+ff/7JI488onOqbHvTo0cPTp06RUxMDC+99JLm9rlz51JUVERGRgbfffedUWuQMIt2w93dneHDhxMUFIS9vT2PPvoogwYN0kyVHTdunGaqbHsyZswYTp06Bdx8G2BhYcEvv/zCk08+iVqtxsrKio4dO9ZbeMMYTDKfWQhdHnjgAbKzswF47LHHeOyxxxq0iYuLIy4uzsSVNd1rr71GVFQU1tbW2Nra8vHHH3Pffffh5eWl2astIiKiRUtxNYXMZzYhc5/P/F7qSRNUYri5Iwa0dQkmIT2zEMDzHxn3/WxzrJ851KD27fI9s0qlIiYmhsDAQMLCwigoKGjrkoRo99plmL/66isqKio4fPgwy5cvJzY2tq1LEqLda5dhzszMZNSoUQAEBAQ0eD8shGioXb5nvnXUj6WlJdXV1XTo8L9y6+5okZ+fb9SF8C9fvoyzs7PRjm8st2vd0Lq1f9IqR2m61qrd98OGtzk5ObFnzx6t7dtlmG8d9aNSqeoFGervaGFs2j4tvx3crnWD1N4c7fI0Ozg4mG+//RaA7OxsBgy4My4tCNES7bJnHjt2LPv27SMoKAi1Wk1ycnJblyREu9cuw2xhYUFSUlJbl6FhqtP51na71g1Se3OYxQgwIUQ7fc8shDCchFkIM3HHh/nMmTPY29sTFham+ffWW2/x1ltvNenxf/31F1u3bjVafcuXL2fEiBEMGTKEoUOHcvToUaP9rLomTpyo2Vq3tcTGxhIWFkafPn1wdXUlLCwMZ2dnJk6c2KDtvHnzOHfuXKPHCwgI4MyZM61aY126XhvNde7cOb7++msAwsLCyM/Pb61SgXb6AZip9e3bt9kv3BMnTrBr1y6effbZ1i2KmwvF7dq1i6ysLBQKBT/++CNTp07l+PHjrf6zTGHVqlUAfPLJJ+Tn57N8+XLS0tK0ftiZmJho4uq0a8lr41YHDhwgPz9f6zTP1nDH98zapKWlaXqLHj16MHLkSObPn88XX3yBv78/ISEhTJgwAZVKxdtvv82BAwc0o9Fak4ODA+fOnWPDhg1cuHABb29vcnJyOHnyJEOHDiUsLIxx48Zx9epV1Go1c+fOxc/PD29vb3bu3Anc7A39/f3x9/dn7dq1AERFRTFr1ixGjhzJQw89xLFjx4Cbq0n6+Pjw6KOPmnRyy6lTp4iIiGDQoEGaecu1PVdcXBzh4eEEBQXx888/s2jRIgYNGsQTTzxBSUmJyWqsS9fvtHZk1p49e4iKigLqv36WL1/O1q1b2bVrF3BzCeFhw4bh7+/P6dOnW1yX9Mzc7AFrl00FmDlzpubr8+fPc+zYMRwdHXnqqad45ZVXGD9+PBs3bqS0tJRFixaRlJRklMsR999/P7t27eK9994jPj4eW1tb3n77bVauXMmGDRvo27cv69evZ8WKFfj6+lJSUkJOTg5Xrlxh9erVWFpaUlhYSHZ2NtXV1YSEhDBs2DDg5ovsww8/5KOPPkKpVBIfH8/atWs5efIkFhYWDBo0qNWfjy4VFRV89dVX1NTU4Orq2mAhAk9PT9auXUtubi4ZGRl8//33XL9+HXd3d6PXdutrY9q0aTp/p9rUff14eXmRn5/PmDFjWL16NaNHj+a5554jLi6O7du38+qrr7aoVgkzDU+l6n7t5OSEo6MjAKtXr2bZsmX861//wtPTkyeeeMKodRUUFGBvb8+GDRsAyM3NJSIigoqKCl544QUAqqqqcHd3p1OnTgQGBgLQpUsXlixZwsqVKwkNDUWhUGBlZUVAQAB5eXkA+Pj4AODi4kJWVha//fYb/fr1o2PHjgD4+fkZ9bnV1b9/f83PvXXYLkDv3r0B+PXXX/H19cXCwgJ7e3uTjAy89bXR2O+0Vt2rvXVfP7eq/YN57733Ulxc3OJa5TRbj7rrNimVSuLi4khPT0etVvPll19iYWGh2XmhtZ04cYK5c+dSWVkJgIeHB507d6ZXr15s3LiRtLQ0VqxYQWRkJJ6ennz//fcAXL16lZEjR+Lp6UlmZiZwM/SHDh3S9Ga3Lvvq7u7OTz/9xI0bN6ipqeGHH34wynPSRt8StLX/B3379iUnJweVSkVZWVmDEJmCrt+pjY0NRUVFAJq3LXVrr/267multZfelZ7ZAH5+fkRGRtKpUyfs7OyIjIykoqKCkydPkpiYyLx581r15z355JP8/PPPDB48GDs7O1QqFStXrsTFxYUpU6ZQXV2NQqFg/fr1uLu7k5qaSkhICNXV1bz55ptERESQlpZGYGAglZWVTJgwgYEDB2r9Wc7Ozrz22msEBQXh7OzM3Xff3arPpTV4e3sTERHB4MGD6d69O127djV5DZGRkVp/pzNmzGD69Ols2bIFDw8PrY8dMGAAb7/9ts7/g5aSEWBCmAk5zRbCTEiYhTATEmYhzISEWQgzIWEWwkxImO9gaWlpdO3albCwMIYMGUJQUJBJNgUXxiHXme9ww4YNY9u2bQBcv36dIUOG4OHhgbe3d9sWJgwmPbPQsLOzY9asWWzfvp2FCxcSEhJCYGAgn3/+OQA5OTkMHjyYYcOGMXHiRKKiojhz5gwDBgwgLCyMFStWaJ0EAmg9nmhd0jOLerp168aqVavw8fEhMzOTiooKAgICeOSRR4iJiWHTpk3069ePRYsWceHCBQCKi4s5evQo1tbWBAQENJgEEhISQmFhYYPjtfWm6eZGwizqOXv2LJMmTWLTpk2a2UJVVVWcOXOGixcv0q9fPwBCQ0M1p+dubm5YW1sD8PPPPzeYBHLy5EmOHj3a4HhyKt+6JMxCo7S0lI8++ogZM2YwdOhQlEolKpWKJUuW0LNnT1xcXMjLy6Nv376aPZWh/mSC3r17s3HjRlxdXcnKyqKoqAhra2utxxOtS8J8hztw4ABhYWGaLYDi4+MZO3YssbGxhIaGcv36dcaOHUunTp14//33mT59OnZ2dlhbW3P//fc3ON4HH3ygdRJIWlpag+OJ1iUTLUSTrVu3jgkTJuDs7MzixYuxtrbmjTfeaOuyxP8nPbNosm7duhEeHo6dnR0ODg58+umnbV2SqEN6ZiHMhFxnFsJMSJiFMBMSZiHMhIRZCDMhYRbCTPw/X7ryXv7Wlv8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 240.945x216 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Set the palette using the name of a palette:\n",
    "sns.set_palette(\"Blues\")\n",
    "\n",
    "#Make sure there aren't duplicates because of demo\n",
    "\n",
    "rel_only_df = rel_limit_df[[\"ptid\",\"relation_mrn\",\"rel_degree\"]]\n",
    "print( len(rel_limit_df) )\n",
    "\n",
    "rel_only_df.drop_duplicates(inplace=True)\n",
    "print( len(rel_only_df) )\n",
    "\n",
    "#Count of individuals by degree, can be case OR control\n",
    "degree_counts_all = rel_limit_df.groupby(\"rel_degree\").agg({\"mrn\":'count'}).reset_index()\n",
    "\n",
    "degree_counts = degree_counts_all[ degree_counts_all[\"rel_degree\"].isin([\"First\",\"Second\",\"Third\",\"Fourth\"]) ]\n",
    "#degree_counts = degree_counts_all\n",
    "\n",
    "#Divide by 2 since reciprical of every relationship is included\n",
    "degree_counts[\"unique_rel_count\"] = degree_counts[\"mrn\"]/2\n",
    "degree_counts.sort_values(by=\"unique_rel_count\", ascending=False, inplace=True)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "fig.set_size_inches(figsize_1colw,3)\n",
    "sns.barplot(x=\"rel_degree\", y=\"unique_rel_count\", \\\n",
    "                 data=degree_counts)\n",
    "sns.despine()\n",
    "\n",
    "ax.set_xlabel(\"Degree\")\n",
    "ax.set_ylabel(\"Number of Relationships\")\n",
    "sns.despine()\n",
    "\n",
    "ax.bar_label(ax.containers[0])\n",
    "\n",
    "save_and_track_figure(fig,\"relationship_degrees.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 168,
   "id": "44a41e99",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:57.185527Z",
     "start_time": "2023-11-03T20:32:35.466354Z"
    }
   },
   "outputs": [],
   "source": [
    "#Map hhid and bid to patients\n",
    "rel_limit_df[\"hhid\"] = rel_limit_df[\"ptid\"].astype(int).map(mrn2hhid)\n",
    "rel_limit_df[\"bid\"] = rel_limit_df[\"ptid\"].astype(int).map(mrn2bid)\n",
    "\n",
    "rel_limit_df[\"rel_hhid\"] = rel_limit_df[\"relation_mrn\"].astype(int).map(mrn2hhid)\n",
    "rel_limit_df[\"rel_bid\"] = rel_limit_df[\"relation_mrn\"].astype(int).map(mrn2bid)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 169,
   "id": "c76caca2",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:57.301049Z",
     "start_time": "2023-11-03T20:32:57.190303Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5197226870666986"
      ]
     },
     "execution_count": 169,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Percent relatives that share hhid (will be the same as bid overlap)\n",
    "len( rel_limit_df[ rel_limit_df[\"hhid\"] == rel_limit_df[\"rel_hhid\"]][\"ptid\"].unique())\\\n",
    "/len(rel_limit_df[\"ptid\"].unique())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 170,
   "id": "e05cdc56",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:57.478212Z",
     "start_time": "2023-11-03T20:32:57.303863Z"
    }
   },
   "outputs": [],
   "source": [
    "share_hhid_df = rel_limit_df[ (rel_limit_df['hhid'] == rel_limit_df['rel_hhid'] ) &\\\n",
    "            (rel_limit_df['mrn'] != rel_limit_df['relation_mrn'] )\n",
    "            \n",
    "            ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 171,
   "id": "5323050b",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:57.590265Z",
     "start_time": "2023-11-03T20:32:57.482383Z"
    }
   },
   "outputs": [],
   "source": [
    "#percent that live together\n",
    "\n",
    "degree_hhid_list = list()\n",
    "for degree in [\"First\",\"Second\",\"Third\",\"Fourth\"]:\n",
    "    \n",
    "    shared_hhid_count = len( share_hhid_df[share_hhid_df[\"rel_degree\"] == degree][[\"mrn\",\"relation_mrn\"]].drop_duplicates() )\n",
    "    all_degree = len( rel_limit_df[rel_limit_df[\"rel_degree\"] == degree][[\"mrn\",\"relation_mrn\"]].drop_duplicates() )\n",
    "    \n",
    "    percent_live_together = shared_hhid_count/all_degree\n",
    "    #must divide by 2 since reciprical\n",
    "    total_live_together = shared_hhid_count/2\n",
    "    total_with_that_rel = all_degree/2\n",
    "    \n",
    "    degree_hhid_list.append([degree,total_live_together,total_with_that_rel,percent_live_together])\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 172,
   "id": "080bb9a6",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:57.642970Z",
     "start_time": "2023-11-03T20:32:57.593432Z"
    }
   },
   "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>Degree</th>\n",
       "      <th>N Lives Together</th>\n",
       "      <th>N</th>\n",
       "      <th>Percent lives together</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>First</td>\n",
       "      <td>3294.0</td>\n",
       "      <td>6124.0</td>\n",
       "      <td>0.537884</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Second</td>\n",
       "      <td>255.0</td>\n",
       "      <td>1171.0</td>\n",
       "      <td>0.217763</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Third</td>\n",
       "      <td>13.0</td>\n",
       "      <td>190.0</td>\n",
       "      <td>0.068421</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Fourth</td>\n",
       "      <td>9.0</td>\n",
       "      <td>133.0</td>\n",
       "      <td>0.067669</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   Degree  N Lives Together       N  Percent lives together\n",
       "0   First            3294.0  6124.0                0.537884\n",
       "1  Second             255.0  1171.0                0.217763\n",
       "2   Third              13.0   190.0                0.068421\n",
       "3  Fourth               9.0   133.0                0.067669"
      ]
     },
     "execution_count": 172,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pd.DataFrame(degree_hhid_list,columns=[\"Degree\",\"N Lives Together\",\"N\",\"Percent lives together\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 173,
   "id": "3ec6791e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:57.737185Z",
     "start_time": "2023-11-03T20:32:57.645720Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.5197226870666986"
      ]
     },
     "execution_count": 173,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#Percent relatives that share hhid\n",
    "len( rel_limit_df[ rel_limit_df[\"bid\"] == rel_limit_df[\"rel_bid\"]][\"ptid\"].unique())\\\n",
    "/len(rel_limit_df[\"ptid\"].unique())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab35f2f0",
   "metadata": {},
   "source": [
    "## Get age by week"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "id": "422ff8ab",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:57.925865Z",
     "start_time": "2023-11-03T20:32:57.740083Z"
    }
   },
   "outputs": [],
   "source": [
    "mrn2age = dict( zip( demo_df[\"ptid\"], demo_df[\"age\"] ) )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "id": "bf9b02e4",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:32:58.068143Z",
     "start_time": "2023-11-03T20:32:57.929564Z"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "regex = r'(\\d+_\\d+_\\d+)'\n",
    "\n",
    "pull_date = lambda trait_name: datetime.strptime(re.findall(regex, trait_name)[0]\\\n",
    "                                                 ,'%Y-%m-%d').date()\n",
    "\n",
    "pull_date_window = lambda trait_name: datetime.strptime(re.findall(regex, trait_name)[0]\\\n",
    "                                                 ,'%m_%d_%Y').date()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "id": "7c9b5281",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-03T20:34:08.964534Z",
     "start_time": "2023-11-03T20:32:58.071529Z"
    }
   },
   "outputs": [],
   "source": [
    "\n",
    "\n",
    "directory = cummulative_directory  # directory of gzip cummulative files\n",
    "\n",
    "# Get a list of all gzip files in the directory\n",
    "gzip_files = [file for file in os.listdir(directory) if file.endswith('.gz')]\n",
    "\n",
    "# # List to store the DataFrames\n",
    "\n",
    "cases = dict()\n",
    "controls = dict()\n",
    "\n",
    "\n",
    "percent_female_list = []\n",
    "# Iterate over each gzip file and read its contents into a DataFrame\n",
    "for file_name in gzip_files:\n",
    "#     print(file_name)\n",
    "    file_path = os.path.join(directory, file_name)\n",
    "    with gzip.open(file_path, 'rt') as file:\n",
    "        df = pd.read_csv(file, delimiter='\\t')  # Adjust the delimiter and other options as per your file format\n",
    "\n",
    "        #Map familiy IDS\n",
    "        df[\"famid\"] = df[\"ptid\"].map(mrn2famid)\n",
    "        df[\"sex\"] = df[\"ptid\"].map(mrn2sex)\n",
    "        \n",
    "        \n",
    "        #For every week, pull out cases and control IDs, and their ages\n",
    "        for pheno in df[\"pheno\"].unique():\n",
    "            \n",
    "            #Limit dataframe to those patients in solar run (ie for w/ proband has at least 1 case and 2 fam members)\n",
    "            df_solar = pull_fam_count(pheno, df, print_res=False, return_patient_limited_df=True)\n",
    "            \n",
    "            week = pull_date_window(pheno)\n",
    "            cases[week] = list(df_solar[ (df_solar[\"value\"] == 1) ][\"ptid\"].map(mrn2age))\n",
    "            controls[week] = list( df_solar[ (df_solar[\"value\"] == 0) ][\"ptid\"].map(mrn2age))\n",
    "            \n",
    "         \n",
    "            percent_female_list.append( (week, round( len( df_solar[ df_solar[\"sex\"] == \"F\"] ) / len(df_solar) * 100 , 2)))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 242,
   "id": "e1fa949e",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-06T23:58:21.487804Z",
     "start_time": "2023-11-06T23:58:18.783594Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAGoCAYAAACTwUbTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAwpklEQVR4nO3deXgUVd728buzYNjCGh0iizqooCLykFEQMoAsgiA6LMoOYgggAdlRYZBNFiUgYISExeCIMzIOAo8CKmICARETkWdQgsTBAIKYAIaIZO16/+BND1lJoKu7Ovl+rotLu/tUnV9Vd3LnnNNdbTMMwxAAALAEL3cXAAAA/otgBgDAQghmAAAshGAGAMBCCGYAACyEYAYAwEIIZtyQ06dPl3mb3Nxc/fzzzyZUg7LgeShaRkaGzp8/7+4yUIERzChk1qxZevLJJzVo0CANGTJEly9f1rp16/TTTz/la5ecnKw333wz331z585VTEyM1qxZU+z+4+Li9OGHH2r79u2Kj4+/7jqffvppRUZGOm5///33uuOOO5SZmXnd+7xw4YI6dOig2NjYa7bdsWOHBg8erCFDhujll1++7j6v16BBg/Ldfu211zR8+HD17dtXH3zwwTW3z3seivPjjz9qxowZpe6/KP369ZMk7d+/X82bN5ckffbZZyW+PiQV+xo6e/bsNbe9Ue+9956+/vprzZo1S0lJSab2VVwfc+fOvea21zr/JT2ed343b96sf//739cuFC7l4+4CYE2LFy9W48aN9dVXX2n16tU6f/68Ll++rDFjxsjLy0uNGjVSvXr1FBsbq82bN+uNN97Q008/raNHjyo4OFjvvvuuEhIS1K5dO2VkZKh+/frq1KmThg0bpjvvvFNffPGF+vTpI19fX61cuVKHDh1SRkaGIiMj9dhjj6lDhw46evSooqOj5e3tLUmaPXu2UlNTZRiGBg8erO+++05hYWGOmteuXatRo0bpvffe05AhQxQZGalvvvlGx44d04wZM5ScnKyEhASlp6fr9ddf10svveQ4lsmTJ0u68ks5JydHt99+uwYOHCg/Pz91795dFy9e1JYtWzR48GD16tVLOTk5WrlypTZv3iybzaZ3331Xly9f1qRJk2S321WtWjW9+OKLGjt2rHx9fTVs2DBVrlxZf/vb3/T7779rzJgx2rt3r5KSkpSRkaHVq1dLkgzD0JgxYxz7CAsL0/Dhw9W+fXtlZWVpxowZGjFihGrVqqVff/3VcewHDx7UhQsXtG7dOhmGoejoaJ0/f15jx47Ndwx79+5V5cqV1a5dOyUmJuqLL76Qj4+P49g+//xzSVLt2rX1zDPPOM7J559/rosXLyo6OlqVKlWSJP3nP//R+PHjlZmZqZCQEO3cuVPTpk1TaGiooqKiJEl33XWXTpw4oU8++URdu3ZVUlKSdu/erWeffVbTp0/XxYsXVbVqVc2ZM0dTpkxRbm6ubrvtNgUFBUmSJk2apL59+6pVq1aSpOXLl2vatGnav39/vnO5YsUKrV+/Xi+//LIGDhyoiIiIfPuaN2+eOnfurNOnT2vZsmUKDQ1VZGSkbDabkpKS9PLLL8vb21sjR45UXFycbrrpJtWtW1dz587VhQsX9NJLL0lSvj5Xr16tnJwcTZo0Sffcc4/CwsK0cOFCjRgxQhMmTNDOnTs1cOBALVq0SJL04IMPqm3btgoPD1dubq6efPJJSdJvv/2mQYMG6fXXX1fdunUlyfHav/r5euCBB/TKK6/ot99+0+zZsyVdCdikpCSFhIRo2LBhWrx4scaMGaN69epJkpKSkvL15+/vr1WrVikjI0Ndu3bVE088obFjx+qNN95wyu8NOAcjZpSoWbNmOnr0qON2cnKymjdvrp49e6pNmzZq166datasqR49emjEiBGOdn379nUEV0Ft2rRx/FKSrozcVq1apccee0zbtm1TjRo1NGPGDDVs2NAxSk9PT9eZM2e0YsUKNWzYUJUqVVLz5s0VHBwsScrJydG3336rsWPH6p133pF0ZZS2cuVKPfroo5KkqKgoValSRd7e3jp48GC+Y8nTtWtXtWvXTv/7v/+rsLAwrV27Vhs3bpQkDR06VL169ZIkpaSk6I9//KNsNpskacCAAfLx8VH//v3VpUsXJSQk6OLFi7p48aJ69Oih++67TxEREapevbpq166tPXv2KDk5WYGBgQoJCXH0n5OTk28fkvTAAw9o5syZSk5O1meffaYePXpo2bJl8vPzc2yXmJjoCC+bzaZnnnlGf//73wsdw+OPP67w8HB99NFH+Z6HoUOH6u6771ZgYKBWrFih1NRUx3TuiRMnVLt2bY0cOdLxR5IkNWjQQK+//rruuOMO2Ww2HT58WLGxsY7nRJI6deqkuLg4paSkaMCAAfr444/1yy+/KD09XV9++aWqVq2qH374Qdu3b9epU6dUrVo1HThwQNKVEG7QoIHjuKQro3h/f/9C57JRo0Y6ffq0Tp06pWPHjhXaV48ePTRt2jRdunRJR48eVaNGjRzPXWRkpBYvXqy1a9dq1apV+c7L888/ryVLlujTTz8t1GdmZqbCw8N1zz33SJI6duyouLg4eXt7Ky4uTmfOnNGHH36ojIwM+fv764svvlBkZKT8/PxUu3Zt7d69W5I0ZMgQTZ482RHKV7v6+YqKitKCBQsUFRWlatWqFWorXfkjauLEiXrllVfk5eVVqL9169YpOjpazz77rOO1wrS99RDMKFF8fLzuvPNOSVfWJGfPnq3AwEBNnDjR8YtNUqFfFHmh4efnJx8fH+Xk5DhGeFdvJ10ZJeax2+2qUqWKJMnX11d2u93R5urt8u7P89FHH+ncuXOaMGGCkpOTdfDgQWVlZeXrz9/fXwsXLlRISIgaNGiQ71gKunr/ef9/9THefPPNOn78uKP2V199Vd98842ioqLUqFEj1apVS9WrV9f8+fOVlpam+fPnKzs7Wy+99JKmTp2q5s2ba/DgwerYsaNmz56ts2fPSpJOnTqVbx+SHOfD29tbXl5ejmO6OiTvueceffnll45zNXXq1CKPoUqVKo7zWvD5K3hO82536dJF/fv316pVq3To0CHH43nPsY+Pj3x8fNS6dWstXLhQvXv3drRp3bq1YmJiVLduXTVv3lyJiYmqWbOmcnNz1bJlSy1cuFB9+/aVzWZTp06dtHDhQkco9unTRzExMY7nMe8cSCp0Lp966iktWrRI999/v3JzcwvtK++569Wrl0JCQjR48OBCx2mz2QqdF39/f/n4+Cg3N7dQnwVfEx07dlR0dLQ6dOigI0eOOI6zb9++WrBggdq3b6/c3FyFhoZq1qxZeuihhyRJU6ZMybckc7Wrn6/s7GxHbWfOnHGc+5ycHBmGobS0NMfrw8vLS15eXoX6y3v86teOr69vkX3DfZjKRpEmT56sqlWrSpJWrVql8PBweXt7680331SVKlX0P//zP6pVq5b27NnjGDFcbdOmTTp+/Lgef/xxtWjRQmFhYfrss89UtWpVBQYGas6cOXrqqackSW3bttW4ceOUk5Oj8PDwItdH/f39dcstt2jSpEmy2Wxq0aJFvsf//ve/a+vWrQoICNDXX3+tlStX6qGHHtL48eOVkJCgtm3bqnv37hoxYoQuXbqkNWvWaP78+Y5jKah///6aOHGiNmzYoP79+ystLS3f497e3ho2bJh69+6tKlWqqHHjxmrYsKF+/fVXvf322zpz5ozsdrtmz56tOnXqqHPnzgoMDNSIESNks9n017/+VTt27NCxY8cUEBCQL4QL7uNqjzzyiEaNGqV9+/YpOzvbcX/z5s21ZcsWDRgwQFlZWRo4cKCCg4NLPIaCz0OzZs0UFRWlyZMnq169eo4R3LfffqtPPvlE3t7eatCggWP748eP65VXXtGvv/6q+++/X3Xq1NHBgwcdf0hIV37p//LLLxo6dKgkKTMzU127dtX999+vN954Q+PGjZOfn5/mzp2rYcOG6fDhw7r99tsVGBiowMBAPfTQQwoPD9eLL77oqPnXX3/V2LFj853Le++9V3v27NG0adNUp06dQvvK06VLF61atUoNGzZ03Ddy5EhNmzZNfn5+CgsLU2ZmplatWqUmTZrkO18F+yzI399fZ86cUXBwsGJjY9WxY0fde++9Gj16tN5//3116tRJo0aN0ksvvSQ/Pz/HqLV169Y6cuSI9u7dqzZt2hTab56QkBBNnTpVWVlZmjlzpuM5mz9/vhITE1W1alU99dRTGjdunAICAuTt7V2ov2effVbDhw9XTk6OHnnkEdntdtWoUaPYPuEeNr7EAuXVa6+9prNnzyojI0Ph4eG66aab3F1SuZWWlqYhQ4ZowYIFRf6h5iznz5/X+++/r9DQ0Ovavl+/fgoNDdUjjzzi5Mo80wcffKD77rvPMSsGayCYAQCwENaYAQCwEIIZAAALIZgBALAQjwjmrl27ursEAABcwiOCOTU11d0lAADgEh4RzAAAVBRcYAQAYEnZ2dk6deqUMjIy3F3KDfPz81P9+vVLdaU1ghkAYEmnTp1S9erVddtttxW6lK8nMQxD586d06lTp3T77bdfsz1T2QAAS8rIyFCdOnU8OpSlK9dhr1OnTqlH/oyYAQCWlRfKHyWdVWau/RqtC7vJ20vdG9/i7LLKrCx/XBDMAADLu55QLmm7jIwMTZo0SZmZmTp37pwWLlyou++++0ZKdBqCGQBQ4axcuVJPPPGEunTpopSUFO3fv1/h4eE6c+aMFi1apE8//VRJSUnKyMjQ6tWrNX36dF28eFFVq1bVlClTNHbsWPn6+mrYsGHq0KGDU2tjjRkAUOEkJiaqVatWkqSAgADdeeedGjRokO6880599dVXSk5OVmBgoEJCQvTtt9/qyy+/VNWqVfXDDz/o1KlTunjxonr06KH77rvP6bURzACACueee+7R/v37JUkpKSlq0aKFzp07p/vvv1+GYWjw4MHq2LGjZs+erdTUVLVs2VILFy5U3759deutt2r+/PlKS0vT/PnznV4bU9kAAMu7ydvrut/8VZTQ0FCNHj1a77zzjn777Te98MIL+uyzz3ThwgX96U9/0q5du3Ts2DEFBASodevW2rBhg8aNGyc/Pz+1b99es2fPVp06ddS5c+cbPbRCPOL7mIOCghQfH+/uMgAALnTkyBE1bdrU3WU4TWmPh6lsAAAshGAGAMBCCGYAAK4hLS1Nly5dcklfvPkLAGB5UXui9HvW72XerkqlKgoNDi10f1kvMPL6669r4MCBaty4cZGPR0dHq379+urUqVOZayyIYPYgV1/SraT37JW2nRn7pG/6dkffZuyTvq3V9/WEcknbrVy5Uk2bNlWrVq104cIFbdy4Ud9//73sdrtmzpypJUuWKDAwUMeOHdOsWbO0bds2nT9/XidPnlStWrX07LPPKiIiQn5+furevbtjvyEhIfrtt9900003af369ddVM1PZHiTvRXutX4ClbWfGPumbvt3Rt6uPx919u+t5dGffzpaYmKghQ4ZIkjp37qwffvhB69at09KlS/Xmm28qMzNTEyZM0NChQ7Vv3z5169ZN7du3l7+/v8LDw3Xw4EGFhYVp7dq12rhxo2O/WVlZql+/viZNmnTdtRHMAIAKp+AFRt577z3HY3b7lc9LV6lSRb6+vrLb7flG9tWqVXO0ubq9JD3//PO68847NX78eOXm5l5XbUxlAwAqnLwLjKSkpOimm27Sjh07FBoaKi8vL7344ouaN29evvZ33HGHNm3aJG9vb0lS//79NXHiRG3YsEH9+/dXWlqaJGnjxo36z3/+o6ZNmzralhXBDACwvCqVqlz3m7+KUrlyZUVHRys+Pl5BQUGSpHbt2jkej46OliS1b9/ecV/Tpk0dbevWrau333670H47deqUb5/Xg2AGAFhe3jurSxt6NxqO7sQaMwDAsi5fvuzyN4Y5m2EYunz5cqnbM2IGAFhSvXr19NNPPyk7O9txX2pqqo4cOXLNbZ3d7kb36evrq3r16pWqH77EwsPYbLZS/fVY2nZm7JO+6Zu+PW+f9O36vovDVDYAABbi1GD+8ssvHe9gS0pKUtu2bRUcHKzRo0c7Puc1e/ZsPfjgg3r44Yd14MABZ3YPAIDHc1owv/rqqwoJCVFGRoYkaeLEiZo3b5727NkjwzC0ZcsWff3114qNjdWXX36pf/zjHxozZoyzugcAoFxwWjD/8Y9/1KZNmxy3ExISHJ8J69atm3bu3Km4uDh16dJFNptNDRs2VE5OjlJSUpxVAgAAHs9pwdy7d2/5+vo6bhuG4biEWfXq1ZWWlqaLFy+qRo0ajjZ59xclKipKQUFBCgoKIrwBABWGaW/+8vL6767T09NVs2ZN+fv7Kz09vdD9RQkNDVV8fLzi4+MVEBBgVpkAAFiKacHcokULxcTESJK2b9+u4OBgtWnTRh9//LHsdrtOnDghu92uunXrmlUCAAAex7QLjISHh2vEiBHKyspS06ZN1adPH3l7eys4OFitW7eW3W5XRESEWd1Lss53xTqrbwBA+VfuLzBSUT+4Tt/0Td/lt28z9knfru+7OFxgBAAACyGYAQCwEI/8EgvWZAEA5ZVHjpjzwphQBgCUNx4ZzAAAlFceOZUNwDnc+fE/lqSAojFiBiqw0i4LlWX5yIx9AhUJwQwAgIVUqKnsq6fOJNde4Yu+Xd+3Gfssb317gvJ2zumbvq+lQo2Yr546c9Y02/W0c0ffZh6PVfs2Y5/lrW9PUN7OOX3T97VUqGAGAMDqCGYAACyEYAYAwEIIZgAALIRgBgDAQghmAAAshGAGAMBCCGYAACyEYAYAwEIIZgAALIRgBgDAQghmAAAspFx8u1R5+iYdAEDFVi5GzOXpm3QAABVbuQhmAADKC4IZAAALIZgBALAQghkAAAshmAEAsBCCGQAACyGYAQCwEIIZAAALIZgBALAQghkAAAshmAEAsBCCGQAACyGYAQCwEIIZAAALIZgBALAQghkAAAshmAEAsBCCGQAACyGYAQCwEIIZAAALIZgBALAQghkAAAshmAEAsBCCGQAACyGYAQCwEIIZAAALIZgBALAQghkAAAshmAEAsBCCGQAACyGYAQCwEIIZAAALIZgBALAQghkAAAshmAEAsBAfM3eenZ2toUOH6scff5S3t7dWr14tHx8fDRs2TDabTffdd58iIiLk5cXfBwAASCaPmLdt26acnBzt27dPM2fO1PTp0zVx4kTNmzdPe/bskWEY2rJli5klAADgUUwN5rvuuks5OTmy2+26ePGifH19lZCQoHbt2kmSunXrpp07d5pZAgAAHsXUqexq1arpxx9/VJMmTZSamqoPP/xQu3fvls1mkyRVr15daWlpRW4bFRWlqKgoSVJKSoqZZQIAYBmmjpiXLl2qRx99VN9//70OHTqkoUOHKisry/F4enq6atasWeS2oaGhio+PV3x8vAICAswsEwAAyzA1mGvVqqUaNWpIkmrXrq3s7Gy1aNFCMTExkqTt27crODjYzBIAAPAopk5lT5gwQcOHD1dwcLCysrI0f/58BQUFacSIEcrKylLTpk3Vp08fM0sAAMCjmL7GvHHjxkL3x8bGmtktAAAeiw8QAwBgIQQzAAAWQjADAGAhBDMAABZCMAMAYCEEMwAAFkIwAwBgIQQzAAAWQjADAGAhBDMAABZCMAMAYCEEMwAAFkIwAwBgIQQzAAAWQjADAGAhBDMAABZCMAMAYCEEMwAAFkIwAwBgIQQzAAAWQjADAGAh5TaYGzZqJJvNJkmy2Wxq2KiRmysCAODafNxdgFlOnjihfyWedtzu3STQjdUAAFA65XbEDACAJyKYVfppb6bHAQBmK7dT2WVR2mlvpscBAGZjxGwCRtYAgOvFiNkEjKwBANeLEbObOXt9uyyjddbWAcB6GDG7mbPXt8syWnf2Phs2aqSTJ05IuhLgDRo21Ink5GL7BwAURjDDaZjCB4Abx1Q23KI8TeG7c6q/YN8sNQCejxEz3KI8TeG7c6agYN+u7h+A8zFiBlyIN9IBuBaCGXChvBFu3r+8N8sVZEaAmz2F78plgYretzuWYsrSN24MU9mABZkxPW72FL4Z+6Rvz+wbN4YRMwAAFkIwAwBgIQQzAMCpytPaujuwxgwAcKrytLZ+9RUNJddc1ZBgBgCgGO540xtT2QAAOIGzpscZMQMA4ATO+pgjI2YAACyEYAYAwEIIZgAALIRgBgDAQghmAAAsxOOC2cpXawEA4EZ53Mel3Pml9AAAmM3jRswAAJRnBDMAABZCMAMAYCEEMwAAFkIwAwBgIQQzAAAWQjADAGAhBDMAABZiejAvWLBArVu3VsuWLbV27VolJSWpbdu2Cg4O1ujRo2W3280uAQAAj2FqMMfExGjfvn3au3evYmNjdfLkSU2cOFHz5s3Tnj17ZBiGtmzZYmYJAAB4FFOD+eOPP1azZs30l7/8RY8//rh69OihhIQEtWvXTpLUrVs37dy508wSAADwKKZeKzs1NVXJycn68MMPdfz4cfXs2VN2u93xJRTVq1dXWlpakdtGRUUpKipKkpSSklJsHz6+Pvm+1KJ+g/o6eeKkk48EAADXMDWY69SpoyZNmqhSpUq6++675efnp5Mn/xua6enpqlmzZpHbhoaGKjQ0VJIUFBRUbB852TlaunOp4/aEThOcUzwAAG5g6lR227ZttWPHDhmGodOnT+vSpUvq2LGjYmJiJEnbt29XcHCwmSUAAOBRTB0x9+jRQ7t379aDDz4ou92uiIgI3X777RoxYoSysrLUtGlT9enTx8wSAADwKKZ/H/Orr75a6L7Y2FizuwUAwCNxgREAACyEYAYAwEIIZgAALIRgBgDAQghmAAAshGAGAMBCKkwwF7x0p81mU4OGDdxcFQAA+Zn+OWarKHjpTqn4y3eW9vrbXKcbAOBsFWbEXBZ5IZ7379TJUzfUjtE6ilPwtVHc66K07czYZ1n6did3Hjd90/eN7jPfdqVqhRtSltE6KpbSfglLWb6sxdn79JQvinHncdM3fd/oPq/GiBkAAAshmC2mPE0v5tXojukmT+m7vKmo0+iAMzGVbTHlaXqxLFP4Zk9LWbXv8qaiTqMDzsSIGQAACylxxHzkyBE1aNBAb731llq2bKmHH37YVXXdIJsMQ5ICHfcYhvT6Z24rCACAUilxxDx//nx9/vnnunjxotasWeOqmuBEZqx3su4HAOYpccRst9v17rvvasaMGZo1a5aLSnIGQzabTf9KPO24p3eTQC3d6caS3MSM9U7W/QDAPCUG8wsvvKDffvtN/v7+mjNnjqtqAgCgwipxKnv58uW69957deutt+ruu+92VU1wE6aoAcD9Shwxe3t7a/z48Y7b69atM7seuBFT1ADgfiUG8/PPP6+vvvpKhmEoLS3NVTUBAFBhlTiV/cILL+j//u//JElfffWVSwoCAKAiKzGY//CHP+iJJ55QjRo1VKtWLVfVBABAhVViMPft21e33HKLdu3apXr16rmqJgAAKqwS15izs7MVHR2tRx55pByPmEt7lTCuJgYAMF+JI+bNmzdLkrp166a1a9e6oh4AACq0YoN527Zt8vb2VmZmpn788UddvnzZlXW5kCGbTdp09LTj3///KO91tgMA4PoVG8yHDx/WyZMn5ePjo9WrV2vBggWurAsAgAqp2GCeOnWqFi1apNjYWJ09e5ZLcgIA4ALFvvkrISFBc+fO1cKFC9W+fXsXlgQAQMVVbDCvX79eGzZsUNWqVV1ZzzV4wjujPaFGAIBVFRvMy5cvd2UdAABA1/i4lPV4wjujPaFGAIBVeVgwAwBQvpV45S+YjfVoAEB+BLNHIMABoKJgKtutWI8GAOTHiNkjGLLZbPpX4mnHPb2bBGrpzoLtyjKydvaXdzCqd73C51xy1fPoKX0DnocRMwAAFsKIuVwp7ci6LG2d3Q7OU/icS656Hj2lb3fO5DCLhOtDMMMN3DmNXp76hvO4cxq9vC0f8Dq/UUxlAyjH3PkGS97cievDiBlu4M5p9PLUN5zHjGl0d/ZdnpYuzBitWxsjZgAALIQRMwDAwswYrbtzbf3aGDEDAGAhjJgBABWMO9fWr40RMwAAFkIwAwBgIQQzAAAWQjADAGAhBDMAABZCMAMAYCEEMwAAFkIwAwBgIQQzAAAWQjADAGAhBDMAABZCMAMAYCEuCeZffvlFDRo0UGJiopKSktS2bVsFBwdr9OjRstvtrigBAACPYHowZ2dna+TIkapcubIkaeLEiZo3b5727NkjwzC0ZcsWs0sAAMBjmB7MkydP1qhRoxQYeOXLoxMSEtSuXTtJUrdu3bRzZ9HfiRUVFaWgoCAFBQUpJSXF7DIBALAEU4M5OjpaAQEBevTRRx33GcaV76yUpOrVqystLa3IbUNDQxUfH6/4+HgFBASYWSYAAJbhY+bO161bJ5vNpp07d+qbb77RkCFD9MsvvzgeT09PV82aNc0sAQAAj2LqiHn37t2KjY1VTEyMHnjgAb399tvq1q2bYmJiJEnbt29XcHCwmSUAAOBRTB0xFyU8PFwjRoxQVlaWmjZtqj59+ri6BAAALMtlwZw3Spak2NhYV3ULAIBH4QIjAABYCMEMAICFEMwAAFgIwQwAgIUQzAAAWAjBDACAhRDMAABYiMsvMOJMvZoEypCkThMc9413VzEAADgBI2YAACzEo0fMmxJPq3eTQC3dudRx34ROE66MogEA8ECMmAEAsBCCGQAACyGYAQCwEIIZAAAL8eg3f5WkQcOG6t0k0N1lAABQJuV2xHwiOVmGceX92Xn/BQDA6srtiBkACirtRYnKcvEiZ++TvitW30XxuGA2Y4q6tPt0djsAAAryuGA+kZwsSbLZbDIMQzabzWX7dHY7AK5V2osSleXiRc7eJ31XrL6LUm7XmAEA8EQeN2IG4Nncue4HeAKC2c3cub7N2joAWA/B7GbuXN9mbR3u4M51P8ATEMwAisV0MuB6BDPcoqJO4bN8AOBaCGa4RUWdwje777z/dxamkwHXI5iB8sRmM/2qRADMRTADFcTV094ljWKZHgfci2AGypP/P+1d1HTy1dPepW3Hu+sB1+PKXwAAWEjFGDEXse4m3eCaWinX8oAbxmsNqFAqRjADcDpP+ZgYa+bwNBUjmItYT5Nu8KMcJazRAU5l0ddaWT6q5c6PibFmDk9TMYK5LJw9bWjGNDpcr7Svi7I830xRAygCwQwAKn9XeqNvz0UwF+TsaUMzptHheqV9XZTl+bboFHVF5c5p9LK2M2OfFa1vKyOYreQGpkuLbFfeMJ0MwEmsPLImmAEAFU5pR9bumEYnmK3kBqZLK8QUKNPJAFzMHdPoXPkLAAALYcRcnpmx1so6LwCYihEzAAAWwoi5PDNjrZV1XgAwFSNmAAAshGAGAMBCCGYAACyEYAYAwEIIZgAALIRgBgDAQghmAAAshGAGAMBCCGYAACyEYAYAwEIIZgAALIRgBgDAQghmAAAshGAGAMBCCGYAACzE1GDOzs7W4MGDFRwcrAcffFBbt25VUlKS2rZtq+DgYI0ePVp2u93MEgAA8CimBvM777yjOnXqaM+ePdqxY4fCwsI0ceJEzZs3T3v27JFhGNqyZYuZJQAA4FFMDea+fftq7ty5kiTDMOTj46OEhAS1a9dOktStWzft3LnTzBIAAPAopgZztWrVVL16daWnp6tPnz6aN2+eDMOQzWaTJFWvXl1paWlFbhsVFaWgoCAFBQUpJSXFzDIBALAM09/8dfLkSXXo0EGDBw/WgAED5OX13y7T09NVs2bNIrcLDQ1VfHy84uPjFRAQYHaZAABYgqnBfPbsWXXp0kWLFi3S8OHDJUktWrRQTEyMJGn79u0KDg42swQAADyKj5k7nz9/vi5cuKC5c+c61pqXLVumcePGKSsrS02bNlWfPn3MLAEAAI9iajAvW7ZMy5YtK3R/bGysmd0CAOCxuMAIAAAWQjADAGAhBDMAABZCMAMAYCEEMwAAFkIwAwBgIQQzAAAWQjADAGAhBDMAABZCMAMAYCEEMwAAFkIwAwBgIQQzAAAWQjADAGAhBDMAABZCMAMAYCEEMwAAFkIwAwBgIQQzAAAWQjADAGAhBDMAABZCMAMAYCEEMwAAFkIwAwBgIQQzAAAWQjADAGAhBDMAABZCMAMAYCE+7i4AAOCZ6jeorwmdJuS7fSPtzGBGjWYfN8EMALguJ0+clCTZbDYZhnHD7cqitKFnRo1mHzfBDACwhLKMMM0Ie6sgmAFUKJ4wtekpfTtbeQ7bsiCYAZQL5Wlq01P6hjkIZgDlAoGC8sLjg9md0y4Aro8nTKsC7uLxwcxfyYDnYVoVKB4XGAEAwEIIZgAALMTjp7IBWAdrwsCNI5gBOA1rwsCNYyobAAALYcQMoERMTwOuRTADKBHT04BrMZUNAICFEMwAAFgIwQwAgIUQzAAAWAhv/gIqKN5tDVgTwQxUULzbGrAmprIBALAQghkAAAthKhsoZ1g7BjwbwQyUM6wdA56NqWwAACzEI0fMNpvN8d/SjggKTu/l3Qd4AqangYrDI4P5eqbnmN6DJ+P1C1QcTGUDAGAhBDMAAFZiuEFubq4xcuRIo1WrVka7du2MY8eOldi+ZcuWJT4uyfGvpPuvvl1S2/oN6he67cx2ZuyTvst330XdLq5dWdo6u527+jbreOjbOn2bfTyu7rskbllj3rx5szIyMvTFF19o//79mjRpkrZs2XLd+zOKWXMreH9x7a71mJnt6Ju+S9uutK/zsrR1djv6pm+z+jZjn+7suyRumcqOi4tT165dJUmtWrVSfHy8O8oAAMBy3BLMFy9eVI0aNRy3vb29lZOTk69NVFSUgoKCFBQUpJSUFFeXCACAW7glmP39/ZWenu64bbfb5eOTf1Y9NDRU8fHxio+PV0BAgKtLBADALdwSzG3atNG2bdskSfv371ezZs3cUQYAAJbjljd//eUvf9Gnn36qhx9+WIZh6K233nJHGQAAWI7NuJ63jLlYUFAQbxADAFQIXGAEAAALIZgBALAQghkAAAshmAEAsBCCGQAACyGYAQCwEIIZAAALIZgBALAQghkAAAvxiCt/1a1bV7fddluRj6WkpFj+Sy48oUbJM+qkRufxhDo9oUbJM+qkRudxRp1169bVjh07inzMI4K5JJ5wuU5PqFHyjDqp0Xk8oU5PqFHyjDqp0XnMrpOpbAAALIRgBgDAQjw+mENDQ91dwjV5Qo2SZ9RJjc7jCXV6Qo2SZ9RJjc5jdp0ev8YMAEB54vEjZgAAyhOCGQAAC/FxV8efffaZ/vGPf+j333/XXXfdpZ9++kl2u12rVq3SihUrlJycrJSUFEVGRuq1117Tr7/+qptvvlnTp0937GPTpk3asWOHY7vo6GgdOnRIdrtdERERjnb//ve/FR4eLkkaMGCAunTporNnz6p79+7XfMu7GXVGRERo7969qlKlihYsWKB69eo52u7atUvvvvuu1qxZo6VLl+r7779XjRo1NG/ePPn4FP10mVFjWFiYMjIytG/fPv3rX/9Ss2bNJEn79u1TZGSkJGnZsmVavHixTp06pYMHD2rBggV67LHHXFbjjBkzdOnSJQUGBurFF190tCtYY1xcnHbu3Kmff/5Z4eHhuvXWW019viVpw4YNOnr0qObMmaOTJ09q0KBBio2NzdcmLCxMNptNzZo10/DhwzVq1Cj5+fmpbdu26tevn0trXLp0ab7tqlWrVqivqVOnateuXU5/TVauXFljx45Vr1691KlTp2Kf74LbFVfjTz/9pA8++EC5ubmKiIhQ5cqVnXIuS1unj4+Pdu3apfXr12v9+vVF9jV16lT985//dPrPTmlr3LRpk/bt26e0tDS9+eabjnNUsMbDhw9r3759On/+vBYvXuy0n53S1vn222/ru+++06lTp/TWW2/ley6v/l05Y8YMl51LKf/PTkREhI4cOaJKlSppyZIl+dpdXeOGDRtKfS5luMnWrVsNu91ufP3114aXl5dhGIaxfv16Y8eOHcaWLVsMwzCMJUuWGDExMcYzzzxjGIZhjBw50khPT3fs4+mnn8633YcffmgYhmE899xzxvnz5x3tYmNjjQsXLhjnz583nnvuOcNutxsTJ040evXq5ZY6Bw8ebIwePdqYMmWKkZub62h37NgxY/HixcbAgQMNwzCMvn37GoZhGO+//77j2FxVY149s2bNytfX0KFDjcuXLxuxsbHGqlWrDMMwjHPnzhnjxo1z+Xns06ePYRiGsXjxYuPw4cPF1jht2jSje/fuRufOnY2MjAzT69y3b58xa9YsY/r06YZhGMbkyZONnj175usnMTHRmDlzpmEYhjFgwADjwoULRlxcnJGdnW0MHjzY5TVevd2BAweK7OvVV191+mvywIEDxrp164xRo0YZn376ab59FHy+S1vj9u3bjZycHGPRokXGwYMHnXIuy1JnamqqMX36dMfPcXF1GoZzf3bKUuOAAQOM8ePHG3PmzCmxxrz9b9q0ydi4caPL63z66aeNp556yujdu3e+dgV/V7ryXBb82cmrefbs2UZiYmKxNZblXLptKvvxxx/X77//ruXLl6tLly6SpHr16unMmTPq2bOnkpKS9N1336lx48a6+eabJUl16tRRamqqYx++vr75tuvevbvWrFmjH374QVWqVHG0+/Of/6xq1arpr3/9q8LCwrRq1SoNGjSoxL+kzazzueeeU0REhBo3bpzvyi+NGzfWpEmTHLcHDx6skJAQ7dmzJ9/+XFGjJC1dulTjxo3L11dubq78/PzytXvjjTc0evRol5/Hrl27avTo0Tp8+HC+dgVr7Nmzp7Zu3apevXpp//79ptfZunVrDR061HH7tddeU/Xq1fP1k5qa6pgpsdlsqlGjhh566CH16tVLf/rTn1xe49XbtWzZssi+hg4d6vTXZMuWLfXMM8/ooYceKrSPgs93aWvs2rWrtm3bpi1btqhRo0ZOOZdlqXPBggWaNm1aiX3lnXtn/uyUpcYTJ05oyZIlyszM1HfffVdsjT179lRKSoo++OAD9ejRw+V1hoWF6b333lPjxo2VnJzsaFfwd6Urz2XBn52ifncWVWNZzqXbgjk1NVVjx47VnDlzdMstt0iSzpw5oz/84Q86cOCAli9frjfeeEMBAQE6f/68pCuXQUtLS1O/fv0UFRUlb2/vfNvFxcUpJCREI0aM0ObNm9WvXz+NHz9ely9f1pgxYzRy5Eg1bdpUn376qSIjI3XgwAH985//dHmdCQkJstlsqlWrllJTU9WvXz+98sorhfrOyMjQmjVrdN9996lBgwYurfHy5csyDEO1atXS6dOnHTVWqVJFmZmZjnaSdPz4cTVp0sTl59Hf318rV67UrbfeKl9f32JrXLJkiby8vBQQEKBz586ZXmdJli9frn79+qlu3br6+eefJUl2u13nz59XUlKStm7dqi+++MLlNV693c8//+w4l1f3dfPNNzv9Nenllf9X0KFDh4p9vktb4969e/X4449r+fLleuutt5xyLktb5++//65jx45pypQpOnDggGJjY4utU3Luz05ZzmWjRo1ks9lUu3btfL+DCtaYlJSkl156ScuXL7/mQMaMOpcuXSpJCggI0E8//VTs70pXnsuC8tqeOXNGlSpVKrbGspxLt31casiQIUpJSVGdOnXUuXNnxcXFyTAMRUZGqkmTJnrwwQdls9n04osvauvWrTp9+rRuvfVWvfDCC459fPTRR9q8ebNju6VLl+r48eO6dOmSwsPDVadOHUnSzJkzFRMTo4YNG6pNmzaOv6oGDRqkd955x+V1rly5Ut99950Mw9Dy5csdf3Hlyatrw4YN2rVrl2w2m6Kiogq9WMys8ZtvvlFMTEyhv0oTEhK0YsUKZWdnKzIyUpmZmZo1a5ZWrFjh8vO4bNkyHT16VP7+/nrttdeKrfFvf/ubDhw4IC8vL0VERMjPz8/UOiXpxx9/1Jo1azRv3rx8z+nVJk2apMuXLysoKEhPP/20QkNDVbduXQUEBGjGjBkuq3Hu3Lm666678m137733FurrySefVGZmptNfk/fee6+io6NVv379fGt5BZ/vorYrqcacnBz99a9/VePGjZ1yLktbZ164FHzOC9bZoUMHp//slLbGt956SwkJCapUqZJef/112Wy2Imtcvny5br75Zvn5+SkkJETt27d3aZ2zZs1SamqqfHx8ijxPeef43LlzLjuXUv6f77Vr1+rrr7+Wn5+f471MRdX45z//udTnks8xAwBgIXxcCgAACyGYAQCwEIIZAAALIZgBALAQghmoIPKuJrZ//341b95c0pUrIa1Zs6bYbWJiYkp8HIDzue2SnABc66677tKJEyf0ySefqGvXrkpKStLu3bt14sQJHTp0SFWrVtWcOXM0ZcoU5ebm6rbbblNQUJCkKx/t6tu3r1q1auXmowDKP0bMQAXRqVMnxcXFKSUlRQMGDNDHH3+s3bt36+TJk6patap++OEHbd++XadOnVK1atV04MABSVcuitKgQQNCGXARghmoIFq3bq2YmBjVrVtXzZs3V2JiomrWrKmWLVtq4cKF6tu3r2w2mzp16qSFCxfqySeflCT16dNHMTExysrKcu8BABUEwQxUEL6+vvrll18cVzHKzMxUSEiILly4oHHjxik+Pl6PPvqodu/erTFjxuj06dOSpMDAQI0ePbrIqxoBcD6u/AUAgIUwYgYAwEIIZgAALIRgBgDAQghmAAAshGAGAMBCCGYAACzk/wERD7soH3DGOAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 496.063x432 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Graph box and whisker plot of ages per week for cases and controls to see why it changes over time\n",
    "weeks = list(cases.keys())[::2]  # Select every other week\n",
    "\n",
    "\n",
    "sorted_weeks = sorted(weeks, key=lambda x: x if isinstance(x, datetime) \\\n",
    "                 else datetime(x.year, x.month, x.day))\n",
    "\n",
    "case_ages = [cases[week] for week in sorted_weeks]\n",
    "control_ages = [controls[week] for week in sorted_weeks]\n",
    "\n",
    "# Create a new figure\n",
    "fig = plt.figure(figsize=(figsize_2colw,6))\n",
    "\n",
    "# Calculate the positions for box plots\n",
    "case_positions = np.arange(len(sorted_weeks)) - 0.2\n",
    "control_positions = np.arange(len(sorted_weeks)) + 0.2\n",
    "\n",
    "# Create the box plot for cases\n",
    "case_box = plt.boxplot(case_ages, positions=case_positions, widths=0.4, patch_artist=True, boxprops=dict(facecolor='lightblue'))\n",
    "\n",
    "# Create the box plot for controls\n",
    "control_box = plt.boxplot(control_ages, positions=control_positions, widths=0.4, patch_artist=True, boxprops=dict(facecolor='darkseagreen'))\n",
    "\n",
    "# Customize the style of the median line\n",
    "for box in case_box['medians']:\n",
    "    box.set(color='yellow', linewidth=2)  # Change median line color and thickness\n",
    "\n",
    "for box in control_box['medians']:\n",
    "    box.set(color='red', linewidth=2)  # Change median line color and thickness\n",
    "\n",
    "# Add x-axis labels and title with increased font size\n",
    "plt.xlabel(\"Week\", fontsize=7)\n",
    "plt.ylabel(\"Years\", fontsize=7)\n",
    "plt.title(\"Distribution of Ages for Cases and Controls by Week (every other week included)\", fontsize=7)\n",
    "\n",
    "# Set x-axis tick labels with every 4th label\n",
    "plt.xticks(range(0, len(weeks), 4), sorted_weeks[::4], fontsize=7)\n",
    "\n",
    "# Create custom legend\n",
    "cases_patch = mpatches.Patch(color='lightblue', label='Cases')\n",
    "controls_patch = mpatches.Patch(color='darkseagreen', label='Controls')\n",
    "plt.legend(handles=[cases_patch, controls_patch], fontsize=7, loc='upper right')\n",
    "\n",
    "\n",
    "# Remove top and right borders\n",
    "sns.despine()\n",
    "\n",
    "plt.tight_layout()\n",
    "\n",
    "# Save the figure using fig instead of plt\n",
    "save_and_track_figure(fig, \"age_dist.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 244,
   "id": "15a0e237",
   "metadata": {
    "ExecuteTime": {
     "end_time": "2023-11-06T23:59:12.294830Z",
     "start_time": "2023-11-06T23:59:12.285349Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "3174"
      ]
     },
     "execution_count": 244,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#sample size\n",
    "\n",
    "len( case_ages[-1] ) + len( control_ages[-1] )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac8041df",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "covid_herit",
   "language": "python",
   "name": "covid_herit"
  },
  "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.8.10"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  },
  "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": 5
}