MemoHR fMRI ROI - Predict memory outcomes

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_betas_ItemSourceDM_rf_new_ST.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)
actdata$RegionAP <- droplevels(actdata$RegionAP)
actdata$Region <- droplevels(actdata$Region)
## Get condition info from beta descriptions
actdata$contrast <- as.character(actdata$contrast)
actdata %>%
  mutate(emotion = factor(str_sub(contrast,29,31))) %>%
  filter(emotion!="Oth") %>%
  mutate(session = factor(str_sub(contrast,26,26)),
         emotion = factor(str_sub(contrast,29,31)),
         memory = str_sub(contrast,32,35)) %>%
  separate(memory,c("memory","tmp"),sep="_")  %>% # this is the line that throws the "too few values" warning - not a real issue
  select(roi_file,subject,contrast,activity,Region,Hemisphere,RegionAP,emotion,session,memory) -> actdata
actdata %>%
  mutate(memory = factor(ifelse(memory=="-SI","R-SI",memory),
                         levels=c("R-SC","R-SI","K","M")),
         recmem = factor(ifelse(memory=="R-SI"|memory=="R-SC","Rec","Non-Rec"),levels=c("Non-Rec","Rec")),
         sourcemem = factor(ifelse(memory=="R-SI","SI",
                                ifelse(memory=="R-SC","SC",NA)),levels=c("SI","SC"))) -> actdata 
str(actdata)

'data.frame':   49668 obs. of  12 variables:
 $ roi_file  : Factor w/ 12 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 1 1 ...
 $ contrast  : chr  "spm_spm:beta (0001) - Sn(1) EmoR-SC_1*bf(1)" "spm_spm:beta (0002) - Sn(1) EmoR-SC_2*bf(1)" "spm_spm:beta (0003) - Sn(1) EmoR-SI_1*bf(1)" "spm_spm:beta (0004) - Sn(1) EmoR-SI_2*bf(1)" ...
 $ activity  : num  2.151 -0.958 9.683 -6.312 -1.97 ...
 $ Region    : Factor w/ 4 levels "Amyg","PHC","PRC",..: 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 ...
 $ RegionAP  : Factor w/ 6 levels "Amyg","PHC","PRC",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ emotion   : Factor w/ 2 levels "Emo","Neu": 1 1 1 1 1 1 1 1 1 1 ...
 $ session   : Factor w/ 6 levels "1","2","3","4",..: 1 1 1 1 1 1 1 1 1 1 ...
 $ memory    : Factor w/ 4 levels "R-SC","R-SI",..: 1 1 2 2 2 3 3 3 3 3 ...
 $ recmem    : Factor w/ 2 levels "Non-Rec","Rec": 2 2 2 2 2 1 1 1 1 1 ...
 $ sourcemem : Factor w/ 2 levels "SI","SC": 2 2 1 1 1 NA NA NA NA NA ...

actdata %>%
  group_by(subject,roi_file) %>%
  mutate(meanact = mean(activity),
         sdact = sd(activity),
         activity.z = (activity - meanact)/sdact,
         activity.z = ifelse(abs(activity.z)>3,NA,activity.z)) -> actdata # remove outliers
## Average over hemispheres
actdata %>%
  group_by(subject,RegionAP,session,emotion,memory,recmem,sourcemem,contrast) %>%
  summarize(activity = mean(activity.z))  -> avgdata 
## Cast RegionAP to columns
avgdata %>%
  spread(RegionAP,activity) -> avgdata.roi
## Define colors
emocolor <- '#d01c8b'
neucolor <- '#4dac26'

2 Basic plots

2.1 Sanity checks

To make sure averaged results look roughly the way that I expect from prior analyses

avgdata %>%
  group_by(subject,memory,emotion,RegionAP) %>%
  summarize(activity = mean(activity)) %>%
  ggplot(aes(x=memory,y=activity,color=emotion)) + geom_boxplot() + facet_grid(.~RegionAP)

To check out distributions

avgdata.roi %>%
  ggplot(aes(x=Amyg,y=PRC)) + facet_grid(.~recmem) + stat_smooth(method="lm",se=FALSE) + geom_point(aes(color=subject)) 

3 Regression predicting recollection outcomes

