Binding fMRI: Analyses of feature encoding and integration

The analyses presented here test subsequent memory effects during encoding of a multi-feature event (from the GitHub repository: https://github.com/memobc/paper-bindingfmri).

Mean single trial beta estimates at i) onset and ii) offset of each visual display (onset N and offset N are 6s apart) were extracted from each of 204 ROIs across the cortex and medial temporal lobe (see ~/data/bindingfmri_wholebrain_ROIs_info.csv). Onsets are used to test the relationship between brain activity and memory for the upcoming event, whereas offsets are used to test the relationship between brain activity and memory for the preceding event.

Each event includes 3 associations encoded in parallel – an object with an associated i) color, ii) sound, and iii) scene location.

Subsequent memory detail is calculated as the number of features (0-3) later correctly recalled. Feature-specific estimates of subsequent memory are also analyzed (correct, 1, vs. incorrect, 0). The relationship between activity (each ROI’s beta series) and subsequent memory is estimated using linear mixed effects (LME) models with the following pipeline:

1. Data are cleaned by removing any trials with outlying betas (z > |4|), calculated within subject per ROI
   
2. Each LME model is initially specified with a maximal random effects structure (including variable slopes for all fixed effects, as well as intercepts, by subject), but to avoid convergence problems and singular fits:
   

1. Random effects correlations are constrained to zero for all models
   
2. The random effect structure of models with > 2 fixed effects is simplified by iteratively removing slopes that do not improve the model fit

3. The significance of the fixed effects (slope > or < 0?) is tested using t-tests with Satterthwaite’s approximation method

Behavioral Data

### Load in behavioral data for all 27 subjects:
behavData <- read.csv('../data/AllSubjects_bindingfMRI-behavior.csv', header = TRUE)
behav_subs <- unique(behavData$SubID)
cat('Number of subjects =',length(behav_subs))

Number of subjects = 27

# ** note. behavioral data is already sorted by subject and encoding trial onset ** #
# Create additional columns for subsequent memory measures of interest -- Sound, Color, & Scene features ----------------------------
behavData$Sound <- 0
behavData$Sound[behavData$SoundCorrect == 1 & behavData$SoundConfidence == 1] <- 1 # sound only sucessful if correct *and* high confidence
# define 'success' thresholds for color and scene trials -> 75% probability within von Mises for 27 subjects aggregate data
cCutOff = 42
sCutOff = 24
# feature accuracy
behavData$Color <- 0
behavData$Color[behavData$ColAbsError <=cCutOff] <- 1 #correct if error (how far from target) less than or equal to threshold
behavData$Scene <- 0
behavData$Scene[behavData$SceAbsError <=sCutOff] <- 1 #correct if error (how far from target) less than or equal to threshold
# overall memory detail (accounting for feature success - 0-3)
behavData$MemoryDetail <- behavData$Sound + behavData$Color + behavData$Scene
#### separate trial-unique conditions based on feature combinations recalled --------------------
## 1 feature
behavData$Color.Only   <- 0
behavData$Sound.Only <- 0
behavData$Scene.Only   <- 0
behavData$Color.Only[behavData$Color == 1 & behavData$Sound == 0 & behavData$Scene == 0]  <- 1
behavData$Sound.Only[behavData$Color == 0 & behavData$Sound == 1 & behavData$Scene == 0]  <- 1
behavData$Scene.Only[behavData$Color == 0 & behavData$Sound == 0 & behavData$Scene == 1]  <- 1
## 2 features
behavData$Color.Scene   <- 0
behavData$Sound.Color <- 0
behavData$Scene.Sound <- 0
behavData$Color.Scene[behavData$Color == 1 & behavData$Sound == 0 & behavData$Scene == 1] <- 1
behavData$Sound.Color[behavData$Color == 1 & behavData$Sound == 1 & behavData$Scene == 0] <- 1
behavData$Scene.Sound[behavData$Color == 0 & behavData$Sound == 1 & behavData$Scene == 1] <- 1
## 3 features
behavData$Color.Scene.Sound  <- 0
behavData$Color.Scene.Sound[behavData$Color == 1 & behavData$Sound == 1 & behavData$Scene == 1] <- 1
### run a check to make sure no trial has been allocated to more than one combination:
nPresent <- rowSums(behavData[,c("Color.Only","Sound.Only","Scene.Only","Color.Scene","Sound.Color","Scene.Sound","Color.Scene.Sound")])
if (max(nPresent) < 1) {
  cat('ERROR: trial(s) assigned to more than one combination\n')
}
# double check encoding data is sorted correctly (by subject, then encoding onset (by run and time))
behavData <- behavData[order(behavData$TotalEncodingTrial),]
nTotal <- nrow(behavData)
# print percentage correct for each feature
feature.accuracy <- behavData[,c("SubID","StudyTrial","Color","Sound","Scene")]
sub.behav <- feature.accuracy %>% gather(Feature,Accuracy,Color:Scene) %>%
                          group_by(SubID, Feature) %>%
                          summarise(M = mean(Accuracy))
kable(sub.behav %>% group_by(Feature) %>%
                          summarise(Memory = mean(M), SE = se(M)),format="pandoc",
                                    caption="Feature Accuracy (proportion of trials labeled as correct)")

Table: Feature Accuracy (proportion of trials labeled as correct)

Feature       Memory          SE
--------  ----------  ----------
Color      0.6604424   0.0240285
Scene      0.6606996   0.0331897
Sound      0.5626543   0.0338430

### print total number of trials that uniquely fall into each possible combination:
cat('\nNone recalled trials =', nrow(subset(behavData,MemoryDetail == 0)),'(',round(((nrow(subset(behavData,MemoryDetail == 0)))/nTotal)*100,2),'%)\n')

None recalled trials = 456 ( 12.03 %)

cat('Color only trials =', sum(behavData$Color.Only),'(',round(((sum(behavData$Color.Only))/nTotal)*100,2),'%)\n')

Color only trials = 384 ( 10.13 %)

cat('Scene only trials =', sum(behavData$Scene.Only),'(',round(((sum(behavData$Scene.Only))/nTotal)*100,2),'%)\n')

Scene only trials = 268 ( 7.07 %)

cat('Sound only trials =', sum(behavData$Sound.Only),'(',round(((sum(behavData$Sound.Only))/nTotal)*100,2),'%)\n')

Sound only trials = 162 ( 4.27 %)

cat('Color & Scene only trials =', sum(behavData$Color.Scene),'(',round(((sum(behavData$Color.Scene))/nTotal)*100,2),'%)\n')

Color & Scene only trials = 542 ( 14.29 %)

cat('Scene & Sound only trials =', sum(behavData$Scene.Sound),'(',round(((sum(behavData$Scene.Sound))/nTotal)*100,2),'%)\n')

Scene & Sound only trials = 397 ( 10.47 %)

cat('Sound & Color only trials =', sum(behavData$Sound.Color),'(',round(((sum(behavData$Sound.Color))/nTotal)*100,2),'%)\n')

Sound & Color only trials = 277 ( 7.3 %)

cat('Color & Scene & Sound trials =', sum(behavData$Color.Scene.Sound),'(',round(((sum(behavData$Color.Scene.Sound))/nTotal)*100,2),'%)\n')

Color & Scene & Sound trials = 1306 ( 34.44 %)

Check for behavioral (memory) trial dependencies

Check if encoding success at trial N is related to encoding success on trial N+1 within subjects

behavData$EncodingEvent <- rep(1:24, nrow(behavData)/24) #event number within run
# overall detail, then each feature separately
features = c('MemoryDetail','Color','Sound','Scene')
for (f in 1:length(features)) {
  my.col <- which(colnames(behavData) %in% features[f])
  correlations <- data.frame(array(0,c(length(behav_subs),2)))
  colnames(correlations) <- c('SubID','Z')
  for (s in 1:length(behav_subs)){
    subData <- subset(behavData, SubID == behav_subs[s])
    correlations$SubID[s] <- behav_subs[s]
    trialn  <- subset(subData, EncodingEvent < 24) #all but last event of run for N vector (1:23) correlated with ...
    trialn1 <- subset(subData, EncodingEvent > 1)  #all but first event of the run for N+1 vector (2:24)
    correlations$Z[s] <- fisherz(cor(trialn[,my.col],trialn1[,my.col]))
  }
  
  R = fisherz2r(mean(correlations$Z))
  SEmean = se(correlations$Z)
  cat('Mean correlation (r) between trial N and N+1 for',features[f],'within subject =',R,', SE =',SEmean,'\n')
}

Mean correlation (r) between trial N and N+1 for MemoryDetail within subject = 0.1005285 , SE = 0.02729053 
Mean correlation (r) between trial N and N+1 for Color within subject = 0.01079016 , SE = 0.01759081 
Mean correlation (r) between trial N and N+1 for Sound within subject = 0.1127844 , SE = 0.03943163 
Mean correlation (r) between trial N and N+1 for Scene within subject = 0.06142598 , SE = 0.02007916 

Get Onset Data

Onsets are used to test the relationship between activity and memory for the upcoming event

############################# single trial activity data ##############################
# get subjects' data
subjects <- list.files(path='../data/single-trial-betas/event-onset/',full.names=TRUE, recursive = TRUE, pattern = 'onset-wholebrain')
allData  = do.call(rbind, lapply(subjects, function(x) {read.csv(x, header = TRUE)}))
allData$SubID=as.factor(allData$SubID)
allData$ROI <- as.factor(allData$ROI)
allData <- allData[,-c(5)] # remove vox numbers - not needed here
rois = levels(allData$ROI)
cat('Number of subjects =', length(unique(allData$SubID)),'\n')

