MemoHR fMRI ROI Activation Data

1 Load & Process data

savepath <- '../../Analysis/ROI_estimates/'
behavpath <- '../../Data/Behavioral/summary/'
## Load data from csv
### Activation estimates
actdata<-read.csv(paste0(savepath,'cons_ashs_ItemSourceDM_rf_new.csv')) # created with output_roi_cons script in Matlab
### ROI metadata
roidata<-read.csv('../../Analysis/ROI_estimates/roi_info.csv') # contains various labeling schemes for ROIs
actdata <- merge(actdata,roidata,by="roi_file")
rm(roidata)
### Clean up data
# first deal with s21, who had to have EmoK subbed in for EmoM (no EmoM trials)
actdata <- unique(actdata) # remove doubled-up EmoK
tmp <-subset(actdata,subject=="s21" & contrast=="EmoK") # get the EmoK for relabel -> these trial types get averaged anyway
tmp$contrast <- "EmoM"
actdata <- rbind(actdata,tmp) # pull back into main dataframe
### Define new variables
actdata$Emotion <- factor(sapply(actdata$contrast,simplify=TRUE,
                                 function(x) ifelse(grepl('Emo',x),'Emo',
                                                    ifelse(grepl('Neu',x),'Neu',NA))))
actdata <- subset(actdata,!is.na(actdata$Emotion)) # get rid of all other contrasts
actdata$RecMem <- factor(sapply(actdata$contrast,simplify=TRUE,
                                      function(x) ifelse(grepl('R',x),'Rec',
                                                         ifelse(grepl('K',x),'Non-Rec',
                                                                ifelse(grepl('M',x),'Non-Rec',NA)))),
                               levels=c("Rec","Non-Rec"))
actdata$SourceMem <- factor(sapply(actdata$contrast,simplify=TRUE,
                                   function(x) ifelse(grepl('SC',x),'Corr',
                                                      ifelse(grepl('SI',x),'Incorr',NA))))
actdata$BothMem <- factor(ifelse(actdata$RecMem=="Rec" & actdata$SourceMem=="Corr","R+S",
                                 ifelse(actdata$RecMem=="Rec" & actdata$SourceMem=="Incorr","R-S",
                                        ifelse(actdata$RecMem=="Non-Rec","NR",NA))),
                          levels=c("R+S","R-S","NR"))
str(actdata)

'data.frame':   4928 obs. of  15 variables:
 $ roi_file     : Factor w/ 28 levels "AMY_L","AMY_R",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ subject      : Factor w/ 22 levels "s04","s05","s06",..: 1 1 1 1 1 1 1 1 2 2 ...
 $ contrast     : Factor w/ 9 levels "ALL","EmoK","EmoM",..: 4 5 2 3 8 9 6 7 4 5 ...
 $ activity     : num  1.582 0.652 1.389 1.901 0.687 ...
 $ numvox       : int  520 520 520 520 520 520 520 520 517 517 ...
 $ Region       : Factor w/ 12 levels "Amyg","BasolatAmyg",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Hemisphere   : Factor w/ 2 levels "L","R": 1 1 1 1 1 1 1 1 1 1 ...
 $ MainRegion   : Factor w/ 3 levels "Amyg","CorticalMTL",..: NA NA NA NA NA NA NA NA NA NA ...
 $ AP           : Factor w/ 3 levels "body","head",..: NA NA NA NA NA NA NA NA NA NA ...
 $ full_filename: Factor w/ 40 levels "rseg_AMY_L.nii",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ RegionAP     : Factor w/ 21 levels "Amyg","BasolatAmyg",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ Emotion      : Factor w/ 2 levels "Emo","Neu": 1 1 1 1 2 2 2 2 1 1 ...
 $ RecMem       : Factor w/ 2 levels "Rec","Non-Rec": 1 1 2 2 1 1 2 2 1 1 ...
 $ SourceMem    : Factor w/ 2 levels "Corr","Incorr": 1 2 NA NA 1 2 NA NA 1 2 ...
 $ BothMem      : Factor w/ 3 levels "R+S","R-S","NR": 1 2 3 3 1 2 3 3 1 2 ...

kable(table(actdata$BothMem,actdata$subject))

	s04	s05	s06	s07	s08	s09	s10	s11	s13	s14	s16	s17	s18	s19	s20	s21	s22	s23	s25	s26	s27	s28