avgdata.roi %>%
  filter(!is.na(recmem)) %>%
  filter(!is.na(Amyg+PRC+PHC+wholeHipp+wholeHipp.head+wholeHipp.body)) -> avgdata.roi.rec # this limits to trials with intact data

3.1 Each region alone

# Set up base models, no emotion
intonly <- glmer(recmem ~ (1|subject/session),data=avgdata.roi.rec,family=binomial)
amyg <- glmer(recmem ~ Amyg + (1|subject/session),data=avgdata.roi.rec,family=binomial)
hipp <- glmer(recmem ~ wholeHipp + (1|subject/session),data=avgdata.roi.rec,family=binomial)
prc <- glmer(recmem ~ PRC + (1|subject/session),data=avgdata.roi.rec,family=binomial)
phc <- glmer(recmem ~ PHC + (1|subject/session),data=avgdata.roi.rec,family=binomial)
ahipp <- glmer(recmem ~ wholeHipp.head + (1|subject/session),data=avgdata.roi.rec,family=binomial)
anova(intonly,amyg)

Data: avgdata.roi.rec
Models:
intonly: recmem ~ (1 | subject/session)
amyg: recmem ~ Amyg + (1 | subject/session)
        Df  AIC  BIC logLik deviance Chisq Chi Df    Pr(>Chisq)    
intonly  3 5176 5195  -2585     5170                               
amyg     4 5140 5165  -2566     5132  38.3      1 0.00000000061 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(intonly,hipp)

Data: avgdata.roi.rec
Models:
intonly: recmem ~ (1 | subject/session)
hipp: recmem ~ wholeHipp + (1 | subject/session)
        Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)    
intonly  3 5176 5195  -2585     5170                            
hipp     4 5160 5185  -2576     5152  18.9      1   0.000014 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(intonly,prc)

Data: avgdata.roi.rec
Models:
intonly: recmem ~ (1 | subject/session)
prc: recmem ~ PRC + (1 | subject/session)
        Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)    
intonly  3 5176 5195  -2585     5170                            
prc      4 5153 5178  -2572     5145  25.6      1 0.00000042 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(intonly,phc)

Data: avgdata.roi.rec
Models:
intonly: recmem ~ (1 | subject/session)
phc: recmem ~ PHC + (1 | subject/session)
        Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)    
intonly  3 5176 5195  -2585     5170                            
phc      4 5163 5188  -2577     5155  15.5      1   0.000084 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# make sure that the amygdala~memory effect is modulated by emotion (which it should be)
amygemo <- glmer(recmem ~ Amyg + emotion + (1|subject/session),data=avgdata.roi.rec,family=binomial)
amygemoxx <- glmer(recmem ~ Amyg*emotion + (1|subject/session),data=avgdata.roi.rec,family=binomial)
anova(amygemo,amygemoxx)

Data: avgdata.roi.rec
Models:
amygemo: recmem ~ Amyg + emotion + (1 | subject/session)
amygemoxx: recmem ~ Amyg * emotion + (1 | subject/session)
          Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)    
amygemo    5 5125 5156  -2557     5115                            
amygemoxx  6 5114 5152  -2551     5102  12.7      1    0.00037 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

# Set up rest of base models, with emotion
hippemo <- glmer(recmem ~ wholeHipp + emotion + (1|subject/session),data=avgdata.roi.rec,family=binomial)
prcemo <- glmer(recmem ~ PRC + emotion + (1|subject/session),data=avgdata.roi.rec,family=binomial)
phcemo <- glmer(recmem ~ PHC + emotion + (1|subject/session),data=avgdata.roi.rec,family=binomial)
ahippemo <- glmer(recmem ~ wholeHipp.head + emotion + (1|subject/session),data=avgdata.roi.rec,family=binomial)
hippemoxx <- glmer(recmem ~ wholeHipp*emotion + (1|subject/session),data=avgdata.roi.rec,family=binomial)
prcemoxx <- glmer(recmem ~ PRC*emotion + (1|subject/session),data=avgdata.roi.rec,family=binomial)
phcemoxx <- glmer(recmem ~ PHC*emotion + (1|subject/session),data=avgdata.roi.rec,family=binomial)
ahippemoxx <- glmer(recmem ~ wholeHipp.head*emotion + (1|subject/session),data=avgdata.roi.rec,family=binomial)
 

