causality_w3.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Proof of propensity score delivering balance\n",
    "\n",
    "$$ P(X=x \\mid \\pi(x) = p, A=1) = \\frac{P(A=1 \\mid \\pi(x)=p, X=x)P(X=x \\mid \\pi(x)=p)}{P(A=1 \\mid \\pi(x) = p)}$$\n",
    "\n",
    "$$ P(A = 1 \\mid \\pi(x)=p) = p $$\n",
    "$$ P(A = 1 \\mid \\pi(x)=p, X=x) = p $$\n",
    "\n",
    "$$ P(X=x \\mid \\pi(x) = p, A=1) = P(X=x \\mid \\pi(x)=p) \\text{ e.g. independent of treatment}$$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width=6, repr.plot.height=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "install.packages(c(\"tableone\", \"Matching\", \"MatchIt\", \"optmatch\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(tableone)\n",
    "library(Matching)\n",
    "library(MatchIt)\n",
    "library(ggplot2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "std_diff <- function(dfr, tcol, gcol='treat'){\n",
    "    xt <- dfr[dfr[[gcol]]==1,tcol]\n",
    "    xc <- dfr[dfr[[gcol]]==0,tcol]\n",
    "    nt = 2; nc = 2;  \n",
    "    #nt = length(xt);nc = length(xc);\n",
    "    (mean(xt) - mean(xc))/(sqrt( ((nt-1)*var(xt) + (nc-1)*var(xc))/(nt+nc-2)))\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data(lalonde)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The outcome is re78 – post-intervention income.\n",
    "\n",
    "The treatment is treat – which is equal to 1 if the subject received the labor training and equal to 0 otherwise.\n",
    "\n",
    "The potential confounding variables are: age, educ, black, hispan, married, nodegree, re74, re75."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### q1\n",
    "Find the standardized differences for all of the confounding variables (pre-matching). What is the standardized difference for married (to nearest hundredth)?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xt = lalonde[lalonde$treat==1,]$married\n",
    "xc = lalonde[lalonde$treat==0,]$married\n",
    "#nt = length(xt);nc = length(xc);\n",
    "nt = 2; nc = 2;\n",
    "(mean(xt) - mean(xc))/(sqrt( ((nt-1)*var(xt) + (nc-1)*var(xc))/(nt+nc-2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> q1: 0.719"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### q2\n",
    "What is the raw (unadjusted) mean of real earnings in 1978 for treated subjects minus the mean of real earnings in 1978 for untreated subjects?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean(lalonde[lalonde$treat==1,]$re78) - mean(lalonde[lalonde$treat==0,]$re78)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">q2: -$635"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### q3\n",
    "What are the minimum and maximum values of the estimated propensity score?\n",
    "\n",
    "Fit a propensity score model. Use a logistic regression model, where the outcome is treatment. Include the 8 confounding variables in the model as predictors, with no interaction terms or non-linear terms (such as squared terms). Obtain the propensity score for each subject."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data(lalonde)\n",
    "fit <- glm(treat ~ age + educ + black + hispan + married + nodegree + re74 + re75,\n",
    "           data=lalonde, family = binomial(link = \"logit\"))\n",
    "prop_hat <- predict(fit, newdata = lalonde, type=\"response\")\n",
    "lalonde$pscore <- prop_hat\n",
    "min(prop_hat); max(prop_hat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> q3: pscore (0.009, 0.853)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```g1 = ggplot(lalonde, aes(x=lalonde$pscore, group=treat, fill=treat)) + \n",
    "    geom_histogram(position=\"dodge\",binwidth=0.025) + theme_bw()\n",
    "multiplot(g1,  cols=2)```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### q4\n",
    "\n",
    "Now carry out propensity score matching using the Match function.\n",
    "\n",
    "Before using the Match function, first do:\n",
    "\n",
    ">set.seed(931139)\n",
    "\n",
    "Setting the seed will ensure that you end up with a matched data set that is the same as the one used to create the solutions.\n",
    "\n",
    "Use options to specify pair matching, without replacement, no caliper.\n",
    "\n",
    "Match on the propensity score itself, not logit of the propensity score. Obtain the standardized differences for the matched data.\n",
    "\n",
    "What is the standardized difference for married?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One alternative would be\n",
    "```\n",
    "m.out <- matchit(treat ~ age + educ + black + hispan + married + nodegree + re74 + re75,\n",
    "                 data=lalonde, method = \"nearest\")\n",
    "p1<- plot(m.out, type='jitter')\n",
    "p2<- plot(m.out, type='hist')\n",
    "```\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set.seed(931139)\n",
    "pmatch <- Match(Tr = lalonde$treat, M=1, X=lalonde$pscore, replace=FALSE, caliper=NaN)\n",
    "matched <- lalonde[unlist(pmatch[c('index.treated', 'index.control')]), ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "std_diff(matched, 'married', 'treat')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> q4: -0.027"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### q5\n",
    "For the propensity score matched data:\n",
    "Which variable has the largest standardized difference?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lapply(c('age', 'nodegree', 're74', 'black'), function(v) std_diff(matched, v, 'treat'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> q5: black, 0.85"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### q6\n",
    "Re-do the matching, but use a caliper this time. Set the caliper=0.1 in the options in the Match function.\n",
    "Again, before running the Match function, set the seed: 931139\n",
    "\n",
    "How many matched pairs are there?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "set.seed(931139)\n",
    "pmatch <- Match(Tr = lalonde$treat, M=1, X=lalonde$pscore, replace=FALSE, caliper=0.1)\n",
    "matched <- lalonde[unlist(pmatch[c('index.treated', 'index.control')]), ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "nrow(matched)/2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> q6: 111"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### q7\n",
    "Use the matched data set (from propensity score matching with caliper=0.1) to carry out the outcome analysis.\n",
    "\n",
    "For the matched data, what is the mean of real earnings in 1978 for treated subjects minus the mean of real earnings in 1978 for untreated subjects?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mean(matched[matched$treat==1,]$re78) - mean(matched[matched$treat==0,]$re78)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> q7: 1246.81"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### q8\n",
    "\n",
    "Use the matched data set (from propensity score matching with caliper=0.1) to carry out the outcome analysis.\n",
    "\n",
    "Carry out a paired t-test for the effect of treatment on earnings. What are the values of the 95% confidence interval?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t.test(matched[matched$treat==1,c('re78')],\n",
    "       matched[matched$treat==0,c('re78')],\n",
    "       conf.level = 0.95, paired=TRUE)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> q8: (-420.03, 2913.64)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tbone <- CreateTableOne(data=matched, strata='treat', test=TRUE, smd=TRUE )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "--- "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Multiple plot function\n",
    "#\n",
    "# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)\n",
    "# - cols:   Number of columns in layout\n",
    "# - layout: A matrix specifying the layout. If present, 'cols' is ignored.\n",
    "#\n",
    "# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),\n",
    "# then plot 1 will go in the upper left, 2 will go in the upper right, and\n",
    "# 3 will go all the way across the bottom.\n",
    "#\n",
    "multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {\n",
    "  library(grid)\n",
    "\n",
    "  # Make a list from the ... arguments and plotlist\n",
    "  plots <- c(list(...), plotlist)\n",
    "\n",
    "  numPlots = length(plots)\n",
    "\n",
    "  # If layout is NULL, then use 'cols' to determine layout\n",
    "  if (is.null(layout)) {\n",
    "    # Make the panel\n",
    "    # ncol: Number of columns of plots\n",
    "    # nrow: Number of rows needed, calculated from # of cols\n",
    "    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),\n",
    "                    ncol = cols, nrow = ceiling(numPlots/cols))\n",
    "  }\n",
    "\n",
    " if (numPlots==1) {\n",
    "    print(plots[[1]])\n",
    "\n",
    "  } else {\n",
    "    # Set up the page\n",
    "    grid.newpage()\n",
    "    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))\n",
    "\n",
    "    # Make each plot, in the correct location\n",
    "    for (i in 1:numPlots) {\n",
    "      # Get the i,j matrix positions of the regions that contain this subplot\n",
    "      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))\n",
    "\n",
    "      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,\n",
    "                                      layout.pos.col = matchidx$col))\n",
    "    }\n",
    "  }\n",
    "}"
   ]
  }
 ],
 "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": "3.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}