Number of subjects = 27 

cat('Removing influential onset betas, |z| >', outlier_value,'\n')

Removing influential onset betas, |z| > 4 

# ******************************************************** #
### add z score of betas within region and subject to remove outliers:
allData <- allData %>%
               group_by(SubID, ROI) %>%
               mutate(Z = zscore(MeanBeta))
# remove betas that are > XSD -- print this number
allData.clean <- allData
allData.clean$MeanBeta[abs(allData.clean$Z) > outlier_value] <- NA  #mark beta as NA if outlier
allData.clean <- allData.clean[,-c(6)] # remove Z values  
##########################################################
#spread data so each of the 204 ROIs is in its own column (to align with behavioral data)
myBetas = data.frame(spread(allData.clean, ROI, MeanBeta))
myBetas$TotalEncodingTrial <- 1:nTotal
# remove all trials where there is at least one NA (at least one of the ROIs is an outlier on that trial) - 
# to keep the trials analyzed the same across regions:
idxBetas <- complete.cases(myBetas)
myBetas.clean <- myBetas[idxBetas,]
cat('Total number of onsets removed =',sum(!idxBetas),'out of',nTotal,'\n\n')

Total number of onsets removed = 108 out of 3792 

##### merge cleaned betas with behavioral data, removing corresponding behavioral trials:
behavData.clean <- behavData[idxBetas,]
testData <- merge.data.frame(myBetas.clean,behavData.clean,by="TotalEncodingTrial")
# remove duplicate columns for subject, run, and study trial
testData <- testData[,-c(2,3,4)]
colnames(testData)[colnames(testData) == 'SubID.y'] <- 'SubID'
colnames(testData)[colnames(testData) == 'Run.y'] <- 'Run'
subjects <- unique(testData$SubID)
onset.data  <- testData  #store onset betas and behavioral data for upcoming trial
cat('Cleaned and merged onset betas and behavioral trials\n')

Cleaned and merged onset betas and behavioral trials

Get Offset Data

Offsets are used to test the relationship between activity and memory for the preceding event

############################# single trial activity data ##############################
# get subjects' data
subjects <- list.files(path='../data/single-trial-betas/event-offset/',full.names=TRUE, recursive = TRUE, pattern = 'offset-wholebrain')
allData  = do.call(rbind, lapply(subjects, function(x) {read.csv(x, header = TRUE)}))
allData$SubID=as.factor(allData$SubID)
allData$ROI <- as.factor(allData$ROI)
allData <- allData[,-c(5)] # remove vox numbers - not needed here
rois = levels(allData$ROI)
cat('Number of subjects =', length(unique(allData$SubID)),'\n')

Number of subjects = 27 

cat('Removing influential offset betas, |z| >', outlier_value,'\n')

Removing influential offset betas, |z| > 4 

# ******************************************************** #
### add z score of betas within region and subject to remove outliers:
allData <- allData %>%
               group_by(SubID, ROI) %>%
               mutate(Z = zscore(MeanBeta))
# remove betas that are > XSD -- print this number
allData.clean <- allData
allData.clean$MeanBeta[abs(allData.clean$Z) > outlier_value] <- NA  #mark beta as NA if outlier
allData.clean <- allData.clean[,-c(6)] # remove Z values  
##########################################################
#spread data so each of the 204 ROIs is in its own column (to align with behavioral data)
myBetas = data.frame(spread(allData.clean, ROI, MeanBeta))
myBetas$TotalEncodingTrial <- 1:nTotal
# remove all trials where there is at least one NA (at least one of the ROIs is an outlier on that trial) - 
# to keep the trials analyzed the same across regions:
idxBetas <- complete.cases(myBetas)
myBetas.clean <- myBetas[idxBetas,]
cat('Total number of offsets removed =',sum(!idxBetas),'out of',nTotal,'\n\n')

Total number of offsets removed = 113 out of 3792 

##### merge cleaned betas with behavioral data, removing corresponding behavioral trials:
behavData.clean <- behavData[idxBetas,]
testData <- merge.data.frame(myBetas.clean,behavData.clean,by="TotalEncodingTrial")
# remove duplicate columns for subject, run, and study trial
testData <- testData[,-c(2,3,4)]
colnames(testData)[colnames(testData) == 'SubID.y'] <- 'SubID'
colnames(testData)[colnames(testData) == 'Run.y'] <- 'Run'
subjects <- unique(testData$SubID)
offset.data  <- testData  #store offset betas and behavioral data for preceding trial
cat('Cleaned and merged offset betas and behavioral trials\n')

Cleaned and merged offset betas and behavioral trials

Correlate onset/offset betas

1. Separated by 1s ITI: First remove any onsets for which offset N-1 is not available (1s apart) within subject and run, and vice versa. Correlation by roi and subject, then take the overall average within subject, and at the group level.
   
2. Separated by 6s event: Correlation between same-trial onsets and offsets.

######### A) #########
## onset N to offset N-1
onset.data.cor <- onset.data
offset.data.cor <- offset.data
trials_keep <- c()
for (s in 1:length(subjects)) {
  for (r in unique(onset.data.cor$Run[onset.data.cor$SubID == subjects[s]])) {
    subData.onset <- subset(onset.data.cor, SubID == subjects[s] & Run == r)
    subData.offset <- subset(offset.data.cor, SubID == subjects[s] & Run == r)
    for (o in 1:nrow(subData.onset)) {
      if ((subData.onset$StudyTrial[o]-1) %in% subData.offset$StudyTrial) {  #if the trial before current onset is present in offset list, keep
        trials_keep <- c(trials_keep, subData.onset$TotalEncodingTrial[o])
      }
    } 
  }
}
onset.data.cor <- onset.data.cor[onset.data.cor$TotalEncodingTrial %in% trials_keep,]
trials_keep <- c()
for (s in 1:length(subjects)) {
  for (r in unique(offset.data.cor$Run[offset.data.cor$SubID == subjects[s]])) {
    subData.onset <- subset(onset.data.cor, SubID == subjects[s] & Run == r)
    subData.offset <- subset(offset.data.cor, SubID == subjects[s] & Run == r)
    for (o in 1:nrow(subData.offset)) {
      if ((subData.offset$StudyTrial[o]+1) %in% subData.onset$StudyTrial) { #if the trial after current offset is present in onset list, keep
        trials_keep <- c(trials_keep, subData.offset$TotalEncodingTrial[o])
      }
    } 
  }
}
offset.data.cor <- offset.data.cor[offset.data.cor$TotalEncodingTrial %in% trials_keep,]
# now run per ROI and subject
correlations <- data.frame(array(0,c((length(subjects)*length(rois)),3)))
colnames(correlations) <- c('SubID','ROI','Z')
row = 0
for (r in 1:length(rois)) {
 for (s in 1:length(subjects)){
   row = row +1
   sub.onset  <- onset.data.cor[onset.data.cor$SubID == subjects[s],colnames(onset.data.cor) == rois[r]]
   sub.offset <- offset.data.cor[offset.data.cor$SubID == subjects[s],colnames(offset.data.cor) == rois[r]]
   correlations$SubID[row] <- subjects[s]
   correlations$ROI[row]   <- rois[r]
   correlations$Z[row]     <- fisherz(cor(sub.onset,sub.offset))
 }  
}
  
sub.summary <- correlations %>% group_by(SubID) %>%
                   summarise(meanZ = mean(Z))
R = fisherz2r(mean(sub.summary$meanZ))
SEmean = se(sub.summary$meanZ)
cat('Mean correlation (r) between onset N and offset N-1 (separated by 1s ITI) within subject and ROI =',R,', SE =',SEmean,'\n')

Mean correlation (r) between onset N and offset N-1 (separated by 1s ITI) within subject and ROI = 0.8802484 , SE = 0.02030568 

######### B) #########
## onset N to offset N
shared.trials <- intersect(onset.data$TotalEncodingTrial,offset.data$TotalEncodingTrial)
onset.data.cor  <- onset.data[onset.data$TotalEncodingTrial %in% shared.trials,]
offset.data.cor <- offset.data[offset.data$TotalEncodingTrial %in% shared.trials,]
# now run per ROI and subject
correlations <- data.frame(array(0,c((length(subjects)*length(rois)),3)))
colnames(correlations) <- c('SubID','ROI','Z')
row = 0
for (r in 1:length(rois)) {
 for (s in 1:length(subjects)){
   row = row +1
   sub.onset  <- onset.data.cor[onset.data.cor$SubID == subjects[s],colnames(onset.data.cor) == rois[r]]
   sub.offset <- offset.data.cor[offset.data.cor$SubID == subjects[s],colnames(offset.data.cor) == rois[r]]
   correlations$SubID[row] <- subjects[s]
   correlations$ROI[row]   <- rois[r]
   correlations$Z[row]     <- fisherz(cor(sub.onset,sub.offset))
 }  
}
  
sub.summary <- correlations %>% group_by(SubID) %>%
                   summarise(meanZ = mean(Z))
R = fisherz2r(mean(sub.summary$meanZ))
SEmean = se(sub.summary$meanZ)
cat('Mean correlation (r) between onset N and offset N (separated by 6s trial) within subject and ROI =',R,', SE =',SEmean,'\n')