R+S	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56
R-S	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56	56
NR	112	112	112	112	112	112	112	112	112	112	112	112	112	112	112	112	112	112	112	112	112	112

## Create dataframe with 'baseline-corrected' scores, i.e., recollection trial types minus non-recollection
nrbaseline <-
  group_by(actdata,subject,Emotion,roi_file) %>%
  summarize(nr.baseline = mean(activity[RecMem=="Non-Rec"],na.rm=TRUE)) # average over Fam & Miss
actdata.rec <-
  filter(actdata,RecMem=="Rec") %>%
  left_join(nrbaseline)  %>%
  mutate(activity = activity - nr.baseline) # values in actdata.rec are now recollection-related activity

Joining, by = c("roi_file", "subject", "Emotion")

  
## Clean up and save
save(actdata,file=paste0(savepath,'MemoHR_AnatROI_Data_ItemSource.Rdata'))

1.1 Report number of voxels in each ROI

vox <- group_by(actdata,RegionAP) %>%
  summarize(meanvox=mean(numvox),
            minvox=min(numvox),
            maxvox=max(numvox))
kable(vox)

RegionAP	meanvox	minvox	maxvox
Amyg	558.7	397	726
BasolatAmyg	312.5	210	492
CA1.body	194.6	166	237
CA23DG.body	140.3	91	202
CentralAmyg	38.3	13	68
CorticalAmyg	116.5	46	179
MedialAmyg	91.4	51	151
PHC	429.2	271	619
PRC	465.7	247	846
subiculum.body	43.7	29	61
wholeHipp	843.0	600	1071
wholeHipp.body	378.6	312	478
wholeHipp.head	425.9	229	683
wholeHipp.tail	38.6	0	91

2 Plots

2.1 Calculate summary values first

summarydata.all <- group_by(actdata,subject,Emotion,RecMem,SourceMem,Region,RegionAP,MainRegion,AP) %>%
  summarize(activity=mean(activity)) # This part averages across hemispheres, and fam/miss are averaged into non-rec
meandata.all <- group_by(summarydata.all,Region,RegionAP,MainRegion,AP,Emotion) %>%
  group_by(Emotion,RecMem,SourceMem,Region,RegionAP,MainRegion,AP) %>%
  summarize(meanval=mean(activity,na.rm=TRUE),sdval=sd(activity,na.rm=TRUE),
            sem=sdval/sqrt(length(activity)),
            semmin=meanval-sem,
            semmax=meanval+sem,
            nsubj=length(activity))
meandata.all$EmoMem <- interaction(meandata.all$RecMem,meandata.all$Emotion)
meandata.all$MemMem <- interaction(meandata.all$RecMem,meandata.all$SourceMem)
meandata.all$MemMem <- factor(meandata.all$MemMem,levels=c("Rec.Corr","Rec.Incorr","Non-Rec"))
meandata.all$MemMem[is.na(meandata.all$MemMem)] <- "Non-Rec"
summarydata.rec <- group_by(actdata.rec,subject,Emotion,RecMem,SourceMem,Region,RegionAP,MainRegion,AP) %>%
  summarize(activity=mean(activity))  # This part averages across hemispheres
meandata.rec <- group_by(summarydata.rec,Region,RegionAP,MainRegion,AP,Emotion) %>%
  group_by(Emotion,RecMem,SourceMem,Region,RegionAP,MainRegion,AP) %>%
  summarize(meanval=mean(activity,na.rm=TRUE),sdval=sd(activity,na.rm=TRUE),
            sem=sdval/sqrt(length(activity)),
            semmin=meanval-sem,
            semmax=meanval+sem,
            nsubj=length(activity))
meandata.rec$EmoMem <- interaction(meandata.rec$RecMem,meandata.rec$Emotion)
meandata.rec$MemMem <- interaction(meandata.rec$RecMem,meandata.rec$SourceMem)
meandata.rec$MemMem <- factor(meandata.rec$MemMem,levels=c("Rec.Corr","Rec.Incorr","Non-Rec"))
meandata.rec$MemMem[is.na(meandata.rec$MemMem)] <- "Non-Rec"

2.2 Plots with all bars

2.2.1 Main regions

