I used the following packages for this analysis: coin, ggplot2, tidyr, ggpubr, and lsmeans.

Hothorn, T., Hornik, K., Van De Wiel, M. A., & Zeileis, A. (2008). Implementing a class of permutation pests: the coin package. Kassambara, A. (2017). ggpubr:“ggplot2” based publication ready plots. R package version 0.1. 6.

Lenth, R. V. (2016). Least-squares means: the R package lsmeans. Journal of statistical software, 69(1), 1-33.

Wickham, H. (2016). tidyr: Easily Tidy Data with ’spread ()’and ’gather ()’Functions. R package version 0.6. 0.

Wickham, H. (2010). ggplot2: elegant graphics for data analysis. J Stat Softw, 35(1), 65-88.

library(coin)
library(ggplot2)
library(tidyr)
library(ggpubr)
library(lsmeans)

I normalized all of the measurement data by dividing the difference between the participant’s brain activation at the blood sugar level of interest and the participant’s brain activation at baseline by the participant’s brain activation at baseline. For example, to calculate the normalized brain activation of the Right Thalamus at a blood sugar level of 45, I computed (Right Thalamus Activation at 45 - Right Thalamus Activation at 90)/Right Thalamus Activation at 90. Participants typically had higher activations at lower blood sugar levels, and thus, many of the normalized brain activation measurements were above 1.0. I saved the normalized activation data in a separate spreadsheet, which I call here at the beginning of my analysis.

df.normalize <-read.csv("C:/Users/wischj/Documents/Hypoglycemia_df_normalize.csv", header = TRUE)
df.normalize$GroupDay <- as.factor(paste(df.normalize$Group,df.normalize$Day, sep = ""))

R’s Coin package allows you to perform permutation tests like the one I did here. I tested for independence between brain region activation (Region) and whether the participant was diabetic or typical, and if typical-whether it was day 1 or day 2 of the test (GroupDay) and the participant’s Glucose Level. I used permutation testing because of the relatively low number of samples and the unknown true distribution of the various brain activation regions.

Nichols, T. E., & Holmes, A. P. (2002). Nonparametric permutation tests for functional neuroimaging: a primer with examples. Human brain mapping, 15(1), 1-25.

IndependenceCheck<-function(Region){
  it<-independence_test(Region + GroupDay~ Glucose.level, data = df.normalize,
                        ytrafo = function(data)
                          trafo(data, numeric_trafo = rank_trafo),
                        teststat = "quadratic")
  return(it)}

I ran an independence test on all of the brain region measurements, and stored these in a list.

IClist<-list()
for (i in 6:88){
  IC<-IndependenceCheck(df.normalize[,i])
  IClist[[i]]<-IC
}

pList<-data.frame(ncol = 4, nrow = 88)
for (i in 6:length(IClist)){
  pList[i, 1]<-names(df.normalize[i])
    pList[i, 2]<-pvalue(IClist[[i]])
    pList[i, 3]<-0
    pList[i, 4]<-0
}

Then I used the Benjamini-Hocberg procedure to account for multiple comparisons. This is not a multiple comparisons correction I was familiar with; however, it was referenced in Lindquest & Mejia, 2015 as an appropriate approach for neuroimaging (and noted as being less conservative than Bonferroni, which seemed appealing to me because of our sample size).

I used a false discovery rate of 5%, but could easly change this if there is a different FDR that is the norm in neuroimaging.

After applying Benjamini-Hocberg to our independence test results, there are four regions that pop out: the left thalamus, the right thalamus, the left pallidum, and the right pallidum. I was happy to see statistical support for your expectation that the thalamus would be significantly affected by your experimental treatment. Does it make sense that the pallidum is also observed to be significantly affected by glucose level?

Lindquist, M. A., & Mejia, A. (2015). Zen and the art of multiple comparisons. Psychosomatic medicine, 77(2), 114.

pList<-pList[6:88,]
colnames(pList)<-c("Region", "pvalue", "Rank", "BH")
pList_sorted<-pList[order(pList$pvalue),]
pList_sorted$Rank<-seq(from=1, to = 83, by = 1)
pList_sorted$BH<-pList_sorted$Rank/length(pList_sorted)*0.05

Then I used MANOVA and ANOVA, as well as a post hoc Tukey test, to further explore the relationships between the thalamus and pallidum and the participant’s status as diabetic or not and their glucose level. The results are all below, but basically, they show that glucose level has a significant effect on brain activation for all four regions. The only participant grouping that suggested a trend for brain activation was the Right Pallidum, where diabetic participants differed from typical day 2 patients (p = 0.0606, so only a trend. Not significant.). These results were somewhat less impressive than the permutation tests results.