3.2 Amyg combined with others

3.2.1 PRC

3.2.1.1 Additive effects of amygdala & PRC

amygprcemo <- glmer(recmem ~ Amyg*emotion + PRC*emotion + (1|subject/session),data=avgdata.roi.rec,family=binomial)
anova(amygemoxx,amygprcemo) # additive effect of PRC

Data: avgdata.roi.rec
Models:
amygemoxx: recmem ~ Amyg * emotion + (1 | subject/session)
amygprcemo: recmem ~ Amyg * emotion + PRC * emotion + (1 | subject/session)
           Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)   
amygemoxx   6 5114 5152  -2551     5102                           
amygprcemo  8 5107 5157  -2546     5091  10.7      2     0.0047 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(prcemoxx,amygprcemo) # additive effect of Amyg

Data: avgdata.roi.rec
Models:
prcemoxx: recmem ~ PRC * emotion + (1 | subject/session)
amygprcemo: recmem ~ Amyg * emotion + PRC * emotion + (1 | subject/session)
           Df  AIC  BIC logLik deviance Chisq Chi Df  Pr(>Chisq)    
prcemoxx    6 5135 5173  -2562     5123                             
amygprcemo  8 5107 5157  -2546     5091  32.3      2 0.000000099 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Conclusion: Subsequent recollection is predicted well by amygdala and its interaction with emotion, but it improves when PRC is added.

3.2.1.2 Interactive effects of amygdala & PRC

amygprcxxemo <- glmer(recmem ~ Amyg*emotion + PRC*emotion + Amyg*PRC + (1|subject/session),data=avgdata.roi.rec,family=binomial)
anova(amygprcemo,amygprcxxemo)

Data: avgdata.roi.rec
Models:
amygprcemo: recmem ~ Amyg * emotion + PRC * emotion + (1 | subject/session)
amygprcxxemo: recmem ~ Amyg * emotion + PRC * emotion + Amyg * PRC + (1 | subject/session)
             Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
amygprcemo    8 5107 5157  -2546     5091                          
amygprcxxemo  9 5104 5160  -2543     5086  5.48      1      0.019 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

amygprcemoxxx <- glmer(recmem ~ Amyg*emotion*PRC + (1|subject/session),data=avgdata.roi.rec,family=binomial)
anova(amygprcxxemo,amygprcemoxxx)

Data: avgdata.roi.rec
Models:
amygprcxxemo: recmem ~ Amyg * emotion + PRC * emotion + Amyg * PRC + (1 | subject/session)
amygprcemoxxx: recmem ~ Amyg * emotion * PRC + (1 | subject/session)
              Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
amygprcxxemo   9 5104 5160  -2543     5086                        
amygprcemoxxx 10 5106 5169  -2543     5086  0.01      1       0.94

Conclusion: Knowing the interaction of amygdala & PRC improves the model. There is not a three-way interaction with emotion (i.e., the amyg-PRC interaction is similar for negative & neutral).

3.2.1.3 Visualization

trellis.par.set(strip.background=list(col="gray90"))

Note: The default device has been opened to honour attempt to modify trellis settings

visreg(amygprcxxemo, "PRC", scale="response",by="Amyg",cond=list(emotion="Emo"),rug=2,
       xlab="PRC activity", ylab="P(subsequent recollection)",
       breaks=c(-1,0,1),strip.names=c("low amygdala activity",'medium amygdala activity','high amygdala activity'),
       line=list(col=c(emocolor)))

trellis.par.set(strip.background=list(col="gray90"))
visreg(amygprcxxemo, "PRC", scale="response",by="Amyg",cond=list(emotion="Neu"),rug=2,
       xlab="PRC activity", ylab="P(subsequent recollection)",
       breaks=c(-1,0,1),strip.names=c("low amygdala activity",'medium amygdala activity','high amygdala activity'),
       line=list(col=c(neucolor)))

# can't seem to change rug color separately from line color - creating another version to include in overlay
trellis.par.set(strip.background=list(col="gray90"))
visreg(amygprcxxemo, "PRC", scale="response",by="Amyg",cond=list(emotion="Neu"),rug=2,
       xlab="PRC activity", ylab="P(subsequent recollection)",
       breaks=c(-1,0,1),strip.names=c("low amygdala activity",'medium amygdala activity','high amygdala activity'),
       line=list(col=c("gray50")))

