invivo_Rstat.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shift function analysis\n",
    "\n",
    "<img src=\"https://upload.wikimedia.org/wikipedia/commons/thumb/1/1b/R_logo.svg/1280px-R_logo.svg.png\" width=\"50px\"/>\n",
    "\n",
    "To run this notebook locally, please make sure that you installed the following libraries by:\n",
    "\n",
    "```r\n",
    "install.packages(tibble)\n",
    "install.packages(png)\n",
    "install.packages(ggplot2)\n",
    "install.packages(R.matlab)\n",
    "install.packages(plotly)\n",
    "install.packages(htmlwidgets)\n",
    "install.packages(tidyverse)\n",
    "install.packages(plyr)\n",
    "install.packages(dplyr)\n",
    "install.packages(gridExtra)\n",
    "install.packages(repr)\n",
    "install.packages(PairedData)\n",
    "devtools::install_github('GRousselet/rogme')\n",
    "install.packages(\"tidyverse\")\n",
    "```\n",
    "\n",
    "Please also set the derivatives directory:\n",
    "\n",
    "```r\n",
    "derivativesDir <- \"/tmp/VENUS/qMRFlow/\"\n",
    "```\n",
    "\n",
    "```{admonition} Note\n",
    "The following code cell loads all the necessary R libraries and defines some helper functions.\n",
    "```\n",
    "\n",
    "```{warning} \n",
    "The data is expected to be found in the `/tmp/VENUS` directory.\n",
    "\n",
    "Please do not execute this notebook first. Python notebooks (`phantom_python` or `invivo_python`) uses repo2data to download the data on the BinderHub instance. Executing either one of them first will make the data available to this notebook as well.\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": [
     "hide-output",
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "R.matlab v3.6.2 (2018-09-26) successfully loaded. See ?R.matlab for help.\n",
      "\n",
      "\n",
      "Attaching package: ‘R.matlab’\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:base’:\n",
      "\n",
      "    getOption, isOpen\n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘plotly’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:ggplot2’:\n",
      "\n",
      "    last_plot\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:stats’:\n",
      "\n",
      "    filter\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:graphics’:\n",
      "\n",
      "    layout\n",
      "\n",
      "\n",
      "── \u001b[1mAttaching packages\u001b[22m ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.1 ──\n",
      "\n",
      "\u001b[32m✔\u001b[39m \u001b[34mtidyr  \u001b[39m 1.1.3     \u001b[32m✔\u001b[39m \u001b[34mdplyr  \u001b[39m 1.0.7\n",
      "\u001b[32m✔\u001b[39m \u001b[34mreadr  \u001b[39m 2.1.0     \u001b[32m✔\u001b[39m \u001b[34mstringr\u001b[39m 1.4.0\n",
      "\u001b[32m✔\u001b[39m \u001b[34mpurrr  \u001b[39m 0.3.4     \u001b[32m✔\u001b[39m \u001b[34mforcats\u001b[39m 0.5.1\n",
      "\n",
      "── \u001b[1mConflicts\u001b[22m ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──\n",
      "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mfilter()\u001b[39m masks \u001b[34mplotly\u001b[39m::filter(), \u001b[34mstats\u001b[39m::filter()\n",
      "\u001b[31m✖\u001b[39m \u001b[34mdplyr\u001b[39m::\u001b[32mlag()\u001b[39m    masks \u001b[34mstats\u001b[39m::lag()\n",
      "\n",
      "------------------------------------------------------------------------------\n",
      "\n",
      "You have loaded plyr after dplyr - this is likely to cause problems.\n",
      "If you need functions from both plyr and dplyr, please load plyr first, then dplyr:\n",
      "library(plyr); library(dplyr)\n",
      "\n",
      "------------------------------------------------------------------------------\n",
      "\n",
      "\n",
      "Attaching package: ‘plyr’\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:dplyr’:\n",
      "\n",
      "    arrange, count, desc, failwith, id, mutate, rename, summarise,\n",
      "    summarize\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:purrr’:\n",
      "\n",
      "    compact\n",
      "\n",
      "\n",
      "The following objects are masked from ‘package:plotly’:\n",
      "\n",
      "    arrange, mutate, rename, summarise\n",
      "\n",
      "\n",
      "\n",
      "Attaching package: ‘gridExtra’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:dplyr’:\n",
      "\n",
      "    combine\n",
      "\n",
      "\n",
      "Loading required package: MASS\n",
      "\n",
      "\n",
      "Attaching package: ‘MASS’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:dplyr’:\n",
      "\n",
      "    select\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:plotly’:\n",
      "\n",
      "    select\n",
      "\n",
      "\n",
      "Loading required package: gld\n",
      "\n",
      "Loading required package: mvtnorm\n",
      "\n",
      "Loading required package: lattice\n",
      "\n",
      "\n",
      "Attaching package: ‘PairedData’\n",
      "\n",
      "\n",
      "The following object is masked from ‘package:base’:\n",
      "\n",
      "    summary\n",
      "\n",
      "\n"
     ]
    }
   ],
   "source": [
    "options(warn = -1)\n",
    "library(rogme)\n",
    "library(tibble)\n",
    "library(ggplot2)\n",
    "library(R.matlab)\n",
    "library(plotly)\n",
    "library(htmlwidgets)\n",
    "library(tidyverse)\n",
    "library(plyr)\n",
    "library(dplyr)\n",
    "library(png)\n",
    "library(grid)\n",
    "library(gridExtra)\n",
    "library(repr)\n",
    "library(PairedData)\n",
    "set.seed(123)\n",
    "\n",
    "percentage_difference <- function(value, value_two) {\n",
    "  pctd = abs(value - value_two) / ((value + value_two) / 2) * 100\n",
    "  return(pctd)\n",
    "} \n",
    "\n",
    "obj2plotly <- function(curObj,acq,plim){\n",
    "\n",
    "if (acq==\"rth\"){\n",
    "color=\"#d62728\"    \n",
    "}else{\n",
    "color=\"#1f77b4\"\n",
    "}\n",
    "    \n",
    "p<-plot_hsf_pb(curObj, viridis_option=\"D\", ind_line_alpha=0.9, ind_line_size=1, gp_line_colour = \"red\", gp_point_colour=color,gp_point_size=4,gp_line_size=1.5) \n",
    "p$labels$colour = \"\"\n",
    "p$data$participants = revalue(p$data$participants,c(\"1\"=\"sub1\", \"2\"=\"sub2\",\"3\"=\"sub3\"))\n",
    "\n",
    "th <- ggplot2::theme_bw() + ggplot2::theme(\n",
    "  legend.title = element_blank(),  \n",
    "  axis.line=element_blank(),\n",
    "    legend.position = \"none\",\n",
    "  axis.ticks = element_blank(),\n",
    "  #axis.text = element_blank(),\n",
    "  axis.title = element_blank(),\n",
    "  panel.border = element_rect(colour = \"black\", fill=NA, size=0.5)\n",
    "  ) \n",
    "\n",
    "k <-p  + th + aes(x = dec, y = diff, colour = participants) + ggtitle(\"\")\n",
    "\n",
    "k$labels$colour = revalue(k$labels$colour,c(\"red\"=\"group\"))\n",
    "\n",
    "axx <- list(\n",
    "  zeroline = FALSE,\n",
    "  tickmode = \"array\",\n",
    "  tickvals = c(\"0.1\",\"0.2\",\"0.3\",\"0.4\",\"0.5\",\"0.6\",\"0.7\",\"0.8\",\"0.9\"),\n",
    "  ticktext = c(\"q1\",\"q2\",\"q3\",\"q4\",\"median\",\"q6\",\"q7\",\"q8\",\"q9\"),\n",
    "  showticklabels= TRUE,\n",
    "  ticksuffix=\" \"\n",
    ")\n",
    "    \n",
    "axy <- list(\n",
    "  zeroline = FALSE,   \n",
    "  range = c(-plim,plim),\n",
    "  overlaying = \"y2\",\n",
    "  showticklabels= FALSE\n",
    "  # The range is fixed per panel (3 MAD from both side) to infer the relative effect size   \n",
    ")\n",
    "\n",
    "axy2 <- list(\n",
    "  range = c(-plim,plim),\n",
    "    showticklabels= TRUE,\n",
    "    ticksuffix=\" \"\n",
    "  # The range is fixed per panel (3 MAD from both side) to infer the relative effect size   \n",
    ") \n",
    "figObj <- ggplotly(k) %>% \n",
    "            layout(height=800,width=800,xaxis=axx,yaxis = axy, yaxis2 = axy2, margin=c(5,5,5,5)) %>% # SWITCH BACK TO 400 x 400\n",
    "            style(traces=1,mode=\"lines+markers\",opacity=1,marker=list(size=10),line = list(shape = \"spline\",color=\"purple\",width=5)) %>%\n",
    "            style(traces=2,mode=\"lines+markers\",opacity=1,marker=list(size=10),line = list(shape = \"spline\",color=\"hotpink\",width=5)) %>% \n",
    "            style(traces=3,mode=\"lines+markers\",opacity=1,marker=list(size=10),line = list(shape = \"spline\",color=\"mediumvioletred\",width=5)) %>% \n",
    "            style(traces=4,opacity=1,line = list(color=\"black\",width=2,dash=\"dash\"),yaxis = 'y2') %>% # zero line\n",
    "            style(traces=5,mode = \"lines+markers\",opacity=0.9,line = list(shape = \"spline\",width=3,color=\"black\"),marker = list(color=color,size=30,line=list(color=\"black\",width=1))) %>% # group markers\n",
    "            style(traces=6,opacity=0,line = list(shape = \"spline\",color=\"#762a83\",width=3))\n",
    "return(figObj)    \n",
    "} #END FUNCTION DEF\n",
    "\n",
    "get_vector <- function(subject, session, metric, nSamples,derivativesDir){\n",
    "  derivativesDir <- paste(derivativesDir,subject ,\"/\",sep = \"\")  \n",
    "  curMat <- readMat(paste(derivativesDir,session,\"/stat/\",subject,\"_\",session,\"_desc-wm_metrics.mat\",sep = \"\"))\n",
    "  curMetric <- curMat[[metric]]\n",
    "  if (metric == \"MTR\"){\n",
    "      curMetric <- curMetric[!sapply(curMetric, is.nan)]\n",
    "      curMetric <- curMetric[curMetric > 35]\n",
    "      curMetric <- curMetric[curMetric < 70]\n",
    "  }else if (metric == \"MTsat\"){\n",
    "      curMetric <- curMetric[!sapply(curMetric, is.nan)]\n",
    "      curMetric <- curMetric[curMetric > 1]\n",
    "      curMetric <- curMetric[curMetric < 8]\n",
    "  }else if (metric == \"T1\"){\n",
    "      curMetric <- curMetric[!sapply(curMetric, is.nan)]\n",
    "      curMetric <- curMetric[curMetric > 0]\n",
    "      curMetric <- curMetric[curMetric < 3]    \n",
    "  } \n",
    "  curMetric <- sample (curMetric, size=nSamples, replace =F)\n",
    "  return(curMetric)  \n",
    "} #END FUNCTION DEF\n",
    "\n",
    "get_df <- function(ses1,ses2,metric,nSamples){\n",
    "df2 <- tibble(participant=as_factor(c(rep(1,nSamples*2),rep(2,nSamples*2),rep(3,nSamples*2))))\n",
    "df2[\"condition\"] = as_factor(rep(c(rep(ses1,nSamples),rep(ses2,nSamples)),3))\n",
    "curComp = c(get_vector(\"sub-invivo1\",ses1,metric,nSamples,derivativesDir),\n",
    "       get_vector(\"sub-invivo1\",ses2,metric,nSamples,derivativesDir),\n",
    "       get_vector(\"sub-invivo2\",ses1,metric,nSamples,derivativesDir),\n",
    "       get_vector(\"sub-invivo2\",ses2,metric,nSamples,derivativesDir),\n",
    "       get_vector(\"sub-invivo3\",ses1,metric,nSamples,derivativesDir),\n",
    "       get_vector(\"sub-invivo3\",ses2,metric,nSamples,derivativesDir))\n",
    "df2[\"metric\"] = curComp\n",
    "return(df2)\n",
    "} #END FUNCTION DEF\n",
    "\n",
    "display_sf <- function(derivativesDir,sub,scn1,scn2,metric,implementation){\n",
    "    \n",
    "    T1DIR = paste(derivativesDir,sub,\"/\",\"T1_SF_PLOTS\",sep=\"\")\n",
    "    MTRDIR = paste(derivativesDir,sub,\"/\",\"MTR_SF_PLOTS\",sep=\"\")\n",
    "    MTSDIR = paste(derivativesDir,sub,\"/\",\"MTS_SF_PLOTS\",sep=\"\")\n",
    "    ses_venus=data.frame(row.names=c(\"G1\",\"S1\",\"S2\") , val=c(\"rth750\",\"rthPRI\",\"rthSKY\"))\n",
    "    ses_native=data.frame(row.names=c(\"G1\",\"S1\",\"S2\") , val=c(\"vendor750\",\"vendorPRI\",\"vendorSKY\"))\n",
    "    \n",
    "    if (metric==\"T1\"){\n",
    "        curDir = T1DIR\n",
    "    }else if (metric==\"MTR\"){\n",
    "        curDir = MTRDIR\n",
    "    }else if (metric==\"MTsat\"){\n",
    "        curDir = MTSDIR\n",
    "    }\n",
    "\n",
    "    sel1 = ses_venus[scn1,]\n",
    "    sel2 = ses_venus[scn2,]\n",
    "    \n",
    "    file1 = paste(curDir,\"/\",sub,\"_\",sel1,\"vs\",sel2,\"_\",metric,\"_rev.png\",sep=\"\")\n",
    "    print(file1)\n",
    "    file2 = paste(curDir,\"/\",sub,\"_\",sel2,\"vs\",sel1,\"_\",metric,\"_rev.png\",sep=\"\")\n",
    "    \n",
    "    if (file.exists(file1)){\n",
    "        plot1 <- readPNG(file2)\n",
    "        plot2 <- readPNG(paste(curDir,\"/\",sub,\"_\",ses_native[scn1,],\"vs\",ses_native[scn2,],\"_\",metric,\"_rev.png\",sep=\"\"))\n",
    "        grid.arrange(rasterGrob(plot2),rasterGrob(plot1),nrow=1,top=textGrob(label=\"NATIVE (blue) vs VENUS (red)\",size=12))\n",
    "    }else if (file.exists(file2)){\n",
    "        plot1 <- readPNG(file2)\n",
    "        plot2 <- readPNG(paste(curDir,\"/\",sub,\"_\",ses_native[scn2,],\"vs\",ses_native[scn1,],\"_\",metric,\"_rev.png\",sep=\"\"))\n",
    "        grid.arrange(rasterGrob(plot2),rasterGrob(plot1),nrow=1,top=textGrob(\"NATIVE (blue) vs VENUS (red)\"))\n",
    "    }else {\n",
    "        print(\"No file\")\n",
    "    }\n",
    "    \n",
    "}\n",
    "\n",
    "display_hsf <- function(derivativesDir,scn1,scn2,metric){\n",
    "    \n",
    "    curDir = paste(derivativesDir,\"HSF\",sep=\"\")\n",
    "    ses_venus=data.frame(row.names=c(\"G1\",\"S1\",\"S2\") , val=c(\"rth750\",\"rthPRI\",\"rthSKY\"))\n",
    "    ses_native=data.frame(row.names=c(\"G1\",\"S1\",\"S2\") , val=c(\"vendor750\",\"vendorPRI\",\"vendorSKY\"))\n",
    "\n",
    "    sel1 = ses_venus[scn1,]\n",
    "    sel2 = ses_venus[scn2,]\n",
    "    \n",
    "    file1 = paste(curDir,\"/\",sel1,\"vs\",sel2,\"_HSF_\",metric,\".png\",sep=\"\")\n",
    "    print(file1)\n",
    "    file2 = paste(curDir,\"/\",sel2,\"vs\",sel1,\"_HSF_\",metric,\".png\",sep=\"\")\n",
    "    \n",
    "    if (file.exists(file1)){\n",
    "        plot1 <- readPNG(file1)\n",
    "        plot2 <- readPNG( paste(curDir,\"/\",ses_native[scn1,],\"vs\",ses_native[scn2,],\"_HSF_\",metric,\".png\",sep=\"\"))\n",
    "        plot3 <- readPNG( paste(curDir,\"/\",ses_venus[scn1,],\"vs\",ses_venus[scn2,],\"_HSFdif_\",metric,\".png\",sep=\"\"))\n",
    "        plot4 <- readPNG( paste(curDir,\"/\",ses_native[scn1,],\"vs\",ses_native[scn2,],\"_HSFdif_\",metric,\".png\",sep=\"\"))\n",
    "        grid.arrange(rasterGrob(plot2),rasterGrob(plot1),rasterGrob(plot4),rasterGrob(plot3),nrow=2,ncol=2,top=textGrob(\"NATIVE (blue) vs VENUS (red)\"))\n",
    "    }else if (file.exists(file2)){\n",
    "        plot1 <- readPNG(file2)\n",
    "        plot2 <- readPNG( paste(curDir,\"/\",ses_native[scn2,],\"vs\",ses_native[scn1,],\"_HSF_\",metric,\".png\",sep=\"\"))\n",
    "        grid.arrange(rasterGrob(plot2),rasterGrob(plot1),nrow=1,top=textGrob(\"NATIVE (blue) vs VENUS (red)\"))\n",
    "    }else {\n",
    "        print(\"No file\")\n",
    "    }\n",
    "    \n",
    "}\n",
    "\n",
    "\n",
    "# derivativesDir <- \"/tmp/VENUS/qMRFlow/\"\n",
    "\n",
    "# subs = list(\"sub-invivo1\",\"sub-invivo2\",\"sub-invivo3\")\n",
    "metrics = list(\"T1\",\"MTR\",\"MTsat\")\n",
    "# sessions <- list(\"ses-rth750rev\",\"ses-rthPRIrev\",\"ses-rthSKYrev\",\"ses-vendor750rev\",\"ses-vendorPRIrev\",\"ses-vendorSKYrev\")\n",
    "\n",
    "options(repr.plot.width=16, repr.plot.height=12)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## `T1`, `MTR` and `MTSat` shift function analysis\n",
    "\n",
    "`````{admonition} Note\n",
    ":class: tip\n",
    "Reproduces Figure 6 from the [article](https://www.biorxiv.org/content/10.1101/2021.12.27.474259v2) for the Participant-3 and includes supplementary analyses for the remaining participants.\n",
    "`````\n",
    "\n",
    "The following code cell performs a percentile bootstrapped comparison of dependent samples using Harrell-Davis quantile estimator. Outputs will be written as `png` and `csv` files to the following directories:\n",
    "\n",
    "* `/tmp/VENUS/qMRFLow/sub-invivo#/T1_SF_PLOTS` for T1 comparisons\n",
    "* `/tmp/VENUS/qMRFLow/sub-invivo#/MTR_SF_PLOTS` for MTR comparisons\n",
    "* `/tmp/VENUS/qMRFLow/sub-invivo#/MTS_SF_PLOTS` for MTsat comparisons\n",
    "\n",
    "```{warning}\n",
    "Executing the following code block takes about 10 - 15 minutes per metric! \n",
    "\n",
    "The outputs are already available in the directories listed above. You can skip the following cell and explore the outputs in the executable Jupyter Lab interface (requires starting a Binder session).\n",
    "\n",
    "To re-generate the outputs, please un-comment the code block (`ctrl` or `cmd` + `/`).\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 43404\n",
      "[1] \"rthSKYrev\"\n",
      "[1] 40294\n",
      "[1] \"rthSKYxa30\"\n",
      "[1] 47699\n",
      "[1] \"vendorSKYrev\"\n",
      "[1] 46346\n",
      "[1] \"vendorSKYxa30\"\n",
      "[1] 5.535536\n",
      "[1] 8.95549\n",
      "[1] 10.72789\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Scale for 'y' is already present. Adding another scale for 'y', which will\n",
      "replace the existing scale.\n",
      "\n",
      "Saving 178 x 178 mm image\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 4.312201\n",
      "[1] 7.570094\n",
      "[1] 7.751622\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Scale for 'y' is already present. Adding another scale for 'y', which will\n",
      "replace the existing scale.\n",
      "\n",
      "Saving 178 x 178 mm image\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 40155\n",
      "[1] \"rthSKYrev\"\n",
      "[1] 37682\n",
      "[1] \"rthSKYxa30\"\n",
      "[1] 43292\n",
      "[1] \"vendorSKYrev\"\n",
      "[1] 42334\n",
      "[1] \"vendorSKYxa30\"\n",
      "[1] 0.3499895\n",
      "[1] 0.6529195\n",
      "[1] 1.637242\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Scale for 'y' is already present. Adding another scale for 'y', which will\n",
      "replace the existing scale.\n",
      "\n",
      "Saving 178 x 178 mm image\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 0.02361323\n",
      "[1] 0.92132\n",
      "[1] 1.566945\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Scale for 'y' is already present. Adding another scale for 'y', which will\n",
      "replace the existing scale.\n",
      "\n",
      "Saving 178 x 178 mm image\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 41516\n",
      "[1] \"rthSKYrev\"\n",
      "[1] 38468\n",
      "[1] \"rthSKYxa30\"\n",
      "[1] 45849\n",
      "[1] \"vendorSKYrev\"\n",
      "[1] 44387\n",
      "[1] \"vendorSKYxa30\"\n",
      "[1] 6.36797\n",
      "[1] 8.511569\n",
      "[1] 9.524855\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Scale for 'y' is already present. Adding another scale for 'y', which will\n",
      "replace the existing scale.\n",
      "\n",
      "Saving 178 x 178 mm image\n",
      "\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] 1.281373\n",
      "[1] 3.326174\n",
      "[1] 4.020143\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Scale for 'y' is already present. Adding another scale for 'y', which will\n",
      "replace the existing scale.\n",
      "\n",
      "Saving 178 x 178 mm image\n",
      "\n"
     ]
    }
   ],
   "source": [
    "subs <- list('sub-invivo1')\n",
    "data_path <- \"/Users/agah/Desktop/neuropoly/ismrm23/ds-sysupgrade\"\n",
    "derivativesDir <- paste(data_path,\"/derivatives/\",\"qMRFlow/\",sep=\"\")\n",
    "\n",
    "metrics = list(\"T1\",\"MTR\",\"MTsat\")\n",
    "for (metric in metrics){\n",
    "for (sub in subs) {\n",
    "\n",
    "T1DIR = paste(derivativesDir,sub,\"/\",\"T1_SF_PLOTS\",sep=\"\")\n",
    "MTRDIR = paste(derivativesDir,sub,\"/\",\"MTR_SF_PLOTS\",sep=\"\")\n",
    "MTSDIR = paste(derivativesDir,sub,\"/\",\"MTS_SF_PLOTS\",sep=\"\")\n",
    "dir.create(T1DIR)\n",
    "dir.create(MTRDIR)\n",
    "dir.create(MTSDIR)\n",
    "\n",
    "sessions = list('ses-rthSKYrev','ses-rthSKYxa30_acq-retest','ses-vendorSKYrev','ses-vendorSKYxa30_acq-retest')\n",
    "sessions2 = list('ses-rthSKYrev','ses-rthSKYxa30','ses-vendorSKYrev','ses-vendorSKYxa30')\n",
    "sessions3 = list('rthSKYrev','rthSKYxa30','vendorSKYrev','vendorSKYxa30')\n",
    "#sessions <- list(\"ses-rth750rev\",\"ses-rthPRIrev\",\"ses-rthSKYrev\",\"ses-vendor750rev\",\"ses-vendorPRIrev\",\"ses-vendorSKYrev\")\n",
    "\n",
    "# Create dataframe with 37k samples per session\n",
    "nSamples = 37000\n",
    "df <- tibble(N = 1:nSamples)\n",
    "for (i in seq(from=1, to=length(sessions), by=1)) {\n",
    "  curMat <- readMat(paste(derivativesDir,sub,\"/\",sessions2[i],\"/stat/\",sub,\"_\",sessions[i],\"_desc-gm_metrics.mat\",sep = \"\"))\n",
    "  curMetric <- curMat[[metric]]\n",
    "  if (metric == \"MTR\"){\n",
    "      OUTDIR = MTRDIR\n",
    "      PLTRANGE = 10.1 \n",
    "      curMetric <- curMetric[curMetric > 35]\n",
    "      curMetric <- curMetric[curMetric < 70]\n",
    "  }else if (metric == \"MTsat\"){\n",
    "      OUTDIR = MTSDIR\n",
    "      PLTRANGE = 1.7\n",
    "      curMetric <- curMetric[curMetric > 1]\n",
    "      curMetric <- curMetric[curMetric < 8]\n",
    "  }else if (metric == \"T1\"){\n",
    "      OUTDIR = T1DIR\n",
    "      PLTRANGE = 0.5\n",
    "      curMetric <- curMetric[curMetric > 0]\n",
    "      curMetric <- curMetric[curMetric < 3]    \n",
    "  }\n",
    "  # Do not choose the same element more than once\n",
    "  print(length(curMetric))\n",
    "  curMetric <- sample (curMetric, size=nSamples, replace =F)          \n",
    "  varName <- toString(sessions3[i])\n",
    "  #varName <- substr(varName,5,nchar(varName)-3)\n",
    "  print(varName)\n",
    "  #df <- mutate(df, !!varName := curMetric)\n",
    "  df[varName] = curMetric  \n",
    "}\n",
    "             \n",
    "neutral = list(\"rth750\",\"rthPRI\",\"rthSKY\")\n",
    "native  = list(\"vendor750\",\"vendorPRI\",\"vendorSKY\")\n",
    "combns <- combn(c(1,2,3),2)\n",
    "\n",
    "for (j in seq(from=1,to=2, by=1)){\n",
    "\n",
    "if (j==1){\n",
    "    sel1 = 'rthSKYrev'\n",
    "    sel2 = 'rthSKYxa30'\n",
    "    }\n",
    "else\n",
    "{\n",
    "    sel1 = 'vendorSKYrev'\n",
    "    sel2 = 'vendorSKYxa30'\n",
    "}\n",
    "\n",
    "paste(sel1,\"vs\",sel2,\" T1\",' in progress',sep=\"\")\n",
    "    \n",
    "grp1 <- df[[sel1]]\n",
    "grp2 <- df[[sel2]]\n",
    "\n",
    "sfInput <- mkt2(grp1,grp2,group_labels=c(sel2,sel1))\n",
    "sf <- shiftdhd_pbci(\n",
    "  data = sfInput,\n",
    "  formula = obs ~ gr,\n",
    "  nboot = 250,\n",
    "  alpha = 0.05/10\n",
    ")\n",
    "\n",
    "sf[[1]]['pctdif'] = percentage_difference(sf[[1]][[sel1]],sf[[1]][[sel2]])\n",
    "write.csv(sf,paste(OUTDIR,\"/\",sub,\"_\",sel1,\"vs\",sel2,\"_\",metric,\"_GM\",'.csv',sep=\"\"))\n",
    "print(min(percentage_difference(sf[[1]][[sel1]],sf[[1]][[sel2]])))\n",
    "print(percentage_difference(sf[[1]][[sel1]],sf[[1]][[sel2]])[5])\n",
    "print(max(percentage_difference(sf[[1]][[sel1]],sf[[1]][[sel2]])))\n",
    "\n",
    "th <- ggplot2::theme_dark() + ggplot2::theme(\n",
    "  panel.border = element_blank(),\n",
    "  panel.background = element_rect(fill = \"black\"),\n",
    "  plot.background = element_rect(fill = \"black\"),\n",
    "  legend.background = element_rect(fill=\"white\",size=0.5),  \n",
    "  legend.text=element_text(color=\"white\",size=12),    \n",
    "  axis.line=element_blank(),\n",
    "  axis.text = element_blank(),\n",
    "  axis.title = element_blank(),\n",
    "  axis.text.x = element_text(colour=\"white\",size=18,face='bold'),\n",
    "  axis.text.y = element_text(colour=\"white\",size=18,face='bold')\n",
    "  ) \n",
    "\n",
    "if (j==1){\n",
    "    lnclr = \"darkgray\"\n",
    "    flclr = \"red\"\n",
    "} else\n",
    "{\n",
    "    lnclr = \"darkgray\"\n",
    "    flclr = \"blue\"\n",
    "}\n",
    "    \n",
    "psf <- plot_sf(\n",
    "  data = sf,\n",
    "  plot_theme = 1,\n",
    "  symb_col = \"black\",\n",
    "  symb_fill = flclr,\n",
    "  symb_size = 8,\n",
    "  symb_shape = 21,\n",
    "  diffint_col = \"white\",\n",
    "  diffint_size = 1.5,\n",
    "  q_line_col = \"black\",\n",
    "  q_line_alpha = 0.8,\n",
    "  q_line_size = 1.8,\n",
    "  theme2_alpha = NULL\n",
    ") \n",
    "\n",
    "curPlot <- psf[[1]] + geom_hline(yintercept = 0,colour=lnclr, linetype = \"longdash\",size=1) + ylim(-PLTRANGE,PLTRANGE)\n",
    "ggsave(\n",
    "  filename = paste(OUTDIR,\"/\",sub,\"_\",sel1,\"vs\",sel2,\"_\",metric,\"_GM\",'.png',sep=\"\"),\n",
    "  plot = curPlot,\n",
    "  device = NULL,\n",
    "  path = NULL,\n",
    "  scale = 1,\n",
    "  width = NA,\n",
    "  height = NA,\n",
    "  units = c(\"mm\"),\n",
    "  dpi = 300,\n",
    "  limitsize = TRUE\n",
    ")\n",
    "}}\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "df <- mutate(df, !!varName := curMetric)  \n",
    "df[varName] = curMetric"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table class=\"dataframe\">\n",
       "<caption>A tibble: 74000 × 2</caption>\n",
       "<thead>\n",
       "\t<tr><th scope=col>gr</th><th scope=col>obs</th></tr>\n",
       "\t<tr><th scope=col>&lt;fct&gt;</th><th scope=col>&lt;dbl&gt;</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "\t<tr><td>vendorSKYrev</td><td>54.25170</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>53.68619</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>54.19176</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>54.55876</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>54.16074</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>53.06513</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>49.47857</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>52.58469</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>51.25350</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>55.64699</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>53.03757</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>55.15184</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>53.14215</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>53.65039</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>52.94985</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>46.17459</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>52.66718</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>53.02032</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>53.09795</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>55.91122</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>56.12119</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>54.44606</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>50.42297</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>53.39558</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>53.36021</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>52.37278</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>55.23981</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>52.04903</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>53.90504</td></tr>\n",
       "\t<tr><td>vendorSKYrev</td><td>55.23962</td></tr>\n",
       "\t<tr><td>⋮</td><td>⋮</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>46.96217</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>48.95001</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>54.44372</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>53.15768</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>53.61669</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>56.04972</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>52.01628</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>55.18143</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>50.95692</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>52.49494</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>58.77942</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>57.86157</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>53.78435</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>57.31303</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>51.61267</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>55.87280</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>54.57355</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>52.56245</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>52.01884</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>53.09877</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>51.38911</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>56.87275</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>55.88947</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>48.54071</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>50.94834</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>53.36089</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>51.46187</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>54.62327</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>56.09983</td></tr>\n",
       "\t<tr><td>vendorSKYxa30</td><td>55.15531</td></tr>\n",
       "</tbody>\n",
       "</table>\n"
      ],
      "text/latex": [
       "A tibble: 74000 × 2\n",
       "\\begin{tabular}{ll}\n",
       " gr & obs\\\\\n",
       " <fct> & <dbl>\\\\\n",
       "\\hline\n",
       "\t vendorSKYrev & 54.25170\\\\\n",
       "\t vendorSKYrev & 53.68619\\\\\n",
       "\t vendorSKYrev & 54.19176\\\\\n",
       "\t vendorSKYrev & 54.55876\\\\\n",
       "\t vendorSKYrev & 54.16074\\\\\n",
       "\t vendorSKYrev & 53.06513\\\\\n",
       "\t vendorSKYrev & 49.47857\\\\\n",
       "\t vendorSKYrev & 52.58469\\\\\n",
       "\t vendorSKYrev & 51.25350\\\\\n",
       "\t vendorSKYrev & 55.64699\\\\\n",
       "\t vendorSKYrev & 53.03757\\\\\n",
       "\t vendorSKYrev & 55.15184\\\\\n",
       "\t vendorSKYrev & 53.14215\\\\\n",
       "\t vendorSKYrev & 53.65039\\\\\n",
       "\t vendorSKYrev & 52.94985\\\\\n",
       "\t vendorSKYrev & 46.17459\\\\\n",
       "\t vendorSKYrev & 52.66718\\\\\n",
       "\t vendorSKYrev & 53.02032\\\\\n",
       "\t vendorSKYrev & 53.09795\\\\\n",
       "\t vendorSKYrev & 55.91122\\\\\n",
       "\t vendorSKYrev & 56.12119\\\\\n",
       "\t vendorSKYrev & 54.44606\\\\\n",
       "\t vendorSKYrev & 50.42297\\\\\n",
       "\t vendorSKYrev & 53.39558\\\\\n",
       "\t vendorSKYrev & 53.36021\\\\\n",
       "\t vendorSKYrev & 52.37278\\\\\n",
       "\t vendorSKYrev & 55.23981\\\\\n",
       "\t vendorSKYrev & 52.04903\\\\\n",
       "\t vendorSKYrev & 53.90504\\\\\n",
       "\t vendorSKYrev & 55.23962\\\\\n",
       "\t ⋮ & ⋮\\\\\n",
       "\t vendorSKYxa30 & 46.96217\\\\\n",
       "\t vendorSKYxa30 & 48.95001\\\\\n",
       "\t vendorSKYxa30 & 54.44372\\\\\n",
       "\t vendorSKYxa30 & 53.15768\\\\\n",
       "\t vendorSKYxa30 & 53.61669\\\\\n",
       "\t vendorSKYxa30 & 56.04972\\\\\n",
       "\t vendorSKYxa30 & 52.01628\\\\\n",
       "\t vendorSKYxa30 & 55.18143\\\\\n",
       "\t vendorSKYxa30 & 50.95692\\\\\n",
       "\t vendorSKYxa30 & 52.49494\\\\\n",
       "\t vendorSKYxa30 & 58.77942\\\\\n",
       "\t vendorSKYxa30 & 57.86157\\\\\n",
       "\t vendorSKYxa30 & 53.78435\\\\\n",
       "\t vendorSKYxa30 & 57.31303\\\\\n",
       "\t vendorSKYxa30 & 51.61267\\\\\n",
       "\t vendorSKYxa30 & 55.87280\\\\\n",
       "\t vendorSKYxa30 & 54.57355\\\\\n",
       "\t vendorSKYxa30 & 52.56245\\\\\n",
       "\t vendorSKYxa30 & 52.01884\\\\\n",
       "\t vendorSKYxa30 & 53.09877\\\\\n",
       "\t vendorSKYxa30 & 51.38911\\\\\n",
       "\t vendorSKYxa30 & 56.87275\\\\\n",
       "\t vendorSKYxa30 & 55.88947\\\\\n",
       "\t vendorSKYxa30 & 48.54071\\\\\n",
       "\t vendorSKYxa30 & 50.94834\\\\\n",
       "\t vendorSKYxa30 & 53.36089\\\\\n",
       "\t vendorSKYxa30 & 51.46187\\\\\n",
       "\t vendorSKYxa30 & 54.62327\\\\\n",
       "\t vendorSKYxa30 & 56.09983\\\\\n",
       "\t vendorSKYxa30 & 55.15531\\\\\n",
       "\\end{tabular}\n"
      ],
      "text/markdown": [
       "\n",
       "A tibble: 74000 × 2\n",
       "\n",
       "| gr &lt;fct&gt; | obs &lt;dbl&gt; |\n",
       "|---|---|\n",
       "| vendorSKYrev | 54.25170 |\n",
       "| vendorSKYrev | 53.68619 |\n",
       "| vendorSKYrev | 54.19176 |\n",
       "| vendorSKYrev | 54.55876 |\n",
       "| vendorSKYrev | 54.16074 |\n",
       "| vendorSKYrev | 53.06513 |\n",
       "| vendorSKYrev | 49.47857 |\n",
       "| vendorSKYrev | 52.58469 |\n",
       "| vendorSKYrev | 51.25350 |\n",
       "| vendorSKYrev | 55.64699 |\n",
       "| vendorSKYrev | 53.03757 |\n",
       "| vendorSKYrev | 55.15184 |\n",
       "| vendorSKYrev | 53.14215 |\n",
       "| vendorSKYrev | 53.65039 |\n",
       "| vendorSKYrev | 52.94985 |\n",
       "| vendorSKYrev | 46.17459 |\n",
       "| vendorSKYrev | 52.66718 |\n",
       "| vendorSKYrev | 53.02032 |\n",
       "| vendorSKYrev | 53.09795 |\n",
       "| vendorSKYrev | 55.91122 |\n",
       "| vendorSKYrev | 56.12119 |\n",
       "| vendorSKYrev | 54.44606 |\n",
       "| vendorSKYrev | 50.42297 |\n",
       "| vendorSKYrev | 53.39558 |\n",
       "| vendorSKYrev | 53.36021 |\n",
       "| vendorSKYrev | 52.37278 |\n",
       "| vendorSKYrev | 55.23981 |\n",
       "| vendorSKYrev | 52.04903 |\n",
       "| vendorSKYrev | 53.90504 |\n",
       "| vendorSKYrev | 55.23962 |\n",
       "| ⋮ | ⋮ |\n",
       "| vendorSKYxa30 | 46.96217 |\n",
       "| vendorSKYxa30 | 48.95001 |\n",
       "| vendorSKYxa30 | 54.44372 |\n",
       "| vendorSKYxa30 | 53.15768 |\n",
       "| vendorSKYxa30 | 53.61669 |\n",
       "| vendorSKYxa30 | 56.04972 |\n",
       "| vendorSKYxa30 | 52.01628 |\n",
       "| vendorSKYxa30 | 55.18143 |\n",
       "| vendorSKYxa30 | 50.95692 |\n",
       "| vendorSKYxa30 | 52.49494 |\n",
       "| vendorSKYxa30 | 58.77942 |\n",
       "| vendorSKYxa30 | 57.86157 |\n",
       "| vendorSKYxa30 | 53.78435 |\n",
       "| vendorSKYxa30 | 57.31303 |\n",
       "| vendorSKYxa30 | 51.61267 |\n",
       "| vendorSKYxa30 | 55.87280 |\n",
       "| vendorSKYxa30 | 54.57355 |\n",
       "| vendorSKYxa30 | 52.56245 |\n",
       "| vendorSKYxa30 | 52.01884 |\n",
       "| vendorSKYxa30 | 53.09877 |\n",
       "| vendorSKYxa30 | 51.38911 |\n",
       "| vendorSKYxa30 | 56.87275 |\n",
       "| vendorSKYxa30 | 55.88947 |\n",
       "| vendorSKYxa30 | 48.54071 |\n",
       "| vendorSKYxa30 | 50.94834 |\n",
       "| vendorSKYxa30 | 53.36089 |\n",
       "| vendorSKYxa30 | 51.46187 |\n",
       "| vendorSKYxa30 | 54.62327 |\n",
       "| vendorSKYxa30 | 56.09983 |\n",
       "| vendorSKYxa30 | 55.15531 |\n",
       "\n"
      ],
      "text/plain": [
       "      gr            obs     \n",
       "1     vendorSKYrev  54.25170\n",
       "2     vendorSKYrev  53.68619\n",
       "3     vendorSKYrev  54.19176\n",
       "4     vendorSKYrev  54.55876\n",
       "5     vendorSKYrev  54.16074\n",
       "6     vendorSKYrev  53.06513\n",
       "7     vendorSKYrev  49.47857\n",
       "8     vendorSKYrev  52.58469\n",
       "9     vendorSKYrev  51.25350\n",
       "10    vendorSKYrev  55.64699\n",
       "11    vendorSKYrev  53.03757\n",
       "12    vendorSKYrev  55.15184\n",
       "13    vendorSKYrev  53.14215\n",
       "14    vendorSKYrev  53.65039\n",
       "15    vendorSKYrev  52.94985\n",
       "16    vendorSKYrev  46.17459\n",
       "17    vendorSKYrev  52.66718\n",
       "18    vendorSKYrev  53.02032\n",
       "19    vendorSKYrev  53.09795\n",
       "20    vendorSKYrev  55.91122\n",
       "21    vendorSKYrev  56.12119\n",
       "22    vendorSKYrev  54.44606\n",
       "23    vendorSKYrev  50.42297\n",
       "24    vendorSKYrev  53.39558\n",
       "25    vendorSKYrev  53.36021\n",
       "26    vendorSKYrev  52.37278\n",
       "27    vendorSKYrev  55.23981\n",
       "28    vendorSKYrev  52.04903\n",
       "29    vendorSKYrev  53.90504\n",
       "30    vendorSKYrev  55.23962\n",
       "⋮     ⋮             ⋮       \n",
       "73971 vendorSKYxa30 46.96217\n",
       "73972 vendorSKYxa30 48.95001\n",
       "73973 vendorSKYxa30 54.44372\n",
       "73974 vendorSKYxa30 53.15768\n",
       "73975 vendorSKYxa30 53.61669\n",
       "73976 vendorSKYxa30 56.04972\n",
       "73977 vendorSKYxa30 52.01628\n",
       "73978 vendorSKYxa30 55.18143\n",
       "73979 vendorSKYxa30 50.95692\n",
       "73980 vendorSKYxa30 52.49494\n",
       "73981 vendorSKYxa30 58.77942\n",
       "73982 vendorSKYxa30 57.86157\n",
       "73983 vendorSKYxa30 53.78435\n",
       "73984 vendorSKYxa30 57.31303\n",
       "73985 vendorSKYxa30 51.61267\n",
       "73986 vendorSKYxa30 55.87280\n",
       "73987 vendorSKYxa30 54.57355\n",
       "73988 vendorSKYxa30 52.56245\n",
       "73989 vendorSKYxa30 52.01884\n",
       "73990 vendorSKYxa30 53.09877\n",
       "73991 vendorSKYxa30 51.38911\n",
       "73992 vendorSKYxa30 56.87275\n",
       "73993 vendorSKYxa30 55.88947\n",
       "73994 vendorSKYxa30 48.54071\n",
       "73995 vendorSKYxa30 50.94834\n",
       "73996 vendorSKYxa30 53.36089\n",
       "73997 vendorSKYxa30 51.46187\n",
       "73998 vendorSKYxa30 54.62327\n",
       "73999 vendorSKYxa30 56.09983\n",
       "74000 vendorSKYxa30 55.15531"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sfInput"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "'rthSKYrev'"
      ],
      "text/latex": [
       "'rthSKYrev'"
      ],
      "text/markdown": [
       "'rthSKYrev'"
      ],
      "text/plain": [
       "[1] \"rthSKYrev\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "neutral = list(\"rthSKYxa30\",\"vendorSKYxa30\")\n",
    "native  = list(\"rthSKYrev\",\"vendorSKYrev\")\n",
    "\n",
    "combns <- combn(c(1,2,3),2)\n",
    "\n",
    "for (j in seq(from=1,to=2, by=1)){\n",
    "for (i in seq(from=1,to=2, by=1)){\n",
    "\n",
    "if (j==1){\n",
    "    sel1 = toString(neutral[combns[,i][1]])\n",
    "    sel2 = toString(neutral[combns[,i][2]])\n",
    "    }\n",
    "else\n",
    "{\n",
    "    sel1 = toString(native[combns[,i][1]])\n",
    "    sel2 = toString(native[combns[,i][2]])\n",
    "}\n",
    "\n",
    "paste(sel1,\"vs\",sel2,\" T1\",' in progress',sep=\"\")\n",
    "}}\n",
    "\n",
    "sel1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explore shift function plots\n",
    "\n",
    "```{admonition} Tip\n",
    ":class: tip\n",
    "You can start an interactive BinderHub session to explore VENUS (red) vs NATIVE (blue) shift functions side by side for a selected participant, metric and scanner pairs. \n",
    "\n",
    "* Participants\n",
    "    * `sub-invivo1`\n",
    "    * `sub-invivo2`\n",
    "    * `sub-invivo3`    \n",
    "* Metrics\n",
    "    * `T1`\n",
    "    * `MTR`\n",
    "    * `MTsat`\n",
    "* Scanners\n",
    "    * `G1`\n",
    "    * `S1`\n",
    "    * `S2`\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "'/Users/agah/Desktop/neuropoly/ismrm23/ds-sysupgrade /derivatives/ qMRFlow/sub-invivo1/ses-rthSKYrev/stat/sub-invivo1_ses-rthSKYrev_desc-wm_metrics.mat'"
      ],
      "text/latex": [
       "'/Users/agah/Desktop/neuropoly/ismrm23/ds-sysupgrade /derivatives/ qMRFlow/sub-invivo1/ses-rthSKYrev/stat/sub-invivo1\\_ses-rthSKYrev\\_desc-wm\\_metrics.mat'"
      ],
      "text/markdown": [
       "'/Users/agah/Desktop/neuropoly/ismrm23/ds-sysupgrade /derivatives/ qMRFlow/sub-invivo1/ses-rthSKYrev/stat/sub-invivo1_ses-rthSKYrev_desc-wm_metrics.mat'"
      ],
      "text/plain": [
       "[1] \"/Users/agah/Desktop/neuropoly/ismrm23/ds-sysupgrade /derivatives/ qMRFlow/sub-invivo1/ses-rthSKYrev/stat/sub-invivo1_ses-rthSKYrev_desc-wm_metrics.mat\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "paste(derivativesDir,sub,\"/\",sessions2[i],\"/stat/\",sub,\"_\",sessions[i],\"_desc-wm_metrics.mat\",sep = \"\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "data_path <- \"/Users/agah/Desktop/neuropoly/ismrm23/ds-sysupgrade\"\n",
    "derivativesDir <- paste(data_path,\"/derivatives/\",\"qMRFlow/\")\n",
    "\n",
    "\n",
    "subject =\"sub-invivo1\" # or sub-invivo2 or sub-invivo3\n",
    "metric =\"T1\" # or MTR or MTsat\n",
    "scanner1 = \"S1\" # or S1 or S2\n",
    "scanner2 = \"S1\" # or G1 or S2\n",
    "\n",
    "\n",
    "display_sf(derivativesDir,subject,scanner1,scanner2,metric,implementation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Perform hierarchical shift function (HSF) analysis\n",
    "\n",
    "`````{admonition} Note\n",
    ":class: tip\n",
    "Reproduces Figure 7 from the [article](https://www.biorxiv.org/content/10.1101/2021.12.27.474259v2).\n",
    "`````\n",
    "\n",
    "The following code cell iterates over all the comparisons and saves interactive HSF figures.\n",
    "\n",
    "* `/tmp/VENUS/HSF` \n",
    "\n",
    "```{warning}\n",
    "Executing the following code block takes about 10 - 15 minutes! \n",
    "\n",
    "The (static) outputs are already available at `/tmp/VENUS/HSF`. You can skip the following cell and explore the outputs in the executable Jupyter Lab interface (requires starting a Binder session).\n",
    "\n",
    "To re-generate the outputs and to save interactive plots, please un-comment the code block (`ctrl` or `cmd` + `/`).\n",
    "\n",
    "You can see interactive visualization from an iteration using\n",
    "\n",
    "```r\n",
    "obj2plotly(out,clr,PLTRANGE)\n",
    "```\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "# HSFDIR = paste(derivativesDir,\"HSF/\",sep=\"\")\n",
    "# #HSFDIR = paste(\"HSF/\",sep=\"\")\n",
    "# dir.create(HSFDIR)\n",
    "\n",
    "# sessions <- list(\"ses-rth750rev\",\"ses-rthPRIrev\",\"ses-rthSKYrev\",\"ses-vendor750rev\",\"ses-vendorPRIrev\",\"ses-vendorSKYrev\")\n",
    "# nSamples = 37000\n",
    "\n",
    "# neutral = list(\"ses-rth750rev\",\"ses-rthPRIrev\",\"ses-rthSKYrev\")\n",
    "# native  = list(\"ses-vendor750rev\",\"ses-vendorPRIrev\",\"ses-vendorSKYrev\")\n",
    "# combns <- combn(c(1,2,3),2)\n",
    "\n",
    "# metrics = list(\"T1\",\"MTR\",\"MTsat\")\n",
    "# for (metric in metrics){\n",
    "\n",
    "# if (metric == \"MTR\"){\n",
    "#       PLTRANGE = 10.3 \n",
    "#   }else if (metric == \"MTsat\"){\n",
    "#       PLTRANGE = 1.8\n",
    "#   }else if (metric == \"T1\"){\n",
    "#       PLTRANGE = 0.6\n",
    "#   }\n",
    "# for (j in seq(from=1,to=2, by=1)){\n",
    "# for (i in seq(from=1,to=3, by=1)){\n",
    "\n",
    "# if (j==1){\n",
    "#     sel1 = toString(neutral[combns[,i][1]])\n",
    "#     sel2 = toString(neutral[combns[,i][2]])\n",
    "#     clr = \"rth\"\n",
    "#     }\n",
    "# else\n",
    "# {\n",
    "#     clr = \"vendor\"\n",
    "#     sel1 = toString(native[combns[,i][1]])\n",
    "#     sel2 = toString(native[combns[,i][2]])\n",
    "# }\n",
    "\n",
    "# print(paste(sel1,\"vs\",sel2,\" T1\",' in progress',sep=\"\"))\n",
    "# cur_df <-get_df(sel1,sel2,metric,nSamples)\n",
    "# out <- hsf_pb(\n",
    "#   data = cur_df,\n",
    "#   formula = metric ~ condition + participant,\n",
    "#   qseq = seq(0.1, 0.9, 0.1),\n",
    "#   tr = 0,\n",
    "#   alpha = 0.05,\n",
    "#   qtype = 8,\n",
    "#   todo = c(1, 2),\n",
    "#   nboot = 750)\n",
    "# htmlwidgets::saveWidget(obj2plotly(out,clr,PLTRANGE), file = paste(substr(sel1,5,nchar(sel1)-3),\"vs\",substr(sel2,5,nchar(sel2)-3),\"_HSF_\",metric,\".html\",sep=\"\"))\n",
    "# }}}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Explore HSF plots\n",
    "\n",
    "```{figure} ../assets/hsf_exp.jpg\n",
    "---\n",
    "height: 600px\n",
    "name: hsfexp\n",
    "---\n",
    "explanation of the hierarchical shift function analysis.\n",
    "```\n",
    "\n",
    "```{admonition} Tip\n",
    ":class: tip\n",
    "You can start an interactive BinderHub session to explore VENUS (red) vs NATIVE (blue) HSF side by side for a selected metric and scanner pairs. \n",
    "\n",
    "* Metrics\n",
    "    * `T1`\n",
    "    * `MTR`\n",
    "    * `MTsat`\n",
    "* Scanners\n",
    "    * `G1`\n",
    "    * `S1`\n",
    "    * `S2`\n",
    "\n",
    "In addition to the HSF plots, respective bootstrapped differences will also be displayed (second row).\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "metric =\"MTsat\" # or MTR or MTsat\n",
    "scanner1 = \"G1\" # or S1 or S2\n",
    "scanner2 = \"S1\" # or G1 or S2\n",
    "\n",
    "\n",
    "display_hsf(\"/tmp/VENUS/\",scanner1,scanner2,metric)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bootstrapped differences at each decile for HSF\n",
    "\n",
    "`````{admonition} Supplementary\n",
    ":class: tip\n",
    "Creates bootstrapped difference plots displayed in the above cell.\n",
    "`````"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "# HSFDIR = paste(derivativesDir,\"HSF/\",sep=\"\")\n",
    "# dir.create(HSFDIR)\n",
    "\n",
    "# sessions <- list(\"ses-rth750rev\",\"ses-rthPRIrev\",\"ses-rthSKYrev\",\"ses-vendor750rev\",\"ses-vendorPRIrev\",\"ses-vendorSKYrev\")\n",
    "# nSamples = 37000\n",
    "\n",
    "# neutral = list(\"ses-rth750rev\",\"ses-rthPRIrev\",\"ses-rthSKYrev\")\n",
    "# native  = list(\"ses-vendor750rev\",\"ses-vendorPRIrev\",\"ses-vendorSKYrev\")\n",
    "# combns <- combn(c(1,2,3),2)\n",
    "\n",
    "# metrics = list(\"T1\",\"MTR\",\"MTsat\")\n",
    "# for (metric in metrics){\n",
    "\n",
    "# if (metric == \"MTR\"){\n",
    "#       PLTRANGE = 10.3 \n",
    "#   }else if (metric == \"MTsat\"){\n",
    "#       PLTRANGE = 1.8\n",
    "#   }else if (metric == \"T1\"){\n",
    "#       PLTRANGE = 0.6\n",
    "#   }\n",
    "# for (j in seq(from=1,to=2, by=1)){\n",
    "# for (i in seq(from=1,to=3, by=1)){\n",
    "\n",
    "# if (j==1){\n",
    "#     sel1 = toString(neutral[combns[,i][1]])\n",
    "#     sel2 = toString(neutral[combns[,i][2]])\n",
    "#     clr = \"#d62728\"\n",
    "#     }\n",
    "# else\n",
    "# {\n",
    "#     clr = \"#1f77b4\"\n",
    "#     sel1 = toString(native[combns[,i][1]])\n",
    "#     sel2 = toString(native[combns[,i][2]])\n",
    "# }\n",
    "\n",
    "# print(paste(sel1,\"vs\",sel2,\" T1\",' in progress',sep=\"\"))\n",
    "    \n",
    "\n",
    "# cur_df <-get_df(sel1,sel2,metric,nSamples)\n",
    "# out <- hsf_pb(\n",
    "#   data = cur_df,\n",
    "#   formula = metric ~ condition + participant,\n",
    "#   todo = c(1, 2),\n",
    "#   nboot = 750)\n",
    "# obj <- plot_hsf_pb_dist(out,fill_colour=clr,point_interv=\"median\")\n",
    "# ggsave(\n",
    "#   filename = paste(substr(sel1,5,nchar(sel1)-3),\"vs\",substr(sel2,5,nchar(sel2)-3),\"_HSFdif_\",metric,\".png\",sep=\"\"),\n",
    "#   plot = obj,\n",
    "#   device = NULL,\n",
    "#   path = NULL,\n",
    "#   scale = 1,\n",
    "#   width = NA,\n",
    "#   height = NA,\n",
    "#   units = c(\"mm\"),\n",
    "#   dpi = 300,\n",
    "#   limitsize = TRUE\n",
    "# )\n",
    "# #orca(obj, file = paste(substr(sel1,5,nchar(sel1)-3),\"vs\",substr(sel2,5,nchar(sel2)-3),\"_HSFdif_\",metric,\".png\",sep=\"\"))    \n",
    "# }}}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Significance test\n",
    "\n",
    "Paired comparison of difference scores between VENUS and NATIVE implementations.\n",
    "\n",
    "`````{admonition} Note\n",
    ":class: tip\n",
    "The values are drawn from the `csv` files available at `/tmp/VENUS/qMRFlow/sub-invivo#/_SF_` to construct the data frame. \n",
    "`````"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "data <- data.frame(participant = rep(1:3, each=12),\n",
    "                   sequence = rep(c(\"native\",\"venus\"), times=18),\n",
    "                   metric = rep(c(\"T1\",\"T1\",\"T1\",\"T1\",\"MTR\",\"MTR\",\"MTR\", \"MTR\",\"MTsat\",\"MTsat\",\"MTsat\", \"MTsat\"),3),\n",
    "                   pctdif = c(31.2,10.9,16.30,1.8,17.7,8.3,16.9,3.6,14.5,6.7,25.7,3.2,\n",
    "                             20.9,6.6,10.8,3.7,15.5,6.4,14.8,8.3,17.7,6.7,25.2,3.2,\n",
    "                            14.2,0.1,14.5,3.6,14.7,6.3,16.3,2.2,21.9,10.3,26.1,0.6))\n",
    "\n",
    "data <- data.frame(participant = rep(1:3, each=12),\n",
    "                   sequence = rep(c(\"native\",\"venus\"), times=18),\n",
    "                   metric = rep(c(\"T1\",\"T1\",\"T1\",\"T1\",\"MTR\",\"MTR\",\"MTR\", \"MTR\",\"MTsat\",\"MTsat\",\"MTsat\", \"MTsat\"),3),\n",
    "                   pctdif = c(31.2,10.9,16.30,1.8,17.7,8.3,16.9,3.6,14.5,6.7,25.7,3.2,\n",
    "                             20.9,6.6,10.8,3.7,15.5,6.4,14.8,8.3,17.7,6.7,25.2,3.2,\n",
    "                            14.2,0.1,14.5,3.6,14.7,6.3,16.3,2.2,21.9,10.3,26.1,0.6))\n",
    "\n",
    "# group_by(data, sequence) %>%\n",
    "#   dplyr::summarise(\n",
    "#     count = n(),\n",
    "#     median = median(pctdif, na.rm = TRUE),\n",
    "#     IQR = IQR(pctdif, na.rm = TRUE)\n",
    "#   )\n",
    "\n",
    "\n",
    "neutral <- data$pctdif[data$sequence==\"venus\"]\n",
    "native <- data$pctdif[data$sequence==\"native\"]\n",
    "pd <- paired(native,neutral)\n",
    "plt <- plot(pd, type = \"profile\") + ggtitle(\"NATIVE vs VENUS inter-vendor percent difference paired comparison\") + theme_gray()\n",
    "# Plotly outputs will not be rendered in Jupyter Book.\n",
    "# ggplotly(plt)\n",
    "plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ]
   },
   "outputs": [],
   "source": [
    "# Across scores\n",
    "res <- wilcox.test(pctdif ~ sequence,data=data, paired=TRUE) \n",
    "p.adjust(res$p.value, method = \"BH\", n = 6*3)\n",
    "print('Across metrics significance')\n",
    "print(res)\n",
    "# T1\n",
    "res <- wilcox.test(subset(data,sequence==\"native\" & metric==\"T1\")$pctdif,subset(data,sequence==\"venus\" & metric==\"T1\")$pctdif,paired=TRUE,alternative = \"greater\")\n",
    "print('T1 significance before correction')\n",
    "print(res)\n",
    "print('T1 significance after correction')\n",
    "p.adjust(res$p.value, method = \"bonferroni\", n = 3)\n",
    "\n",
    "\n",
    "# MTR\n",
    "res <- wilcox.test(subset(data,sequence==\"native\" & metric==\"MTR\")$pctdif,subset(data,sequence==\"venus\" & metric==\"MTR\")$pctdif,paired=TRUE,alternative = \"greater\")\n",
    "print('MTR significance before correction')\n",
    "print(res)\n",
    "print('MTR significance after correction')\n",
    "p.adjust(res$p.value, method = \"bonferroni\", n = 3)\n",
    "\n",
    "# MTsat\n",
    "res <- wilcox.test(subset(data,sequence==\"native\" & metric==\"MTsat\")$pctdif,subset(data,sequence==\"venus\" & metric==\"MTsat\")$pctdif,paired=TRUE,alternative = \"greater\")\n",
    "print('MTsat significance before correction')\n",
    "print(res)\n",
    "print('MTsat significance after correction')\n",
    "p.adjust(res$p.value, method = \"bonferroni\", n = 3)"
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Edit 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.1.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}