Mean correlation (r) between onset N and offset N (separated by 6s trial) within subject and ROI = 0.1379871 , SE = 0.01122709 

Whole brain analyses

Note - each chunk will take a few minutes to run.

Memory detail

Where across the brain does activity at linearly scale with memory quantity (0,1,2,3) in the upcoming trial (using onset) or the preceding trial (using offset)?

curFile = 'memorydetail_effects_onsetoffset.csv'
### set up data frame for mixed effect model results:
mixed.data <- data.frame(array(0,c(length(rois)*2,8)))
colnames(mixed.data) <- c('ROI','Phase','Feature','Beta','SE','T','Df','Pvalue')
row = 0
### loop over onset (~ upcoming event memory) and offset (~ preceding event memory)
for (phase in 1:2) {
  
  # gather data so behavioral variables are replicated over each ROI (now as a single column factor)
  if (phase == 1) {
    model.data <- gather(onset.data, ROI, Beta, rois)
    time = 'Onset'
  } else if (phase == 2) {
    model.data <- gather(offset.data, ROI, Beta, rois)
    time = 'Offset'
  }
  
  # ******** run linear mixed effects model within each ROI and store results **********
  for (r in 1:length(rois)) {
    
    thisROI <- as.character(rois[r])
    
    roiData <- subset(model.data, ROI == thisROI)
    roiData$MemoryDetail <- scale(roiData$MemoryDetail, scale = FALSE) # mean-center memory predictor
    
    # run lmer and get fixed effect stats
    full.model <- suppressMessages(lmer(Beta ~ MemoryDetail + 
                                          (1 + MemoryDetail||SubID), data = roiData))
    stats <- summary(full.model)$coefficients
    nIV <- nrow(stats) - 1 # number of fixed effects, discounting intercept
  
    ### add fixed effects and p values to data frame
    mixed.data$ROI[(row+1):(row+nIV)]      <- thisROI
    mixed.data$Phase[(row+1):(row+nIV)]    <- time
    mixed.data$Feature[(row+1):(row+nIV)]  <- row.names(stats)[2:(nIV+1)]
    mixed.data[(row+1):(row+nIV),c("Beta","SE","T","Df","Pvalue")] <- stats[2:(nIV+1),
                                                                           c("Estimate","Std. Error","t value","df","Pr(>|t|)")]
    row = row + nIV
  }
}
#### FDR-correct all p-values for multiple comparisons
mixed.data$Pfdr <- 1
p.values <- mixed.data$Pvalue
mixed.data$Pfdr <- p.adjust(p.values, method = "fdr")
## save out parameter estimates and p values as csv file for creating a .nii map of results
write.csv(mixed.data, file = curFile, row.names = FALSE)
cat('Results saved in',curFile)

Results saved in memorydetail_effects_onsetoffset.csv

Control analysis 1

This analysis is the same as above, but we are checking that the results look similar if we explicitly include memory detail on trial N-1 or N+1 as a covariate when testing the relationship between onset/offset N and memory detail on trial N.

curFile = 'memorydetail_effects_covariate.csv'
### set up data frame for mixed effect model results:
mixed.data <- data.frame(array(0,c(length(rois)*2,8)))
colnames(mixed.data) <- c('ROI','Phase','Feature','Beta','SE','T','Df','Pvalue')
row = 0
### loop over onset and offset
for (phase in 1:2) {
  
  trials.n <- c()
  trials.n1 <- c()
    
  # gather data so behavioral variables are replicated over each ROI (now as a single column factor)
  if (phase == 1) {
    model.data <- onset.data
    time = 'Onset'
    # now loop through and only retain trial N if trial N-1 (immediately adjacent to onset N) is present (within subject and run)
    for (n in 2:nrow(model.data)) {  #can't include first trial in N vector
      if (model.data$EncodingEvent[n-1] == model.data$EncodingEvent[n]-1) {
           trials.n  <- c(trials.n, n)
           trials.n1 <- c(trials.n1, n-1)
      }
    }
  } else if (phase == 2) {
    model.data <- offset.data
    time = 'Offset'
    # now loop through and only retain trial N if trial N+1 (immediately adjacent to offset N) is present (within subject and run)
    for (n in 1:(nrow(model.data)-1)) {  #can't include final trial N vector
      if (model.data$EncodingEvent[n+1] == model.data$EncodingEvent[n]+1) {
           trials.n  <- c(trials.n, n)
           trials.n1 <- c(trials.n1, n+1)
      }
    }
  }
  
  # format data for MemoryDetail and MemoryDetailN1 predictors
  new.data    <- model.data[trials.n,]
  new.data.n1 <- model.data[trials.n1,]
  new.data$MemoryDetailN1  <- new.data.n1$MemoryDetail
  model.data <- gather(new.data, ROI, Beta, rois)
  
  # ******** run linear mixed effects within each ROI and store results **********
  for (r in 1:length(rois)) {
    
    thisROI <- as.character(rois[r])
    
    roiData <- subset(model.data, ROI == thisROI)
    roiData$MemoryDetail   <- scale(roiData$MemoryDetail, scale = FALSE) # mean-center predictor
    roiData$MemoryDetailN1 <- scale(roiData$MemoryDetailN1, scale = FALSE) # mean-center predictor
   
    # run lmer and get fixed effect stats
    full.model <- suppressMessages(lmer(Beta ~ MemoryDetail + MemoryDetailN1 +
                                          (1 + MemoryDetail + MemoryDetailN1||SubID), data = roiData))
    stats <- summary(full.model)$coefficients
    nIV <- nrow(stats) - 2 # number of fixed effect predictors - but don't need to store covariate results or intercept
  
    ### add fixed effects and p values to data matrix
    mixed.data$ROI[(row+1):(row+nIV)]      <- thisROI
    mixed.data$Phase[(row+1):(row+nIV)]    <- time
    mixed.data$Feature[(row+1):(row+nIV)]  <- row.names(stats)[2:(nIV+1)]
    mixed.data[(row+1):(row+nIV),c("Beta","SE","T","Df","Pvalue")] <- stats[2:(nIV+1),
                                                                           c("Estimate","Std. Error","t value","df","Pr(>|t|)")]
    row = row + nIV
  }
}
#### FDR-correct all p-values for multiple comparisons
mixed.data$Pfdr <- 1
p.values <- mixed.data$Pvalue
mixed.data$Pfdr <- p.adjust(p.values, method = "fdr")
## save out parameter estimates and p values as csv file for creating a .nii map of results
write.csv(mixed.data, file = curFile, row.names = FALSE)
cat('Results saved in',curFile)

Results saved in memorydetail_effects_covariate.csv

Control analysis 2

Instead of modeling time points at the beginning and end of each event, this analysis models each event (from onset to offset) as a 6s epoch.

############################# single trial activity data ##############################
# get subjects' data
subjects <- list.files(path='../data/single-trial-betas/6s-event/',full.names=TRUE, recursive = TRUE, pattern = 'epoch-wholebrain')
allData  = do.call(rbind, lapply(subjects, function(x) {read.csv(x, header = TRUE)}))
allData$SubID=as.factor(allData$SubID)
allData$ROI <- as.factor(allData$ROI)
allData <- allData[,-c(5)] # remove vox numbers - not needed here
rois = levels(allData$ROI)
cat('Number of subjects =', length(unique(allData$SubID)),'\n')

Number of subjects = 27 

cat('Removing influential event betas, |z| >', outlier_value,'\n')

Removing influential event betas, |z| > 4 

# ******************************************************** #
### add z score of betas within region and subject to remove outliers:
allData <- allData %>%
               group_by(SubID, ROI) %>%
               mutate(Z = zscore(MeanBeta))
# remove betas that are > XSD -- print this number
allData.clean <- allData
allData.clean$MeanBeta[abs(allData.clean$Z) > outlier_value] <- NA  #mark beta as NA if outlier
allData.clean <- allData.clean[,-c(6)] # remove Z values  
##########################################################
#spread data so each of the 204 ROIs is in its own column (to align with behavioral data)
myBetas = data.frame(spread(allData.clean, ROI, MeanBeta))
myBetas$TotalEncodingTrial <- 1:nTotal
# remove all trials where there is at least one NA (at least one of the ROIs is an outlier on that trial) - 
# to keep the trials analyzed the same across regions:
idxBetas <- complete.cases(myBetas)
myBetas.clean <- myBetas[idxBetas,]
cat('Total number of events removed =',sum(!idxBetas),'out of',nTotal,'\n\n')

Total number of events removed = 59 out of 3792 

##### merge cleaned betas with behavioral data, removing corresponding behavioral trials:
behavData.clean <- behavData[idxBetas,]
testData <- merge.data.frame(myBetas.clean,behavData.clean,by="TotalEncodingTrial")
# remove duplicate columns for subject, run, and study trial
testData <- testData[,-c(2,3,4)]
colnames(testData)[colnames(testData) == 'SubID.y'] <- 'SubID'
colnames(testData)[colnames(testData) == 'Run.y'] <- 'Run'
subjects <- unique(testData$SubID)
epoch.data  <- testData  #store event betas and behavioral data for current trial
cat('Cleaned and merged event betas and behavioral trials\n')

Cleaned and merged event betas and behavioral trials

