Part 2. Modelling the behavioral data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Part 2. Modelling the behavioral data\n",
    "\n",
    "In Part 1, we had a look at options for preprocessing fMRI data. Now, let's turn to the other essential ingredient of model-based fMRI-analyses: modelling behavioral data.\n",
    "\n",
    "There are Python packages that allows you to easily fit decision-making models such as the drift diffusion model (DDM). Most notably, [`hddm`](http://ski.clps.brown.edu/hddm_docs/) developed at Brown University is a widely used package.\n",
    "\n",
    "However, `EMC` (as used in the tutorials earlier this week) has some advantages. For example, `EMC` offers much more models (including a multitude of LBA-model variants, racing diffusion models) and it offers a better sampler. So you might want to use `EMC` to fit a model, and then use `Nipype` to create your fMRI-processing pipeline.\n",
    "\n",
    "So, we're just going to use `EMC` to fit the behavioral data for now. Note that _this_ notebook you're looking at right now (`Part 2. Modelling the behavioral data`) is actually _not_ an interactive Python notebook ('.ipynb'), but it is an interactive **R** notebook ('.irnb')! That means that instead of running Python code, we can run R-code in here.\n",
    "\n",
    "For example, to draw some random numbers from a normal distribution (4,1), we can use this R-function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rnorm(10, 4, 1)\n",
    "\n",
    "# note further that the help-function works as usual:\n",
    "# ?rnorm\n",
    "\n",
    "# this increases the plotting window size - this line only needs to be run once\n",
    "options(repr.plot.width = 18, repr.plot.height = 10)\n",
    "\n",
    "# and plotting is done in-line\n",
    "plot(1:10, 1:10, main='In line plot')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### So let's use EMC to fit the behavioral data\n",
    "\n",
    "First, set-up some directories, so R knows where to find EMC, the data, and where to store the results of these analyses (the 'res(ults)Dir')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## run this first once to set-up some thing for this server -- can forget about this\n",
    "library(RhpcBLASctl)\n",
    "omp_set_num_threads(1)\n",
    "blas_set_num_threads(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "projDir <- '/home/neuro/model_based_fmri'\n",
    "dataDir <- file.path('/data', 'bids')\n",
    "resDir <- file.path('/data', 'behavior_fits')\n",
    "dir.create(resDir, showWarnings=FALSE)\n",
    "\n",
    "# Load EMC\n",
    "library(EMC2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Load the data, and reformat to a EMC-approved form"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load data\n",
    "dat <- read.csv(file.path(dataDir, 'behavior.tsv'), sep='\\t')\n",
    "\n",
    "# Remove unnecessary columns, rename according to EMC's wishes\n",
    "dat <- dat[,c('pp', 'cond', 'stimulus', 'response', 'RT')]\n",
    "colnames(dat) <- c('subjects', 'E', 'S', 'R', 'rt')  # S = stimulus, E = emphasis, R = response, rt = response time\n",
    "\n",
    "# Remove all extremely fast responses (<.15s)\n",
    "dat <- dat[dat$rt>.15,]\n",
    "\n",
    "# Ensure the following columns are of type factor\n",
    "dat$subjects <- factor(dat$subjects)\n",
    "dat$E <- factor(dat$E, levels=c('spd', 'acc'))\n",
    "dat$S <- ifelse(dat$S=='stimright', 'right', 'left')\n",
    "dat$S <- factor(dat$S, levels=c('left', 'right'))\n",
    "dat$R <- factor(dat$R, levels=c('left', 'right'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### The next step is to set up the model specification"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_defective_density(dat,factors=c('E', \"S\"),layout=c(2,2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Emat <- matrix(c(-1,1),ncol=1)\n",
    "dimnames(Emat) <- list(NULL,\"acc-spd\")     # estimate the contrast between ACC-SPD trials\n",
    "\n",
    "# Average rate = intercept, and rate d = difference (match-mismatch) contrast\n",
    "ADmat <- matrix(c(-1/2,1/2),ncol=1,dimnames=list(NULL,\"d\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Next, we fit a series of models\n",
    "Unfortunately the behavioral data in this dataset is rather limited in trial numbers (only 80).\n",
    "As a result, models with more complexity than 1 parameter varying between SAT conditions, do not sample very well, rendering the resulting parameters rather unreliable.\n",
    "Hence, here we just do a simple model comparison that compares three models that allow B, v, and t0 to vary between conditions respectively.\n",
    "In the model that only allows B to vary between conditions, we further constrain `B_{spd}` to 0.05 and only estimate the contrast `B_Eacc-spd`; freely estimating `B_{spd}` also led to sampling issues.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# convenience function to fit a model and return the samples and posterior predictives\n",
    "fit_model <- function(data, design, fileName) {\n",
    "    samplers <- make_emc(data,design,type=\"standard\",rt_resolution=.02)\n",
    "    if(file.exists(fileName)) {\n",
    "        samplers <- EMC2:::loadRData(fileName)\n",
    "        return(samplers)\n",
    "    }\n",
    "    samplers <- fit(samplers, fileName=fileName, iter=1000, cores_per_chain=8, cores_for_chains=3, verbose=TRUE, verboseProgress=TRUE, max_tries = 20)\n",
    "}\n",
    "\n",
    "# Only B affected by E\n",
    "design_B <- design(\n",
    "  factors=list(subjects=levels(dat$subjects),S=levels(dat$S),E=levels(dat$E)),\n",
    "  Rlevels=levels(dat$R),matchfun=function(d)d$S==d$lR,\n",
    "  contrasts=list(lM=ADmat,lR=ADmat,S=ADmat,E=Emat),\n",
    "  formula=list(v~lM,sv~1,B~E,A~1,t0~1),       # <-- only B affected by E\n",
    "  constants=c(sv=log(1), B=log(0.05)),\n",
    "  model=LBA)\n",
    "fit_model(dat, design_B, \"samples/sJoN_B.RData\")\n",
    "\n",
    "# Only v affected by E\n",
    "design_v <- design(\n",
    "  factors=list(subjects=levels(dat$subjects),S=levels(dat$S),E=levels(dat$E)),\n",
    "  Rlevels=levels(dat$R),matchfun=function(d)d$S==d$lR,\n",
    "  contrasts=list(lM=ADmat,lR=ADmat,S=ADmat,E=Emat),\n",
    "  formula=list(v~E+lM,sv~1,B~1,A~1,t0~1),     # <-- only v affected by E\n",
    "  constants=c(sv=log(1)),\n",
    "  model=LBA)\n",
    "fit_model(dat, design_v, \"samples/sJoN_v.RData\")\n",
    "\n",
    "\n",
    "# only t0 affected by E\n",
    "design_t0 <- design(\n",
    "  factors=list(subjects=levels(dat$subjects),S=levels(dat$S),E=levels(dat$E)),\n",
    "  Rlevels=levels(dat$R),matchfun=function(d)d$S==d$lR,\n",
    "  contrasts=list(lM=ADmat,lR=ADmat,S=ADmat,E=Emat),\n",
    "  formula=list(v~lM,sv~1,B~1,A~1,t0~E),       # <-- only t0 affected by E\n",
    "  constants=c(sv=log(1)),\n",
    "  model=LBA)\n",
    "fit_model(dat, design_t0, \"samples/sJoN_t0.RData\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Sampling can take a bit of time, so I ran these models in advance, so you can just load the samplers & posterior predictives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sJoN_B <- EMC2:::loadRData('./samples/sJoN_B.RData')\n",
    "sJoN_v <- EMC2:::loadRData('./samples/sJoN_v.RData')\n",
    "sJoN_t0 <- EMC2:::loadRData('./samples/sJoN_t0.RData')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model comparison\n",
    "Let's start by doing a very basic model comparison based on information criteria (BPIC and DIC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "compare(list(B=sJoN_B,\n",
    "             v=sJoN_v,\n",
    "             t0=sJoN_t0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "Of these three models, the threshold-only model clearly wins both according to the DIC and BPIC."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Did sampling go well?\n",
    "Let's then assume the threshold-only model is the best model. What do the chains of the corresponding sampler look like?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width = 18, repr.plot.height = 20)  # here we adjust the plotting area in jupyter; you can adjust the width and height if this doesn't fit your screen well.\n",
    "plot(sJoN_B, selection='LL', layout=c(4,5))    # this actually plots; i.c. the posterior log likelihoods per subject"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The LL plots look like \"big hairy caterpillars\", as intended, so that's good. Next we can check the posteriors of the parameters at the subject level (alpha)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width = 18, repr.plot.height = 5)  # here we adjust the plotting area in jupyter; you can adjust the width and height if this doesn't fit your screen well.\n",
    "plot(sJoN_B, selection='alpha', layout=c(1,5), subject='197')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generally, `t0`, `v_lMd`, and `B_Eacc-spd` seem fine. `v` seems to sample a bit worse, the shapes suggest autocorrelations. What about the population level?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width = 18, repr.plot.height = 5)  # here we adjust the plotting area in jupyter; you can adjust the width and height if this doesn't fit your screen well.\n",
    "plot(sJoN_B, selection='mu', layout=c(1,5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That looks a bit similar for `v` as in the subject level. \n",
    "\n",
    "So let's check out how autocorrelated the chains arse. How many effective samples do we have?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "check(emc=sJoN_B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "so we have at least 350 effective samples for v (group-level); which isn't too much but good enough; and Rhat is all close to 1.01, so mixing is sufficient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fits\n",
    "Next, let's visualize how well this model actually fits the data. E.g., are there any obvious misfits? Does the model generally describe the data well? Can it capture the difference between SAT conditions?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ppJoN_B <- predict(sJoN_B, n_cores=19)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "options(repr.plot.width = 18, repr.plot.height = 10)                                    # change plotting area size a bit again\n",
    "plot_fit(dat,ppJoN_B,layout=c(2,2),factors=c(\"E\",\"S\"),lpos=\"right\",xlim=c(.25,1.5))     # aggregate over subjects and plot 'overall' CDFs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The overall (across-subject) fit seem to look reasonable (but not great). The model seems to capture the median/mean RTs and accuracy okay. Note that individual fits look a lot worse:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": []
   },
   "outputs": [],
   "source": [
    "options(repr.plot.width = 18, repr.plot.height = 5)                              #  \n",
    "plot_fit(dat, ppJoN_B, layout=c(1,2), factors=c(\"subjects\", \"E\", \"S\"))           # plot by subject"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given the limited amount of behavioral data, this is to be expected. Most importantly for the purposes of the current tutorial is that the SAT manipulation seems to be captured reasonably well by the threshold contrast `B_Eacc-spd`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# What about the parameters? \n",
    "Is `B_Eacc-spd` 'significant' at the population level? Let's first just visualize the posteriors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_pars(sJoN_B,layout=c(1,6),selection=\"mu\") \n",
    "# Note that the red lines are the priors, so we can see here that the posteriors really updated (ie we didn't just sample from the prior)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also calculate the proportion of density above/below 0 for the `B_Eacc-spd` parameter"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "credible(sJoN_B, \"B_Eacc-spd\", selection='mu', mu=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And get the Savage-Dickey ratio at 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hypothesis(sJoN_B,\"B_Eacc-spd\", selection='mu')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So pretty much the entire posterior lies above 0, indicating that there is an effect, and the Bayes Factor for that effect is extremely high."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# In sum\n",
    "we have a model that: \n",
    "- wins simple model comparison, \n",
    "- shows reasonably good sampling properties, \n",
    "- fits the behavioral data reasonably well, and \n",
    "- shows a between-condition contrast in the threshold parameter, captured by `B_Eacc-spd`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For the model-based fMRI analyses, we need the individual subject threshold-difference values (`B_Eacc-spd`), so let's extract these and save"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "parameters_df = parameters(sJoN_B, mapped=TRUE, selection=\"alpha\")\n",
    "\n",
    "median_parameters = aggregate(.~subjects, parameters_df, median)   # this calculates all parameters' median by subject\n",
    "median_parameters\n",
    "\n",
    "# and save\n",
    "write.csv(median_parameters, file=file.path(resDir, 'parameters_per_subject_lba_B_emc2_group2.csv'))"
   ]
  },
  {
   "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": "4.4.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}