# All regions
plotcolors <- c('#d01c8b','#f1b6da',"gray70",'#4dac26','#b8e186',"gray70")
meandata.all %>%
  filter(RegionAP=="Amyg" | RegionAP=="wholeHipp") %>%
  ungroup() %>%
  mutate(RegionAP = factor(RegionAP,levels=c("Amyg","wholeHipp")),
         RegionLabel = factor(ifelse(RegionAP=="wholeHipp","Hipp","Amyg"),
                              levels=c("Amyg","Hipp")),
         MemoryCond = factor(ifelse(MemMem=="Rec.Corr","Rec + Src",ifelse(MemMem=="Rec.Incorr","Rec - Src","Non-Rec")),
                             levels=c("Rec + Src","Rec - Src","Non-Rec")),
         EmoMem = interaction(MemoryCond,Emotion)) %>%
  ggplot(.,aes(x=MemoryCond,y=meanval)) + geom_bar(stat="identity",position="dodge",aes(fill=EmoMem)) +
  geom_errorbar(aes(ymin=meanval-sem,ymax=meanval+sem,width=.25)) + theme_minimal(18) +
  theme(axis.title.x=element_blank()) + facet_grid(.~RegionLabel+Emotion) +
  xlab("Mean activity (au)") + scale_fill_manual(values=plotcolors) +  guides(color=FALSE, fill=FALSE) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

meandata.all %>%
  filter(RegionAP=="PHC" | RegionAP=="PRC") %>%
  ungroup() %>%
  mutate(RegionAP = factor(RegionAP,levels=c("PRC","PHC")),
         MemoryCond = factor(ifelse(MemMem=="Rec.Corr","Rec + Src",ifelse(MemMem=="Rec.Incorr","Rec - Src","Non-Rec")),
                             levels=c("Rec + Src","Rec - Src","Non-Rec")),
         EmoMem = interaction(MemoryCond,Emotion)) %>%
  ggplot(.,aes(x=MemoryCond,y=meanval)) + geom_bar(stat="identity",position="dodge",aes(fill=EmoMem)) +
  geom_errorbar(aes(ymin=meanval-sem,ymax=meanval+sem,width=.25)) + theme_minimal(18) +
  theme(axis.text.x=element_blank(),axis.title.x=element_blank()) + facet_grid(.~RegionAP+Emotion) +
  xlab("Mean activity (au)") + scale_fill_manual(values=plotcolors) +  guides(color=FALSE, fill=FALSE) + 
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

2.3 Plots with bars subtracting non-recollection

2.3.1 Single Plots

plotcolors <- c('#d01c8b','#f1b6da','#4dac26','#b8e186')
singlePlotROI <- function(roi) {
  meandata.rec %>%
  filter(RegionAP==roi) %>%
  ungroup() %>%
  mutate(MemoryCond = factor(ifelse(MemMem=="Rec.Corr","Rec + Src",ifelse(MemMem=="Rec.Incorr","Rec - Src","Non-Rec")),
                             levels=c("Rec + Src","Rec - Src","Non-Rec")),
         EmoMem = interaction(MemoryCond,Emotion)) %>%
  ggplot(.,aes(x=MemoryCond,y=meanval)) + geom_bar(stat="identity",position="dodge",aes(fill=EmoMem)) +
  geom_errorbar(aes(ymin=meanval-sem,ymax=meanval+sem,width=.25)) + theme_minimal(18) +
  theme(axis.text.x=element_blank(),axis.title.x=element_blank()) + facet_grid(.~Emotion,switch="x") +
  ylab("Recollection-related activity") + scale_fill_manual(values=plotcolors) + guides(color=FALSE) +
    ggtitle(roi) 
}
singlePlotROI_fixedLim <- function(roi,limits) { # fix limits for some plots for better visualization
  singlePlotROI(roi) + ylim(limits)
}
# Plot all regions
singlePlotROI_fixedLim('Amyg',c(-.25,1))

singlePlotROI_fixedLim('wholeHipp',c(-.25,1))

singlePlotROI('PRC')

singlePlotROI('PHC')

3 Stats

3.1 All regions

statdata1 <- subset(actdata.rec,RegionAP=="wholeHipp" | Region=="Amyg"  | Region=="PRC"  | Region=="PHC") 
summary(droplevels(statdata1$Region)) # Note that hemisphere is still included but ezANOVA will average them

     Amyg       PHC       PRC wholeHipp 
      176       176       176       176 