curFile = 'memorydetail_effects_eventepoch.csv'
### set up data frame for mixed effect model results:
mixed.data <- data.frame(array(0,c(length(rois),8)))
colnames(mixed.data) <- c('ROI','Phase','Feature','Beta','SE','T','Df','Pvalue')
row = 0
model.data <- gather(epoch.data, ROI, Beta, rois)
time = 'Epoch'
# ******** run linear mixed effects within each ROI and store results **********
for (r in 1:length(rois)) {
  
  thisROI <- as.character(rois[r])
  
  roiData <- subset(model.data, ROI == thisROI)
  roiData$MemoryDetail <- scale(roiData$MemoryDetail, scale = FALSE) # mean-center predictor
  
  # run lmer and get fixed effect stats
  full.model <- suppressMessages(lmer(Beta ~ MemoryDetail + 
                                        (1 + MemoryDetail||SubID), data = roiData))
  stats <- summary(full.model)$coefficients
  nIV <- nrow(stats) - 1 # number of fixed effects, discounting intercept
    
  ### add fixed effects and p values to data matrix
  mixed.data$ROI[(row+1):(row+nIV)]      <- thisROI
  mixed.data$Phase[(row+1):(row+nIV)]    <- time
  mixed.data$Feature[(row+1):(row+nIV)]  <- row.names(stats)[2:(nIV+1)]
  mixed.data[(row+1):(row+nIV),c("Beta","SE","T","Df","Pvalue")] <- stats[2:(nIV+1),
                                                                         c("Estimate","Std. Error","t value","df","Pr(>|t|)")]
  row = row + nIV
}
#### FDR-correct all p-values for multiple comparisons
mixed.data$Pfdr <- 1
p.values <- mixed.data$Pvalue
mixed.data$Pfdr <- p.adjust(p.values, method = "fdr")
## save out parameter estimates and p values as csv file for creating a .nii map of results
write.csv(mixed.data, file = curFile, row.names = FALSE)
cat('Results saved in',curFile)

Results saved in memorydetail_effects_eventepoch.csv

Control analysis 3

Here, we demonstrate the similarity of our onset- and offset-based analyses to a model where we disregard any distinction between correlated offset N and onset N+1 signals, modeling them together as a single 1s ITI. ITI activity is predicted by subsequent memory for the upcoming or preceding event.

############################# single trial activity data ##############################
# get subjects' data
subjects <- list.files(path='../data/single-trial-betas/1s-ITI/',full.names=TRUE, recursive = TRUE, pattern = 'ITI-wholebrain')
allData  = do.call(rbind, lapply(subjects, function(x) {read.csv(x, header = TRUE)}))
allData$SubID=as.factor(allData$SubID)
allData$ROI <- as.factor(allData$ROI)
allData <- allData[,-c(5)] # remove vox numbers - not needed here
rois = levels(allData$ROI)
cat('Number of subjects =', length(unique(allData$SubID)),'\n')

Number of subjects = 27 

cat('Removing influential ITI betas, |z| >', outlier_value,'\n')

Removing influential ITI betas, |z| > 4 

# ******************************************************** #
### add z score of betas within region and subject to remove outlier boundaries:
allData <- allData %>%
               group_by(SubID, ROI) %>%
               mutate(Z = zscore(MeanBeta))
# remove betas that are > XSD -- print this number
allData.clean <- allData
allData.clean$MeanBeta[abs(allData.clean$Z) > outlier_value] <- NA  #mark beta as NA if outlier
allData.clean <- allData.clean[,-c(6)] # remove Z values
##########################################################
#spread data so each of the 204 ROIs is in its own column
myBetas = data.frame(spread(allData.clean, ROI, MeanBeta))
myBetas$TotalEvent <- 1:nrow(myBetas)
forwardBehav <- behavData[behavData$EncodingEvent > 1,]  #if predicting forward from ITI, we can use all but behavioral trial 1 (first trial per run)
forwardBehav$TotalEvent <- 1:nrow(forwardBehav)
backwardBehav <- behavData[behavData$EncodingEvent < 24,] #if predicting backward from ITI, we can use all but behavioral trial 24 (final trial per run)
backwardBehav$TotalEvent <- 1:nrow(backwardBehav)
### CLEAN DATA and merge with behavior
# remove all trials where there is at least one NA (at least one of the ROIs is an outlier on that trial):
idxBetas <- complete.cases(myBetas)
myBetas.clean <- myBetas[idxBetas,]
forwardBehav.clean <- forwardBehav[idxBetas,]
backwardBehav.clean <- backwardBehav[idxBetas,]
cat('Total number of ITIs removed =',nrow(myBetas)-nrow(myBetas.clean),'out of',nrow(myBetas),'\n\n')

Total number of ITIs removed = 98 out of 3634 

forward.data <- merge.data.frame(myBetas.clean,forwardBehav.clean,by="TotalEvent")
backward.data <- merge.data.frame(myBetas.clean,backwardBehav.clean,by="TotalEvent")
# remove duplicate columns for subject, study trial, and run ----------
forward.data <- forward.data[,-c(2,3,4)]
colnames(forward.data)[colnames(forward.data) == 'SubID.y'] <- 'SubID'
colnames(forward.data)[colnames(forward.data) == 'Run.y'] <- 'Run'
backward.data <- backward.data[,-c(2,3,4)]
colnames(backward.data)[colnames(backward.data) == 'SubID.y'] <- 'SubID'
colnames(backward.data)[colnames(backward.data) == 'Run.y'] <- 'Run'
subjects <- unique(forward.data$SubID)
cat('Cleaned and merged ITI betas and behavioral trials\n')

Cleaned and merged ITI betas and behavioral trials

curFile = 'memorydetail_effects_ITI.csv'
### set up data frame for mixed effect model results:
mixed.data <- data.frame(array(0,c(length(rois)*2,8)))
colnames(mixed.data) <- c('ROI','Phase','Feature','Beta','SE','T','Df','Pvalue')
row = 0
### loop over onset and offset
for (phase in 1:2) {
  
  # gather data so behavioral variables are replicated over each ROI (now as a single column factor)
  if (phase == 1) {
    model.data <- gather(forward.data, ROI, Beta, rois)
    time = 'Upcoming'
  } else if (phase == 2) {
    model.data <- gather(backward.data, ROI, Beta, rois)
    time = 'Preceding'
  }
  
  # ******** run linear mixed effects within each ROI and store results **********
  for (r in 1:length(rois)) {
    
    thisROI <- as.character(rois[r])
    
    roiData <- subset(model.data, ROI == thisROI)
    roiData$MemoryDetail <- scale(roiData$MemoryDetail, scale = FALSE) # mean-center predictor
    
    # run lmer and get fixed effect stats
    full.model <- suppressMessages(lmer(Beta ~ MemoryDetail + 
                                          (1 + MemoryDetail||SubID), data = roiData))
    stats <- summary(full.model)$coefficients
    nIV <- nrow(stats) - 1 # number of fixed effects, discounting intercept
    
    ### add fixed effects and p values to data matrix
    mixed.data$ROI[(row+1):(row+nIV)]      <- thisROI
    mixed.data$Phase[(row+1):(row+nIV)]    <- time
    mixed.data$Feature[(row+1):(row+nIV)]  <- row.names(stats)[2:(nIV+1)]
    mixed.data[(row+1):(row+nIV),c("Beta","SE","T","Df","Pvalue")] <- stats[2:(nIV+1),
                                                                             c("Estimate","Std. Error","t value","df","Pr(>|t|)")]
    row = row + nIV
  }
}
#### FDR-correct p-values for multiple comparisons
mixed.data$Pfdr <- 1
p.values <- mixed.data$Pvalue
mixed.data$Pfdr <- p.adjust(p.values, method = "fdr")
## save out parameter estimates and p values as csv file for creating a .nii map of results
write.csv(mixed.data, file = curFile, row.names = FALSE)
cat('Results saved in',curFile)

Results saved in memorydetail_effects_ITI.csv

Unique feature encoding

After accounting for covariance with the other two features in the model, where is onset/offset activity uniquely predicted by successful color, scene, or sound encoding?

curFile = 'feature_effects_onsetoffset.csv'
### set up data frame for mixed effect model results:
mixed.data <- data.frame(array(0,c(length(rois)*6,8)))  # *6 (3 features x 2 times)
colnames(mixed.data) <- c('ROI','Phase','Feature','Beta','SE','T','Df','Pvalue')
row = 0
### loop over onset and offset
for (phase in 1:2) {
  
  # gather data so behavioral variables are replicated over each ROI (now as a single column factor)
  if (phase == 1) {
    model.data <- gather(onset.data, ROI, Beta, rois)
    time = 'Onset'
  } else if (phase == 2) {
    model.data <- gather(offset.data, ROI, Beta, rois)
    time = 'Offset'
  }
  
  # ******** run linear mixed effects within each ROI and store results for each feature fixed effect **********
  for (r in 1:length(rois)) {
    
    thisROI <- as.character(rois[r])
    roiData <- subset(model.data, ROI == thisROI)
    roiData$Sound  <- scale(roiData$Sound, scale = FALSE) # mean-center
    roiData$Color  <- scale(roiData$Color, scale = FALSE) # mean-center
    roiData$Scene  <- scale(roiData$Scene, scale = FALSE) # mean-center
    
    # run lmer, simplify random effect structure as > 2 fixed effects, and get fixed effect stats
    full.model <- suppressMessages(suppressWarnings(lmer(Beta ~ Sound + Color + Scene + 
                                                           (1 + Sound + Color + Scene||SubID), data = roiData)))
    full.model <- suppressMessages(suppressWarnings(get_model(step(full.model, 
                                                                   reduce.random = TRUE, reduce.fixed = FALSE))))
    stats <- summary(full.model)$coefficients
    nIV <- nrow(stats) - 1 # number of fixed effects, discounting intercept
    
    ### add fixed effects and p values to data matrix
    mixed.data$ROI[(row+1):(row+nIV)]      <- thisROI
    mixed.data$Phase[(row+1):(row+nIV)]    <- time
    mixed.data$Feature[(row+1):(row+nIV)]  <- row.names(stats)[2:(nIV+1)]
    mixed.data[(row+1):(row+nIV),c("Beta","SE","T","Df","Pvalue")] <- stats[2:(nIV+1),
                                                                           c("Estimate","Std. Error","t value","df","Pr(>|t|)")]
    row = row + nIV
  }
}
#### FDR-correct p-values for multiple comparisons
mixed.data$Pfdr <- 1
p.values <- mixed.data$Pvalue
mixed.data$Pfdr <- p.adjust(p.values, method = "fdr")
## save out parameter estimates and p values as csv file for creating a .nii map of results
write.csv(mixed.data, file = curFile, row.names = FALSE)
cat('Results saved in',curFile)

