Itoga_simulation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "import statsmodels.api as sm\n",
    "\n",
    "from src.IterativeFitting import IterativeFitting as IF\n",
    "from src.CorrFuncs import covariance_matrix, trend_est"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reading in data\n",
    "study_df = pd.read_excel(\"SBPvPAD_data.xlsx\")\n",
    "\n",
    "# Taking care of non-log values in dataframe\n",
    "study_df.iloc[0,5:8] = np.log(study_df.iloc[0,5:8].to_numpy().astype(np.float64))\n",
    "study_df.iloc[5,5:8] = np.log(study_df.iloc[5,5:8].to_numpy().astype(np.float64))\n",
    "\n",
    "# Creating Itoga-specific dataframe\n",
    "study_df_i = study_df.loc[study_df[\"Author\"] == \"Itoga\"]\n",
    "study_df_i_bref = study_df_i.iloc[[0]]\n",
    "study_df_i_aref = study_df_i.iloc[1:]\n",
    "study_df_i = pd.concat([study_df_i_aref,study_df_i_bref],ignore_index=True)\n",
    "\n",
    "# Create exposure levels relative to reference exposure\n",
    "x_i = study_df_i[\"dose\"].to_numpy()[1:] - study_df_i[\"dose\"].to_numpy()[0]\n",
    "\n",
    "# Get log-odds and corresponding variance estimates\n",
    "L_i = study_df_i[\"logOR\"].to_numpy()[1:]\n",
    "v_i = study_df_i[\"std_error\"].to_numpy()[1:]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 4.8,  9.2, 13.9, 19.7, -5.9])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Author</th>\n",
       "      <th>dose</th>\n",
       "      <th>cases</th>\n",
       "      <th>controls</th>\n",
       "      <th>subjects</th>\n",
       "      <th>logOR</th>\n",
       "      <th>low</th>\n",
       "      <th>high</th>\n",
       "      <th>std_error</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>129.8</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>2738.333333</td>\n",
       "      <td>2988.0</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.0000</td>\n",
       "      <td>0.000</td>\n",
       "      <td>0.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>134.6</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>6048.333333</td>\n",
       "      <td>6298.0</td>\n",
       "      <td>-0.062</td>\n",
       "      <td>-0.1500</td>\n",
       "      <td>-0.023</td>\n",
       "      <td>0.032071</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>139.0</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>9341.333333</td>\n",
       "      <td>9591.0</td>\n",
       "      <td>0.077</td>\n",
       "      <td>0.0004</td>\n",
       "      <td>0.150</td>\n",
       "      <td>0.037778</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>143.7</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>8543.333333</td>\n",
       "      <td>8793.0</td>\n",
       "      <td>0.110</td>\n",
       "      <td>0.0220</td>\n",
       "      <td>0.210</td>\n",
       "      <td>0.047475</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>149.5</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>4401.333333</td>\n",
       "      <td>4651.0</td>\n",
       "      <td>0.190</td>\n",
       "      <td>0.1100</td>\n",
       "      <td>0.280</td>\n",
       "      <td>0.042929</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>Itoga</td>\n",
       "      <td>123.9</td>\n",
       "      <td>249.666667</td>\n",
       "      <td>786.333333</td>\n",
       "      <td>1036.0</td>\n",
       "      <td>0.230</td>\n",
       "      <td>0.1600</td>\n",
       "      <td>0.310</td>\n",
       "      <td>0.037879</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  Author   dose       cases     controls  subjects  logOR     low   high  \\\n",
       "0  Itoga  129.8  249.666667  2738.333333    2988.0  0.000  0.0000  0.000   \n",
       "1  Itoga  134.6  249.666667  6048.333333    6298.0 -0.062 -0.1500 -0.023   \n",
       "2  Itoga  139.0  249.666667  9341.333333    9591.0  0.077  0.0004  0.150   \n",
       "3  Itoga  143.7  249.666667  8543.333333    8793.0  0.110  0.0220  0.210   \n",
       "4  Itoga  149.5  249.666667  4401.333333    4651.0  0.190  0.1100  0.280   \n",
       "5  Itoga  123.9  249.666667   786.333333    1036.0  0.230  0.1600  0.310   \n",
       "\n",
       "   std_error  \n",
       "0   0.000000  \n",
       "1   0.032071  \n",
       "2   0.037778  \n",
       "3   0.047475  \n",
       "4   0.042929  \n",
       "5   0.037879  "
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "study_df_i"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initializing simulation slope and intercept parameters\n",
    "beta0 = -3.6289\n",
    "beta1 = 0.0246\n",
    "\n",
    "# Generating plausible x values to throw into probability generator\n",
    "xs = np.random.uniform(low=0,high=20,size=40000)\n",
    "\n",
    "# Function to generate probabilities of being a case v non-case\n",
    "p = lambda x: (np.exp(beta0 + beta1*x))/(1 + np.exp(beta0 + beta1*x))\n",
    "\n",
    "# Actually calculating probabilities on exposures as defined\n",
    "px = p(xs)\n",
    "\n",
    "# Actually assigning to case or not\n",
    "outcomes = np.array([np.random.binomial(n=1,p=p,size=1)[0] for p in px])\n",
    "\n",
    "# Constructing and sorting dataframe of outcomes and exposure\n",
    "df = np.stack([outcomes,xs],axis=1)\n",
    "df = df[np.argsort(df[:, 1])]\n",
    "\n",
    "# Observations at each category level\n",
    "C1 = df[df[:,1] < x_i[0]]\n",
    "C2 = df[np.logical_and(df[:,1] >= x_i[0], df[:,1] < x_i[1])]\n",
    "C3 = df[np.logical_and(df[:,1] >= x_i[1], df[:,1] < x_i[2])]\n",
    "C4 = df[np.logical_and(df[:,1] >= x_i[2], df[:,1] < x_i[3])]\n",
    "C5 = df[np.logical_and(df[:,1] >= x_i[3], df[:,1] < x_i[4])]\n",
    "C6 = df[df[:,1] >= x_i[-1]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(40000, 2)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Getting numbers of cases vs non-cases in each category\n",
    "cases1 = np.sum(C1[:,0])\n",
    "noncases1 = C1[:,0].shape[0] - cases1\n",
    "\n",
    "cases2 = np.sum(C2[:,0])\n",
    "noncases2 = C2[:,0].shape[0] - cases2\n",
    "\n",
    "cases3 = np.sum(C3[:,0])\n",
    "noncases3 = C3[:,0].shape[0] - cases3\n",
    "\n",
    "cases4 = np.sum(C4[:,0])\n",
    "noncases4 = C4[:,0].shape[0] - cases4\n",
    "\n",
    "cases5 = np.sum(C5[:,0])\n",
    "noncases5 = C5[:,0].shape[0] - cases5\n",
    "\n",
    "cases6 = np.sum(C6[:,0])\n",
    "noncases6 = C6[:,0].shape[0] - cases6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cases5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/p0/hl6rzzgd2v5802nvrkllqms40000gp/T/ipykernel_93809/2871705025.py:2: RuntimeWarning: invalid value encountered in scalar divide\n",
      "  crude = lambda x,y: (x/y) / (cases1/noncases1)\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "all input arrays must have the same shape",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[26], line 26\u001b[0m\n\u001b[1;32m     23\u001b[0m cats \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mappend(cats,in_cat5)\n\u001b[1;32m     24\u001b[0m cats \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mappend(cats,in_cat6)\n\u001b[0;32m---> 26\u001b[0m cats_out_df \u001b[38;5;241m=\u001b[39m \u001b[43mnp\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mstack\u001b[49m\u001b[43m(\u001b[49m\u001b[43m[\u001b[49m\u001b[43mdf\u001b[49m\u001b[43m[\u001b[49m\u001b[43m:\u001b[49m\u001b[43m,\u001b[49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43mcats\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43maxis\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m1\u001b[39;49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/numpy/core/shape_base.py:449\u001b[0m, in \u001b[0;36mstack\u001b[0;34m(arrays, axis, out, dtype, casting)\u001b[0m\n\u001b[1;32m    447\u001b[0m shapes \u001b[38;5;241m=\u001b[39m {arr\u001b[38;5;241m.\u001b[39mshape \u001b[38;5;28;01mfor\u001b[39;00m arr \u001b[38;5;129;01min\u001b[39;00m arrays}\n\u001b[1;32m    448\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(shapes) \u001b[38;5;241m!=\u001b[39m \u001b[38;5;241m1\u001b[39m:\n\u001b[0;32m--> 449\u001b[0m     \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mall input arrays must have the same shape\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m    451\u001b[0m result_ndim \u001b[38;5;241m=\u001b[39m arrays[\u001b[38;5;241m0\u001b[39m]\u001b[38;5;241m.\u001b[39mndim \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m\n\u001b[1;32m    452\u001b[0m axis \u001b[38;5;241m=\u001b[39m normalize_axis_index(axis, result_ndim)\n",
      "\u001b[0;31mValueError\u001b[0m: all input arrays must have the same shape"
     ]
    }
   ],
   "source": [
    "# Creating a function to get crude OR estimates\n",
    "crude = lambda x,y: (x/y) / (cases1/noncases1)\n",
    "\n",
    "# Actually calculating the crude OR estimates\n",
    "crudeor1 = crude(cases1,noncases1)\n",
    "crudeor2 = crude(cases2,noncases2)\n",
    "crudeor3 = crude(cases3,noncases3)\n",
    "crudeor4 = crude(cases4,noncases4)\n",
    "crudeor5 = crude(cases5,noncases5)\n",
    "crudeor6 = crude(cases6,noncases6)\n",
    "\n",
    "# Getting categories defined\n",
    "in_cat1 = np.zeros(C1.shape[0]) + 2\n",
    "in_cat2 = np.zeros(C2.shape[0]) + 3\n",
    "in_cat3 = np.zeros(C3.shape[0]) + 4\n",
    "in_cat4 = np.zeros(C4.shape[0]) + 5\n",
    "in_cat5 = np.zeros(C5.shape[0]) + 6\n",
    "in_cat6 = np.zeros(C6.shape[0]) + 7\n",
    "\n",
    "cats = np.append(in_cat1,in_cat2)\n",
    "cats = np.append(cats,in_cat3)\n",
    "cats = np.append(cats,in_cat4)\n",
    "cats = np.append(cats,in_cat5)\n",
    "cats = np.append(cats,in_cat6)\n",
    "\n",
    "cats_out_df = np.stack([df[:,0],cats], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(79376,)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cats.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Getting subjects and total number of cases\n",
    "N = np.array([cases1+noncases1,cases2+noncases2,cases3+noncases3,cases4+noncases4,cases5+noncases5])\n",
    "M1 = cases1 + cases2 + cases3 + cases4 + cases5\n",
    "\n",
    "# Initialization\n",
    "A0 = M1*N[1:]/(N.sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we will use the log-odds and variance estimates from Itoga directly to construct the covariance matrix for the adjusted method, and we will use them again to estimate the slope coefficient on the standard, non-correlation corrected method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/correlationCorrection/lib/python3.12/site-packages/cvxpy/reductions/solvers/solving_chain.py:354: FutureWarning: \n",
      "    You specified your problem should be solved by ECOS. Starting in\n",
      "    CXVPY 1.6.0, ECOS will no longer be installed by default with CVXPY.\n",
      "    Please either add an explicit dependency on ECOS or switch to our new\n",
      "    default solver, Clarabel, by either not specifying a solver argument\n",
      "    or specifying ``solver=cp.CLARABEL``.\n",
      "    \n",
      "  warnings.warn(ECOS_DEP_DEPRECATION_MSG, FutureWarning)\n"
     ]
    }
   ],
   "source": [
    "it_fit_ex = IF(L_i,A0,N,M1)\n",
    "A, B, a0, b0 = it_fit_ex.convexProgram()\n",
    "\n",
    "C = covariance_matrix(A,B,a0,b0,v_i**2)\n",
    "inv_C = np.linalg.inv(C)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recovering slope estimate for corrected correlation\n",
    "vb_star = 1/(np.dot(x_i,np.dot(inv_C,x_i)))\n",
    "b_star = vb_star*(np.dot(x_i,np.dot(inv_C,L_i)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.014094088127677437"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b_star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recovering standard slope estimate\n",
    "vb = 1/(np.dot(x_i,np.dot(np.linalg.inv(np.diag(v_i**2)),x_i)))\n",
    "b = vb*(np.dot(x_i,np.dot(np.linalg.inv(np.diag(v_i**2)),L_i)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0077151064467757356"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "correlationCorrection",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}