anova_mod1<-ezANOVA(statdata1,dv=activity,wid=subject,within=list(BothMem,Emotion,Region),type=3,detailed=TRUE)
print(anova_mod1)

$ANOVA
                  Effect DFn DFd    SSn  SSd      F         p p<.05     ges
1            (Intercept)   1  21 70.806 57.0 26.070 0.0000467     * 0.22872
2                BothMem   1  21  8.005 19.0  8.835 0.0072650     * 0.03244
3                Emotion   1  21  1.515 36.8  0.865 0.3629812       0.00631
4                 Region   3  63  6.031 45.4  2.790 0.0476819     * 0.02464
5        BothMem:Emotion   1  21  0.283 18.5  0.322 0.5766273       0.00118
6         BothMem:Region   3  63  6.078 16.4  7.791 0.0001671     * 0.02483
7         Emotion:Region   3  63  7.185 31.5  4.782 0.0045822     * 0.02921
8 BothMem:Emotion:Region   3  63  0.297 14.1  0.442 0.7234704       0.00124

$`Mauchly's Test for Sphericity`
                  Effect     W      p p<.05
4                 Region 0.517 0.0236     *
6         BothMem:Region 0.515 0.0229     *
7         Emotion:Region 0.645 0.1243      
8 BothMem:Emotion:Region 0.432 0.0055     *

$`Sphericity Corrections`
                  Effect   GGe   p[GG] p[GG]<.05   HFe    p[HF] p[HF]<.05
4                 Region 0.721 0.06791           0.807 0.060881          
6         BothMem:Region 0.699 0.00108         * 0.778 0.000661         *
7         Emotion:Region 0.818 0.00822         * 0.935 0.005648         *
8 BothMem:Emotion:Region 0.655 0.64192           0.722 0.660713          

rm(statdata1,anova_mod1)

4 Follow-up tests

4.1 Amygdala vs PRC

statdata1 <- subset(actdata.rec,RegionAP=="PRC" | Region=="Amyg") 
summary(droplevels(statdata1$Region)) # Note that hemisphere is still included but ezANOVA will average them

Amyg  PRC 
 176  176 

anova_mod1<-ezANOVA(statdata1,dv=activity,wid=subject,within=list(BothMem,Emotion,Region),type=3,detailed=TRUE)
print(anova_mod1$ANOVA)

statdata2 <- subset(actdata.rec,RegionAP=="PRC" )
anova_mod2<-ezANOVA(statdata2,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod2$ANOVA)

statdata3 <- subset(actdata.rec,Region=="Amyg")
anova_mod3<-ezANOVA(statdata3,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod3$ANOVA)

rm(statdata1,statdata2,statdata3,anova_mod1,anova_mod2,anova_mod3)

4.2 PHC vs whole hippocampus

statdata1 <- subset(actdata.rec,RegionAP=="wholeHipp" | Region=="PHC") 
summary(droplevels(statdata1$Region)) # Note that hemisphere is still included but ezANOVA will average them

      PHC wholeHipp 
      176       176 

anova_mod1<-ezANOVA(statdata1,dv=activity,wid=subject,within=list(BothMem,Emotion,Region),type=3,detailed=TRUE)
print(anova_mod1$ANOVA)

statdata2 <- subset(actdata.rec,RegionAP=="wholeHipp" )
anova_mod2<-ezANOVA(statdata2,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod2$ANOVA)

statdata3 <- subset(actdata.rec,Region=="PHC")
anova_mod3<-ezANOVA(statdata3,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod3$ANOVA)

rm(statdata1,statdata2,statdata3,anova_mod1,anova_mod2,anova_mod3)

4.3 Amygdala vs whole hippocampus

statdata1 <- subset(actdata.rec,RegionAP=="wholeHipp" | Region=="Amyg") 
summary(droplevels(statdata1$Region)) # Note that hemisphere is still included but ezANOVA will average them

     Amyg wholeHipp 
      176       176 

anova_mod1<-ezANOVA(statdata1,dv=activity,wid=subject,within=list(BothMem,Emotion,Region),type=3,detailed=TRUE)
print(anova_mod1$ANOVA)

statdata2 <- subset(actdata.rec,RegionAP=="wholeHipp" )
anova_mod2<-ezANOVA(statdata2,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod2$ANOVA)

statdata3 <- subset(actdata.rec,Region=="Amyg")
anova_mod3<-ezANOVA(statdata3,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod3$ANOVA)