Results saved in feature_effects_onsetoffset.csv

Control analysis - event epoch

Like for overall memory detail, here we are predicting activity through the entire 6s trial (onset to offset) with subsequent memory for each feature.

curFile = 'feature_effects_eventepoch.csv'
### set up data frame for mixed effect model results:
mixed.data <- data.frame(array(0,c(length(rois)*3,8)))  # *3 features
colnames(mixed.data) <- c('ROI','Phase','Feature','Beta','SE','T','Df','Pvalue')
row = 0
model.data <- gather(epoch.data, ROI, Beta, rois)
time = 'Epoch'
# ******** run linear mixed effects within each ROI and store results for each feature fixed effect **********
for (r in 1:length(rois)) {
  
  thisROI <- as.character(rois[r])
  
  roiData <- subset(model.data, ROI == thisROI)
  roiData$Sound  <- scale(roiData$Sound, scale = FALSE) # mean-center
  roiData$Color  <- scale(roiData$Color, scale = FALSE) # mean-center
  roiData$Scene  <- scale(roiData$Scene, scale = FALSE) # mean-center
  
  # run lmer, simplify random effect structure as > 2 fixed effects, and get fixed effect stats
  full.model <- suppressMessages(suppressWarnings(lmer(Beta ~ Sound + Color + Scene + 
                                                         (1 + Sound + Color + Scene||SubID), data = roiData)))
  full.model <- suppressMessages(suppressWarnings(get_model(step(full.model, 
                                                                 reduce.random = TRUE, reduce.fixed = FALSE))))
  stats <- summary(full.model)$coefficients
  nIV <- nrow(stats) - 1 # number of fixed effects, discounting intercept
  
  ### add fixed effects and p values to data matrix
  mixed.data$ROI[(row+1):(row+nIV)]      <- thisROI
  mixed.data$Phase[(row+1):(row+nIV)]    <- time
  mixed.data$Feature[(row+1):(row+nIV)]  <- row.names(stats)[2:(nIV+1)]
  mixed.data[(row+1):(row+nIV),c("Beta","SE","T","Df","Pvalue")] <- stats[2:(nIV+1),
                                                                          c("Estimate","Std. Error","t value","df","Pr(>|t|)")]
  row = row + nIV
}
#### FDR-correct all p-values for multiple comparisons
mixed.data$Pfdr <- 1
p.values <- mixed.data$Pvalue
mixed.data$Pfdr <- p.adjust(p.values, method = "fdr")
## save out parameter estimates and p values as csv file for creating a .nii map of results
write.csv(mixed.data, file = curFile, row.names = FALSE)
cat('Results saved in',curFile)

Results saved in feature_effects_eventepoch.csv

ROI visualization

Define regions and data

Note, here we are predicting memory on trial N from onset N and offset N (in the same model) to test and contrast unique variance of fixed effects

############ define rois ##############
# 1. Significant offset memory detail response (activity relates to memory for preceding event)
HIPP    <- c('LH_HIPP')
# 3. Significant onset memory detail response (activity relates to memory for upcoming event)
IFG   <- c('LH_ContA_PFCl_2','LH_DefaultB_PFCv_4','LH_ContA_PFCl_1')
# 2. Significant feature onset effects (activity relates to memory for upcoming event)
VTC   <- c('RH_DorsAttnA_TempOcc_1','RH_LimbicA_TempPole_4','RH_VisCent_ExStr_2')  # Color
AUD   <- c('RH_SomMotB_Aud_1','LH_SomMotB_Aud_1','RH_SomMotB_Aud_2','LH_SomMotB_Aud_2')  # Sound
MTL   <- c('RH_DefaultC_PHC_1','LH_DefaultC_PHC_1','RH_DefaultC_Rsp_1','LH_DefaultC_Rsp_1')  # Scene
myRois <- c(HIPP,IFG,VTC,AUD,MTL)
newRois <- c('HIPP','IFG','VTC','AUD','MTL')
mycolors <- c("#7c67a2","#ff9000","#F64F59","#1D976C","#009FFF")
cat('Regions analyzed:\n\n',
    'HIPP:',HIPP,'\n',
    'IFG:',IFG,'\n',
    'VTC:',VTC,'\n',
    'AUD:',AUD,'\n',
    'MTL:',MTL,'\n')

Regions analyzed:

 HIPP: LH_HIPP 
 IFG: LH_ContA_PFCl_2 LH_DefaultB_PFCv_4 LH_ContA_PFCl_1 
 VTC: RH_DorsAttnA_TempOcc_1 RH_LimbicA_TempPole_4 RH_VisCent_ExStr_2 
 AUD: RH_SomMotB_Aud_1 LH_SomMotB_Aud_1 RH_SomMotB_Aud_2 LH_SomMotB_Aud_2 
 MTL: RH_DefaultC_PHC_1 LH_DefaultC_PHC_1 RH_DefaultC_Rsp_1 LH_DefaultC_Rsp_1 

### make sure each trial to-be-modeled has an associated onset and offset beta
matching_trials <- intersect(onset.data$TotalEncodingTrial,offset.data$TotalEncodingTrial)
onset.data.roi  <- onset.data[onset.data$TotalEncodingTrial %in% matching_trials,]
offset.data.roi <- offset.data[offset.data$TotalEncodingTrial %in% matching_trials,]
behav.roi <- onset.data.roi[,206:244] #grab behavioral data for all trials with valid onsets & offsets
########## average betas within these new ROIs ###########
##### onset data for rois
onset.data.roi <- gather(onset.data.roi, ROI, Beta, myRois)
nleft <- length(rois) - length(myRois) + 1
onset.data.roi <- onset.data.roi[,-c(1:nleft)] #remove other regions
### allocate 'new' ROI
onset.data.roi$NewROI[onset.data.roi$ROI %in% HIPP] <- 'HIPP'
onset.data.roi$NewROI[onset.data.roi$ROI %in% IFG] <- 'IFG'
onset.data.roi$NewROI[onset.data.roi$ROI %in% VTC] <- 'VTC'
onset.data.roi$NewROI[onset.data.roi$ROI %in% AUD] <- 'AUD'
onset.data.roi$NewROI[onset.data.roi$ROI %in% MTL] <- 'MTL'
### now summarise activity within these new ROIs per trial
model.data <- onset.data.roi
grouped.data <- model.data %>% group_by(SubID,StudyTrial,NewROI) %>%
                                 summarise(NewBeta = mean(Beta))
grouped.data <- spread(grouped.data, NewROI, NewBeta)
myROIData <- merge(grouped.data,behav.roi)
model.data <- gather(myROIData, ROI, Beta, AUD:VTC)
model.data.onset <- model.data
##### offset data for rois
offset.data.roi <- gather(offset.data.roi, ROI, Beta, myRois)
nleft <- length(rois) - length(myRois) + 1
offset.data.roi <- offset.data.roi[,-c(1:nleft)]
### allocate 'new' ROI
offset.data.roi$NewROI[offset.data.roi$ROI %in% HIPP] <- 'HIPP'
offset.data.roi$NewROI[offset.data.roi$ROI %in% IFG] <- 'IFG'
offset.data.roi$NewROI[offset.data.roi$ROI %in% VTC] <- 'VTC'
offset.data.roi$NewROI[offset.data.roi$ROI %in% AUD] <- 'AUD'
offset.data.roi$NewROI[offset.data.roi$ROI %in% MTL] <- 'MTL'
### now summarise within these new ROIs
model.data <- offset.data.roi
grouped.data <- model.data %>% group_by(SubID,StudyTrial,NewROI) %>%
                                 summarise(NewBeta = mean(Beta))
grouped.data <- spread(grouped.data, NewROI, NewBeta)
myROIData <- merge(grouped.data,behav.roi)
model.data <- gather(myROIData, ROI, Beta, AUD:VTC)
model.data.offset <- model.data
###### combined onset and offset data
model.data.onset$Phase <- 'Onset'
model.data.offset$Phase <- 'Offset'
model.data <- rbind(model.data.onset,model.data.offset)

Onset/Offset subsequent memory

Now we are predicting memory detail or feature memory on trial N with onset N and offset N regressors to get unique variance and to run offset - onset contrast per ROI

