mice_med_in_R.ipynb

{
  "cells": [
    {
      "metadata": {
        "_uuid": "bad66ce8a7db3f5f31f92260d324fe02f8bf9062",
        "_cell_guid": "0c7c1a5a-6be1-403e-ba44-63bde753002c"
      },
      "cell_type": "markdown",
      "source": "## Missing Data: Using Multiple Imputation to Conduct a Mediation in R\n\nMultiple imputation is an extremely helpful and powerful tool when you have missing data. As a child development researcher, my data is particularly prone to missingness. Parents might not want to return to the lab, children get sleepy or fussy- it happens.\n\nHowever, there are ways to successfully and accurately reduce errors and bias caused by missing data. I found an article by Enders (2012) by to be extremely helpful in explaining the benefits of imputation as opposed to case-wise deletion- which is the most common practice in developmental research. \n\nBefore going forward, I used van Buuren's article to help me create the code, think about multiple imputation on a theoretical level, and also decide which parameters to use. I highly recommend you read through his article before pursuing multiple imputation!\n\nFor a recent project, I decided to use multiple imputation. I had rates of missingness ranging from 0% to 23%. Many sources (e.g., Little & Rubin, 2002; Bodner, 2008; White, Royston, & Wood, 2011) suggest creating as many datasets as the average rate of missing. So, in my case, I had an average of 10% missing, so I created 10 imputed datasets. \n\n**Note, if you want your results to be consistent (for sharing or publishing purposes, make sure to set your seed when imputing!)\n\nUsing the 'mice' package in R was very easy, and I had little trouble generating imputed datasets or pooling them. I also found visualizing the imputed data and comparing it to the original data using the 'VIM' package to be helpful. Using 'mice' to run linear regression models was also fairly simple.\n\nWhere I ran into trouble was using 'mice' and 'lavaan' to run a mediation analysis using my imputed data sets. Here is how I solved it- I hope it helps!"
    },
    {
      "metadata": {
        "_uuid": "da797d8d2bfc6c9f1b334c6288bf9e2b9861254e",
        "_cell_guid": "6d413b1a-fc61-4ccc-a55f-7cbc69d61cdc"
      },
      "cell_type": "markdown",
      "source": "First, load your libraries (download if needed). Next, read in a dataset ('airquality') and create some fake missing data. Then, have a look at the data:"
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "ba3b096a43d290bed0febd6980638fbd4c64455d",
        "trusted": false,
        "_cell_guid": "311feb50-0c5b-4633-8140-ab8f18d33a48"
      },
      "cell_type": "code",
      "source": "#Load Packages\nlibrary(VIM)\nlibrary(mice)\nlibrary(lattice)\n\n# Use the built-in airquality dataset\ndata <- airquality\n\n#create missing data\ndata[80:81,3] <- rep(NA, 2)\ndata[4:15,3] <- rep(NA,12)\ndata[1:5,2] <- rep(NA, 5)\n\n# Removing categorical variables\ndata <- data[-c(5,6)]\nsummary(data)",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "_uuid": "e35213c23fded3efc378ce96104cba1e8025fbb2",
        "_cell_guid": "7990d5e4-6061-48e2-913b-bcb0c56481a3"
      },
      "cell_type": "markdown",
      "source": "Next, let's look at how much data is missing for each variable. Be wary of missing data patterns higher than 5%."
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "e8dd1ff0c1aa374efe599d8acfcbfc4fa132a7d2",
        "trusted": false,
        "_cell_guid": "2a518b45-b825-4570-be41-890a5396fd60"
      },
      "cell_type": "code",
      "source": "# Look for missing > 5% variables\npMiss <- function(x){sum(is.na(x))/length(x)*100}\n\n# Check each column\napply(data,2,pMiss)\n\n# Check each row\napply(data,1,pMiss)",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "_kg_hide-output": true,
        "_uuid": "51d49624b0d2f09e4961ec71a35fd17bcca04612",
        "_cell_guid": "a838edfe-8cdf-4986-ad9b-780a9f1f85bc"
      },
      "cell_type": "markdown",
      "source": "Now let's visualize our missing data:"
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "c1d21ffbd3d2ef7f7c44b25ced563c1880f711b8",
        "trusted": false,
        "_cell_guid": "cd4cfcff-0e02-4087-be19-0b429bf11278"
      },
      "cell_type": "code",
      "source": "# Missing data pattern\nmd.pattern(data)",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "bc39c01ac8e31f1f8a2073ea710e51f7da241cb5",
        "trusted": false,
        "_cell_guid": "a0eaf789-8637-4e7f-8257-8f00ace68605"
      },
      "cell_type": "code",
      "source": "# Plot of missing data pattern\naggr_plot <- aggr(data, col=c('navyblue','red'), numbers=TRUE, sortVars=TRUE, labels=names(data), cex.axis=.7, gap=3, ylab=c(\"Histogram of missing data\",\"Pattern\"))",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "b3580ad3968a7a229c4da4b3361a9626e341a9fa",
        "trusted": false,
        "_cell_guid": "333507b2-b619-4b1e-a674-0982b4131e96"
      },
      "cell_type": "code",
      "source": "# Box plot\nmarginplot(data[c(1,2)])",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "_uuid": "b0790dd01307f317504a5ab58124cedbf14f1b54",
        "_cell_guid": "0dc0163d-b6b4-4b6b-b5dd-748cc4bfba01"
      },
      "cell_type": "markdown",
      "source": "Alright, now time to impute!!"
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "b900a8a726e272971f9db0605dd22d0e992f74ba",
        "trusted": false,
        "_cell_guid": "bda25703-ef16-4b82-9bff-4665f4a04fc7"
      },
      "cell_type": "code",
      "source": "# Impute missing data using mice\n#about 10% average missing data, so maxit= 10\ntempData <- mice(data,m=5,maxit=10,meth='pmm',seed=500)\nsummary(tempData)",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "a1be5573755f22ece99ba25db778357e490a485a",
        "trusted": false,
        "_cell_guid": "47ff2cab-a2fa-40f5-ae2a-befe6affc740"
      },
      "cell_type": "code",
      "source": "# Get imputed data (for the Ozone variable)\ntempData$imp$Ozone",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "872d3bb150e49c5bf67db4b03daadcfc9a4c2d30",
        "trusted": false,
        "_cell_guid": "75cef5d9-b2f1-4518-8684-97dbaad939c2"
      },
      "cell_type": "code",
      "source": "# Possible imputation models provided by mice() are\nmethods(mice)",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "9f6ebb72ab7b3c274494b5121ad79c2e7cd3dcf1",
        "trusted": false,
        "_cell_guid": "a9ed797f-c2b1-4cfe-8ea5-05a38f40291a"
      },
      "cell_type": "code",
      "source": "# What imputation method did we use?\ntempData$meth",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "3ea7a0b67048abdf7c720a6daf0b5d2ccff5dd98",
        "trusted": false,
        "_cell_guid": "e1cffbf1-f07d-465f-a613-253c4b37c72d"
      },
      "cell_type": "code",
      "source": "# Get completed datasets (observed and imputed)\ncompletedData <- complete(tempData,1)\nsummary(completedData)",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "_uuid": "fcf0f37eb128289dd094f39e0a3ffcc5609981f7",
        "_cell_guid": "f0ce19a5-9e65-443e-9a0b-63571328c1a0"
      },
      "cell_type": "markdown",
      "source": " Now, let's have a look at how the new imputed data (in red) looks compared to out original data (in blue):\n"
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "a9c4a5878d48111ae7d5b848e643dc8713eff07f",
        "trusted": false,
        "_cell_guid": "b7c47da0-1371-44d5-88dd-e9f5bb8b74a1"
      },
      "cell_type": "code",
      "source": "# Plots of imputed vs. orginal data\n# Scatterplot Ozone vs all\nxyplot(tempData,Ozone ~ Wind+Temp+Solar.R,pch=18,cex=1)",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "17859c56fbcd4a08644d6918b380e24af579ca97",
        "trusted": false,
        "_cell_guid": "646100a4-1e6d-4391-a7a2-459a4af5d5e2"
      },
      "cell_type": "code",
      "source": "# Density plot original vs imputed dataset\ndensityplot(tempData)",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "4fd4d348b93fc316e63d5bc22e6d7c09101a5036",
        "trusted": false,
        "_cell_guid": "bf5893f3-5147-4737-bd34-9aac8a7c166c"
      },
      "cell_type": "code",
      "source": "# Another take on the density: stripplot()\nstripplot(tempData, pch = 20, cex = 1.2)",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "_kg_hide-output": true,
        "_uuid": "9e5bf79678ced858eaf1be35e54195660d675f2f",
        "_cell_guid": "6c7c21ae-47c7-41c9-8c80-c17a277fd463"
      },
      "cell_type": "markdown",
      "source": "See the van Buren article for guidelines on how to measure goodness of fit for imputed data."
    },
    {
      "metadata": {
        "_uuid": "397c8a98e6851923c0872fc8b8c3850faf2f797e",
        "_cell_guid": "11ed90f6-cb4f-447f-8374-7191e1950b21"
      },
      "cell_type": "markdown",
      "source": "Finally- create a data frame from the imputed datasets, create your mediation model, extract and pool your parameters and check your final mediation model!!"
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "0b228fa70236e07b719b2dddfd6bff04e19e959f",
        "trusted": false,
        "_cell_guid": "9730e42e-30dc-4955-b44d-4e23c1d34606"
      },
      "cell_type": "code",
      "source": "# IMPUTE\n# create imputed dataframe\nimp1 <- miceadds::datlist_create(tempData)",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "48b69a712ec75afce4d640b7e9adcc5e3d229a9c",
        "trusted": false,
        "_cell_guid": "7d905362-84cb-43cf-8fbc-2d3c07a5230d"
      },
      "cell_type": "code",
      "source": "#create correlation table\ncorr_mice = miceadds::micombine.cor(mi.res=tempData )\ncorr_mice",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "c9c6f8b075778aaac92f0e98568b26cd3e09c760",
        "trusted": false,
        "_cell_guid": "a03d53e9-7e1d-4318-85c1-f8ab98827ff5"
      },
      "cell_type": "code",
      "source": "# look at correlation coefficients\nattr(corr_mice, \"r_matrix\") \n# look at p values\npvals = (corr_mice$p) ",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "8dab899182fc5a645eb83032bb4af4159e6ff6ce",
        "trusted": false,
        "_cell_guid": "a715bc70-8944-4a61-a27f-bb9e60736596"
      },
      "cell_type": "code",
      "source": "# Mediation\n##Create your mediation model\nmediation <- ' \n# direct effect\nTemp ~ cprime*Ozone\n# mediator \nSolar.R ~ a*Ozone\nTemp ~ b*Solar.R\n# indirect effect\nab := a*b\ntotal := cprime + (a*b)\ndirect:= cprime\n'",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "257d4c43702bfc7cdbd3a7383554d06ec83e7b14",
        "trusted": false,
        "_cell_guid": "47456874-c6b1-4721-9daf-a977c48a02d0"
      },
      "cell_type": "code",
      "source": "# analysis based on all imputed datasets\nmod6b <- lapply( imp1 , FUN = function(data){\n  res <- lavaan::sem(mediation , data = data )\n  return(res)\n} )",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "5a8eef34cf44f48c30ddf6e94b1dc46b940a6c73",
        "trusted": false,
        "_cell_guid": "a3248981-1ea6-43e2-8311-d8079a5decad"
      },
      "cell_type": "code",
      "source": "# extract all parameters\nqhat <- lapply( mod6b , FUN = function(ll){\n  h1 <- lavaan::parameterEstimates(ll)\n  parnames <- paste0( h1$lhs , h1$op , h1$rhs )\n  v1 <- h1$est\n  names(v1) <- parnames\n  return(v1)\n} )\nse <- lapply( mod6b , FUN = function(ll){\n  h1 <- lavaan::parameterEstimates(ll)\n  parnames <- paste0( h1$lhs , h1$op , h1$rhs )\n  v1 <- h1$se\n  names(v1) <- parnames\n  return(v1)\n} )",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "a0ee2f64a284859eee3d518828098d3736464090",
        "trusted": false,
        "_cell_guid": "cf72298c-3421-450c-9a3b-78899be3775a"
      },
      "cell_type": "code",
      "source": "# use mitml for mediation\nse2 <- lapply( se , FUN = function(ss){ ss^2 } ) # input variances\nresults <- mitml::testEstimates(qhat=qhat, uhat=se2)",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "collapsed": true,
        "_uuid": "183dee31f2d2685ac1a2a5d5585f2541938018d3",
        "trusted": false,
        "_cell_guid": "c21a30bd-06bb-4de1-abbb-e80d86c30867"
      },
      "cell_type": "code",
      "source": "#look at your results!\nresults",
      "execution_count": null,
      "outputs": []
    },
    {
      "metadata": {
        "_uuid": "87ea83d997842cb11a2dc376ff0564ccead296a1",
        "_cell_guid": "23eade16-0f7e-4eac-8b85-e7f7aa861c01"
      },
      "cell_type": "markdown",
      "source": "We did it!\n\nFor the full code, see documentation here on GitHub. \n\nHappy Imputing!"
    },
    {
      "metadata": {
        "_uuid": "6ea99caa830918fb9cd0f2d000b710f4523cee25",
        "_cell_guid": "ae2c00a3-0cd4-4941-b43d-20aaca02b613"
      },
      "cell_type": "markdown",
      "source": "### References:\n1. Bodner, T. E. (2008). What improves with increased missing data imputations? Structural Equation Modeling, 15(4), 651-675. doi:http://dx.doi.org/10.1080/10705510802339072\n2. Enders, C. K. (2013), Dealing With Missing Data in Developmental Research. Child Dev Perspectives, 7, 27–31. doi:10.1111/cdep.12008\n3. Little, R. J., & Rubin, D. B. (2002). Single imputation methods. Statistical Analysis with Missing Data, Second Edition, 59-74. doi:http://dx.doi.org/10.1002/9781119013563\n4. van Buuren, S., & Groothuis-Oudshoorn, K. (2011). mice: Multivariate Imputation by Chained Equations in R. Journal of Statistical Sotware, 45(3), 1-67. doi:http://dx.doi.org/10.18637/jss.v045.i03\n5. White, I. R., Royston, P., & Wood, A. M. (2011). Multiple imputation using chained equations: issues and guidance for practice. Statistics in Medicine, 30(4), 377-399. doi:http://dx.doi.org/10.1002/sim.4067"
    },
    {
      "metadata": {
        "_uuid": "a88908d31fa0c46dc37cb462a8ecfa1a6c458bb8",
        "trusted": false,
        "_cell_guid": "57b4d887-2ea3-46c5-beba-3e2f2194c2be"
      },
      "cell_type": "code",
      "source": "",
      "execution_count": null,
      "outputs": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "R",
      "language": "R",
      "name": "ir"
    },
    "language_info": {
      "mimetype": "text/x-r-source",
      "version": "3.4.2",
      "name": "R",
      "pygments_lexer": "r",
      "file_extension": ".r",
      "codemirror_mode": "r"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 1
}