rm(statdata1,statdata2,statdata3,anova_mod1,anova_mod2,anova_mod3)

4.4 Subregions within amygdala

statdata1 <- subset(actdata.rec,MainRegion=="Amyg")
summary(droplevels(statdata1$Region))

 BasolatAmyg  CentralAmyg CorticalAmyg   MedialAmyg 
         176          176          176          176 

anova_mod1<-ezANOVA(statdata1,dv=activity,wid=subject,within=list(BothMem,Emotion,Region),type=3,detailed=TRUE)
print(anova_mod1)

$ANOVA
                  Effect DFn DFd      SSn   SSd       F        p p<.05      ges
1            (Intercept)   1  21 128.3956 125.0 21.5705 0.000139     * 0.151510
2                BothMem   1  21   0.0438  48.3  0.0191 0.891486       0.000061
3                Emotion   1  21  31.9045 108.7  6.1660 0.021548     * 0.042486
4                 Region   3  63  27.5015 174.9  3.3028 0.025861     * 0.036838
5        BothMem:Emotion   1  21   0.3365  36.5  0.1937 0.664322       0.000468
6         BothMem:Region   3  63   0.4582  54.1  0.1778 0.911030       0.000637
7         Emotion:Region   3  63   3.6340 131.8  0.5789 0.631021       0.005028
8 BothMem:Emotion:Region   3  63   1.5538  39.8  0.8190 0.488233       0.002156

$`Mauchly's Test for Sphericity`
                  Effect     W       p p<.05
4                 Region 0.681 0.18170      
6         BothMem:Region 0.421 0.00446     *
7         Emotion:Region 0.794 0.47438      
8 BothMem:Emotion:Region 0.774 0.41049      

$`Sphericity Corrections`
                  Effect   GGe  p[GG] p[GG]<.05   HFe  p[HF] p[HF]<.05
4                 Region 0.842 0.0341         * 0.967 0.0274         *
6         BothMem:Region 0.625 0.8244           0.685 0.8431          
7         Emotion:Region 0.869 0.6081           1.004 0.6310          
8 BothMem:Emotion:Region 0.880 0.4756           1.018 0.4882          

statdata2 <- subset(actdata.rec,Region=="BasolatAmyg")
anova_mod2<-ezANOVA(statdata2,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod2$ANOVA)

statdata3 <- subset(actdata.rec,Region=="CentralAmyg")
anova_mod3<-ezANOVA(statdata3,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod3$ANOVA)

statdata4 <- subset(actdata.rec,Region=="CorticalAmyg")
anova_mod4<-ezANOVA(statdata4,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod4$ANOVA)

statdata5 <- subset(actdata.rec,Region=="MedialAmyg")
anova_mod5<-ezANOVA(statdata5,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod5$ANOVA)

rm(statdata1,statdata2,statdata3,statdata4,statdata5,anova_mod1,anova_mod2,anova_mod3,anova_mod4,anova_mod5)

4.5 Posterior vs anterior hippocampus

statdata1 <- subset(actdata.rec,(RegionAP=="wholeHipp.head" | RegionAP=="wholeHipp.body"))
summary(droplevels(statdata1$RegionAP))

wholeHipp.body wholeHipp.head 
           176            176 

anova_mod1<-ezANOVA(statdata1,dv=activity,wid=subject,within=list(BothMem,Emotion,RegionAP),type=3,detailed=TRUE)
print(anova_mod1$ANOVA)

statdata2 <- subset(actdata.rec,(RegionAP=="wholeHipp.head") )
anova_mod2<-ezANOVA(statdata2,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod2$ANOVA)

statdata3 <- subset(actdata.rec,(RegionAP=="wholeHipp.body" ) )
anova_mod3<-ezANOVA(statdata3,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod3$ANOVA)

rm(statdata1,statdata2,statdata3,anova_mod1,anova_mod2,anova_mod3)

4.6 Subregions within posterior hippocampus

statdata1 <- subset(actdata.rec,MainRegion=="Hipp")
summary(droplevels(statdata1$Region))

      CA1    CA23DG subiculum 
      176       176       176 

anova_mod<-ezANOVA(statdata1,dv=activity,wid=subject,within=list(BothMem,Emotion,Region),type=3,detailed=TRUE)
print(anova_mod)

$ANOVA
                  Effect DFn DFd     SSn   SSd       F       p p<.05       ges