####################################
mixed.data <- data.frame(array(0,c(length(newRois)*2,7)))
colnames(mixed.data) <- c('ROI','Feature','Beta','SE','T','Df','Pvalue')
row = 0
for (r in 1:length(newRois)) {
  
  this.roi <- newRois[r]
  
  roiData <- subset(model.data, ROI == this.roi)
  roiData <- spread(roiData, Phase, Beta)  #make onset and offset activity different columns
  roiData$Onset  <- scale(roiData$Onset, scale = TRUE) # z-score for direct comparison of betas
  roiData$Offset <- scale(roiData$Offset, scale = TRUE) # z-score for direct comparison of betas
  
  if (this.roi == 'HIPP' | this.roi == 'IFG') {  # Memory Detail
       roiData$MemoryDetail <- scale(roiData$MemoryDetail, scale = TRUE) # z-score for direct visual comparison across features
       # run lmer, simplify random effect structure, and get fixed effect stats
       full.model <- suppressMessages(lmer(MemoryDetail ~ Onset + Offset + 
                                                     (1 + Onset + Offset||SubID), data = roiData))
       
  } else if (this.roi == 'VTC') {  # Color
       roiData$Color <- scale(roiData$Color, scale = TRUE) # z-score for direct visual comparison across features
       # run lmer, simplify random effect structure, and get fixed effect stats
       full.model <- suppressMessages(lmer(Color ~ Onset + Offset + 
                                              (1 + Onset + Offset||SubID), data = roiData))
  } else if (this.roi == 'AUD') {  # Sound
       roiData$Sound <- scale(roiData$Sound, scale = TRUE) # z-score for direct visual comparison across features
       # run lmer, simplify random effect structure, and get fixed effect stats
       full.model <- suppressMessages(lmer(Sound ~ Onset + Offset + 
                                              (1 + Onset + Offset||SubID), data = roiData))
  } else if (this.roi == 'MTL') {  # Scene
       roiData$Scene <- scale(roiData$Scene, scale = TRUE) # z-score for direct visual comparison across features
       # run lmer, simplify random effect structure, and get fixed effect stats
       full.model <- suppressMessages(lmer(Scene ~ Onset + Offset + 
                                              (1 + Onset + Offset||SubID), data = roiData))
  }
  stats <- summary(full.model)$coefficients
  nIV <- nrow(stats) - 1 # number of fixed effects, discounting intercept
    
  # one-sample t-test (first regressor is intercept) and effect size
  con <- contest(full.model, c(0,-1,1), joint = FALSE, ddf = "Satterthwaite")
  print(kable(con, 
              format="pandoc", caption=paste(this.roi,': offset - onset',sep="")))
  ### add fixed effects and p values to data matrix
  mixed.data$ROI[(row+1):(row+nIV)]      <- as.character(this.roi)
  mixed.data$Feature[(row+1):(row+nIV)]  <- row.names(stats)[2:(nIV+1)]
  mixed.data[(row+1):(row+nIV),c("Beta","SE","T","Df","Pvalue")] <- stats[2:(nIV+1),
                                                                         c("Estimate","Std. Error","t value","df","Pr(>|t|)")]
  row = row + nIV
}

Table: HIPP: offset - onset

  Estimate   Std. Error         df    t value       lower       upper    Pr(>|t|)
----------  -----------  ---------  ---------  ----------  ----------  ----------
 0.0778628    0.0314208   52.53739   2.478067   0.0148277   0.1408979   0.0164591


Table: IFG: offset - onset

  Estimate   Std. Error         df     t value        lower        upper    Pr(>|t|)
----------  -----------  ---------  ----------  -----------  -----------  ----------
 -0.069933    0.0289432   64.63027   -2.416221   -0.1277428   -0.0121233   0.0185171


Table: VTC: offset - onset

   Estimate   Std. Error         df     t value        lower       upper    Pr(>|t|)
-----------  -----------  ---------  ----------  -----------  ----------  ----------
 -0.0342065    0.0302195   63.75277   -1.131933   -0.0945814   0.0261684   0.2619021


Table: AUD: offset - onset

   Estimate   Std. Error         df     t value        lower        upper    Pr(>|t|)
-----------  -----------  ---------  ----------  -----------  -----------  ----------
 -0.1526207    0.0316429   49.82872   -4.823219   -0.2161827   -0.0890586   0.0000137


Table: MTL: offset - onset

   Estimate   Std. Error         df     t value        lower        upper   Pr(>|t|)
-----------  -----------  ---------  ----------  -----------  -----------  ---------
 -0.1021936    0.0298283   64.05106   -3.426061   -0.1617816   -0.0426056   0.001074

# order factors for plotting
mixed.data$ROI     <- factor(mixed.data$ROI, levels=newRois)
mixed.data$Feature <- factor(mixed.data$Feature, levels=c('Onset','Offset'))
### print full output:
print(kable(mixed.data,format="pandoc",caption = 'Memory onset and offset fixed effects for ROIs'))

Table: Memory onset and offset fixed effects for ROIs

ROI    Feature          Beta          SE            T         Df      Pvalue
-----  --------  -----------  ----------  -----------  ---------  ----------
HIPP   Onset      -0.0065245   0.0244801   -0.2665232   24.62175   0.7920557
HIPP   Offset      0.0713383   0.0183772    3.8818908   25.48367   0.0006532
IFG    Onset       0.0916657   0.0208859    4.3888752   28.30980   0.0001441
IFG    Offset      0.0217326   0.0184680    1.1767737   28.28564   0.2490942
VTC    Onset       0.0696376   0.0192857    3.6108400   26.60833   0.0012466
VTC    Offset      0.0354312   0.0211580    1.6746018   25.15619   0.1064018
AUD    Onset       0.1247330   0.0254266    4.9056194   26.17756   0.0000424
AUD    Offset     -0.0278876   0.0167154   -1.6683770   26.59410   0.1069743
MTL    Onset       0.1036571   0.0189948    5.4571405   28.65378   0.0000074
MTL    Offset      0.0014635   0.0211917    0.0690607   25.07645   0.9454891

###### plot
ggplot(mixed.data, aes(x=ROI, y=Beta, fill=Feature)) +
     geom_hline(yintercept=0,size=1) +
     geom_bar(stat = "identity", alpha = 0.75, size = 1, width = 0.75,
              color = "black", position = position_dodge(1)) +
     scale_fill_manual(values = c("#BFBFBF","#404040")) +
     ggtitle("Unique effects of pre- and post-event activity\non subsequent memory") +
     geom_errorbar(aes(ymin = Beta - SE, ymax = Beta + SE), width = 0.4, size = 1,
                   color = "black", position = position_dodge(1)) +
     theme(plot.title = element_text(hjust = 0.5, size=28, margin=margin(0,0,20,0)),
          axis.line.x = element_line(color = "black", size = 1), axis.line.y = element_line(color = "black", size = 1),
          axis.text = element_text(size = 22), axis.title.y = element_text(size = 26),
          axis.title.x = element_blank(), panel.background = element_blank(),
          text = element_text(family="Helvetica"), legend.title = element_blank(),
          plot.margin = margin(1, 1, 1, 1, "cm"))

#ggsave('Rois.png',plot=last_plot(),width=11,height=5,units="in",dpi=500) 

Modulation of onset/offset relationship by memory detail

Here, we are testing if the relationship between onset activity of IFG/VTC/AUD/PHC-RSC and HIPP offset activity is modulated by subsequent memory detail (strongest synchrony if more features recalled?)

####################################
mixed.data <- data.frame(array(0,c(4*3,7)))
colnames(mixed.data) <- c('ROI','Predictor','Beta','SE','T','Df','Pvalue')
row = 0
new.df <- model.data[model.data$ROI == 'HIPP' & model.data$Phase == 'Offset',c("SubID","Beta","MemoryDetail")]
colnames(new.df)[2] <- 'HIPP'
new.df$IFG <- model.data[model.data$ROI == 'IFG' & model.data$Phase == 'Onset',c("Beta")]
new.df$VTC <- model.data[model.data$ROI == 'VTC' & model.data$Phase == 'Onset',c("Beta")]
new.df$AUD <- model.data[model.data$ROI == 'AUD' & model.data$Phase == 'Onset',c("Beta")]
new.df$MTL <- model.data[model.data$ROI == 'MTL' & model.data$Phase == 'Onset',c("Beta")]
new.df[,2:7] <- scale(new.df[,2:7], scale = TRUE) # z score all variables
for (r in 1:(length(newRois)-1)) {
  
  if (r == 1) {  # IFG
       full.model <- suppressMessages(lmer(HIPP ~ IFG*MemoryDetail + 
                                             (1 + IFG*MemoryDetail||SubID), data = new.df))
  } else if (r == 2) {  # VTC
       full.model <- suppressMessages(lmer(HIPP ~ VTC*MemoryDetail + 
                                             (1 + VTC*MemoryDetail||SubID), data = new.df))
  } else if (r == 3) {  # AUD
       full.model <- suppressMessages(lmer(HIPP ~ AUD*MemoryDetail + 
                                             (1 + AUD*MemoryDetail||SubID), data = new.df))
  } else if (r == 4) {  # MTL
       full.model <- suppressMessages(lmer(HIPP ~ MTL*MemoryDetail + 
                                             (1 + MTL*MemoryDetail||SubID), data = new.df))
  }
  stats <- summary(full.model)$coefficients
  nIV <- nrow(stats) - 1 #removing intercept
  ### add fixed effects and p values to data matrix
  mixed.data$ROI[(row+1):(row+nIV)]        <- as.character(newRois[r+1])
  mixed.data$Predictor[(row+1):(row+nIV)]  <- row.names(stats)[2:(nIV+1)]
  mixed.data[(row+1):(row+nIV),c("Beta","SE","T","Df","Pvalue")] <- stats[2:(nIV+1),
                                                                          c("Estimate","Std. Error","t value","df","Pr(>|t|)")]
  row = row + nIV
}
### print full output:
print(kable(mixed.data,format="pandoc",caption = 'Hipp Offset ~ [ROI] Onset * Memory Detail'))