3.2.2 Hippocampus

amyghipp <- glmer(recmem ~ Amyg + wholeHipp + (1|subject/session),data=avgdata.roi.rec,family=binomial)
amyghippxx <- glmer(recmem ~ Amyg*wholeHipp  + (1|subject/session),data=avgdata.roi.rec,family=binomial)
anova(amyghipp,amyghippxx)

Data: avgdata.roi.rec
Models:
amyghipp: recmem ~ Amyg + wholeHipp + (1 | subject/session)
amyghippxx: recmem ~ Amyg * wholeHipp + (1 | subject/session)
           Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
amyghipp    5 5136 5167  -2563     5126                        
amyghippxx  6 5138 5176  -2563     5126  0.15      1        0.7

amyghippemo <- glmer(recmem ~ Amyg*emotion + wholeHipp*emotion + (1|subject/session),data=avgdata.roi.rec,family=binomial)
amyghippemoxx <- glmer(recmem ~ Amyg*wholeHipp*emotion + (1|subject/session),data=avgdata.roi.rec,family=binomial)
anova(amyghippemo,amyghippemoxx)

Data: avgdata.roi.rec
Models:
amyghippemo: recmem ~ Amyg * emotion + wholeHipp * emotion + (1 | subject/session)
amyghippemoxx: recmem ~ Amyg * wholeHipp * emotion + (1 | subject/session)
              Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
amyghippemo    8 5110 5161  -2547     5094                        
amyghippemoxx 10 5113 5176  -2547     5093  0.93      2       0.63

visreg(amyghippxx, "wholeHipp", scale="response",by="Amyg",rug=2,xlab="Hippocampal activity", ylab="P(subsequent recollection)",
       breaks=c(-1,0,1),strip.names=c("low amygdala activity",'medium amygdala activity','high amygdala activity'))

# visreg(amyghippxx, "Amyg", scale="response",by="wholeHipp",rug=2,xlab="Amygdala activity", ylab="P(subsequent recollection)",
#        breaks=c(-1,0,1))

The amygdala-hippocampal interaction is not significant, with or without emotion.

4 Regression predicting source memory

avgdata.roi %>%
  filter(!is.na(sourcemem))  %>%
  filter(!is.na(Amyg+PRC+PHC+wholeHipp+wholeHipp.head+wholeHipp.body)) -> avgdata.roi.src # this limits to only recollection trials with intact data

4.1 Each region alone

intonly <- glmer(sourcemem ~  (1|subject/session),data=avgdata.roi.src,family=binomial)
amyg <- glmer(sourcemem ~ Amyg + (1|subject/session),data=avgdata.roi.src,family=binomial)
hipp <- glmer(sourcemem ~ wholeHipp + (1|subject/session),data=avgdata.roi.src,family=binomial)
prc <- glmer(sourcemem ~ PRC + (1|subject/session),data=avgdata.roi.src,family=binomial)
phc <- glmer(sourcemem ~ PHC + (1|subject/session),data=avgdata.roi.src,family=binomial)
anova(intonly,amyg)

Data: avgdata.roi.src
Models:
intonly: sourcemem ~ (1 | subject/session)
amyg: sourcemem ~ Amyg + (1 | subject/session)
        Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
intonly  3 3106 3123  -1550     3100                        
amyg     4 3108 3130  -1550     3100   0.3      1       0.58

anova(intonly,hipp)

Data: avgdata.roi.src
Models:
intonly: sourcemem ~ (1 | subject/session)
hipp: sourcemem ~ wholeHipp + (1 | subject/session)
        Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)   
intonly  3 3106 3123  -1550     3100                           
hipp     4 3097 3120  -1545     3089  10.5      1     0.0012 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

anova(intonly,prc)

Data: avgdata.roi.src
Models:
intonly: sourcemem ~ (1 | subject/session)
prc: sourcemem ~ PRC + (1 | subject/session)
        Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
intonly  3 3106 3123  -1550     3100                        
prc      4 3106 3129  -1549     3098  2.05      1       0.15

anova(intonly,phc)

Data: avgdata.roi.src
Models:
intonly: sourcemem ~ (1 | subject/session)
phc: sourcemem ~ PHC + (1 | subject/session)
        Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)    
