IL181 Graphical Model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# $$\\text{Structural Causal Model}$$\n",
    "# $$\\text{of Instrumental Variables}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\text{Introduction}$\n",
    "Curve-fitting promises that regularities in nature embody the concept of causality. As much as highly correlated events, such as cloudiness and rain, might lead us to think that the first causes the second, it doesn't tell us much about the mechanism by which this sequence of events occurs. Furthermore, mere associations -as strong as they might be- wouldn't offer us an explanation as to why witnessing clouds doesn't necessarily mean that it will rain.\n",
    "Pearl et al. propose a framework to investigate these mechanisms based on statistical dependencies between variables. In other words, we can only discuss causality if a pair of variables X, Y are conditionally dependent. In such instances, we enumerate five possibilities:\n",
    "- X causes Y\n",
    "- Y causes X\n",
    "- X and Y share a common cause\n",
    "- X and Y are statistically dependent after conditioning on a variable Z that has both X and Y as causes (e.g., If we know that all Minervans are either smart or tall, selection bias occur when we condition on a person being a Minervan, our knowledge of their cognitive ability will infrom us about their height despite both being unconditionally independent).\n",
    "- A combination of the above\n",
    "\n",
    "The path to constructing a causal graph from observational data is known as **causal search** or **causal discovery**. Although a single dataset can have many potential causal graphs, our understanding of the phenomenon -even partially- would help us decide the most appropriate graph.\n",
    "\n",
    "In this report, we first dive into the structural causal model concept, discuss its utility, and finally present an overview of instrumental variables as an application."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\text{Structural Causal Models}$\n",
    "Pearl proposed a mechanism that would allow data scientists to craft conclusions that transcend associations and be able to pose queries on data rather than passively observing them. The idea is to construct a causal graph by which the data is believed to be generated, then expose it to observational data to deduce whether the links between variables are correctly drawn.\n",
    "A causal graph is a set of nodes representing random variables (covariates) linked together with directed edges. Each edge reflects whether a variable is believed to cause another. In other words, if X has a directed edge towards Y, then our knowledge about X informs us about the state of the variable Y.\n",
    "Before diving into examples of causal graphs, it's essential to showcase the ladder of causality proposed by Pearl. The ladder illustrates the levels of our understanding of a given phenomenon and the type of questions we can pose at each level.\n",
    "\n",
    "1. **Associations:** refer to the pattern we observe in data with little to no information about the process by which they're generated. The sort of questions we can pose are merely about correlation of events e.g., what would observing clouds tell me about whether it would rain? what would a survey inform us about an election?\n",
    "\n",
    "2. **Interventions:** At this stage we established our knowledge of the causal graph behind the data, and the question revolve around tweaking levers and observing how the model behaves e.g., if increase the temperature in the room, what would happen to the chocolat? how rising the interest rate affect the economic performance? Generally, it refers to the expected outcome after intervening in the system. As explained in the introduction, the causal search requires us to temper with the levers in the system to arrive at the most appropriate causal structure.\n",
    "\n",
    "3. **Counterfactuals:** The top rung of the ladder include questions that about a world different from the one that generated the observations. Counterfactuals are investigating what would happen had we not intervened in a system. \n",
    "\n",
    "Structural causal graphs can take many shapes, but the main building blocks are:\n",
    "1. **Chains:** Are a sequence of three variables or more with edges that have the same direction\n",
    "![chain.png](https://raw.githubusercontent.com/Tahahaha7/Causal_Inference_Tutorial/master/Chain.png)\n",
    "2. **Forks:** Also known as common cause structure when a variable has two directed edges towards two other variables\n",
    "![fork.png](https://raw.githubusercontent.com/Tahahaha7/Causal_Inference_Tutorial/master/Forks.png)\n",
    "3. **Colliders:** Are nodes that have two incoming directed edges from two different nodes\n",
    "![collider.png](https://raw.githubusercontent.com/Tahahaha7/Causal_Inference_Tutorial/master/Colliders.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\text{Instrumental Variables (Theoritical Background)}$\n",
    "The path between cause and effect can be paved with confounders, the challenge is to find an unbiased estimate of the causal effect without having to measure/control all variables in the system. Instrumental variables are common in social sciences studies as a strategy for finding Direct Causal Effect between two measureable variables $X \\rightarrow Y$ without necessary measuring the confounder variable H.\n",
    "\n",
    "Given the graph in following figure I will provide a theoretical proof that IV lead to an unbiased estimate of causal effect. However, the premise of IV rests on the following assumptions:\n",
    "- The variables in the system have a linear function relationship, so is to say that the regression coefficient of Y on X is the same as the \n",
    "- The IV affect the outcome variable $Y$ only through the cause $X$ and doesn't have a direct effect to $Y$\n",
    "- The IV doesn't affect the confounder variable $H$\n",
    "- The effect of IV on $X$ is statistically significant\n",
    "\n",
    "![graph.jpg](https://raw.githubusercontent.com/Tahahaha7/Causal_Inference_Tutorial/master/graph.jpg)\n",
    "\n",
    "Assuming that the exogenous variables $N_{X}$ and $N_{Y}$ are Independent and identically distributed random variables and that the functional relationship between the variables is linear.\n",
    "\n",
    "$$X=\\alpha\\cdot I+\\beta\\cdot H+N_x$$\n",
    "\n",
    "$$Y=\\delta\\cdot X+\\sigma\\cdot H+N_y$$\n",
    "\n",
    "The objective is to find the direct causal effect $\\delta$ of $X$ on $Y$. First, we substitute $X$ into the equation of $Y$, it yeilds the following equation:\n",
    "\n",
    "$$Y=\\delta\\left[\\alpha\\cdot I+\\beta\\cdot H+N_x\\right]+\\sigma\\cdot H+N_y$$\n",
    "\n",
    "We then factor out the remaining variables:\n",
    "$$Y=\\left[\\delta\\cdot\\beta+\\sigma\\right]\\cdot H+\\delta\\cdot\\left[\\alpha\\cdot I\\right]+\\gamma\\cdot N_x+N_y$$\n",
    "\n",
    "Notice that $I$ is independent of $H$ which relates to our previous assumption that the Instrumental Variable doesn't affect the confounder variable of interest. The methodology then uses Two-Stage Least Squares to estimate the regression coefficient of I on X which is $\\alpha$, then the fitted values $\\alpha \\cdot I$ will be regressed on Y to find $\\delta$, our direct causal effect of $X$ on $Y$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## $\\text{Instrumental Variables (Application)}$\n",
    "\n",
    "One of the prominent research questions in labor economics is estimating how wage is determined by education and work experience. Card (1993) studies the economic return to schooling and uses proximity to college as an instrumental variable to account for confounders.\n",
    "The dataset is a survey of high school graduates containing their wages, education, and other socio-economic metrics. The standard approach is to regress wages on education. However, education is not randomly assigned across the surveyed population, and we might fall for selection bias."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "    gender        ethnicity        score       fcollege   mcollege    home     \n",
       " male  :2139   other   :3050   Min.   :28.95   no :3753   no :4088   no : 852  \n",
       " female:2600   afam    : 786   1st Qu.:43.92   yes: 986   yes: 651   yes:3887  \n",
       "               hispanic: 903   Median :51.19                                   \n",
       "                               Mean   :50.89                                   \n",
       "                               3rd Qu.:57.77                                   \n",
       "                               Max.   :72.81                                   \n",
       " urban          unemp             wage           distance         tuition      \n",
       " no :3635   Min.   : 1.400   Min.   : 6.590   Min.   : 0.000   Min.   :0.2575  \n",
       " yes:1104   1st Qu.: 5.900   1st Qu.: 8.850   1st Qu.: 0.400   1st Qu.:0.4850  \n",
       "            Median : 7.100   Median : 9.680   Median : 1.000   Median :0.8245  \n",
       "            Mean   : 7.597   Mean   : 9.501   Mean   : 1.803   Mean   :0.8146  \n",
       "            3rd Qu.: 8.900   3rd Qu.:10.150   3rd Qu.: 2.500   3rd Qu.:1.1270  \n",
       "            Max.   :24.900   Max.   :12.960   Max.   :20.000   Max.   :1.4042  \n",
       "   education      income       region    \n",
       " Min.   :12.00   low :3374   other:3796  \n",
       " 1st Qu.:12.00   high:1365   west : 943  \n",
       " Median :13.00                           \n",
       " Mean   :13.81                           \n",
       " 3rd Qu.:16.00                           \n",
       " Max.   :18.00                           "
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "library(dagitty)\n",
    "library(lavaan)\n",
    "library(AER)\n",
    "data(\"CollegeDistance\")\n",
    "summary(CollegeDistance)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0gAAANICAMAAADKOT/pAAAAM1BMVEUAAABNTU1oaGh8fHyM\njIyampqnp6eysrK9vb2+vr7Hx8fQ0NDZ2dnh4eHp6enw8PD////ojgWfAAAACXBIWXMAABJ0\nAAASdAHeZh94AAAdgklEQVR4nO3biVYjuRJF0TSGgmb0/39tOzM9gg0eQtK9obPXW1WugqYk\nZRxP8IYVgLsNrRcAZEBIQABCAgIQEhCAkIAAhAQEICQgACEBAQgJCEBIQABCAgIQEhCAkIAA\nhAQEICQgACEBAQgJCEBIQABCAgIQEhCAkIAAhAQEICQgACEBAQgJCEBIQABCAgIQEhCAkIAA\nhAQEICQgACEBAQgJCEBIQABCAgIQEhCAkIAAhAQEICQgACEBAQgJCEBIQABC6sHH02J4evvl\nE76ehuHfZV9rYGRO4VQ68G+YPHye/YzH9YefL/tihHQSp5Lf87BYPxp9rX87W9IwnI/sx6fG\nrCoZTiW9z21AT8PTuc+5og5COolTSe/f9knb1+PL+NvLw/Aw3Vg/DD0Oi+fxxuj4Q9tfd580\nfqXF+oXU/KH1Zy42n/n1MDyuVm/LYVj+9jIsOUJKbzl8HP9xtFyNDSyG6bXRNqTDD61W25A2\nn7T58OP0oceDz3wc36d4mb/GS/XdqSCk9I6fi70Oi4/Vx2J4HT+w/FoX8LD9nKMPbf/L/Sdt\nP7z+y7fxL7+Ww9vm46vVYqz1dfpafSKk9I5Dehynfyxh/MD79sObh5mDD23/y/0nPU633uab\nYzxf41O6+ePr3zp+WjcipPSOQ9r8aW5k9xf7X75/6Putzc1h+3Rw87f/1k/xPo6eQnaGkNJ7\n3L1GevsqFtLqeXwpdf799fQIKb3n7bt27+NLmKCQdl99f/Pt3wOvkZDX7vtIy/FNte0Locef\nIX3/0PtxSPOH3/c3V7v//OQf+tLvzvvxNP1kw/jtoNWZt+Z+vGv3sG7ua3kc0tv+XbvpM1cv\n+xrX/8Er79ohueXhz9qd+mbR5ub+Qy+7bxkdfNL0zaOn6eb8meMj3Sak1/mfeK+/ORGE1IPX\ndQPL180fXhY/fnxhm8PuQ+N7B0/fX0iNf3nwkw3D0+f+v5x/sqHfjggJiEBIQABCAgIQEhCA\nkIAAhAQEICQgACEBAQgJCEBIQABCAgIQEhCAkIAAhIRv/vvlTziHkPANId2CkHDsezmUdBFC\nwjFCugkh4Rgh3YSQcORnN5R0CULCEUK6DSHhCCHdhpBwhJBuQ0g4dKoaSroAIeEQId2IkHCI\nkG5ESDhwuhlK+hsh4QAh3YqQcICQbkVI2DtXDCX9iZCwR0g3IyTsEdLNCAk753uhpL8QEnYI\n6XaEhB1Cuh0hYYeQbkdI2PqtFkr6AyFhi5DuQEjYIqQ7EBI2fm+Fkn5HSNggpHsQEjYI6R6E\nhNlfpVDSrwgJM0K6CyFhRkh3ISRM/u6Ekn5DSJgQ0n0ICRNCug8hYUJI9yEkjC6phJJ+QUgY\nEdKdCAkjQroTIWF1aSOUdB4hYUVI9yMkrAjpfoSEywuhpLMICYQUgJBASAEICVf0QUnnEBII\nKQAhgZACEBIIKQAh4Zo6KOkMQgIhBSAkEFIAQuredW1Q0mmE1D1CikBI3SOkCITUu2vLoKST\nCKl3hBSCkHpHSCEIqXPXd0FJpxBS5wgpBiF1jpBiEFLnCCkGIfXtlioo6QRC6hshBSGkvhFS\nEELq2m1NUNJPhNQ1QopCSF0jpCiE1LNbi6CkHwipZ4QUhpB6RkhhCKlnhBSGkDp2ew+U9B0h\ndYyQ4hBSxwgpDiH1654aKOkbQuoXIQUipH4RUiBC6tZ9LVDSMULqFiFFIqRuEVIkQurVvSVQ\n0hFC6hUhhSKkXhFSKELqFSGFIqRO3d8BJR0ipE4RUixC6hQhxSKkPkVUQEkHCKlPhBSMkPpE\nSMEIqUsxDVDSHiF1iZCiEVKXCCkaIfUoqgBK2iGkHhFSOELqESGFI6QeEVI4QupQ3PxT0hYh\ndYiQ4hFShwgpHiH1J3L6KWmDkPpDSAUQUn8IqQBC6k7s7FPSjJC6Q0glEFJ3CKkEQupN9ORT\n0oSQekNIRRBSbwipCELqDSEVQUidiZ97ShoRUmcIqQxC6gwhlUFIfSkx9ZS0IqTeEFIhhNQX\nQiqEkLpSZuYpiZA6Q0ilEFJXCKkUQupJqYmnJELqCiEVQ0g9IaRiCKknhFQMIXWk3LxTEiF1\nhJDKIaSOEFI5hNSPktPefUmE1A9CKoiQ+kFIBRFSN8rOeu8lEVI3CKkkQuoGIZVESL0oPemd\nl0RIvSCkogipF4RUFCH1gpCKIqROlJ/zvksipE4QUlmE1AlCKouQ+lBjyrsuiZD6QEiFEVIf\nCKkwQupCnRnvuSRC6gIhlUZIXSCk0gipB7UmvOOSCKkHhFQcIfWAkIojpB4QUnGE1IF6891v\nSYTUAUIqj5A6QEjlEVJ+Nae725IIKT9CqoCQ8iOkCggpvbqz3WtJhJQeIdVASOkRUg2ElF3t\nye60JELKjpCqIKTsCKkKQsqOkKogpOTqz3WfJRFScoRUByElR0h1EFJuLaa6y5IIKTdCqoSQ\nciOkSggptTYz3WNJhJQaIdVCSKkRUi2ElFmrie6wJELKjJCqIaTMCKkaQsqMkKohpMTazXN/\nJRFSYoRUDyElRkj1EFJeLae5u5IIKS9CqoiQ8iKkiggprbaz3FtJhJQWIdVESGkRUk2ElBYh\n1URIWbWe5Nb/fmWElFXrQW7971dGSFm1HuTW/35lhJRU+zluv4KaCCmp9mPcfgU1EVJS7ce4\n/QpqIqScFKZYYQ3VEFJOCkOssIZqCCknhSFWWEM1hJSSxgxrrKIOQkpJY4Q1VlEHIaWkMcIa\nq6iDkFLSGGGNVdRBSBmpTLDKOiogpIxUBlhlHRUQUkYqA6yyjgoIKSGd+dVZSWmElJDO+Oqs\npDRCSkhnfHVWUhoh5aM0vUprKYqQ8lEaXqW1FEVI+SgNr9JaiiKkdLRmV2s15RBSOlqjq7Wa\ncggpHa3R1VpNOYSUjtboaq2mHELKRm1y1dZTCCFloza4ausphJCyURtctfUUQkjJ6M2t3opK\nIKRk9MZWb0UlEFIyemOrt6ISCCkXxalVXFM4QspFcWgV1xSOkHJRHFrFNYUjpFQ0Z1ZzVbEI\nKRXNkdVcVSxCSkVzZDVXFYuQUtEcWc1VxSKkTFQnVnVdgQgpE9WBVV1XIELKRHVgVdcViJAS\n0Z1X3ZVFIaREdMdVd2VRCCkR3XHVXVkUQspDeVqV1xaCkPJQHlbltYUgpDyUh1V5bSEIKQ3t\nWdVe3f0IKQ3tUdVe3f0IKQ3tUdVe3f0IKQ3tUdVe3f0IKQv1SVVf350IKQv1QVVf350IKQv1\nQVVf350IKQn9OdVf4T0IKQn9MdVf4T0IKQn9MdVf4T0IKQeHKXVY480IKQeHIXVY480IKQeH\nIXVY480IKQWPGfVY5W0IKQWPEfVY5W0IKQWPEfVY5W0IKQWPEfVY5W0IKQOXCXVZ5w0IKQOX\nAXVZ5w0IKQOXAXVZ5w0IKQGf+fRZ6bUIKQGf8fRZ6bUIKQGf8fRZ6bUIyZ/TdDqt9SqE5M9p\nOJ3WehVC8uc0nE5rvQoh2fOaTa/VXo6Q7HmNptdqL0dI9rxG02u1lyMke16j6bXayxGSO7fJ\ndFvvhQjJndtguq33QoTkzm0w3dZ7IUIy5zeXfiu+BCGZ8xtLvxVfgpDM+Y2l34ovQUjeHKfS\ncc1/IiRvjkPpuOY/EZI3x6F0XPOfCMma50x6rvp3hGTNcyQ9V/07QrLmOZKeq/4dIVnzHEnP\nVf+OkJy5TqTrun9BSM5cB9J13b8gJGeuA+m67l8QkjHfefRd+TmEZMx3HH1Xfg4hGfMdR9+V\nn0NIvpyn0XntJxGSL+dhdF77SYTky3kYndd+EiHZ8p5F79X/REi2vEfRe/U/EZIt71H0Xv1P\nhGTLexS9V/8TIblyn0T39X9DSK7cB9F9/d8Qkiv3QXRf/zeEZMp/Dv13cIiQTPmPof8ODhGS\nKf8x9N/BIULylGEKM+xhh5A8ZRjCDHvYISRPGYYwwx52CMlThiHMsIcdQrKUYwZz7GJGSJZy\njGCOXcwIyVKOEcyxixkhOcoygVn2sSIkT1kGMMs+VoTkKcsAZtnHipAs5Zm/PDshJEN5xi/P\nTgjJUJ7xy7MTQvKTZ/oS7YWQ/KQZvlWivRCSnzTDt0q0F0Lyk2b4Von2Qkh20szeJMtuCMlO\nltGbZdkNIdnJMnqzLLshJDdZJm8ryX4IyU2SwdtJsh9CcpNk8HaS7IeQzCSZuwM5dkRIZnKM\n3aEcOyIkMznG7lCOHRGSlxxTdyzFngjJS4qh+ybFngjJS4qh+ybFngjJS4qh+ybFngjJSoqZ\n+yHDrgjJSoaR+ynDrgjJSoaR+ynDrgjJSYaJOyXBvgjJSYKBOynBvgjJSYKBOynBvgjJSIJ5\nO8N/Z4RkxH/czvHfGSEZ8R+3c/x3Rkg+/KftPPu9EZIP+2H7hf3eCMmH/bD9wn5vhOTDfth+\nYb83QrJhP2u/ct8dIdlwH7Xfue+OkGy4j9rv3HdHSC7cJ+0v5vsjJBfmg/Yn8/0RkgvzQfuT\n+f4IyYT5nF3Ae4eEZMJ7zC7hvUNCMuE9Zpfw3iEhefCesstY75GQPFgP2YWs90hIHqyH7ELW\neyQkD9ZDdiHrPRKSBesZu5jzLgnJgvOIXc55l4RkwXnELue8S0Jy4Dxh1zDeJyE5MB6wqxjv\nk5AcGA/YVYz3SUgGjOfrSr47JSQDvuN1Ld+dEpIB3/G6lu9OCUmf73Rdz3avhKTPdrhuYLtX\nQtJnO1w3sN0rIemzHa4b2O6VkOTZztZNXHdLSPJcR+s2rrslJHmuo3Ub190SkjrXybqV6X4J\nSZ3pYN3MdL+EpM50sG5mul9CEmc6V3fw3DEhifMcq3t47piQxHmO1T08d0xI2jyn6j6WeyYk\nbZZDdSfLPROSNsuhupPlnglJm+VQ3clyz4QkzXKm7ua4a0KS5jhS93PcNSFJcxyp+znumpCU\nOU5UBMN9E5Iyw4EKYbhvQlJmOFAhDPdNSMIM5ymI384JSZjfOEXx2zkhCfMbpyh+OyckXX7T\nFMdu74Sky26YAtntnZB02Q1TILu9E5Iuu2EKZLd3QpJlN0uh3HZPSLLcRimW2+4JSZbbKMVy\n2z0hqXKbpGhm+yckVWaDFM5s/4SkymyQwpntn5BEmc1RAV4nQEiivMaoBK8TICRRXmNUgtcJ\nEJImrykqw+oMCEmT1RAVYnUGhKTJaogKsToDQtJkNUSFWJ0BIUmymqFinE6BkCQ5jVA5TqdA\nSJKcRqgcp1MgJEVOE1SS0TkQkiKjASrK6BwISZHRABVldA6EJMhofgrzOQlCEuQzPqX5nAQh\nCfIZn9J8ToKQBPmMT2k+J0FIenympzybsyAkPTbDU4HNWRCSHpvhqcDmLAhJjs3sVOFyGoQk\nx2V06nA5DUKS4zI6dbicBiGpcZmcWkzOg5DUmAxONSbnQUhqTAanGpPzICQxJnNTkceJEJIY\nj7GpyeNECEmMx9jU5HEihCTGY2xq8jgRQtLiMTV1WZwJIWmxGJrKLM6EkLRYDE1lFmdCSFIs\nZqY6h1MhJCkOI1Ofw6kQkhSHkanP4VQISYnDxLRgcC6EpMRgYJowOBdCUmIwME0YnAshCTGY\nl0b0T4aQhOiPSyv6J0NIQvTHpRX9kyEkIfrj0or+yRCSDv1paUf+bAhJh/ywNCR/NoSkQ35Y\nGpI/G0KSIT8rTamfDiHJUB+VttRPh5BkqI9KW+qnQ0gq1CelNfHzISQV4oPSnPj5EJIK8UFp\nTvx8CEmE+JwI0D4hQhKhPSYKtE+IkERoj4kC7RMiJBHaY6JA+4QISYP2lGiQPiNC0iA9JCKk\nz4iQNEgPiQjpMyIkCdIzIkP5lAhJgvKI6FA+JUKSoDwiOpRPiZAUKE+IEuFzIiQFwgMiRfic\nCEmB8IBIET4nQhIgPB9idE+KkATojoca3ZMiJAG646FG96QISYDueKjRPSlCak93OvTInhUh\ntSc7HIJkz4qQ2pMdDkGyZ0VIzcnOhiTV0yKk5lRHQ5PqaRFSc6qjoUn1tAipNdXJUCV6XoTU\nmuhgyBI9L0JqTXQwZImeFyE1JjoXwjRPjJAa0xwLZZonRkiNaY6FMs0TI6TGNMdCmeaJEVJb\nmlOhTfLMCKktyaEQJ3lmhNSW5FCIkzwzQmpKcibkKZ4aITWlOBL6FE+NkJpSHAl9iqdGSC0p\nToQDwXMjpJYEB8KC4LkRUkuCA2FB8NwIqSHBeTChd3KE1JDeOLjQOzlCakhvHFzonRwhNaQ3\nDi70To6Q2tGbBh9yZ0dI7cgNgxG5syOkduSGwYjc2RFSM3KzYEXt9AipGbVR8KJ2eoTUjNoo\neFE7PUJqRW0S3IidHyG1IjYIdsTOj5BaERsEO2LnR0iNiM2BIa0TJKRGtMbAkdYJElIjWmPg\nSOsECakRrTFwpHWChNSG1hR4kjpDQmpDaghMSZ0hIbUhNQSmpM6QkJqQmgFbSqdISE0ojYAv\npVMkpCaURsCX0ikSUgtKE+BM6BwJqQWhAbAmdI6E1ILQAFgTOkdCakDo+pvTOUlCakDn8rvT\nOUlCakDn8rvTOUlCakDn8rvTOUlCqk/n6vuTOUtCqk/m4icgc5aEVJ/MxU9A5iwJqTqZa5+C\nymkSUnUqlz4HldMkpOpULn0OKqdJSLWpXPksRM6TkGoTufBpiJwnIdUmcuHTEDlPQqpN5MKn\nIXKehFSZyHVPRONECakyjcueicaJElJlGpc9E40TJaS6NK56LhJnSkh1SVz0ZCTOlJDqkrjo\nyUicKSFVJXHN01E4VUKqSuGS56NwqoRUlcIlz0fhVAmpJoUrnpHAuRJSTQIXPCWBcyWkmgQu\neEoC50pINQlc8JQEzpWQKhK43km1P1lCqqj95c6q/ckSUkXtL3dW7U+WkOppf7Xzan62hFRP\n84udWPOzJaR6ml/sxJqfLSFV0/xap9b6dAmpmtaXOrfWp0tI1bS+1Lm1Pl1CqqX1lc6u8fkS\nUi2EVBYhdYKQyiKkThBSWYTUBzoqre0JE1IlhFQaIXWBkEojpB7QUXlNz5iQ6iCk8gipA4RU\nHiHlR0c1tDxlQqqCkGogpPQIqQZCyo6O6mh4zoRUAyHVQUjJEVIdhJQcIdVBSLnRUS3tTpqQ\nKiCkWggpNUKqhZAyo6N6mp01IZVHSPUQUmKEVA8h5UVHNbU6bUIqjpBqIqS0CKkmQsqKjupq\ndN6EVBoh1UVISRFSXYSUFCHVRUg50VFtbU6ckAojpNoIKSVCqo2QMqKj+pqcOSGVRUj1EVJC\nhFQfIeVDRy20OHVCKoqQWiCkdAipBULKho7aaHDuhFQSIbVBSMkQUhuElAwhtUFIudBRK/VP\nnpAKIqRWCCkVQmqFkDKho3aqnz0hlUNI7RBSIoTUDiHlQUct1T59QiqGkFoipDQIqSVCyoKO\n2qp8/oRUCiG1RUhJEFJbhJQEIbVFSDnQUWt1rwAhFUJIrRFSCoTUGiFlQEftVb0GhFQGIbVH\nSAkQUnuE5I+OFNS8CoRUBCEpICR7hKSAkNzRkYaK14GQSiAkDYRkjpA0EJI5QtJASN7oSEW9\nK0FIBRCSCkKyRkgqCMkZHemodi0IKR4h6SAkY4Skg5B80ZGSWleDkMIRkhJCskVISgjJFR1p\nqXQ9CCkaIWkhJFOEpIWQTBGSFkLyREdq6lwRQgpGSGoIyRIhqSEkR3Skp8o1IaRYhKSHkAwR\nkh5C8kNHimpcFUIKRUiKCMkOISkiJDuEpIiQ3NCRpgrXhZAiEZImQjJDSJoIyQsdqSp/ZQgp\nECGpIiQrhKSKkJzQka7i14aQ4hCSLkIyQki6CMkHHSkrfXUIKQwhKSMkG4SkjJBsEJIyQnJB\nR9oKXx9CikJI2gjJBCFpIyQPdKSu7BUipCCEpI6QLBCSOkJyQEf6il4jQopBSPoIyQAh6SMk\nfXTkoORVIqQQhOSAkOQRkgNCkkdIDghJHR15KHidCCkCIXkgJHGE5IGQtNGRi3JXipACEJIL\nQpJGSC4ISRkd+Sh2rQjpfoTkg5CEEZIPQtJFR05KXS1CuhshOSEkWYTkhJBkEZITQlJFR14K\nXS9CuhcheSEkUYTkhZA00ZGbMleMkO5ESG4ISRIhuSEkLf8d/AonRa4cId2KkFwRkpb/dr/A\nSpELR0g3+2/Fw5Kb/fUiJBm7y0FHPnYXjJBk/Lcqc0lQUqmrRki3+2/7PxgpdNUI6Xb/0ZGj\n/zZXLhYh3e6/EhcExRW5boR0u/XloCNH/xGSFjoyVeDCEdId6MgVIQGSCAkIQEhAAEICAhAS\nEICQgACEBAQgJCAAIQEBCAkIQEhAAEICAhASEOCakIbJ4unzr08783dv1ywMuN9i+Bp/+1rP\n7ebGotC/dH1I65R+L+lcSA88+qGyx+F9/O1tPbXT3fj78FjoX7oupPHXr+Xw77Z/SiSk8/9f\nlGHn4K/+/oJvl34iyjl9VV+Gl/G3p+HfPLSbPxdwfUg3PzzKjNrZ/3/kbSHNj7Qyu+vVyav6\nPjyNv62f4c1D+zQ/QhVwQ0jT7+vnnA/Tw+TLw/AwVr6cljgufP7w6nlYPK9W/4bprmD952lA\nv4aH6Utsf2/kTEonYrigDxISceKqznP2uZ7Ux2F8QfIwvlR6e1y/PNk8q/q3WM/nfAXXk7yY\nH6/elsOwvPIV/a2PSMPwOCWynPpYjosdk18svjYhPY9//zZ9+N8+pO2T1tfh+bp1RjtZEiGZ\n+3lVF+PFGZ/QzU/qxiF9np90TCVN8/k0XcHH7SSvP3Vy3ZPA60P6XM5lLMe3QV6HxcfqYzG8\njv/883qNr9sHrPWHXza/LrZ/txpf9k2PtU/DH2/9FXfqQek4id2d1XDw3G39l8tp6bu7tc1T\nwd292vT4vP7j5+P0kIyaflzVx+FjNT0ajY9Kq4/xl2Ec0tfNNM7jO4w317O6fv3/Nsb3MX7C\ndc+ZbnrX7mu8/T6v821azxjycniZnuzNY/U+3fpc7SOaR216cF21fWY3+5nSUUjTndXjt5CW\n2/3v79YOQ9o9Pq//uBhvUlJt367q8/jAMj9ZGuY7+43No9A8vsN4c5zLr7m0G75Rc+v3kTYz\nd/jb53AczuGvB/fZ417ev03Yf42c3OD8XsPr/s5qv43X8W7rac5nd7d2+OHd4/Pmwfj7/UWr\njXbl8MDHZ0DzOw7j+wxPcyKfb8/Lg0u3mc7dtV+/rn/8+LgijOlrXPO5w4/bR79t3mL8PaSv\n+Wlq62d2q9OPSPuQ5hdzb8ebmf7y4F3LbyEdPD5vH5Kr7AR7367q5/piPG8uy/P6GcM4d8vd\nRT4d0up58fd3S78LDOmyR6R1bm+rh/bP7P56jXR0xse3Zt/u1g4+/P0AUM3Jp+sPm2dtD/P1\neFq/jH37/BnS4X/19u+h6GukH7e398GP0+2XzcuDX0P6GJYf7V87/Pmu3V8hfb9bIyQBJ67q\nekI3STwMb7vX8KvPn6+Rvr0wuvLa3RnSwauC1/UTu/ml3dmQPjcbWrR+ZnfB95H+COnH3Roh\nNXfyqv4blpu77edhufme5vvqY34ycfCu3TTJq+ntsod5nKs+Iu3fp/pazG+NfJ4N6WGYX1ys\n7wDaPrP75Scb9rfne6j33WamW8vda6SDu7Wfr5EeCam+M1f1ddi8wTxewbGQ8b2Eyfi3y/0L\no/nm+Mrodf/xy90b0uplMX/n5Gnzkw3LsyG9P8whfTV+V/i3n7Xb3z64s3pYP85+TXdhL+N7\ncf/md+12d2v7V4ZH79p9/4Io6txV/dj/wNf61vRe3NP6nv99fpo3fV/wff6Ml/U9/fz/bJh+\nsuHKnyVqcKXfhtbP7M45+lm7x903vV9231Hafx/p4G5tfqQ9uFfbvlAkJBPTFbv3a9z/Ja61\nLPYTuPc6Cml8E3TzY1jrW0+7t/iHx+luYH+3Nj/Sbu7VFrufbNj/ClXTc72vxxv/7wzHX+r+\nL3HlPzhE9A8E2Px8SsT/2a96SIti/9cq4Fov6yfjDwGPR02e2gH5EBIQgJCAAIQEBCAkIAAh\nAQEICQhASEAAQgICEBIQgJCAAIQEBCAkIAAhAQEICQhASEAAQgICEBIQgJCAAIQEBCAkIAAh\nAQEICQhASEAAQgICEBIQgJCAAIQEBCAkIAAhAQEICQhASEAAQgICEBIQgJCAAIQEBCAkIAAh\nAQEICQjwPyIxqg9VGJ3OAAAAAElFTkSuQmCC",
      "text/plain": [
       "plot without title"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 420,
       "width": 420
      },
      "text/plain": {
       "height": 420,
       "width": 420
      }
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "g <- dagitty('dag {\n",
    "    Proximity [pos=\"0,1\"]\n",
    "    Education [pos=\"1,1\"]\n",
    "    Wages [pos=\"3,1\"]\n",
    "    Confounders [pos=\"2,0\"]\n",
    "    \n",
    "    Proximity -> Education -> Wages\n",
    "    Education <- Confounders -> Wages}')\n",
    "plot(g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Cnfn _||_ Prxm\n",
       "Prxm _||_ Wags | Cnfn, Edct"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "impliedConditionalIndependencies(g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{ Confounders }\n",
       "{ Confounders, Proximity }"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "adjustmentSets(g, \"Education\", \"Wages\", type=\"all\" )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "lm(formula = log(wage) ~ education, data = CollegeDistance)\n",
       "\n",
       "Residuals:\n",
       "     Min       1Q   Median       3Q      Max \n",
       "-0.36208 -0.06271  0.03255  0.07794  0.32436 \n",
       "\n",
       "Coefficients:\n",
       "            Estimate Std. Error t value Pr(>|t|)    \n",
       "(Intercept) 2.213220   0.016228 136.379   <2e-16 ***\n",
       "education   0.002024   0.001166   1.737   0.0825 .  \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "Residual standard error: 0.1435 on 4737 degrees of freedom\n",
       "Multiple R-squared:  0.0006364,\tAdjusted R-squared:  0.0004254 \n",
       "F-statistic: 3.016 on 1 and 4737 DF,  p-value: 0.08249\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Regressing wage on education produces biased estimates\n",
    "wage_model_1 <- lm(log(wage) ~ education, data = CollegeDistance)\n",
    "summary(wage_model_1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "lm(formula = log(wage) ~ unemp + ethnicity + gender + urban + \n",
       "    education, data = CollegeDistance)\n",
       "\n",
       "Residuals:\n",
       "     Min       1Q   Median       3Q      Max \n",
       "-0.39998 -0.08223  0.02833  0.09486  0.37945 \n",
       "\n",
       "Coefficients:\n",
       "                    Estimate Std. Error t value Pr(>|t|)    \n",
       "(Intercept)        2.1519999  0.0168512 127.706   <2e-16 ***\n",
       "unemp              0.0135938  0.0007203  18.874   <2e-16 ***\n",
       "ethnicityafam     -0.0619139  0.0055990 -11.058   <2e-16 ***\n",
       "ethnicityhispanic -0.0535204  0.0052237 -10.246   <2e-16 ***\n",
       "genderfemale      -0.0091150  0.0039785  -2.291   0.0220 *  \n",
       "urbanyes           0.0089393  0.0048005   1.862   0.0626 .  \n",
       "education          0.0006723  0.0011121   0.605   0.5455    \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "Residual standard error: 0.1361 on 4732 degrees of freedom\n",
       "Multiple R-squared:  0.1026,\tAdjusted R-squared:  0.1015 \n",
       "F-statistic:  90.2 on 6 and 4732 DF,  p-value: < 2.2e-16\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Regressing wage on education and other metrics\n",
    "wage_model_2 <- lm(log(wage) ~ unemp + ethnicity + gender + urban + education, data = CollegeDistance)\n",
    "summary(wage_model_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the simple model **wage_model_1**, the effect of education is estimated to be around 0.002. For a model including other covariates **wage_model_2**, the estimate for education's impact is approximately 0.0006. We notice that after including other variables, the effect of education was downgraded by three folds.\n",
    "\n",
    "Proximity to schools can serve as an instrumental variable because it affects whether a person would attend college, but it doesn't affect wages directly (only through college attendance). Card (1993) argues that schools' geographical proximity offers us a lever to mitigate the selection bias among high school graduates.\n",
    "\n",
    "The rationale is to regress education on school proximity, then regress wages on the fitted values to capture an unbiased, direct effect of education on wages. We will experiment with the built-in function for instrumental variables in R using the simple model (only education) and the general model (using other covariates)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] -0.09318309\n"
     ]
    }
   ],
   "source": [
    "# compute the correlation between proximity and education\n",
    "print(cor(CollegeDistance$distance, CollegeDistance$education))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 74,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 0.008683088\n"
     ]
    }
   ],
   "source": [
    "# Perform the first stage regression and compute the fraction of explained variation\n",
    "R2 <- summary(lm(education ~ distance, data = CollegeDistance))$r.squared\n",
    "print(R2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "ivreg(formula = log(wage) ~ education | distance, data = CollegeDistance)\n",
       "\n",
       "Residuals:\n",
       "     Min       1Q   Median       3Q      Max \n",
       "-0.36022 -0.06094  0.03149  0.07747  0.32330 \n",
       "\n",
       "Coefficients:\n",
       "            Estimate Std. Error t value Pr(>|t|)    \n",
       "(Intercept) 2.221281   0.172731  12.860   <2e-16 ***\n",
       "education   0.001441   0.012509   0.115    0.908    \n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "Residual standard error: 0.1435 on 4737 degrees of freedom\n",
       "Multiple R-Squared: 0.0005835,\tAdjusted R-squared: 0.0003725 \n",
       "Wald test: 0.01326 on 1 and 4737 DF,  p-value: 0.9083 \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Estimate the IV regression of log(wage) on education using distance as the instrument\n",
    "wage_iv1 <- ivreg(log(wage) ~ education | distance, data = CollegeDistance)\n",
    "summary(wage_iv1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Call:\n",
       "ivreg(formula = log(wage) ~ unemp + ethnicity + gender + urban + \n",
       "    education | . - education + distance, data = CollegeDistance)\n",
       "\n",
       "Residuals:\n",
       "       Min         1Q     Median         3Q        Max \n",
       "-0.5885016 -0.1191974 -0.0001799  0.1452146  0.4576460 \n",
       "\n",
       "Coefficients:\n",
       "                    Estimate Std. Error t value Pr(>|t|)    \n",
       "(Intercept)        1.2171787  0.2018797   6.029 1.77e-09 ***\n",
       "unemp              0.0142234  0.0009648  14.743  < 2e-16 ***\n",
       "ethnicityafam     -0.0277621  0.0104342  -2.661  0.00782 ** \n",
       "ethnicityhispanic -0.0335043  0.0081520  -4.110 4.02e-05 ***\n",
       "genderfemale      -0.0076101  0.0052865  -1.440  0.15007    \n",
       "urbanyes           0.0064494  0.0063892   1.009  0.31283    \n",
       "education          0.0673242  0.0143812   4.681 2.93e-06 ***\n",
       "---\n",
       "Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1\n",
       "\n",
       "Residual standard error: 0.1805 on 4732 degrees of freedom\n",
       "Multiple R-Squared: -0.5786,\tAdjusted R-squared: -0.5806 \n",
       "Wald test: 54.89 on 6 and 4732 DF,  p-value: < 2.2e-16 \n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "wage_iv2 <- ivreg(log(wage) ~ unemp + ethnicity + gender + urban + education |\n",
    "                      . - education + distance, data = CollegeDistance)\n",
    "summary(wage_iv2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The estimate of education in the IV model is statistically significant compared to the biased model. The effect of education on wages is higher than any other predictor (0.0673***), which complies with our perception that education would have a crucial role in wage estimate predictions."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "4.0.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}