Table: Hipp Offset ~ [ROI] Onset * Memory Detail

ROI   Predictor                 Beta          SE            T           Df      Pvalue
----  -----------------  -----------  ----------  -----------  -----------  ----------
IFG   IFG                  0.1496640   0.0223252    6.7038143     22.37516   0.0000009
IFG   MemoryDetail         0.0679213   0.0231550    2.9333270     22.82346   0.0075071
IFG   IFG:MemoryDetail    -0.0080134   0.0178003   -0.4501845     24.11499   0.6565985
VTC   VTC                  0.0959211   0.0241919    3.9650043     21.03917   0.0007043
VTC   MemoryDetail         0.0728899   0.0235342    3.0971874     23.11032   0.0050628
VTC   VTC:MemoryDetail    -0.0062203   0.0178057   -0.3493454     15.99040   0.7313882
AUD   AUD                  0.1133499   0.0226320    5.0083941     23.36767   0.0000437
AUD   MemoryDetail         0.0724787   0.0233296    3.1067301     23.42577   0.0048983
AUD   AUD:MemoryDetail     0.0110593   0.0168390    0.6567632   1422.24848   0.5114394
MTL   MTL                  0.1210409   0.0258302    4.6860189     24.47365   0.0000882
MTL   MemoryDetail         0.0700037   0.0219351    3.1914049     21.66850   0.0042720
MTL   MTL:MemoryDetail     0.0045381   0.0191412    0.2370848     19.51902   0.8150615

Hippocampal feature combinations

Which specific feature combinations is hippocampus sensitive to? Here, we are predicting L_HIPP onset or offset activity with binary regressors capturing each possible combination of features recalled.

my.plots = list()
plot.max = c(0.6, 1)
c = 0
for (this.roi in c('HIPP')) {
####################################
mixed.data <- data.frame(array(0,c(7*2,7)))
colnames(mixed.data) <- c('Feature','Beta','SE','T','Df','Pvalue','Phase')
row=0 
### loop over onset and offset
for (phase in 1:2) {
  
  # gather data so behavioral variables are replicated over each ROI (now as a single column factor)
  if (phase == 1) {
    roiData <- subset(model.data.onset, ROI == this.roi)
    time = 'Onset'
  } else if (phase == 2) {
    roiData <- subset(model.data.offset, ROI == this.roi)
    time = 'Offset'
  }
  
  # NOTE. not mean-centering here as zero is meaningful -- excluding trials in other conditions. Implicit baseline = none recalled
  full.model <- suppressMessages(suppressWarnings(lmer(Beta ~ Color.Only + Sound.Only + Scene.Only +
                                                              Color.Scene + Sound.Color + Scene.Sound +
                                                              Color.Scene.Sound +
                                                         (1 + Color.Only + Sound.Only + Scene.Only +
                                                              Color.Scene + Sound.Color + Scene.Sound +
                                                              Color.Scene.Sound||SubID), 
                                                       data = roiData)))
  full.model <- suppressMessages(suppressWarnings(get_model(step(full.model, 
                                                                 reduce.random = TRUE, reduce.fixed = FALSE))))
  stats <- summary(full.model)$coefficients
  nIV <- nrow(stats) - 1 # number of fixed effects, discounting intercept
  
  # test if scene binding is > recalling scene alone (first position is intercept)
  con <- contest(full.model, c(0,0,0,-1,(1/3),0,(1/3),(1/3)), joint = FALSE, ddf = "Satterthwaite")
  print(kable(con,
              format="pandoc", caption=paste(this.roi,': Scene Binding vs. Scene Alone at ',time,sep="")))
  # test if scene binding is > recalling color and sound without scene (first position is intercept)
  con <- contest(full.model, c(0,0,0,0,(1/3),-1,(1/3),(1/3)), joint = FALSE, ddf = "Satterthwaite")
  print(kable(con,
              format="pandoc", caption=paste(this.roi,': Scene Binding vs. Color&Sound at ',time,sep="")))
  
  ### add fixed effects and p values to data matrix    
  mixed.data$Phase[(row+1):(row+nIV)]    <- time
  mixed.data$Feature[(row+1):(row+nIV)]  <- row.names(stats)[2:(nIV+1)]
  mixed.data[(row+1):(row+nIV),c("Beta","SE","T","Df","Pvalue")] <- stats[2:(nIV+1),
                                                                          c("Estimate","Std. Error","t value","df","Pr(>|t|)")]
  row = row + nIV
}
##### PLOT
order <- c("Color.Only","Scene.Only","Sound.Only","Sound.Color","Color.Scene","Scene.Sound","Color.Scene.Sound")
# colors for venn diagram space feature combinations
barcolors <- c("#f76c67","#67aaf7","#0bca81","#f6a84f","#9c67f7","#67f7f0","#7F7F7F")
mixed.data$Feature <- factor(mixed.data$Feature, levels=order)
mixed.data$Phase <- factor(mixed.data$Phase, levels=c('Onset','Offset'))
c = c+1
my.plots[[c]] <- ggplot(mixed.data, aes(x=Feature, y=Beta, fill=Feature)) +
  facet_grid(. ~ Phase) +
  geom_hline(yintercept=0,size=1) +
  geom_bar(stat = "identity", alpha = 0.80, size = 1, width = 0.75,
           color = "black", position = position_dodge(1)) +
  scale_fill_manual(values = barcolors) +
  scale_y_continuous(limits = c(-0.3,plot.max[c])) +
  geom_errorbar(aes(ymin = Beta - SE, ymax = Beta + SE), width = 0.4, size = 1,
                color = "black", position = position_dodge(1)) +
  ggtitle(paste(this.roi,'encoding activity vs. none recalled',sep=" ")) +
  theme(plot.title = element_text(hjust = 0.5, size=28, margin=margin(0,0,20,0)),
        axis.line.x = element_line(color = "black", size = 1), axis.line.y = element_line(color = "black", size = 1),
        axis.text.x = element_text(size = 16, angle = 45, hjust=1, vjust=1),
        axis.text.y = element_text(size = 22), axis.title.y = element_text(size = 26),
        axis.title.x = element_blank(), panel.background = element_blank(),
        legend.position="none", strip.text.x = element_text(size=28), text = element_text(family="Helvetica"),
        plot.margin = margin(1, 1, 1, 1, "cm"))
### print full output:
print(kable(mixed.data,format="pandoc",caption = paste('Remembered features (vs. none recalled) for',this.roi,sep = " ")))
}

Table: HIPP: Scene Binding vs. Scene Alone at Onset

  Estimate   Std. Error         df     t value        lower       upper    Pr(>|t|)
----------  -----------  ---------  ----------  -----------  ----------  ----------
 0.0490484    0.1119956   2594.081   0.4379491   -0.1705615   0.2686583   0.6614597


Table: HIPP: Scene Binding vs. Color&Sound at Onset

  Estimate   Std. Error         df     t value        lower       upper    Pr(>|t|)
----------  -----------  ---------  ----------  -----------  ----------  ----------
 0.0309742    0.1129433   2469.901   0.2742455   -0.1904991   0.2524474   0.7839189


Table: HIPP: Scene Binding vs. Scene Alone at Offset

 Estimate   Std. Error         df    t value        lower       upper    Pr(>|t|)
---------  -----------  ---------  ---------  -----------  ----------  ----------
 0.239022    0.1476316   22.43739   1.619043   -0.0668017   0.5448456   0.1194101


Table: HIPP: Scene Binding vs. Color&Sound at Offset

  Estimate   Std. Error         df    t value       lower       upper    Pr(>|t|)
----------  -----------  ---------  ---------  ----------  ----------  ----------
 0.2328095    0.1167075   1285.221   1.994812   0.0038514   0.4617677   0.0462743


Table: Remembered features (vs. none recalled) for HIPP

Feature                    Beta          SE            T           Df      Pvalue  Phase  
------------------  -----------  ----------  -----------  -----------  ----------  -------
Color.Only           -0.1415691   0.1165886   -1.2142619   3555.61457   0.2247284  Onset  
Sound.Only           -0.0489743   0.1539740   -0.3180684   3554.33722   0.7504517  Onset  
Scene.Only           -0.0811641   0.1291471   -0.6284619   3559.88947   0.5297417  Onset  
Color.Scene           0.0131775   0.1103558    0.1194093   3548.08313   0.9049579  Onset  
Sound.Color          -0.0630898   0.1304795   -0.4835231   3574.89334   0.6287540  Onset  
Scene.Sound          -0.0796042   0.1186943   -0.6706652   3540.50522   0.5024776  Onset  
Color.Scene.Sound    -0.0299203   0.1098568   -0.2723576     91.60632   0.7859600  Onset  
Color.Only            0.0875331   0.1178921    0.7424850   3537.85669   0.4578428  Offset 
Sound.Only            0.0388375   0.1556466    0.2495238   3529.99510   0.8029701  Offset 
Scene.Only            0.1604685   0.1598229    1.0040394     30.36143   0.3232950  Offset 
Color.Scene           0.3601762   0.1322209    2.7240494     51.89518   0.0087661  Offset 
Sound.Color           0.1666809   0.1319913    1.2628168   3562.50568   0.2067377  Offset 
Scene.Sound           0.4503095   0.1201740    3.7471461   3525.03630   0.0001817  Offset 
Color.Scene.Sound     0.3879857   0.1128184    3.4390288     90.31313   0.0008853  Offset 