1            (Intercept)   1  21  9.9296 113.3 1.83974 0.18939       0.0216324
2                BothMem   1  21 14.2476  30.6 9.76278 0.00512     * 0.0307503
3                Emotion   1  21  0.1385  83.5 0.03484 0.85372       0.0003084
4                 Region   2  42  5.9476  77.0 1.62135 0.20976       0.0130708
5        BothMem:Emotion   1  21  0.5762  29.4 0.41086 0.52847       0.0012814
6         BothMem:Region   2  42  0.2287  28.4 0.16916 0.84495       0.0005089
7         Emotion:Region   2  42  0.0276  58.2 0.00994 0.99011       0.0000614
8 BothMem:Emotion:Region   2  42  0.6745  28.5 0.49711 0.61182       0.0014997

$`Mauchly's Test for Sphericity`
                  Effect     W      p p<.05
4                 Region 0.717 0.0361     *
6         BothMem:Region 0.715 0.0351     *
7         Emotion:Region 0.832 0.1585      
8 BothMem:Emotion:Region 0.631 0.0101     *

$`Sphericity Corrections`
                  Effect   GGe p[GG] p[GG]<.05   HFe p[HF] p[HF]<.05
4                 Region 0.780 0.215           0.831 0.214          
6         BothMem:Region 0.778 0.791           0.829 0.805          
7         Emotion:Region 0.856 0.982           0.925 0.987          
8 BothMem:Emotion:Region 0.731 0.555           0.772 0.565          

statdata2 <- subset(actdata.rec,Region=="CA1")
anova_mod2<-ezANOVA(statdata2,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod2$ANOVA)

statdata3 <- subset(actdata.rec,Region=="CA23DG")
anova_mod3<-ezANOVA(statdata3,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod3$ANOVA)

statdata4 <- subset(actdata.rec,Region=="subiculum")
anova_mod4<-ezANOVA(statdata4,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod4$ANOVA)

rm(statdata1,anova_mod,statdata2,statdata3,statdata4,anova_mod2,anova_mod3,anova_mod4)

4.7 PRC vs PHC

statdata1 <- subset(actdata.rec,(Region=="PRC" | Region=="PHC") )
anova_mod<-ezANOVA(statdata1,dv=activity,wid=subject,within=list(BothMem,Emotion,Region),type=3,detailed=TRUE)
print(anova_mod$ANOVA)

statdata2 <- subset(actdata.rec,(Region=="PRC") )
anova_mod2<-ezANOVA(statdata2,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod2$ANOVA)