##                Df  Pillai approx F num Df den Df    Pr(>F)    
## Glucose.level   1 0.30893  12.1816      4    109 3.201e-08 ***
## GroupDay        2 0.14575   2.1616      8    220   0.03148 *  
## Residuals     112                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##  Response Left.Thalamus.Proper :
##                Df Sum Sq Mean Sq F value    Pr(>F)    
## Glucose.level   1 2.3765 2.37654  37.195 1.501e-08 ***
## Residuals     114 7.2840 0.06389                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Right.Thalamus.Proper :
##                Df Sum Sq Mean Sq F value    Pr(>F)    
## Glucose.level   1 2.3843 2.38430  35.131 3.351e-08 ***
## Residuals     114 7.7370 0.06787                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Left.Pallidum :
##                Df Sum Sq Mean Sq F value    Pr(>F)    
## Glucose.level   1 1.3345 1.33455  36.771 1.768e-08 ***
## Residuals     114 4.1375 0.03629                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##  Response Right.Pallidum :
##                Df Sum Sq Mean Sq F value    Pr(>F)    
## Glucose.level   1 2.5917 2.59167  39.158 7.062e-09 ***
## Residuals     114 7.5450 0.06618                      
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##  Response Left.Thalamus.Proper :
##              Df Sum Sq  Mean Sq F value Pr(>F)
## GroupDay      2 0.1439 0.071927  0.8541 0.4284
## Residuals   113 9.5167 0.084218               
## 
##  Response Right.Thalamus.Proper :
##              Df Sum Sq Mean Sq F value Pr(>F)
## GroupDay      2 0.1649 0.08245  0.9358 0.3953
## Residuals   113 9.9564 0.08811               
## 
##  Response Left.Pallidum :
##              Df Sum Sq  Mean Sq F value Pr(>F)
## GroupDay      2 0.0457 0.022833  0.4755 0.6228
## Residuals   113 5.4264 0.048021               
## 
##  Response Right.Pallidum :
##              Df Sum Sq  Mean Sq F value Pr(>F)
## GroupDay      2 0.2568 0.128413  1.4687 0.2346
## Residuals   113 9.8798 0.087432
## $lsmeans
##  GroupDay   lsmean         SE  df lower.CL upper.CL
##  01       1.170614 0.04079004 112 1.089794 1.251434
##  02       1.113507 0.04134940 112 1.031579 1.195436
##  11       1.216040 0.03929767 112 1.138177 1.293903
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast    estimate         SE  df t.ratio p.value
##  01 - 02   0.05710649 0.05807081 112   0.983  0.5889
##  01 - 11  -0.04542597 0.05665720 112  -0.802  0.7027
##  02 - 11  -0.10253246 0.05707697 112  -1.796  0.1755
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
## $lsmeans
##  GroupDay   lsmean         SE  df lower.CL upper.CL
##  01       1.203124 0.04207312 112 1.119761 1.286486
##  02       1.111700 0.04265007 112 1.027194 1.196206
##  11       1.200735 0.04053380 112 1.120422 1.281048
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast     estimate         SE  df t.ratio p.value
##  01 - 02   0.091423912 0.05989747 112   1.526  0.2825
##  01 - 11   0.002388909 0.05843939 112   0.041  0.9991
##  02 - 11  -0.089035003 0.05887236 112  -1.512  0.2890
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
## $lsmeans
##  GroupDay   lsmean         SE  df lower.CL upper.CL
##  01       1.150600 0.03102916 112 1.089120 1.212080
##  02       1.105505 0.03145467 112 1.043182 1.167829
##  11       1.118605 0.02989391 112 1.059374 1.177836
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast    estimate         SE  df t.ratio p.value
##  01 - 02   0.04509471 0.04417472 112   1.021  0.5653
##  01 - 11   0.03199546 0.04309938 112   0.742  0.7388
##  02 - 11  -0.01309925 0.04341870 112  -0.302  0.9511
## 
## P value adjustment: tukey method for comparing a family of 3 estimates
## $lsmeans
##  GroupDay   lsmean         SE  df lower.CL upper.CL
##  01       1.164513 0.04113168 112 1.083016 1.246010
##  02       1.112332 0.04169572 112 1.029717 1.194947
##  11       1.244436 0.03962681 112 1.165921 1.322952
## 
## Confidence level used: 0.95 
## 
## $contrasts
##  contrast    estimate         SE  df t.ratio p.value
##  01 - 02   0.05218110 0.05855719 112   0.891  0.6470
##  01 - 11  -0.07992341 0.05713174 112  -1.399  0.3449
##  02 - 11  -0.13210451 0.05755502 112  -2.295  0.0606
## 
## P value adjustment: tukey method for comparing a family of 3 estimates

Finally, I put together a graph that compares the four brain regions we’re interested in with the patient groups and glucose levels. Note that the glucose level of 90 is not shown on the plot because everything has been normalized relative to it. Instead, I included a dashed line at 1.0 to enable easy comparison between the baseline brain activation and brain activation at each of the other glucose levels.

CreatePlot<-function(Region, Title){
  p<-ggplot(df.normalize.significant, aes(x=GroupDay, y = Region, fill = Glucose.level))+geom_boxplot()+
    geom_hline(yintercept=1, linetype="dashed", color="black") + ggtitle(Title) + ylab("Normalized Deviation from Baseline") +
    scale_x_discrete(labels=c("1" = "Typical,\n Day 1", "2" = "Typical, \n Day 2",
                              "3" = "Diabetic")) + xlab("") + theme(plot.title = element_text(hjust = 0.5))
  return(p)
}
grid.arrange(p1, p2, p3, p4, ncol = 2, nrow = 2)