ggarrange(plotlist = my.plots, ncol = c, nrow = 1)

Anterior vs. posterior

##### load in hippocampal voxel data for all subjects:
subjects <- list.files(path='../data/single-trial-betas/event-offset/',full.names=TRUE, recursive = TRUE, pattern = 'offset-LH_HIPP')
allData  = do.call(rbind, lapply(subjects, function(x) {read.csv(x, header = TRUE)}))
allData$SubID=as.factor(allData$SubID)
allData$Segment <- as.factor(allData$Segment)
### add z score of betas within region and subject to remove outliers:
allData <- allData %>%
               group_by(SubID, Segment) %>%
               mutate(Z = zscore(Beta))
# mark betas that are > XSD
allData$Beta[abs(allData$Z) > outlier_value] <- NA  #mark beta as NA if outlier
allData <- allData[,-c(6)] # remove Z values  
# merge with behavioral data
myBetas = spread(allData, Segment, Beta)
myBetas$TotalEncodingTrial <- 1:nTotal
segment.info <- merge(myBetas, behavData, by="TotalEncodingTrial")
segment.info <- segment.info[,-c(2,3,4)]
colnames(segment.info)[colnames(segment.info) == 'SubID.y'] <- 'SubID'
colnames(segment.info)[colnames(segment.info) == 'Run.y'] <- 'Run'
# remove all trials where there is at least one NA (at least one of the ROIs is an outlier on that trial):
idxBetas     <- complete.cases(myBetas)
segment.info <- segment.info[idxBetas,]
segment.info <- gather(segment.info, Segment, Beta, Anterior:Posterior)
  
####################################
mixed.data <- data.frame(array(0,c(7*2,7)))
colnames(mixed.data) <- c('Segment','Feature','Beta','SE','T','Df','Pvalue')
row=0 
for (v in c('Anterior','Posterior')) {
  
  this.data <- subset(segment.info, Segment == v)
  
  # NOTE. not mean-centering here as zero is meaningful -- excluding trials in other conditions. Implicit baseline = none recalled
  full.model <- suppressMessages(suppressWarnings(lmer(Beta ~ Color.Only + Sound.Only + Scene.Only +
                                                              Color.Scene + Sound.Color + Scene.Sound +
                                                              Color.Scene.Sound +
                                                         (1 + Color.Only + Sound.Only + Scene.Only +
                                                              Color.Scene + Sound.Color + Scene.Sound +
                                                              Color.Scene.Sound||SubID),
                                                       data = this.data)))
  full.model <- suppressMessages(suppressWarnings(get_model(step(full.model, 
                                                                 reduce.random = TRUE, reduce.fixed = FALSE))))
  stats <- summary(full.model)$coefficients
  nIV <- nrow(stats) - 1 # number of fixed effects, discounting intercept
  
  # binding: >=2 features vs. 1
  con <- contest(full.model, c(0,-(1/3),-(1/3),-(1/3),(1/4),(1/4),(1/4),(1/4)), joint = FALSE, ddf = "Satterthwaite")
  print(kable(con,
              format="pandoc", caption=paste(this.roi,': Binding (>1 feature vs. 1 feature): ',v,sep="")))
  # test if scene binding is > recalling color and sound without scene (first position is intercept)
  con <- contest(full.model, c(0,0,0,0,(1/3),-1,(1/3),(1/3)), joint = FALSE, ddf = "Satterthwaite")
  print(kable(con,
              format="pandoc", caption=paste(this.roi,': Scene Binding vs. Color&Sound: ',v,sep="")))
  
  ### add fixed effects and p values to data matrix    
  mixed.data$Segment[(row+1):(row+nIV)]  <- v
  mixed.data$Feature[(row+1):(row+nIV)]  <- row.names(stats)[2:(nIV+1)]
  mixed.data[(row+1):(row+nIV),c("Beta","SE","T","Df","Pvalue")] <- stats[2:(nIV+1),
                                                                            c("Estimate","Std. Error","t value","df","Pr(>|t|)")]
  row = row + nIV
}

Table: HIPP: Binding (>1 feature vs. 1 feature): Anterior

  Estimate   Std. Error         df    t value       lower     upper   Pr(>|t|)
----------  -----------  ---------  ---------  ----------  --------  ---------
 0.2518564    0.0977626   2264.288   2.576204   0.0601427   0.44357   0.010052


Table: HIPP: Scene Binding vs. Color&Sound: Anterior

  Estimate   Std. Error         df    t value       lower       upper    Pr(>|t|)
----------  -----------  ---------  ---------  ----------  ----------  ----------
 0.3674827    0.1431639   2476.641   2.566867   0.0867494   0.6482161   0.0103205


Table: HIPP: Binding (>1 feature vs. 1 feature): Posterior

  Estimate   Std. Error         df    t value       lower       upper    Pr(>|t|)
----------  -----------  ---------  ---------  ----------  ----------  ----------
 0.2296148     0.095161   798.2388   2.412909   0.0428194   0.4164102   0.0160502


Table: HIPP: Scene Binding vs. Color&Sound: Posterior

  Estimate   Std. Error         df     t value        lower       upper   Pr(>|t|)
----------  -----------  ---------  ----------  -----------  ----------  ---------
 0.0399494    0.1383117   1092.224   0.2888357   -0.2314374   0.3113361   0.772762

### print full output:
print(kable(mixed.data,format="pandoc",caption = paste('Remembered feature combinations (vs. none recalled) for Anterior vs. Posterior',this.roi,sep = " ")))

Table: Remembered feature combinations (vs. none recalled) for Anterior vs. Posterior HIPP

Segment     Feature                    Beta          SE            T           Df      Pvalue
----------  ------------------  -----------  ----------  -----------  -----------  ----------
Anterior    Color.Only            0.0252935   0.1481130    0.1707715   3751.46478   0.8644127
Anterior    Sound.Only           -0.0729386   0.1950732   -0.3739039   3750.31295   0.7084969
Anterior    Scene.Only            0.1793859   0.1645030    1.0904717   3754.80165   0.2755754
Anterior    Color.Scene           0.4046226   0.1392027    2.9067152   3754.48405   0.0036738
Anterior    Sound.Color           0.0201579   0.1650469    0.1221343   3772.33454   0.9027991
Anterior    Scene.Sound           0.3882230   0.1502669    2.5835574   3743.48907   0.0098163
Anterior    Color.Scene.Sound     0.3700763   0.1420303    2.6056152     87.75649   0.0107723
Posterior   Color.Only           -0.0072212   0.1408325   -0.0512752   3747.72789   0.9591090
Posterior   Sound.Only            0.0901479   0.1854740    0.4860405   3741.94425   0.6269669
Posterior   Scene.Only            0.1133564   0.1563476    0.7250281   3745.43220   0.4684801
Posterior   Color.Scene           0.1661307   0.1700307    0.9770630     39.05305   0.3345519
Posterior   Sound.Color           0.2650805   0.1564853    1.6939645   3761.39199   0.0903548
Posterior   Scene.Sound           0.4083049   0.1419232    2.8769414   3762.88419   0.0040381
Posterior   Color.Scene.Sound     0.3406540   0.1160245    2.9360522   3646.37646   0.0033450

##### PLOT
order <- c("Color.Only","Scene.Only","Sound.Only","Sound.Color","Color.Scene","Scene.Sound","Color.Scene.Sound")
# colors for venn diagram space feature combinations
barcolors <- c("#f76c67","#67aaf7","#0bca81","#f6a84f","#9c67f7","#67f7f0","#7F7F7F")
mixed.data$Feature <- factor(mixed.data$Feature, levels=order)
mixed.data$Segment <- as.factor(mixed.data$Segment)
##### PLOT
ggplot(mixed.data, aes(x=Feature, y=Beta, fill=Feature)) +
  facet_grid(. ~ Segment) +
  geom_hline(yintercept=0,size=1) +
  scale_fill_manual(values = barcolors) +
  scale_y_continuous(limits = c(-0.35,0.65)) +
  geom_bar(stat = "identity", alpha = 0.80, size = 1, width = 0.75,
           color = "black", position = position_dodge(1)) +
  geom_errorbar(aes(ymin = Beta - SE, ymax = Beta + SE), width = 0.4, size = 1,
                color = "black", position = position_dodge(1)) +
  ggtitle("Anterior/posterior hippocampal offset encoding activity") +
  theme(plot.title = element_text(hjust = 0.5, size=28, margin=margin(0,0,20,0)),
        axis.line.x = element_line(color = "black", size = 1), axis.line.y = element_line(color = "black", size = 1),
        axis.text.x = element_text(size = 16, angle = 45, hjust=1, vjust=1),
        axis.text.y = element_text(size = 22), axis.title.y = element_text(size = 26),
        axis.title.x = element_blank(), panel.background = element_blank(),
        legend.position="none", strip.text.x = element_text(size=28), text = element_text(family="Helvetica"),
        plot.margin = margin(1, 1, 1, 1, "cm"))