statdata3 <- subset(actdata.rec,Region=="PHC" )
anova_mod3<-ezANOVA(statdata3,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
print(anova_mod3$ANOVA)

rm(statdata1,statdata2,statdata3,anova_mod1,anova_mod2,anova_mod3)

5 ROI maps

plotROImap <- function(values,colornames) {
  
  ids <- factor(c("Amyg","PRC","PHC","wholeHipp.head","wholeHipp.body","subiculum.body","CA1.body","CA23DG.body","CorticalAmyg","CentralAmyg","MedialAmyg","BasolatAmyg"))
  nid <- length(ids)
  
  # Define polygon coordinates (in this case, 4 each for a rectangle)
  positions <- data.frame(
    id = rep(ids, each = 4),
    order = rep(c(1,2,3,4),nid),
    x = c(1,2,2,1,
          1.5,3,3,1.5,
          3,4,4,3,
          2,3,3,2,
          3,4,4,3,
          3,4,4,3,
          3,4,4,3,
          3,4,4,3,
          1,2,2,1,
          1,2,2,1,
          1,2,2,1,
          1,2,2,1),
    y = c(2,2,3,3,
          1,1,2,2,
          1,1,2,2,
          2,2,3,3,
          2,2,3,3,
          2,2,2.3333,2.3333,
          2.3333,2.3333,2.6666,2.6666,
          2.6666,2.6666,3,3,
          2,2,2.25,2.25,
          2.25,2.25,2.5,2.5,
          2.5,2.5,2.75,2.75,
          2.75,2.75,3,3)
  )
  
  # Calculate midpoint of each poly for label
  positions %>%
    group_by(id) %>%
    summarize(label_x = mean(x),
              label_y = mean(y)) %>%
    left_join(values) %>%
    filter(!is.na(value) & !is.null(value) & value!="NULL") -> values
  
  # Merge id & label info with positions
  datapoly <- merge(values, positions, by=c("id")) %>%
    arrange(id,order) %>%
    ungroup()
  
  # Prettier labels
  values$id[values$id=="wholeHipp.head"] <- "Hipp.ant"
  values$id[values$id=="wholeHipp.body"] <- "Hipp.post"
  
  values$id[values$id=="BasolatAmyg"] <- "Amyg.basolateral"
  values$id[values$id=="MedialAmyg"] <- "Amyg.basomedial"
  values$id[values$id=="CentralAmyg"] <- "Amyg.centromedial"
  values$id[values$id=="CorticalAmyg"] <- "Amyg.cortical"
  # Plot polygons with values coding fill color
  p <- ggplot(datapoly, aes(x=x, y=y)) + geom_polygon(colour="black",aes(fill=value, group=id)) +
    theme_minimal(20) + 
    theme(axis.text = element_blank(),axis.ticks = element_blank()) +
    theme(panel.grid.major = element_blank(),panel.grid.minor=element_blank()) + 
    xlab('') + ylab('') + 
    scale_fill_gradient(low=colornames[1],high=colornames[2],limits=c(0,NA),na.value="gray60",name="effect size") + 
    geom_text(data=values,size=8,aes(label=id,x=label_x,y=label_y))
  return(p)
}

5.1 Get stats summary info for plots

getAnovaEffectSize <- function(df,effectrow) {
  df <- data.frame(df) # bc ezANOVA hates grouped dataframes
  anova_mod <- ezANOVA(df,dv=activity,wid=subject,within=list(BothMem,Emotion),type=3,detailed=TRUE)
  effectsize <- anova_mod$ANOVA[9][[1]][effectrow]
  return(effectsize)
}
actdata.rec %>%
  filter(RegionAP %in% c("Amyg","wholeHipp.head","wholeHipp.body","PHC","PRC")) %>%
  group_by(RegionAP) %>%
  do(memory = getAnovaEffectSize(.,1),
     source = getAnovaEffectSize(.,2),
     emotion = getAnovaEffectSize(.,3)) %>%
  unnest() -> effectsizes.roi
actdata.rec %>%
  filter(RegionAP %in% c("CorticalAmyg","CentralAmyg","MedialAmyg","BasolatAmyg","CA1.body","CA23DG.body","subiculum.body","wholeHipp.head","PHC","PRC")) %>%
  group_by(RegionAP) %>%
  do(memory = getAnovaEffectSize(.,1),
     source = getAnovaEffectSize(.,2),
     emotion = getAnovaEffectSize(.,3)) %>%
  unnest() -> effectsizes.roi.subreg
print(effectsizes.roi)

print(effectsizes.roi.subreg)

5.2 Major subdivisions

map.mem <- data.frame(id=as.character(effectsizes.roi$RegionAP),value=as.numeric(effectsizes.roi$memory))
plotROImap(map.mem,c('white','orange')) + ggtitle("Recollection-related activity")

map.emo <- data.frame(id=as.character(effectsizes.roi$RegionAP),value=as.numeric(effectsizes.roi$emotion))
plotROImap(map.emo,c('white','tomato')) + ggtitle("Effect of emotion on recollection-related activity")

map.src <- data.frame(id=as.character(effectsizes.roi$RegionAP),value=as.numeric(effectsizes.roi$source))
plotROImap(map.src,c('white','skyblue1')) + ggtitle("Effect of context encoding on recollection-related activity")

5.3 With subregions

map.mem2 <- data.frame(id=as.character(effectsizes.roi.subreg$RegionAP),value=as.numeric(effectsizes.roi.subreg$memory))
plotROImap(map.mem2,c('white','orange')) + ggtitle("Recollection-related activity")

map.emo2 <- data.frame(id=as.character(effectsizes.roi.subreg$RegionAP),value=as.numeric(effectsizes.roi.subreg$emotion))
plotROImap(map.emo2,c('white','tomato')) + ggtitle("Effect of emotion on recollection-related activity")

map.src2 <- data.frame(id=as.character(effectsizes.roi.subreg$RegionAP),value=as.numeric(effectsizes.roi.subreg$source))
plotROImap(map.src2,c('white','skyblue1')) + ggtitle("Effect of context encoding on recollection-related activity")