TrajectoryPredictionAnalysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Five Maternal Characteristics Determine the 12-month Course of Postpartum Depression Symptoms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Objective: \n",
    "Postpartum depression is a heterogeneous disorder in phenotype and etiology. Characterizing the longitudinal course of depressive symptoms over the first year after birth and identifying variables that predict distinct symptom trajectories will expedite efficient mental health treatment planning.  \n",
    "\n",
    "The purpose was to \n",
    "1. Determine 12-month trajectories of postpartum depressive symptoms\n",
    "2. Identify characteristics that predict the trajectories\n",
    "3. Provide a computational algorithm that predicts specific trajectory membership. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Analysis Methods:\n",
    "\n",
    "The STATA/SAS plugin traj was used for estimating group-based trajectory models of SIGH-ADS29 total score during the first year postpartum. The number of trajectories and the shape of the pattern of change for each over time were determined from the comparison of Bayesian Information Criterion (BIC) and visual fit across all possible models.\n",
    "\n",
    "After selecting a final model, we used the newly-created trajectory group membership variable to examine differences for demographic and other potential predictors across these categories.  We employed Pearson’s Chi-square test or Fisher’s exact test for categorical variables and one-way ANOVA for continuous variables. \n",
    "\n",
    "We created a multinomial logistic regression model with the trajectory group membership as the outcome to assess the capacity of baseline variables to predict one-year postpartum depression trajectories. Each demographic predictor was tested in an unadjusted multinomial logistic regression model. Predictors that were individually statistically significant at the 10% level were assessed in the multivariable model. Development of a final predictive model for trajectory membership employed 10-fold cross validation techniques using a penalized multinomial regression model through the R package CARET (Classification and Regression Training). Variables that remained statistical significant at the 10% level were kept in the reduced model. \n",
    "    \n",
    "Missing covariates were handled by multiple imputation using the R package MICE (Multiple Imputation by Chained Equations). The model results were compared with results from complete case analysis. \n",
    "    \n",
    "Post hoc analyses compared the two most severe trajectories with respect to the use of health services at baseline, month 6 and month 12 via a series of generalized linear mixed effect models.  In each model the use of health services was regressed on fixed trajectory groups, time, and the trajectory-by-time interaction terms (assuming a 10% level of significance for interaction terms) with an additional random intercept term. Analyses were performed using STATA 27 and R 3.2.0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 1: Determine the trajectory patterns -STATA traj \n",
    "\n",
    "Tried 2,3 and 4 groups --> 3 groups gave the best BIC\n",
    "\n",
    "At least one group has <5% of the sample if group number >4-->stopped testing at 4 groups\n",
    "\n",
    "Then, determine the appropriate form of the groups\n",
    "intercept only, linear, quadratic and cubic\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/amysheep/NUTrajectoryPredictionAnalysis/master/Graph-213_06142017.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 1 results: Three distinct trajectories of depressive symptoms were identified: 1) gradual remission (50.4%), 2) partial improvement (41.8%), and 3) chronic, non-remission (7.8%). \n",
    "\n",
    "PS: this is an unsupervised learning/clustering problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 2: Which factors are significantly associated with the trajactory patterns?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### some data cleaning using tidyverse\n",
    "\n",
    "library(caret)\n",
    "library(tidyverse)\n",
    "library(nnet)\n",
    "library(broom)\n",
    "\n",
    "data=read.csv(\"../Data/dataforsastraj06142017_out.csv\",h=T)\n",
    "diag=read.csv(\"../Data/diag.csv\")\n",
    "\n",
    "\n",
    "# merge comorbidity data from diag to create Anxiety Yes No variable\n",
    "colnames(diag)[1]=\"id\"\n",
    "diag1=filter(diag,id%in%data$id)%>%select(id,preg,category,short_di)\n",
    "anxid=filter(diag1,short_di%in%c(\"Social Phobia\",\"Posttraumatic Stress Disorder\",\"Panic Disorder With Agoraphobia\",\"Panic Disorder Without Agoraphobia\",\"Obsessive-Compulsive Disorder\",\"Generalized Anxiety Disorder\"),category==\"Anxiety Disorders\")%>%select(id,preg)\n",
    "\n",
    "anxid=unique(anxid)\n",
    "anxid$anx=1\n",
    "dataall=left_join(data,anxid,by=c(\"id\",\"preg\"))%>%mutate(anx=ifelse(is.na(anx),0,anx))\n",
    "\n",
    "# also recoded some of the rare categories in categorical variables to eliminate near-zero-variance issue\n",
    "# code not show here\n",
    "\n",
    "### Create Table 1\n",
    "\n",
    "library(compareGroups)\n",
    "\n",
    "t=compareGroups(X_traj_Group~.,data=t1data,include.miss=F)\n",
    "t1=createTable(t,show.all = T,hide=c(hispanic=\"0\",harm01=\"Never\",marital01=\"Single or divorced\",edu12=\"Some college or less\"),hide.no = \"no\",digits=2)\n",
    "t1 "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 2 results:\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/amysheep/NUTrajectoryPredictionAnalysis/master/Tab1.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 3: Complete case multinomial logistic regression (n=417)\n",
    "\n",
    "--Full model with p<0.1 factors from table 1 (10-fold cross validation R package `caret`-classification and regression training) \n",
    "    Note: Variables are tested for collinearity, near-zero-variance, outliers etc. before modelling.\n",
    "\n",
    "--Penalized multinomial regression model (tuning parameter-weight decay:decay is the regularization parameter to avoid over-fitting)\n",
    "\n",
    "--Use Cross validation to find the optimal tuning parameter that gives the best accuracy \n",
    "\n",
    "--Use upsampling within resampling for class imbalance (group 3 is way too low compared to 1 and 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# specify type of resampling\n",
    "train_control<- trainControl(method=\"cv\", number=10,savePredictions = TRUE,sampling=\"up\")\n",
    "\n",
    "# fitting multinomial logistic regression -multinom function from nnet package\n",
    "set.seed(06132017)\n",
    "model<- train(X_traj_Group~edu12+anx+physab+adltph+parity+n_cmdx+epdstota+gasscore00, \n",
    "              data=t1data, trControl=train_control, method=\"multinom\",\n",
    "              na.action=na.omit,savePredictions = TRUE,\n",
    "              tuneGrid=expand.grid(.decay=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)),trace=F)\n",
    "# optimized model\n",
    "finalmodel1 <- summary(model$finalModel)\n",
    "\n",
    "# there are a few variables had p>0.1 in the above full model\n",
    "# remove them and create a reduced model \n",
    "\n",
    "set.seed(061320171)\n",
    "model2<- train(X_traj_Group~adltph+anx+parity+n_cmdx+gasscore00,\n",
    "               data=t1data, trControl=train_control, method=\"multinom\",\n",
    "               na.action=na.omit,savePredictions = TRUE,\n",
    "               tuneGrid=expand.grid(.decay=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)),trace=F)\n",
    "\n",
    "\n",
    "# use relevel to get 2vs3 comparisons\n",
    "t1data$re_group=relevel(t1data$X_traj_Group,ref=\"2\")\n",
    "model2a<- train(re_group~adltph+anx+parity+n_cmdx+gasscore00, \n",
    "                data=t1data, trControl=train_control, \n",
    "                method=\"multinom\",na.action=na.omit,savePredictions = TRUE,\n",
    "                tuneGrid=expand.grid(.decay=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)),trace=F)\n",
    "\n",
    "# use Broom package to clean up output table and calculate ORs and 95%CIs \n",
    "\n",
    "library(broom)\n",
    "\n",
    "tidy(model2$finalModel)->m2\n",
    "tidy(model2a$finalModel)->m2a\n",
    "\n",
    "table_func <- function(m){\n",
    "  or <- round(exp(m$estimate),3)\n",
    "  or_low <- round(exp(m$estimate-1.96*m$std.error),3)\n",
    "  or_high <- round(exp(m$estimate+1.96*m$std.error),3)\n",
    "  group <- m$y.level\n",
    "  var <- m$term\n",
    "  p <- round(m$p.value,3)\n",
    "  out <- data.frame(cbind(group,var,or,or_low,or_high,p))\n",
    "  return(out)\n",
    "}\n",
    "\n",
    "table_func(m2)[c(2:6,8:12),]\n",
    "table_func(m2a)[c(2:6,8:12),]\n",
    "\n",
    "# confusion matrix\n",
    "\n",
    "filter(model2$pred,decay==0.7)%>%select(obs,pred)->cm\n",
    "confusionMatrix(data=cm$pred,reference=cm$obs)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 3 results :\n",
    "Accuracy was used to select the optimal model using the largest value.\n",
    "The final value used for the model was decay = 0.7.\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/amysheep/NUTrajectoryPredictionAnalysis/master/fig2.png\" style=\"width:500px;height:600px;\">\n",
    "\n",
    "Note: Could also split the data into training/testing sets before model fitting. Did not show code here to save space.\n",
    "\n",
    "The accuracy isn't that high."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 4: Model with multiple imputation (5 different sets of imputations) for missing covariates (R package `MICE`-multiple imputation by chained equations)\n",
    "\n",
    "A. Inspecting the missing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "library(mice)\n",
    "mdpt <- md.pattern(select(t1data,edu12,anx,physab,adltph,parity,n_cmdx,epdstota,gasscore00))\n",
    "colnames(mdpt)[9] <- \"Number of missing Vars\"\n",
    "rownames(mdpt)[7] <- \"Total Missing\"\n",
    "rownames(mdpt) <- paste0(\"N=\",rownames(mdpt))\n",
    "mdpt"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/amysheep/NUTrajectoryPredictionAnalysis/master/fig3.png\" style=\"width:600px;height:600px;\">\n",
    "\n",
    "-- There are six different patterns (n=417 complete data; n=57 missing number of chronic medical conditions; n=14 missing gas;\n",
    "n=9 missing child and adult physical abuse; n=8 missing gas and chronic med conditions; n=2 missing child and adult physical abuse and chronic med conditions )\n",
    "\n",
    "-- N miss=11 for child and adult physical abuse; N miss=22 for GAS; N=67 for number of chronic medical conditions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "B. Create 5 imputated datasets "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "imp <- mice(select(t1data,X_traj_Group,edu12,anx,physab,adltph,parity,n_cmdx,epdstota,gasscore00),\n",
    "            me=c(\"\",\"\",\"\",\"logreg\",\"logreg\",\"\",\"pmm\",\"\",\"pmm\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "C. Diagnostic checking \n",
    "(whether imputated values are plausible e.g.scores and counts are positive numbers)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "stripplot(imp,n_cmdx+adltph~.imp,pch=20,cex=1.2,layout=c(1,2))\n",
    "stripplot(imp,physab+gasscore00~.imp,pch=20,cex=1.2,layout=c(1,2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src=\"https://raw.githubusercontent.com/amysheep/NUTrajectoryPredictionAnalysis/master/fig4.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "D. Analysis of imputed data (Full model with 8 predictors)\n",
    "\n",
    "-- Within each imputed dataset we repeat fitting the prediction model with cross validation \n",
    "\n",
    "-- The final estimates are pooled from the 5 sets of estimates\n",
    "\n",
    "-- EPDS is also a significant factor in addition to the 5 found in the data without imputation\n",
    "\n",
    "-- Final reduced model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fit1=train(X_traj_Group~edu12+anx+physab+adltph+parity+n_cmdx+epdstota+gasscore00, trControl=train_control, method=\"multinom\",savePredictions = TRUE,tuneGrid=expand.grid(.decay=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)),trace=F,data=mice::complete(imp))\n",
    "fit2=train(X_traj_Group~edu12+anx+physab+adltph+parity+n_cmdx+epdstota+gasscore00, trControl=train_control, method=\"multinom\",savePredictions = TRUE,tuneGrid=expand.grid(.decay=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)),trace=F,data=mice::complete(imp,2))\n",
    "fit3=train(X_traj_Group~edu12+anx+physab+adltph+parity+n_cmdx+epdstota+gasscore00, trControl=train_control, method=\"multinom\",savePredictions = TRUE,tuneGrid=expand.grid(.decay=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)),trace=F,data=mice::complete(imp,3))\n",
    "fit4=train(X_traj_Group~edu12+anx+physab+adltph+parity+n_cmdx+epdstota+gasscore00, trControl=train_control, method=\"multinom\",savePredictions = TRUE,tuneGrid=expand.grid(.decay=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)),trace=F,data=mice::complete(imp,4))\n",
    "fit5=train(X_traj_Group~edu12+anx+physab+adltph+parity+n_cmdx+epdstota+gasscore00, trControl=train_control, method=\"multinom\",savePredictions = TRUE,tuneGrid=expand.grid(.decay=c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9)),trace=F,data=mice::complete(imp,5))\n",
    "#The as.mira() function takes the results of repeated complete-data analysis stored as a list, and turns it into a mira object that can be pooled.\n",
    "fit=as.mira(list((fit1$finalModel),(fit2$finalModel),(fit3$finalModel),(fit4$finalModel),(fit5$finalModel)))\n",
    "\n",
    "# calculate pooled estimates\n",
    "beta1=summary(fit1$finalModel)$coefficients\n",
    "beta2=summary(fit2$finalModel)$coefficients\n",
    "beta3=summary(fit3$finalModel)$coefficients\n",
    "beta4=summary(fit4$finalModel)$coefficients\n",
    "beta5=summary(fit5$finalModel)$coefficients\n",
    "\n",
    "betaall=rbind(beta1,beta2,beta3,beta4,beta5)\n",
    "beta2vs1=betaall[rownames(betaall)==2,]\n",
    "beta3vs1=betaall[rownames(betaall)==3,]\n",
    "\n",
    "# pooled beta estimate\n",
    "grandm2vs1=apply(beta2vs1,2,mean)\n",
    "grandm3vs1=apply(beta3vs1,2,mean)\n",
    "\n",
    "# pooled se of beta \n",
    "se1=summary(fit1$finalModel)$standard.errors\n",
    "se2=summary(fit2$finalModel)$standard.errors\n",
    "se3=summary(fit3$finalModel)$standard.errors\n",
    "se4=summary(fit4$finalModel)$standard.errors\n",
    "se5=summary(fit5$finalModel)$standard.errors\n",
    "\n",
    "seall=rbind(se1,se2,se3,se4,se5)\n",
    "se2vs1=seall[rownames(seall)==2,]\n",
    "se3vs1=seall[rownames(seall)==3,]\n",
    "\n",
    "withinse2vs1=apply(se2vs1,2,mean)\n",
    "withinse3vs1=apply(se3vs1,2,mean)\n",
    "\n",
    "# calculate p-value\n",
    "between2vs1=NULL\n",
    "for (i in 1:9){\n",
    "  t=sum((beta2vs1[,i]-grandm2vs1[i])^2)/(5-1)\n",
    "  between2vs1=rbind(between2vs1,t)\n",
    "}\n",
    "\n",
    "between3vs1=NULL\n",
    "for (i in 1:9){\n",
    "  t=sum((beta3vs1[,i]-grandm3vs1[i])^2)/(5-1)\n",
    "  between3vs1=rbind(between3vs1,t)\n",
    "}\n",
    "\n",
    "grandvar2vs1 <- withinse2vs1^2 + ((1 + (1/5)) * between2vs1)\n",
    "grandse2vs1 <- sqrt(grandvar2vs1)\n",
    "#grandse2vs1\n",
    "\n",
    "grandvar3vs1 <- withinse3vs1^2 + ((1 + (1/5)) * between3vs1)\n",
    "grandse3vs1 <- sqrt(grandvar3vs1)\n",
    "#grandse3vs1\n",
    "\n",
    "z2vs1=grandm2vs1/grandse2vs1\n",
    "p2vs1=(1-pnorm(abs(z2vs1),0,1))*2\n",
    "p2vs1=round(p2vs1,3)\n",
    "\n",
    "z3vs1=grandm3vs1/grandse3vs1\n",
    "p3vs1=(1-pnorm(abs(z3vs1),0,1))*2\n",
    "p3vs1=round(p3vs1,3)\n",
    "\n",
    "rownames(p2vs1)=rownames(p3vs1)=colnames(beta1)\n",
    "colnames(p2vs1)=\"p 2vs1\"\n",
    "colnames(p3vs1)=\"p 3vs1\"\n",
    "or2vs1=round(exp(grandm2vs1),3)\n",
    "or3vs1=round(exp(grandm3vs1),3)\n",
    "library(pander)\n",
    "pander(cbind(or2vs1,or3vs1,p2vs1,p3vs1))\n",
    "\n",
    "# remove p>0.1 and repeat the process to get the reduced model\n",
    "\n",
    "# code not show to save space\n",
    "\n",
    "# calculate accuracy from imputed data\n",
    "\n",
    "accudata=rbind(mice::complete(imp),mice::complete(imp,2),mice::complete(imp,3),mice::complete(imp,4),mice::complete(imp,5))\n",
    "accudata1=select(accudata,X_traj_Group,anx,adltph,parity,n_cmdx,epdstota,gasscore00)\n",
    "accudata1=mutate(accudata1,b2x=3.40757026+0.44856672*as.numeric(as.character(anx))+0.62986931*ifelse(adltph==\"Yes\",1,0)+0.13116528*parity+0.15869229*n_cmdx+0.05023554*epdstota-0.08887876*gasscore00,b3x=2.0240581+0.6739823*as.numeric(as.character(anx))+0.3173449*ifelse(adltph==\"Yes\",1,0)+0.5495376*parity+0.1418330*n_cmdx+0.1513239*epdstota-0.1014388*gasscore00,p1=1/(1+exp(b2x)+exp(b3x)),p2=exp(b2x)/(1+exp(b2x)+exp(b3x)),p3=1-p1-p2)\n",
    "\n",
    "pred=apply(select(accudata1,p1,p2,p3),1,which.max)\n",
    "table(accudata1$X_traj_Group,pred)/5\n",
    "\n",
    "accuracy=(172.2+81.8+20.6)/507\n",
    "accuracy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Step 4 results:\n",
    "<img src=\"https://raw.githubusercontent.com/amysheep/NUTrajectoryPredictionAnalysis/master/fig5.png\">\n",
    "\n",
    "Overall accuracy improved, especially for group 3 (20.6/(7+6.4+20.6)=60.58%)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Other analyse (not show here):\n",
    "\n",
    "Used the group status as predictor and run linear mixed effect model to see how the three groups different in utilization of health services at month 0, 6 and 12. \n",
    "\n",
    "Women in the partial improvement trajectory were more likely to consult with their obstetrician or a mental health clinician earlier than women with chronic severe depression despite having similar initial severity. \n",
    "\n",
    "Women with postpartum depression prefer to receive treatment from their obstetrician, which may account for the greater number of visits to the obstetrician rather than mental health specialists. \n",
    "\n",
    "Women with chronic severe depression delayed contact until 12 months of worsening symptoms. This observation supports the integration of mental health into obstetrical care to improve access to care. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Built a simple Shiny app to calculate the predicted group for clinical use\n",
    "\n",
    "<img src=\"https://raw.githubusercontent.com/amysheep/NUTrajectoryPredictionAnalysis/master/fig6.png\" style=\"width:300px;height:500px;\">"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}