SST: ROI analyses
- Is hippocampal activity during planning reduced under stress?
- On a given probe trial, does hippocampal activity during the planning period predict P(shortcut)?
- Across subjects, does mean hippocampal activity during assigned period predict proportion of shortcuts?
- Is CCN activity during planning reduced for stress?
- On a given probe trial, does cognitive control activity predict P(shortcut)?
- Across subjects, does mean cognitive control activity during assigned period predict proportion of shortcuts?
- Does hippocampal activity when arriving at the goal on rep 1 predict p(shortcut) on rep 2?
# Define some functions to get data:
get_behav_data = function(run_type='shortcut', path_chosen='shortcut'){
# run_type (str): either shortcut, or habit
# path_chosen (str): one of the codes from route coders
# Load in onset info (generated from nb SST_behav_fmri_generateData)
######################################################################
filepath = '/Volumes/group/awagner/sgagnon/SST/nav_data/analysis/SST_event_onsets.csv'
d_onsets = read.csv(filepath)
d_onsets = d_onsets %>%
filter(type == get("run_type"),
command != 'SHOCK') %>%
mutate(env = as.numeric(gsub("env", "", env))); dim(d_onsets)[1]/4
# Load in route coding info
######################################################################
filepath = '/Volumes/group/awagner/sgagnon/SST/nav_data/analysis/fmri_route_coding.csv'
d_code = read.csv(filepath)
d_code = d_code %>%
filter(run_type == get("run_type")) %>% # just look at probe trials
filter(!subid %in% c('sst39', 'sst01', 'sst02')) %>% # bad coding
mutate(run_type = factor(run_type),
shortcut = ifelse(final_code == path_chosen, 1, 0),
subid = as.numeric(gsub("sst", "", subid))) # convert subid to numeric
# dim(d_code)
# str(d_code)
# length(unique(d_code$subid))
# str(d_onsets)
# str(d_code)
# Merge w/onsets
d_behav = inner_join(d_onsets, d_code) %>%
mutate(run = scan) %>% # need the scan number for matching w/the neural data
select(subid, env, rep, final_code, shortcut, scan, run,
cond, onset, trial, command)
# dim(d_behav)
# length(unique(d_behav$subid))
# Add in group info
filepath = '/Volumes/group/awagner/sgagnon/SST/scripts/subjects_groups.csv'
d_group = read.csv(filepath) %>%
mutate(subid = as.numeric(gsub("sst", "", subid)))
d_behav = d_behav %>% left_join(d_group)
# str(d_behav)
# dim(d_behav %>% filter(is.na(final_code)))
# length(unique(d_behav$subid))
# str(d_behav)
# dlook = d_behav %>% group_by(subid, rep, command) %>% summarise(n=n())
# dim(dlook %>% filter(n != 12))
return(d_behav)
}
get_pe_brain_data = function(roi_name, exp, condition){
filepath = paste0('/Volumes/group/awagner/sgagnon/SST/analysis/', exp,
'/group/roi/pe_', roi_name, '.csv')
d = read.csv(filepath)
d = d %>% filter(cond == condition)
return(d)
}
merge_behav_brain_data = function(d_behav, roi_name, event_type='ASSIGNED'){
# Load in brain data
######################################################################
basepath = '/Volumes/group/awagner/sgagnon/SST/analysis/mvpa-runtype/group/roi/extractraw_integrated_mvpa-runtype_'
d_behav = d_behav %>% filter(command == get("event_type"))
d = data.frame()
for (hemi in c('lh', 'rh')){
roi = paste0(hemi,'-', roi_name)
roipath = paste0(basepath, roi, '.csv')
droi = read.csv(roipath) %>%
mutate(hemi = hemi)
d = d %>%
bind_rows(droi) %>%
mutate(mask = factor(mask)) %>%
mutate(hemi = factor(hemi))
}
# Average over hemis
d = d %>%
select(-hemi, -mask) %>%
group_by(subid, run, onset, condition) %>%
summarise_each(funs(mean_)) %>%
ungroup()
# length(unique(d$subid))
# Merge dataframes
dm = d %>%
mutate(subid = as.numeric(gsub("sst", "", subid)),
cond = condition) %>%
right_join(d_behav) %>%
filter(rep < 3) %>% # for one subj w/an accidental repeat
mutate(subid = factor(subid),
rep = factor(rep)) %>%
rowwise() %>%
mutate(mean_activity = mean_(c(X4, X6, X8, X10, X12))) %>%
ungroup() # turn off rowwise
# str(dm)
# print(length(unique(dm$subid)))
dlook = dm %>% group_by(subid, rep) %>%
summarise(n=n()) %>% ungroup() %>%
complete(subid, rep, fill = list(n = 0))
print('Missing trials:')
print(dim(dlook %>% filter(n != 12)))
if (dim(dlook %>% filter(n != 12))[1] == 0){
print('Have all trials!')
}else{
print(dlook %>% filter(n != 12))
print(dlook %>% filter(n != 12) %>% select(subid) %>% unlist())
}
dnull = dm %>% filter(is.na(integrated_activity)) # identify trials w/non complete TRs (ie scanner ended)
dnull = dm %>% filter(is.na(mean_activity)) # identify trials w/no samples
dm = dm %>% filter(!is.na(mean_activity))
# str(dm)
dm = dm %>% filter(command == event_type)
return(dm)
}
zscore <- function(vec){
(vec-mean(vec))/sd(vec)
}Is hippocampal activity during planning reduced under stress?
Calculate from extracted PEs:
dm = get_pe_brain_data(roi_name = 'hippocampus', exp = 'mvpa-runtype',
condition = 'ASSIGNED_shortcut')
dm = dm %>% filter(hemi %in% c('rh', 'lh')) %>%
mutate(hemi = factor(hemi),
cond = factor(cond))
str(dm)## 'data.frame': 72 obs. of 11 variables:
## $ X : int 72 73 180 181 288 289 396 397 504 505 ...
## $ cond : Factor w/ 1 level "ASSIGNED_shortcut": 1 1 1 1 1 1 1 1 1 1 ...
## $ hemi : Factor w/ 2 levels "lh","rh": 1 2 1 2 1 2 1 2 1 2 ...
## $ mask_vox : num 172 194 182 207 170 176 204 213 202 214 ...
## $ regspace : Factor w/ 1 level "epi": 1 1 1 1 1 1 1 1 1 1 ...
## $ roi : Factor w/ 1 level "hippocampus": 1 1 1 1 1 1 1 1 1 1 ...
## $ smoothing: Factor w/ 1 level "unsmoothed": 1 1 1 1 1 1 1 1 1 1 ...
## $ subid : Factor w/ 36 levels "sst03","sst04",..: 1 1 2 2 3 3 4 4 5 5 ...
## $ value : num 53.85 -2.28 -61.96 -13.67 -97.4 ...
## $ group : Factor w/ 2 levels "control","stress": 1 1 1 1 1 1 1 1 1 1 ...
## $ remove : logi NA NA NA NA NA NA ...contrasts(dm$group) = c(1,-1); contrasts(dm$group)## [,1]
## control 1
## stress -1contrasts(dm$hemi) = c(1,-1); contrasts(dm$hemi)## [,1]
## lh 1
## rh -1summary(lmer(value ~ group * hemi + (1|subid), data=dm))## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: value ~ group * hemi + (1 | subid)
## Data: dm
##
## REML criterion at convergence: 965.1
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9274 -0.4616 0.0003 0.4165 2.0578
##
## Random effects:
## Groups Name Variance Std.Dev.
## subid (Intercept) 167763 409.6
## Residual 12667 112.5
## Number of obs: 72, groups: subid, 36
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 67.135 69.649 34.000 0.964 0.342
## group1 32.245 69.649 34.000 0.463 0.646
## hemi1 -6.035 13.285 34.010 -0.454 0.653
## group1:hemi1 -14.996 13.285 34.010 -1.129 0.267
##
## Correlation of Fixed Effects:
## (Intr) group1 hemi1
## group1 -0.056
## hemi1 0.000 0.000
## group1:hem1 0.000 0.000 -0.056Calculate from extracted “raw” activity:
d_behav = get_behav_data()
dm = merge_behav_brain_data(d_behav, roi_name = 'hippocampus', event_type = 'ASSIGNED')## [1] "Missing trials:"
## [1] 1 3
## # A tibble: 1 × 3
## subid rep n
## <fctr> <fctr> <dbl>
## 1 20 2 9
## subid
## 20
## 40 Levels: 3 4 5 6 7 9 10 11 12 13 14 15 16 18 19 20 21 22 23 24 26 ... 57str(dm)## Classes 'tbl_df', 'tbl' and 'data.frame': 957 obs. of 23 variables:
## $ subid : Factor w/ 40 levels "3","4","5","6",..: 7 7 7 7 7 7 7 7 7 7 ...
## $ run : num 5 5 5 6 6 6 7 7 7 8 ...
## $ onset : num 11.3 66.3 123.4 11.3 74.3 ...
## $ condition : Factor w/ 9 levels "ARRIVED_habit",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ X0 : num -0.2308 -0.3501 0.1544 -0.0108 -0.1125 ...
## $ X2 : num -0.03899 -0.16166 0.16421 0.34324 -0.00806 ...
## $ X4 : num 0.0587 0.0235 0.5368 0.7856 0.0205 ...
## $ X6 : num 0.221 0.209 0.485 0.601 0.219 ...
## $ X8 : num 0.2493 0.0759 0.4389 0.1872 0.2014 ...
## $ X10 : num -0.244 -0.169 0.301 0.357 0.333 ...
## $ X12 : num -0.1129 -0.347 0.1667 0.3178 -0.0315 ...
## $ integrated_activity: num 0.252 -0.112 1.853 2.09 0.543 ...
## $ cond : chr "ASSIGNED_shortcut" "ASSIGNED_shortcut" "ASSIGNED_shortcut" "ASSIGNED_shortcut" ...
## $ env : num 4 5 6 1 2 3 7 8 9 10 ...
## $ rep : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ final_code : Factor w/ 4 levels "backtrack","habit",..: 4 4 1 1 1 4 2 2 4 4 ...
## $ shortcut : num 1 1 0 0 0 1 0 0 1 1 ...
## $ scan : num 5 5 5 6 6 6 7 7 7 8 ...
## $ trial : num 1 2 3 1 2 3 1 2 3 1 ...
## $ command : Factor w/ 5 levels "ARRIVED","ASSIGNED",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ group : Factor w/ 2 levels "control","stress": 1 1 1 1 1 1 1 1 1 1 ...
## $ remove : int NA NA NA NA NA NA NA NA NA NA ...
## $ mean_activity : num 0.0344 -0.0415 0.3858 0.4496 0.1484 ...contrasts(dm$group) = c(1,-1); contrasts(dm$group)## [,1]
## control 1
## stress -1contrasts(dm$rep) = c(1,-1); contrasts(dm$rep)## [,1]
## 1 1
## 2 -1summary(lmer(mean_activity ~ group * rep + (1+rep|subid), data=dm))## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: mean_activity ~ group * rep + (1 + rep | subid)
## Data: dm
##
## REML criterion at convergence: -286.3
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.4618 -0.4819 0.0936 0.6041 4.4734
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subid (Intercept) 4.743e-03 0.068866
## rep1 6.681e-05 0.008174 0.94
## Residual 3.990e-02 0.199754
## Number of obs: 957, groups: subid, 40
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.037744 0.012660 38.020000 -2.981 0.00498 **
## group1 0.007001 0.012660 38.020000 0.553 0.58350
## rep1 -0.011281 0.006586 38.160000 -1.713 0.09486 .
## group1:rep1 0.006929 0.006586 38.160000 1.052 0.29941
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) group1 rep1
## group1 -0.001
## rep1 0.157 0.002
## group1:rep1 0.002 0.157 -0.003dfwc = dm %>% group_by(subid, group, rep) %>%
summarise(mean_activity = mean(mean_activity)) %>%
group_by(group, rep) %>%
summarise(se = std.error(mean_activity),
mean_activity = mean(mean_activity), n())
dfwc## Source: local data frame [4 x 5]
## Groups: group [?]
##
## group rep se mean_activity `n()`
## <fctr> <fctr> <dbl> <dbl> <int>
## 1 control 1 0.02030501 -0.03509578 20
## 2 control 2 0.02042922 -0.02639111 20
## 3 stress 1 0.02251425 -0.06295473 20
## 4 stress 2 0.01697973 -0.02649052 20p2 = ggplot(dfwc, aes(x=rep, y=mean_activity, group=group, color=group)) +
geom_line(position=pos_dodge, size=1.5) +
geom_errorbar(width=0, aes(ymin=mean_activity-se,
ymax=mean_activity+se), size=1.5, position=pos_dodge) +
geom_point(size=5, position=pos_dodge) +
ylab('Hippocampal activity')+
xlab('Repetition')+
scale_color_manual(values=palette, guide = guide_legend(title = NULL))
p2
On a given probe trial, does hippocampal activity during the planning period predict P(shortcut)?
d_behav = get_behav_data()
dm = merge_behav_brain_data(d_behav, roi_name = 'hippocampus', event_type = 'ASSIGNED')## [1] "Missing trials:"
## [1] 1 3
## # A tibble: 1 × 3
## subid rep n
## <fctr> <fctr> <dbl>
## 1 20 2 9
## subid
## 20
## 40 Levels: 3 4 5 6 7 9 10 11 12 13 14 15 16 18 19 20 21 22 23 24 26 ... 57# dm %>% filter(complete.cases(.))
# dm %>% group_by(subid) %>%
# summarise(mean = mean(mean_activity)) %>% collect %>% .[["mean"]] %>% hist()
# z-score activity within subject
dm = dm %>% group_by(subid) %>%
mutate(mean_activity_z = zscore(mean_activity))
# dm %>% group_by(subid) %>%
# summarise(mean = mean(mean_activity_z)) %>% collect %>% .[["mean"]] %>% hist()
contrasts(dm$group) = c(1,-1); contrasts(dm$group)## [,1]
## control 1
## stress -1contrasts(dm$rep) = c(1,-1); contrasts(dm$rep)## [,1]
## 1 1
## 2 -1res = glmer(shortcut ~ mean_activity_z*rep*group +
(1+mean_activity_z*rep | subid),
control=glmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=2e5)),
family=binomial, data=dm)
summary(res)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## shortcut ~ mean_activity_z * rep * group + (1 + mean_activity_z *
## rep | subid)
## Data: dm
## Control:
## glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik deviance df.resid
## 1266.7 1354.2 -615.3 1230.7 939
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2718 -0.8368 -0.3939 0.8794 2.9448
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subid (Intercept) 0.57734 0.7598
## mean_activity_z 0.04768 0.2184 -0.34
## rep1 0.08690 0.2948 0.80 -0.39
## mean_activity_z:rep1 0.05703 0.2388 -0.87 0.23 -0.97
## Number of obs: 957, groups: subid, 40
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.23313 0.14048 -1.660 0.09701 .
## mean_activity_z 0.21716 0.08593 2.527 0.01150 *
## rep1 -0.24051 0.08648 -2.781 0.00541 **
## group1 0.15651 0.14025 1.116 0.26445
## mean_activity_z:rep1 0.02396 0.08853 0.271 0.78663
## mean_activity_z:group1 -0.13030 0.08516 -1.530 0.12599
## rep1:group1 0.11501 0.08616 1.335 0.18195
## mean_activity_z:rep1:group1 -0.06341 0.08749 -0.725 0.46859
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mn_ct_ rep1 group1 mn_ctvty_z:r1 mn_ctvty_z:g1
## men_ctvty_z -0.163
## rep1 0.412 -0.105
## group1 -0.014 0.025 -0.017
## mn_ctvty_z:r1 -0.336 0.173 -0.294 0.002
## mn_ctvty_z:g1 0.025 -0.064 -0.001 -0.163 -0.007
## rep1:group1 -0.018 0.007 -0.038 0.414 0.044 -0.109
## mn_ctv_:1:1 -0.004 -0.001 0.040 -0.337 -0.063 0.153
## rp1:g1
## men_ctvty_z
## rep1
## group1
## mn_ctvty_z:r1
## mn_ctvty_z:g1
## rep1:group1
## mn_ctv_:1:1 -0.295What if we control for habit performance?
d_behav_habit = get_behav_data(run_type = 'habit', path_chosen = 'habit')## Joining, by = c("env", "rep", "subid")## Joining, by = "subid"d_behav_habit = d_behav_habit %>%
mutate(subid = factor(subid)) %>%
group_by(subid) %>%
summarise(prop_habit = mean(shortcut))
str(d_behav_habit)## Classes 'tbl_df', 'tbl' and 'data.frame': 40 obs. of 2 variables:
## $ subid : Factor w/ 40 levels "3","4","5","6",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ prop_habit: num 0.833 1 0.833 0.833 1 ...dm = dm %>% left_join(d_behav_habit, by='subid')
res = glmer(shortcut ~ mean_activity_z*rep*group + scale(prop_habit)+
(1+mean_activity_z*rep | subid),
control=glmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=2e5)),
family=binomial, data=dm)
summary(res)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: shortcut ~ mean_activity_z * rep * group + scale(prop_habit) +
## (1 + mean_activity_z * rep | subid)
## Data: dm
## Control:
## glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik deviance df.resid
## 1255.1 1347.5 -608.6 1217.1 938
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0327 -0.8252 -0.4031 0.8584 2.8772
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subid (Intercept) 0.32351 0.5688
## mean_activity_z 0.04197 0.2049 -0.15
## rep1 0.09241 0.3040 0.80 -0.25
## mean_activity_z:rep1 0.04295 0.2072 -0.78 0.09 -0.99
## Number of obs: 957, groups: subid, 40
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.23116 0.11554 -2.001 0.045429 *
## mean_activity_z 0.20939 0.08473 2.471 0.013461 *
## rep1 -0.23399 0.08712 -2.686 0.007234 **
## group1 0.20847 0.11628 1.793 0.072986 .
## scale(prop_habit) 0.44449 0.11442 3.885 0.000102 ***
## mean_activity_z:rep1 0.01459 0.08602 0.170 0.865331
## mean_activity_z:group1 -0.12511 0.08430 -1.484 0.137766
## rep1:group1 0.11217 0.08701 1.289 0.197349
## mean_activity_z:rep1:group1 -0.08021 0.08515 -0.942 0.346255
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mn_ct_ rep1 group1 scl(_) mn_ctvty_z:r1
## men_ctvty_z -0.094
## rep1 0.395 -0.061
## group1 -0.019 0.032 -0.019
## scl(prp_hb) -0.016 -0.001 0.022 0.123
## mn_ctvty_z:r1 -0.242 0.144 -0.269 -0.001 -0.049
## mn_ctvty_z:g1 0.032 -0.068 0.000 -0.094 -0.006 -0.021
## rep1:group1 -0.022 0.007 -0.032 0.392 -0.001 0.048
## mn_ctv_:1:1 0.000 -0.016 0.041 -0.242 -0.033 -0.070
## mn_ctvty_z:g1 rp1:g1
## men_ctvty_z
## rep1
## group1
## scl(prp_hb)
## mn_ctvty_z:r1
## mn_ctvty_z:g1
## rep1:group1 -0.060
## mn_ctv_:1:1 0.121 -0.271# plot
d_f = dm %>%
group_by(subid) %>%
mutate(med = ntile(mean_activity_z, 2))
# look at mean activity across tertiles
# glimpse(d_f)
data = d_f %>%
group_by(subid, med) %>%
summarise(mean = mean_(mean_activity_z), n())
# glimpse(data)
ggplot(data, aes(x=med, y=mean, group=subid, color=subid)) +
geom_point(alpha=.5) +
geom_line(alpha=.5) +
theme(legend.position="none")
dlook = d_f %>%
group_by(subid, rep, med) %>%
summarise(n=n()) %>%
ungroup() %>%
complete(subid, rep, med, fill = list(n = 0))
dim(dlook)[1]## [1] 160print('Low trial counts per bin:')## [1] "Low trial counts per bin:"dlook$n## [1] 5 7 7 5 5 7 7 5 6 6 6 6 5 7 7 5 5 7 7 5 8 4 4 8 3 9 9 3 6 6 6 6 6 6 6
## [36] 6 6 6 6 6 9 3 3 9 6 6 6 6 9 3 3 9 5 7 7 5 7 5 5 7 6 6 5 4 8 4 4 8 5 7
## [71] 7 5 6 6 6 6 3 9 9 3 6 6 6 6 8 4 4 8 7 5 5 7 7 5 5 7 6 6 6 6 7 5 5 7 8
## [106] 4 4 8 7 5 5 7 5 7 7 5 5 7 7 5 5 7 7 5 7 5 5 7 4 8 8 4 8 4 4 8 7 5 5 7
## [141] 8 4 4 8 7 5 5 7 7 5 5 7 7 5 5 7 6 6 6 6d_avg = d_f %>%
group_by(subid, group, med, rep) %>%
summarise(shortcut = mean_(shortcut))
dfwc = summarySEwithin(d_avg, measurevar="shortcut", withinvars=c("med", "rep"),
betweenvars = c('group'),
idvar="subid", na.rm=TRUE, conf.interval=.95)
dfwc## group med rep N shortcut shortcut_norm sd se
## 1 control 1 1 20 0.4804762 0.4436334 0.2135020 0.04774050
## 2 control 1 2 20 0.4967857 0.4599430 0.1890770 0.04227890
## 3 control 2 1 20 0.4697222 0.4328795 0.2141335 0.04788170
## 4 control 2 2 20 0.5220833 0.4852406 0.2354155 0.05264050
## 5 stress 1 1 20 0.3245635 0.3614062 0.2372516 0.05305107
## 6 stress 1 2 20 0.3680357 0.4048785 0.2947582 0.06590993
## 7 stress 2 1 20 0.4082738 0.4451166 0.1954066 0.04369423
## 8 stress 2 2 20 0.5734524 0.6102951 0.1918526 0.04289955
## ci
## 1 0.09992202
## 2 0.08849076
## 3 0.10021755
## 4 0.11017784
## 5 0.11103716
## 6 0.13795106
## 7 0.09145308
## 8 0.08978979p2 = ggplot(dfwc, aes(x=med, y=shortcut, group=group, color=group)) +
geom_line(position=pos_dodge, size=1.5) +
geom_errorbar(width=0, aes(ymin=shortcut-se,
ymax=shortcut+se), size=1.5, position=pos_dodge) +
geom_point(size=5, position=pos_dodge) +
facet_grid(.~rep)+
xlab('Hippocampal activity (median split)')+
ylab('p Shortcut')+
scale_color_manual(values=palette, guide = guide_legend(title = NULL))
p2
Collapse over rep (hipp by tertile)
# plot
d_f = dm %>%
group_by(subid) %>%
mutate(tertile = ntile(mean_activity_z, 3))
# look at mean activity across tertiles
# glimpse(d_f)
data = d_f %>%
group_by(subid, tertile) %>%
summarise(mean = mean_(mean_activity_z), n())
d_avg = d_f %>%
group_by(subid, group, tertile) %>%
summarise(shortcut = mean_(shortcut))
dfwc = summarySEwithin(d_avg, measurevar="shortcut", withinvars=c("tertile"),
betweenvars = c('group'),
idvar="subid", na.rm=TRUE, conf.interval=.95)## Automatically converting the following non-factors to factors: tertiledfwc## group tertile N shortcut shortcut_norm sd se
## 1 control 1 20 0.4562500 0.4261905 0.1335867 0.02987088
## 2 control 2 20 0.5125000 0.4824405 0.1417263 0.03169097
## 3 control 3 20 0.4875000 0.4574405 0.2064443 0.04616235
## 4 stress 1 20 0.3589286 0.3889881 0.1864232 0.04168549
## 5 stress 2 20 0.3892857 0.4193452 0.1764992 0.03946642
## 6 stress 3 20 0.5276786 0.5577381 0.1882496 0.04209390
## ci
## 1 0.06252048
## 2 0.06632996
## 3 0.09661890
## 4 0.08724873
## 5 0.08260416
## 6 0.08810354p2 = ggplot(dfwc, aes(x=tertile, y=shortcut, group=group, color=group)) +
geom_line(position=pos_dodge, size=1.5) +
geom_errorbar(width=0, aes(ymin=shortcut-se,
ymax=shortcut+se), size=1.5, position=pos_dodge) +
geom_point(size=5, position=pos_dodge) +
xlab('Hippocampal activity (tertiles)')+
ylab('p Shortcut')+
scale_color_manual(values=palette, guide = guide_legend(title = NULL))
p2
What about just for the first iteration through the environments?
# what about just for rep 1?
res = glmer(shortcut ~ mean_activity_z*group +
(1+mean_activity_z | subid),
control=glmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=2e5)),
family=binomial, data=dm %>% filter(rep == 1))
summary(res)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## shortcut ~ mean_activity_z * group + (1 + mean_activity_z | subid)
## Data: dm %>% filter(rep == 1)
## Control:
## glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik deviance df.resid
## 609.5 638.8 -297.8 595.5 473
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.9521 -0.7522 -0.4030 0.8121 3.1614
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subid (Intercept) 0.9902 0.9951
## mean_activity_z 0.1227 0.3503 -1.00
## Number of obs: 480, groups: subid, 40
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4796 0.1908 -2.514 0.0119 *
## mean_activity_z 0.2436 0.1318 1.848 0.0646 .
## group1 0.2794 0.1903 1.468 0.1420
## mean_activity_z:group1 -0.2132 0.1290 -1.652 0.0984 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mn_ct_ group1
## men_ctvty_z -0.412
## group1 -0.035 0.030
## mn_ctvty_:1 0.024 -0.072 -0.415Across subjects, does mean hippocampal activity during assigned period predict proportion of shortcuts?
mean_hipp = dm %>%
group_by(subid, group, rep) %>%
summarise(mean_activity=mean_(mean_activity)) %>%
ungroup()
# str(mean_hipp)
prop_shortcut = d_behav %>%
group_by(subid, group, rep) %>%
summarise(shortcut = mean(shortcut), n()) %>%
ungroup() %>%
mutate(subid = factor(subid),
rep = factor(rep))
# means
prop_shortcut %>%
group_by(group, rep) %>%
summarise(mean(shortcut), std.error(shortcut), n())## Source: local data frame [4 x 5]
## Groups: group [?]
##
## group rep `mean(shortcut)` `std.error(shortcut)` `n()`
## <fctr> <fctr> <dbl> <dbl> <int>
## 1 control 1 0.4583333 0.05322220 20
## 2 control 2 0.5125000 0.05031879 20
## 3 stress 1 0.3541667 0.05667279 20
## 4 stress 2 0.4972222 0.03377509 20# glimpse(prop_shortcut)
# str(prop_shortcut)
d_gp = left_join(mean_hipp, prop_shortcut)
# str(d_gp)
# averaged over rep
d_avg = d_gp %>%
select(-rep, -group) %>%
group_by(subid) %>%
summarise_each(funs(mean)) %>%
left_join(d_gp %>% select(subid, group) %>% distinct())
# d_avg
hist(d_avg$shortcut)
summary(res <- lm(shortcut ~ scale(mean_activity), data=d_avg))##
## Call:
## lm(formula = shortcut ~ scale(mean_activity), data = d_avg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.41916 -0.13211 -0.04819 0.15972 0.41897
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.45556 0.03070 14.839 <2e-16 ***
## scale(mean_activity) 0.01045 0.03109 0.336 0.739
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1942 on 38 degrees of freedom
## Multiple R-squared: 0.002963, Adjusted R-squared: -0.02328
## F-statistic: 0.1129 on 1 and 38 DF, p-value: 0.7387hist(res$residuals)
with(d_avg %>% mutate(shortcut_log = logit(shortcut)), plot(shortcut, shortcut_log))
summary(lm(shortcut ~ scale(mean_activity), data=d_avg %>% mutate(shortcut = logit(shortcut))))##
## Call:
## lm(formula = shortcut ~ scale(mean_activity), data = d_avg %>%
## mutate(shortcut = logit(shortcut)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.9386 -0.5198 -0.1542 0.6884 2.1635
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.21965 0.14993 -1.465 0.151
## scale(mean_activity) 0.04515 0.15184 0.297 0.768
##
## Residual standard error: 0.9482 on 38 degrees of freedom
## Multiple R-squared: 0.002321, Adjusted R-squared: -0.02393
## F-statistic: 0.08841 on 1 and 38 DF, p-value: 0.7678library("lmtest")
bptest(res)##
## studentized Breusch-Pagan test
##
## data: res
## BP = 0.41525, df = 1, p-value = 0.5193ggplot(d_avg, aes(mean_activity, shortcut, color=group)) +
geom_smooth(method='lm', colour="darkgray") +
geom_point() +
scale_color_manual(values=c('dodgerblue', 'orange')) +
ylab('Proportion shortcut') +
xlab('Hippocampal activity')
Across subjects, there’s no relationship between mean hippocampal activity during planning and the proportion of shortcuts taken. Note that we used OLS here, even though DV is a proportion, but errors look normally distributed, and not heteroskedastic.
What about just for rep 1?
# just rep 1?
d_avg = d_gp %>%
filter(rep == 1) %>%
select(-group, -rep) %>%
group_by(subid) %>%
summarise_each(funs(mean)) %>%
left_join(d_gp %>% select(subid, group) %>% distinct())
# d_avg
summary(lm(shortcut ~ scale(mean_activity), data=d_avg))##
## Call:
## lm(formula = shortcut ~ scale(mean_activity), data = d_avg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.34980 -0.21951 -0.05618 0.17305 0.45546
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.40625 0.03945 10.298 1.5e-12 ***
## scale(mean_activity) 0.03199 0.03995 0.801 0.428
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2495 on 38 degrees of freedom
## Multiple R-squared: 0.01659, Adjusted R-squared: -0.00929
## F-statistic: 0.641 on 1 and 38 DF, p-value: 0.4283# logit transform proportion
summary(lm(shortcut ~ scale(mean_activity), data=d_avg %>% mutate(shortcut = logit(shortcut))))##
## Call:
## lm(formula = shortcut ~ scale(mean_activity), data = d_avg %>%
## mutate(shortcut = logit(shortcut)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.0464 -0.9702 -0.1478 0.8280 2.2995
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5162 0.2002 -2.579 0.0139 *
## scale(mean_activity) 0.1960 0.2027 0.967 0.3398
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.266 on 38 degrees of freedom
## Multiple R-squared: 0.02401, Adjusted R-squared: -0.001677
## F-statistic: 0.9347 on 1 and 38 DF, p-value: 0.3398Still no relationship…
Is CCN activity during planning reduced for stress?
Calculate from extracted PEs:
dm = get_pe_brain_data(roi_name = 'frontoparietal', exp = 'mvpa-runtype',
condition = 'ASSIGNED_shortcut')
str(dm)## 'data.frame': 72 obs. of 11 variables:
## $ X : int 80 81 188 189 296 297 404 405 512 513 ...
## $ cond : Factor w/ 9 levels "ARRIVED_habit",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ hemi : Factor w/ 2 levels "lh","rh": 1 2 1 2 1 2 1 2 1 2 ...
## $ mask_vox : num 1403 1678 1887 2263 1606 ...
## $ regspace : Factor w/ 1 level "epi": 1 1 1 1 1 1 1 1 1 1 ...
## $ roi : Factor w/ 1 level "frontoparietal": 1 1 1 1 1 1 1 1 1 1 ...
## $ smoothing: Factor w/ 1 level "unsmoothed": 1 1 1 1 1 1 1 1 1 1 ...
## $ subid : Factor w/ 36 levels "sst03","sst04",..: 1 1 2 2 3 3 4 4 5 5 ...
## $ value : num 1231 931 470 1028 581 ...
## $ group : Factor w/ 2 levels "control","stress": 1 1 1 1 1 1 1 1 1 1 ...
## $ remove : logi NA NA NA NA NA NA ...contrasts(dm$group) = c(1,-1); contrasts(dm$group)## [,1]
## control 1
## stress -1contrasts(dm$hemi) = c(1,-1); contrasts(dm$hemi)## [,1]
## lh 1
## rh -1summary(lmer(value ~ group * hemi + (1|subid), data=dm))## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: value ~ group * hemi + (1 | subid)
## Data: dm
##
## REML criterion at convergence: 1083.8
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.72831 -0.53833 0.06846 0.49613 1.88057
##
## Random effects:
## Groups Name Variance Std.Dev.
## subid (Intercept) 404370 635.9
## Residual 150966 388.5
## Number of obs: 72, groups: subid, 36
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 1105.347 115.631 34.000 9.559 3.65e-11 ***
## group1 221.165 115.631 34.000 1.913 0.0642 .
## hemi1 4.945 45.861 33.790 0.108 0.9148
## group1:hemi1 -76.649 45.861 33.790 -1.671 0.1039
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) group1 hemi1
## group1 -0.056
## hemi1 0.000 0.000
## group1:hem1 0.000 0.000 -0.056Calculate from extracted “raw” activity:
d_behav = get_behav_data()
dm = merge_behav_brain_data(d_behav, roi_name = 'frontoparietal', event_type = 'ASSIGNED')## [1] "Missing trials:"
## [1] 1 3
## # A tibble: 1 × 3
## subid rep n
## <fctr> <fctr> <dbl>
## 1 20 2 9
## subid
## 20
## 40 Levels: 3 4 5 6 7 9 10 11 12 13 14 15 16 18 19 20 21 22 23 24 26 ... 57str(dm)## Classes 'tbl_df', 'tbl' and 'data.frame': 957 obs. of 23 variables:
## $ subid : Factor w/ 40 levels "3","4","5","6",..: 7 7 7 7 7 7 7 7 7 7 ...
## $ run : num 5 5 5 6 6 6 7 7 7 8 ...
## $ onset : num 11.3 66.3 123.4 11.3 74.3 ...
## $ condition : Factor w/ 9 levels "ARRIVED_habit",..: 4 4 4 4 4 4 4 4 4 4 ...
## $ X0 : num 0.2384 -0.3217 -0.0352 -0.1949 -0.0825 ...
## $ X2 : num 0.3704 -0.22 0.1595 -0.0274 -0.0024 ...
## $ X4 : num 0.695 0.147 0.207 0.421 0.179 ...
## $ X6 : num 0.597 0.258 0.31 0.61 0.219 ...
## $ X8 : num 0.41 0.0632 0.4866 0.2376 0.316 ...
## $ X10 : num 0.1615 -0.2543 0.7298 0.0105 0.3404 ...
## $ X12 : num 0.1687 -0.3773 0.6727 -0.0331 0.3578 ...
## $ integrated_activity: num 2.2721 -0.0399 1.5103 1.1493 0.8405 ...
## $ cond : chr "ASSIGNED_shortcut" "ASSIGNED_shortcut" "ASSIGNED_shortcut" "ASSIGNED_shortcut" ...
## $ env : num 4 5 6 1 2 3 7 8 9 10 ...
## $ rep : Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
## $ final_code : Factor w/ 4 levels "backtrack","habit",..: 4 4 1 1 1 4 2 2 4 4 ...
## $ shortcut : num 1 1 0 0 0 1 0 0 1 1 ...
## $ scan : num 5 5 5 6 6 6 7 7 7 8 ...
## $ trial : num 1 2 3 1 2 3 1 2 3 1 ...
## $ command : Factor w/ 5 levels "ARRIVED","ASSIGNED",..: 2 2 2 2 2 2 2 2 2 2 ...
## $ group : Factor w/ 2 levels "control","stress": 1 1 1 1 1 1 1 1 1 1 ...
## $ remove : int NA NA NA NA NA NA NA NA NA NA ...
## $ mean_activity : num 0.4064 -0.0327 0.4812 0.2493 0.2824 ...contrasts(dm$group) = c(1,-1); contrasts(dm$group)## [,1]
## control 1
## stress -1contrasts(dm$rep) = c(1,-1); contrasts(dm$rep)## [,1]
## 1 1
## 2 -1summary(lmer(mean_activity ~ group * rep + (1+rep|subid), data=dm))## Linear mixed model fit by REML t-tests use Satterthwaite approximations
## to degrees of freedom [lmerMod]
## Formula: mean_activity ~ group * rep + (1 + rep | subid)
## Data: dm
##
## REML criterion at convergence: 320.7
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.6222 -0.5529 0.0860 0.6221 3.0325
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subid (Intercept) 4.862e-03 0.069728
## rep1 5.952e-05 0.007715 1.00
## Residual 7.674e-02 0.277028
## Number of obs: 957, groups: subid, 40
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) -0.012782 0.014204 38.100000 -0.900 0.3738
## group1 0.037053 0.014204 38.100000 2.609 0.0129 *
## rep1 -0.011055 0.009039 424.100000 -1.223 0.2220
## group1:rep1 0.008102 0.009039 424.100000 0.896 0.3706
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) group1 rep1
## group1 -0.001
## rep1 0.103 0.002
## group1:rep1 0.002 0.103 -0.003dfwc = dm %>% group_by(subid, group, rep) %>%
summarise(mean_activity = mean(mean_activity)) %>%
group_by(group, rep) %>%
summarise(se = std.error(mean_activity),
mean_activity = mean(mean_activity), n())
dfwc## Source: local data frame [4 x 5]
## Groups: group [?]
##
## group rep se mean_activity `n()`
## <fctr> <fctr> <dbl> <dbl> <int>
## 1 control 1 0.02382853 0.02131719 20
## 2 control 2 0.02314123 0.02722456 20
## 3 stress 1 0.02464823 -0.06899270 20
## 4 stress 2 0.01961900 -0.03084760 20p2 = ggplot(dfwc, aes(x=rep, y=mean_activity, group=group, color=group)) +
geom_line(position=pos_dodge, size=1.5) +
geom_errorbar(width=0, aes(ymin=mean_activity-se,
ymax=mean_activity+se), size=1.5, position=pos_dodge) +
geom_point(size=5, position=pos_dodge) +
ylab('CCN activity')+
xlab('Repetition')+
scale_color_manual(values=palette, guide = guide_legend(title = NULL))
p2
On a given probe trial, does cognitive control activity predict P(shortcut)?
# z-score activity within subject
dm = dm %>% group_by(subid) %>%
mutate(mean_activity_z = zscore(mean_activity))
contrasts(dm$group) = c(1,-1); contrasts(dm$group)## [,1]
## control 1
## stress -1contrasts(dm$rep) = c(1,-1); contrasts(dm$rep)## [,1]
## 1 1
## 2 -1res = glmer(shortcut ~ mean_activity_z*rep*group +
(1+mean_activity_z*rep | subid),
control=glmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=2e5)),
family=binomial, data=dm)
summary(res)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## shortcut ~ mean_activity_z * rep * group + (1 + mean_activity_z *
## rep | subid)
## Data: dm
## Control:
## glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik deviance df.resid
## 1266.4 1353.9 -615.2 1230.4 939
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.1572 -0.8390 -0.3957 0.8622 3.0425
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subid (Intercept) 0.61313 0.7830
## mean_activity_z 0.06047 0.2459 -0.68
## rep1 0.08220 0.2867 0.87 -0.84
## mean_activity_z:rep1 0.03022 0.1739 -0.75 0.15 -0.64
## Number of obs: 957, groups: subid, 40
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.232169 0.143627 -1.616 0.10599
## mean_activity_z 0.118484 0.086043 1.377 0.16850
## rep1 -0.245317 0.085499 -2.869 0.00411 **
## group1 0.160942 0.143658 1.120 0.26258
## mean_activity_z:rep1 -0.001644 0.081976 -0.020 0.98400
## mean_activity_z:group1 -0.210459 0.086752 -2.426 0.01527 *
## rep1:group1 0.125741 0.085722 1.467 0.14242
## mean_activity_z:rep1:group1 -0.136397 0.082675 -1.650 0.09898 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mn_ct_ rep1 group1 mn_ctvty_z:r1 mn_ctvty_z:g1
## men_ctvty_z -0.307
## rep1 0.449 -0.239
## group1 -0.014 0.020 -0.020
## mn_ctvty_z:r1 -0.240 0.184 -0.182 0.003
## mn_ctvty_z:g1 0.023 -0.038 0.010 -0.311 -0.012
## rep1:group1 -0.021 0.009 -0.038 0.449 0.033 -0.247
## mn_ctv_:1:1 0.005 -0.009 0.036 -0.244 -0.032 0.187
## rp1:g1
## men_ctvty_z
## rep1
## group1
## mn_ctvty_z:r1
## mn_ctvty_z:g1
## rep1:group1
## mn_ctv_:1:1 -0.191What if we control for habit performance?
d_behav_habit = get_behav_data(run_type = 'habit', path_chosen = 'habit')## Joining, by = c("env", "rep", "subid")## Joining, by = "subid"d_behav_habit = d_behav_habit %>%
mutate(subid = factor(subid)) %>%
group_by(subid) %>%
summarise(prop_habit = mean(shortcut))
str(d_behav_habit)## Classes 'tbl_df', 'tbl' and 'data.frame': 40 obs. of 2 variables:
## $ subid : Factor w/ 40 levels "3","4","5","6",..: 1 2 3 4 5 6 7 8 9 10 ...
## $ prop_habit: num 0.833 1 0.833 0.833 1 ...dm = dm %>% left_join(d_behav_habit, by='subid')
res = glmer(shortcut ~ mean_activity_z*rep*group + scale(prop_habit)+
(1+mean_activity_z*rep | subid),
control=glmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=2e5)),
family=binomial, data=dm)
summary(res)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: shortcut ~ mean_activity_z * rep * group + scale(prop_habit) +
## (1 + mean_activity_z * rep | subid)
## Data: dm
## Control:
## glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik deviance df.resid
## 1252.4 1344.9 -607.2 1214.4 938
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.0095 -0.8278 -0.4045 0.8569 3.2936
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subid (Intercept) 0.34881 0.5906
## mean_activity_z 0.05189 0.2278 -0.57
## rep1 0.08364 0.2892 0.89 -0.88
## mean_activity_z:rep1 0.03486 0.1867 -0.91 0.19 -0.63
## Number of obs: 957, groups: subid, 40
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.229747 0.118269 -1.943 0.05207 .
## mean_activity_z 0.111794 0.084343 1.325 0.18502
## rep1 -0.242271 0.085594 -2.830 0.00465 **
## group1 0.214253 0.118714 1.805 0.07111 .
## scale(prop_habit) 0.469350 0.110673 4.241 2.23e-05 ***
## mean_activity_z:rep1 -0.001813 0.082236 -0.022 0.98241
## mean_activity_z:group1 -0.213452 0.085173 -2.506 0.01221 *
## rep1:group1 0.124145 0.085800 1.447 0.14792
## mean_activity_z:rep1:group1 -0.155570 0.083340 -1.867 0.06195 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mn_ct_ rep1 group1 scl(_) mn_ctvty_z:r1
## men_ctvty_z -0.235
## rep1 0.434 -0.240
## group1 -0.018 0.024 -0.026
## scl(prp_hb) -0.013 -0.012 0.006 0.097
## mn_ctvty_z:r1 -0.288 0.185 -0.179 0.001 -0.017
## mn_ctvty_z:g1 0.030 -0.034 0.011 -0.239 -0.035 -0.019
## rep1:group1 -0.028 0.006 -0.031 0.430 0.003 0.036
## mn_ctv_:1:1 0.009 -0.018 0.037 -0.294 -0.029 -0.028
## mn_ctvty_z:g1 rp1:g1
## men_ctvty_z
## rep1
## group1
## scl(prp_hb)
## mn_ctvty_z:r1
## mn_ctvty_z:g1
## rep1:group1 -0.250
## mn_ctv_:1:1 0.198 -0.184# plot
d_f = dm %>%
group_by(subid) %>%
mutate(med = ntile(mean_activity_z, 2))
# look at mean activity across medians
# glimpse(d_f)
data = d_f %>%
group_by(subid, med) %>%
summarise(mean = mean_(mean_activity_z), n())
# glimpse(data)
ggplot(data, aes(x=med, y=mean, group=subid, color=subid)) +
geom_point(alpha=.5) +
geom_line(alpha=.5) +
theme(legend.position="none")
dlook = d_f %>%
group_by(subid, rep, med) %>%
summarise(n=n()) %>%
ungroup() %>%
complete(subid, rep, med, fill = list(n = 0))
dim(dlook)[1]## [1] 160print('Low trial counts per bin:')## [1] "Low trial counts per bin:"dlook$n## [1] 5 7 7 5 5 7 7 5 6 6 6 6 7 5 5 7 4 8 8 4 8 4 4 8 4 8 8 4 5 7 7 5 5 7 7
## [36] 5 5 7 7 5 7 5 5 7 7 5 5 7 7 5 5 7 5 7 7 5 6 6 6 6 6 6 5 4 7 5 5 7 6 6
## [71] 6 6 9 3 3 9 4 8 8 4 6 6 6 6 7 5 5 7 6 6 6 6 8 4 4 8 7 5 5 7 7 5 5 7 7
## [106] 5 5 7 8 4 4 8 7 5 5 7 6 6 6 6 6 6 6 6 5 7 7 5 5 7 7 5 7 5 5 7 7 5 5 7
## [141] 5 7 7 5 7 5 5 7 6 6 6 6 5 7 7 5 7 5 5 7d_avg = d_f %>%
group_by(subid, group, med, rep) %>%
summarise(shortcut = mean_(shortcut))
dfwc = summarySEwithin(d_avg, measurevar="shortcut", withinvars=c("med", "rep"),
betweenvars = c('group'),
idvar="subid", na.rm=TRUE, conf.interval=.95)
dfwc## group med rep N shortcut shortcut_norm sd se
## 1 control 1 1 20 0.4984127 0.4646553 0.2270181 0.05076280
## 2 control 1 2 20 0.5082738 0.4745164 0.1679167 0.03754730
## 3 control 2 1 20 0.4367262 0.4029687 0.1618391 0.03618831
## 4 control 2 2 20 0.5240278 0.4902703 0.2306900 0.05158386
## 5 stress 1 1 20 0.3045238 0.3382812 0.2178555 0.04871397
## 6 stress 1 2 20 0.4711905 0.5049479 0.2545811 0.05692607
## 7 stress 2 1 20 0.3958333 0.4295908 0.2017525 0.04511324
## 8 stress 2 2 20 0.5258333 0.5595908 0.2090096 0.04673598
## ci
## 1 0.10624776
## 2 0.07858741
## 3 0.07574301
## 4 0.10796625
## 5 0.10195951
## 6 0.11914762
## 7 0.09442309
## 8 0.09781952p2 = ggplot(dfwc, aes(x=med, y=shortcut, group=group, color=group)) +
geom_line(position=pos_dodge, size=1.5) +
geom_errorbar(width=0, aes(ymin=shortcut-se,
ymax=shortcut+se), size=1.5, position=pos_dodge) +
geom_point(size=5, position=pos_dodge) +
facet_grid(.~rep)+
xlab('CCN activity (median split)')+
ylab('p Shortcut')+
scale_color_manual(values=palette, guide = guide_legend(title = NULL))
p2
Collapse over rep (hipp by tertile)
# plot
d_f = dm %>%
group_by(subid) %>%
mutate(tertile = ntile(mean_activity_z, 3))
# look at mean activity across tertiles
# glimpse(d_f)
data = d_f %>%
group_by(subid, tertile) %>%
summarise(mean = mean_(mean_activity_z), n())
d_avg = d_f %>%
group_by(subid, group, tertile) %>%
summarise(shortcut = mean_(shortcut))
dfwc = summarySEwithin(d_avg, measurevar="shortcut", withinvars=c("tertile"),
betweenvars = c('group'),
idvar="subid", na.rm=TRUE, conf.interval=.95)## Automatically converting the following non-factors to factors: tertiledfwc## group tertile N shortcut shortcut_norm sd se
## 1 control 1 20 0.4875000 0.4574405 0.1763692 0.03943735
## 2 control 2 20 0.5312500 0.5011905 0.1774731 0.03968420
## 3 control 3 20 0.4375000 0.4074405 0.1889745 0.04225598
## 4 stress 1 20 0.3642857 0.3943452 0.1468519 0.03283709
## 5 stress 2 20 0.4196429 0.4497024 0.1741539 0.03894200
## 6 stress 3 20 0.4919643 0.5220238 0.1947783 0.04355376
## ci
## 1 0.08254331
## 2 0.08305998
## 3 0.08844277
## 4 0.06872882
## 5 0.08150655
## 6 0.09115906p2 = ggplot(dfwc, aes(x=tertile, y=shortcut, group=group, color=group)) +
geom_line(position=pos_dodge, size=1.5) +
geom_errorbar(width=0, aes(ymin=shortcut-se,
ymax=shortcut+se), size=1.5, position=pos_dodge) +
geom_point(size=5, position=pos_dodge) +
xlab('CCN activity (tertiles)')+
ylab('p Shortcut')+
scale_color_manual(values=palette, guide = guide_legend(title = NULL))
p2
What about just for rep 1?
# what about just for rep 1?
res = glmer(shortcut ~ mean_activity_z*group +
(1+mean_activity_z | subid),
control=glmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=2e5)),
family=binomial(link=logit), data=dm %>% filter(rep == 1))
summary(res)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## shortcut ~ mean_activity_z * group + (1 + mean_activity_z | subid)
## Data: dm %>% filter(rep == 1)
## Control:
## glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik deviance df.resid
## 606.2 635.4 -296.1 592.2 473
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.7421 -0.7360 -0.4196 0.8472 3.3787
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subid (Intercept) 1.0482 1.0238
## mean_activity_z 0.1043 0.3229 -1.00
## Number of obs: 480, groups: subid, 40
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4819 0.1949 -2.472 0.01342 *
## mean_activity_z 0.1229 0.1291 0.952 0.34127
## group1 0.2943 0.1950 1.509 0.13125
## mean_activity_z:group1 -0.3583 0.1300 -2.756 0.00585 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mn_ct_ group1
## men_ctvty_z -0.390
## group1 -0.035 0.027
## mn_ctvty_:1 0.029 -0.058 -0.397However, just looking at rep 1, we do see a marginal interaction by group (p < 0.1), and a marginal effect of group (p < 0.1).
Across subjects, does mean cognitive control activity during assigned period predict proportion of shortcuts?
mean_cc = dm %>%
group_by(subid, group, rep) %>%
summarise(mean_activity=mean_(mean_activity)) %>%
ungroup()
prop_shortcut = d_behav %>%
group_by(subid, group, rep) %>%
summarise(shortcut = mean(shortcut)) %>%
ungroup() %>%
mutate(subid = factor(subid),
rep = factor(rep))
# glimpse(prop_shortcut)
prop_shortcut %>%
group_by(group, rep) %>%
summarise(mean(shortcut), n())## Source: local data frame [4 x 4]
## Groups: group [?]
##
## group rep `mean(shortcut)` `n()`
## <fctr> <fctr> <dbl> <int>
## 1 control 1 0.4583333 20
## 2 control 2 0.5125000 20
## 3 stress 1 0.3541667 20
## 4 stress 2 0.4972222 20d_gp = left_join(mean_cc, prop_shortcut)
# str(d_gp)
# averaged over rep
d_avg = d_gp %>%
select(-rep, -group) %>%
group_by(subid) %>%
summarise_each(funs(mean)) %>%
left_join(d_gp %>% select(subid, group) %>% distinct())
# d_avg
hist(d_avg$shortcut)
res = lm(shortcut ~ scale(mean_activity), data=d_avg)
hist(res$residuals)
summary(res)##
## Call:
## lm(formula = shortcut ~ scale(mean_activity), data = d_avg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35099 -0.13119 -0.01809 0.14326 0.30781
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.45556 0.02740 16.624 <2e-16 ***
## scale(mean_activity) 0.08703 0.02775 3.136 0.0033 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1733 on 38 degrees of freedom
## Multiple R-squared: 0.2056, Adjusted R-squared: 0.1847
## F-statistic: 9.833 on 1 and 38 DF, p-value: 0.003301# plot(res)
# transform proportions to logit
with(d_avg %>% mutate(shortcut_log = logit(shortcut)), plot(shortcut, shortcut_log))
summary(lm(shortcut ~ scale(mean_activity), data=d_avg %>% mutate(shortcut = logit(shortcut))))##
## Call:
## lm(formula = shortcut ~ scale(mean_activity), data = d_avg %>%
## mutate(shortcut = logit(shortcut)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.60013 -0.55584 -0.00725 0.60729 1.44403
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2196 0.1328 -1.654 0.10637
## scale(mean_activity) 0.4368 0.1345 3.248 0.00243 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.8399 on 38 degrees of freedom
## Multiple R-squared: 0.2173, Adjusted R-squared: 0.1967
## F-statistic: 10.55 on 1 and 38 DF, p-value: 0.002433# does this model exhibit heteroskedasticity?
# studentized Breusch and Pagan (1979) test of Koenker (1981)
bptest(res)##
## studentized Breusch-Pagan test
##
## data: res
## BP = 1.066, df = 1, p-value = 0.3019# but still try tobit regression (since DV is bounded)
# Tobit regression coefficients are interpreted in the similar manner to OLS regression coefficients; however, the linear effect is on the uncensored latent variable, not the observed outcome (http://stats.idre.ucla.edu/r/dae/tobit-models/)
library(VGAM)
summary(m <- vglm(shortcut ~ scale(mean_activity),
tobit(Upper = 1,
Lower = 0),
data=d_avg))##
## Call:
## vglm(formula = shortcut ~ scale(mean_activity), family = tobit(Upper = 1,
## Lower = 0), data = d_avg)
##
##
## Pearson residuals:
## Min 1Q Median 3Q Max
## mu -2.080 -0.7770 -0.1076 0.8484 1.822
## loge(sd) -0.733 -0.4905 -0.2025 0.2886 2.371
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept):1 0.45556 0.02673 17.046 < 2e-16 ***
## (Intercept):2 -1.77830 0.11276 -15.771 < 2e-16 ***
## scale(mean_activity) 0.08703 0.02710 3.211 0.00132 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Number of linear predictors: 2
##
## Names of linear predictors: mu, loge(sd)
##
## Log-likelihood: 14.3743 on 77 degrees of freedom
##
## Number of iterations: 5# intercept1 = interecept/constant for the model
# intercept2 = exp(int2) = sq(resid variance from OLS)
# what about beta regression?
library(betareg)
summary(m <- betareg(shortcut ~ scale(mean_activity), data=d_avg))##
## Call:
## betareg(formula = shortcut ~ scale(mean_activity), data = d_avg)
##
## Standardized weighted residuals 2:
## Min 1Q Median 3Q Max
## -3.2639 -0.6897 -0.0218 0.7837 1.8976
##
## Coefficients (mean model with logit link):
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1897 0.1105 -1.718 0.08583 .
## scale(mean_activity) 0.3674 0.1147 3.202 0.00136 **
##
## Phi coefficients (precision model with identity link):
## Estimate Std. Error z value Pr(>|z|)
## (phi) 7.368 1.552 4.746 2.07e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Type of estimator: ML (maximum likelihood)
## Log-likelihood: 15.26 on 3 Df
## Pseudo R-squared: 0.2173
## Number of iterations: 14 (BFGS) + 2 (Fisher scoring)ggplot(d_avg, aes(mean_activity, shortcut, color=group)) +
geom_smooth(method='lm', colour="darkgray") +
geom_point() +
scale_color_manual(values=c('dodgerblue', 'orange')) +
ylab('Proportion shortcut') +
xlab('Cognitive control activity')
Yes, it looks like across subjects cognitive control network activity does predict the proportion of shortcuts taken (\(R^2\) = 0.22, \(t\)(29) = 3.11, \(p\) < 0.01).
What about just rep 1?
# just rep 1?
d_avg = d_gp %>%
filter(rep == 1) %>%
select(-group, -rep) %>%
group_by(subid) %>%
summarise_each(funs(mean)) %>%
left_join(d_gp %>% select(subid, group) %>% distinct())
# d_avg
summary(lm(shortcut ~ scale(mean_activity), data=d_avg))##
## Call:
## lm(formula = shortcut ~ scale(mean_activity), data = d_avg)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.36920 -0.19534 -0.02094 0.15811 0.44590
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.40625 0.03579 11.352 8.98e-14 ***
## scale(mean_activity) 0.10844 0.03624 2.992 0.00485 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2263 on 38 degrees of freedom
## Multiple R-squared: 0.1907, Adjusted R-squared: 0.1694
## F-statistic: 8.952 on 1 and 38 DF, p-value: 0.004847ggplot(d_avg, aes(mean_activity, shortcut, color=group)) +
geom_smooth(method='lm', colour="darkgray") +
geom_point() +
scale_color_manual(values=c('dodgerblue', 'orange')) +
ylab('Proportion shortcut') +
xlab('Cognitive control activity')
Same thing holds just looking at rep 1.
Does hippocampal activity when arriving at the goal on rep 1 predict p(shortcut) on rep 2?
It could be the case that hippocampal activity evoked by goal-arrival reflects replay of the previous trial, and might increase the probability of taking a shortcut on the next rep.
d_behav = get_behav_data()
dm = merge_behav_brain_data(d_behav, roi_name = 'hippocampus', event_type = 'ARRIVED')## [1] "Missing trials:"
## [1] 1 3
## # A tibble: 1 × 3
## subid rep n
## <fctr> <fctr> <dbl>
## 1 20 2 9
## subid
## 20
## 40 Levels: 3 4 5 6 7 9 10 11 12 13 14 15 16 18 19 20 21 22 23 24 26 ... 57dm = dm %>% group_by(subid) %>%
mutate(mean_activity_z = zscore(mean_activity))
dm_rep1 = dm %>% filter(rep == 1) %>%
select(subid, group, env, mean_activity_z, shortcut, final_code) %>%
rename(shortcut1=shortcut, code_1=final_code)
dm_rep2 = dm %>% filter(rep == 2) %>%
select(subid, group, env, shortcut, final_code) %>%
rename(shortcut2=shortcut, code_2=final_code)
# dim(dm_rep1)
# dim(dm_rep2)
dmm = left_join(dm_rep1, dm_rep2)
res = glmer(shortcut2 ~ mean_activity_z * scale(shortcut1)+
(1 + mean_activity_z * scale(shortcut1)|subid),
data=dmm, family='binomial',
control=glmerControl(optimizer="bobyqa",
optCtrl=list(maxfun=2e5)))
summary(res)## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## shortcut2 ~ mean_activity_z * scale(shortcut1) + (1 + mean_activity_z *
## scale(shortcut1) | subid)
## Data: dmm
## Control:
## glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 2e+05))
##
## AIC BIC logLik deviance df.resid
## 655.1 713.5 -313.6 627.1 463
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6625 -0.8234 0.4461 0.8111 1.4634
##
## Random effects:
## Groups Name Variance Std.Dev. Corr
## subid (Intercept) 0.12675 0.3560
## mean_activity_z 0.15754 0.3969 0.14
## scale(shortcut1) 0.08924 0.2987 0.04 0.26
## mean_activity_z:scale(shortcut1) 0.03567 0.1889 0.60 0.04
##
##
##
##
## -0.76
## Number of obs: 477, groups: subid, 40
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.02718 0.11816 0.230 0.818
## mean_activity_z -0.06057 0.13568 -0.446 0.655
## scale(shortcut1) 0.52551 0.11785 4.459 8.22e-06 ***
## mean_activity_z:scale(shortcut1) -0.04802 0.13095 -0.367 0.714
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mn_ct_ scl(1)
## men_ctvty_z -0.045
## scl(shrtc1) 0.056 0.123
## mn_ctv_:(1) 0.152 0.072 -0.096Whether or not subjects took a shortcut on the first repetition predicts whether they’ll take a shortcut on the 2nd rep (z=4.1, p < 0.001). However, hippocampal activity when arriving at the goal on rep 1 does not predict taking a shortcut on the following rep through that town (controlling for strategy on rep 1).
# visualize
d_avg = dmm %>% group_by(subid, group, shortcut1) %>%
summarise(mean = mean(shortcut2), n())
dfwc = summarySEwithin(d_avg, measurevar="mean", withinvars=c("shortcut1"),
betweenvars = c('group'),
idvar="subid", na.rm=TRUE, conf.interval=.95)
dfwc## group shortcut1 N mean mean_norm sd se ci
## 1 control 0 19 0.4077409 0.4340937 0.2546880 0.05842943 0.1227557
## 2 control 1 20 0.6087500 0.6246619 0.2498140 0.05586010 0.1169165
## 3 stress 0 19 0.4353156 0.4137645 0.2701817 0.06198393 0.1302234
## 4 stress 1 19 0.6714286 0.6498775 0.2701817 0.06198393 0.1302234ggplot(dfwc, aes(x=shortcut1, y=mean, group=group, fill=group)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(width=0, aes(ymin=mean-se,
ymax=mean+se), size=1.5, position=position_dodge(.9)) +
xlab('Shortcut on rep 1')+
ylab('p Shortcut on rep 2')+
scale_fill_manual(values=palette, guide = guide_legend(title = NULL))
# visualize, splitting by coding
d_avg = dmm %>% group_by(subid, group, code_1) %>%
summarise(mean = mean(shortcut2), n())
dfwc = summarySEwithin(d_avg, measurevar="mean", withinvars=c("code_1"),
betweenvars = c('group'),
idvar="subid", na.rm=TRUE, conf.interval=.95)
dfwc## group code_1 N mean mean_norm sd se ci
## 1 control backtrack 13 0.2179487 0.3013192 0.4106303 0.11388836 0.2481414
## 2 control habit 19 0.4728070 0.4875275 0.3240339 0.07433846 0.1561793
## 3 control other 16 0.4062500 0.4097474 0.3983224 0.09958059 0.2122510
## 4 control shortcut 20 0.6087500 0.6177596 0.2980373 0.06664317 0.1394857
## 5 stress backtrack 12 0.4166667 0.4338360 0.3791579 0.10945347 0.2409055
## 6 stress habit 19 0.3971178 0.3625745 0.2659713 0.06101801 0.1281941
## 7 stress other 17 0.4656863 0.4366838 0.3551032 0.08612517 0.1825772
## 8 stress shortcut 19 0.6714286 0.6368853 0.3326948 0.07632543 0.1603538ggplot(dfwc, aes(x=code_1, y=mean, group=group, fill=group)) +
geom_bar(position=position_dodge(), stat="identity") +
geom_errorbar(width=0, aes(ymin=mean-se,
ymax=mean+se), size=1.5, position=position_dodge(.9)) +
xlab('Shortcut on rep 1')+
ylab('p Shortcut on rep 2')+
scale_fill_manual(values=palette, guide = guide_legend(title = NULL))