intonly  3 3106 3123  -1550     3100                            
phc      4 3087 3110  -1540     3079  20.4      1  0.0000064 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

4.2 Amyg combined with others

4.2.1 Hippocampus

amyghipp <- glmer(sourcemem ~ Amyg + wholeHipp + (1|subject/session),data=avgdata.roi.src,family=binomial)
amyghippxx <- glmer(sourcemem ~ Amyg*wholeHipp  + (1|subject/session),data=avgdata.roi.src,family=binomial)
anova(amyghipp,amyghippxx)

Data: avgdata.roi.src
Models:
amyghipp: sourcemem ~ Amyg + wholeHipp + (1 | subject/session)
amyghippxx: sourcemem ~ Amyg * wholeHipp + (1 | subject/session)
           Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
amyghipp    5 3099 3128  -1545     3089                          
amyghippxx  6 3097 3131  -1542     3085   4.3      1      0.038 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

amyghippemo <- glmer(sourcemem ~ Amyg*emotion + wholeHipp*emotion + (1|subject/session),data=avgdata.roi.src,family=binomial)
amyghippemo2 <- glmer(sourcemem ~ Amyg*emotion + wholeHipp*emotion + Amyg*wholeHipp  + (1|subject/session),data=avgdata.roi.src,family=binomial)
anova(amyghippemo,amyghippemo2)

Data: avgdata.roi.src
Models:
amyghippemo: sourcemem ~ Amyg * emotion + wholeHipp * emotion + (1 | subject/session)
amyghippemo2: sourcemem ~ Amyg * emotion + wholeHipp * emotion + Amyg * wholeHipp + 
amyghippemo2:     (1 | subject/session)
             Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)  
amyghippemo   8 3095 3140  -1539     3079                          
amyghippemo2  9 3092 3143  -1537     3074  4.75      1      0.029 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

amyghippemoxx <- glmer(sourcemem ~ Amyg*wholeHipp*emotion  + (1|subject/session),data=avgdata.roi.src,family=binomial)
anova(amyghippemo2,amyghippemoxx) # does adding the 3-way interaction add anything beyond the simple amyg-hipp intx?

Data: avgdata.roi.src
Models:
amyghippemo2: sourcemem ~ Amyg * emotion + wholeHipp * emotion + Amyg * wholeHipp + 
amyghippemo2:     (1 | subject/session)
amyghippemoxx: sourcemem ~ Amyg * wholeHipp * emotion + (1 | subject/session)
              Df  AIC  BIC logLik deviance Chisq Chi Df Pr(>Chisq)
amyghippemo2   9 3092 3143  -1537     3074                        
amyghippemoxx 10 3092 3150  -1536     3072  1.67      1        0.2

Conclusion: Amygdala & hippocampus interact to support subsequent source memory, but there’s not a three-way interaction with emotion.

4.2.2 Visualization

trellis.par.set(strip.background=list(col="gray90"))

Note: The default device has been opened to honour attempt to modify trellis settings

visreg(amyghippemoxx, "wholeHipp", scale="response",by="Amyg",cond=list(emotion="Emo"),rug=2,
       xlab="Hippocampal activity", ylab="P(subsequent source)",
       breaks=c(-1,0,1),strip.names=c("low amygdala activity",'medium amygdala activity','high amygdala activity'),
       line=list(col=emocolor))

trellis.par.set(strip.background=list(col="gray90"))
visreg(amyghippemoxx, "wholeHipp", scale="response",by="Amyg",cond=list(emotion="Neu"),rug=2,
       xlab="Hippocampal activity", ylab="P(subsequent source)",
       breaks=c(-1,0,1),strip.names=c("low amygdala activity",'medium amygdala activity','high amygdala activity'),
       line=list(col=neucolor))

# can't seem to change rug color separately from line color - creating another version to include in overlay
trellis.par.set(strip.background=list(col="gray90"))
visreg(amyghippemoxx, "wholeHipp", scale="response",by="Amyg",cond=list(emotion="Neu"),rug=2,
       xlab="Hippocampal activity", ylab="P(subsequent source)",
       breaks=c(-1,0,1),strip.names=c("low amygdala activity",'medium amygdala activity','high amygdala activity'),
       line=list(col="gray50"))