{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Shift function analysis\n",
    "\n",
    "\n",
    ":::{tip} Living Preprint\n",
    "If you are running this notebook on Binder provided by NeuroLibre, you don't have to install anything or download data.\n",
    ":::\n",
    "\n",
    "::::{warning} Local run\n",
    "If running 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 after downloading the data:\n",
    "\n",
    "```r\n",
    "derivativesDir <- \".../qMRFlow/\"\n",
    "```\n",
    "\n",
    ":::{important} 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 `.../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",
    ":::\n",
    "::::"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-output",
     "hide-input"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "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 <- \"../data/venus-data/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",
    "* `.../qMRFLow/sub-invivo#/T1_SF_PLOTS` for T1 comparisons\n",
    "* `.../qMRFLow/sub-invivo#/MTR_SF_PLOTS` for MTR comparisons\n",
    "* `.../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": null,
   "metadata": {
    "tags": [
     "hide-input"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "# metrics = list(\"MTR\",\"T1\",\"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-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,\"/\",sessions[i],\"/stat/\",sub,\"_\",sessions[i],\"_desc-wm_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(sessions[i])\n",
    "#   varName <- substr(varName,5,nchar(varName)-3)\n",
    "#   print(varName)\n",
    "#   df <- mutate(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",
    "# 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",
    "#     }\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",
    "# grp1 <- df[[sel1]]\n",
    "# grp2 <- df[[sel2]]\n",
    "\n",
    "# sfInput <- mkt2(grp1,grp2,group_labels=c(sel1,sel2))\n",
    "# sf <- shiftdhd_pbci(\n",
    "#   data = sfInput,\n",
    "#   formula = obs ~ gr,\n",
    "#   nboot = 250,\n",
    "#   alpha = 0.05/30\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,\"_rev\",'.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,\"_rev\",'.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": "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": null,
   "metadata": {
    "tags": [
     "hide-input"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#| label: shiftexample\n",
    "\n",
    "subject =\"sub-invivo3\" # or sub-invivo2 or sub-invivo3\n",
    "metric =\"MTsat\" # or MTR or MTsat\n",
    "scanner1 = \"S1\" # or S1 or S2\n",
    "scanner2 = \"G1\" # 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",
    "* `.../HSF` \n",
    "\n",
    "```{warning}\n",
    "Executing the following code block takes about 10 - 15 minutes! \n",
    "\n",
    "The (static) outputs are already available at `.../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"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "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"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#| label: hsfexample\n",
    "\n",
    "metric =\"MTsat\" # or MTR or T1\n",
    "scanner1 = \"G1\" # or S1 or S2\n",
    "scanner2 = \"S1\" # or G1 or S2\n",
    "\n",
    "\n",
    "display_hsf(\"../data/venus-data/\",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"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "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 `.../qMRFlow/sub-invivo#/_SF_` to construct the data frame. \n",
    "`````"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": [
     "hide-input"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#| label: pairedcomp\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",
    "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"
    ],
    "vscode": {
     "languageId": "r"
    }
   },
   "outputs": [],
   "source": [
    "#| label: sigtest\n",
    "\n",
    "# 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.2.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}