#import and clean data
# get all filepaths for RMET behavioral data
data_files<-Sys.glob(path=paste0("/Volumes/External/WTC_resilience_imaging/data/BIDS_data/sub-*/func/*_task-RMET_dir-AP_run-0*_events.tsv"))
data_names<-gsub("/Volumes/External/WTC_resilience_imaging/data/BIDS_data/","",data_files)
data_names<-gsub("^sub-.*/func/", "", data_names)
data_names<-gsub("*_task-RMET_dir-AP_run-01_events.tsv", ".1", data_names)
data_names<-gsub("*_task-RMET_dir-AP_run-02_events.tsv", ".2", data_names)
data_names<-gsub("sub-", "sub", data_names)
for(i in 1:length(data_files)) { # Head of for-loop
assign(paste0(noquote(data_names[i])), # Read and store data frames
read.table(data_files[i],sep = "\t", header = TRUE))
}
## Warning in read.table(data_files[i], sep = "\t", header = TRUE):
## incomplete final line found by readTableHeader on '/Volumes/External/
## WTC_resilience_imaging/data/BIDS_data/sub-061/func/sub-061_task-RMET_dir-
## AP_run-01_events.tsv'
## Warning in read.table(data_files[i], sep = "\t", header = TRUE):
## incomplete final line found by readTableHeader on '/Volumes/External/
## WTC_resilience_imaging/data/BIDS_data/sub-061/func/sub-061_task-RMET_dir-
## AP_run-02_events.tsv'
rm(sub061.1,sub061.2)
new<-combine(sub002.1, sub002.2, sub003.1, sub003.2, sub004.1, sub004.2, sub005.1,
sub005.2, sub006.1, sub006.2, sub007.1, sub007.2, sub008.1, sub008.2, sub009.1,
sub009.2, sub010.1, sub010.2, sub011.1, sub011.2, sub012.1, sub012.2, sub013.1,
sub013.2, sub014.1, sub014.2, sub015.1, sub015.2, sub016.1, sub016.2, sub017.1,
sub017.2, sub018.1, sub018.2, sub019.1, sub019.2, sub020.1, sub020.2, sub021.1,
sub021.2, sub022.1, sub022.2, sub023.1, sub023.2, sub024.1, sub024.2, sub025.1, sub025.2, sub026.1,
sub026.2, sub027.1, sub027.2, sub028.1, sub028.2, sub029.1, sub029.2, sub030.1,
sub030.2, sub031.1, sub031.2, sub032.1, sub032.2, sub033.1, sub033.2, sub034.1,
sub034.2, sub035.1, sub035.2, sub036.1, sub036.2, sub037.1, sub037.2, sub038.1,
sub038.2, sub039.1, sub039.2, sub040.1, sub040.2, sub041.1, sub041.2, sub042.1,
sub042.2, sub043.1, sub043.2, sub044.1, sub044.2, sub045.1, sub045.2, sub046.1,
sub046.2, sub047.1, sub047.2, sub048.1, sub048.2, sub049.1, sub049.2, sub050.1,
sub050.2, sub051.1, sub051.2, sub052.1, sub052.2, sub053.1, sub053.2, sub054.1,
sub054.2, sub055.1, sub055.2, sub056.1, sub056.2, sub057.1, sub057.2, sub058.1,
sub058.2, sub059.1, sub059.2, sub060.1, sub060.2, sub062.1, sub062.2)
# add variables for subject ID and run
new1 <- new %>%
filter(str_detect(source, ".1$")) %>% mutate(run = "run 1",
bids_id = gsub(".1$","",source),
bids_id = gsub("sub","sub-", bids_id))
new2 <- new %>%
filter(str_detect(source, ".2$")) %>% mutate(run = "run 2",
bids_id = gsub(".2$","",source),
bids_id = gsub("sub","sub-", bids_id))
# define words for social trials
social = c("cautious",
"regretful",
"sceptical",
"anticipating",
"accusing",
"contemplative",
"thoughtful",
"doubtful",
"decisive",
"affectionate",
"impatient",
"grateful",
"ashamed",
"bewildered",
"guilty",
"baffled",
"puzzled",
"indecisive",
"playful",
"upset",
"desire",
"insisting",
"worried",
"fantasizing",
"uneasy",
"despondent",
"preoccupied",
"arrogant",
"guilty",
"confused",
"imploring",
"curious",
"hostile",
"incredulous",
"shy")
# combine run1 and run2 dfs and fix variables
# "resp_speed" is actually the trial-level RT (in seconds) variable
# # *not* "resp_rt"
# from .py:
# # trials.addData('resp_speed', key_resp_trial.rt)
# # trials.addData('resp_rt', key_resp_trial.rt+timer)
dat1<-bind_rows(new1,new2) %>% mutate(trialtype = as.factor(ifelse(!Awords %in% social, "non-social", "social")),
missed = as.factor(ifelse(is.na(resp_speed),1,0)),
resp_RT = resp_speed*1000) %>%
select(c("bids_id","run","trialtype","resp_RT","resp_accuracy","missed"))
hist(dat1$resp_RT)
min(na.omit(dat1$resp_RT))
## [1] 510.096
max(na.omit(dat1$resp_RT))
## [1] 2982.781
# all trials >500ms, so not really concerned about misfire responses (usually <200-300ms)
rmetLong <- dat1 %>%
group_by(bids_id,run,trialtype) %>% filter(missed==0 & resp_accuracy == 1) %>%
summarize(min = min(na.omit(resp_RT)),
max = max(na.omit(resp_RT)),
meanRT = mean(na.omit(resp_RT)),
medianRT=median(na.omit(resp_RT)),
sd = sd(na.omit(resp_RT)))
## `summarise()` has grouped output by 'bids_id', 'run'. You can override using the `.groups` argument.
# accuracy = n correct / n total trials per subj/run/condition
rmetAcc <- dat1 %>%
group_by(bids_id,run,trialtype) %>%
summarize(accuracy=sum(resp_accuracy==1)/(sum(resp_accuracy==1,resp_accuracy==0)))
## `summarise()` has grouped output by 'bids_id', 'run'. You can override using the `.groups` argument.
rmetLong<-left_join(rmetLong,rmetAcc, by=c("bids_id","run","trialtype"))
# if we want to get accuracy counting missed trials as errors...
#acc <- dat1 %>%
# group_by(bids_id,run,trialtype) %>%
# summarize(accuracy=sum(resp_accuracy==1)/(sum(resp_accuracy==1,resp_accuracy==0)))
#left_join(rmetLong,acc, by=c("bids_id","run","trialtype"))
hist(rmetLong$accuracy)
# QUESTION: what to do with people who have accuracy <50%?
# for group variable
grp<-read.table("/Volumes/External/WTC_resilience_imaging/data/BIDS_data/participantscurrent.tsv", sep = "\t", header = TRUE)
grp <- grp %>% dplyr::filter(!BIDS_ID=="") %>% dplyr::rename(bids_id = BIDS_ID, group = group.factor.y) %>%
dplyr::select(c(bids_id, group, age, gender_f, handed_L,traditional_responder.factor)) %>% mutate(group = factor(group, ordered = FALSE),
group_2grp = ifelse(group=="Symptomatic (PTSD)", "PTSD","Non-PTSD"),
group_2grp = as.factor(group_2grp),
group_2grp = relevel(group_2grp, ref="Non-PTSD"))
rmetDat<-left_join(grp,rmetLong, by="bids_id") %>% filter(!is.na(run)) # get rid of subjects w/o behavioral data
write.csv(rmetDat,"/Volumes/External/WTC_resilience_imaging/data/behav_data/behav_data_RMET_113021.csv")
saveRDS(rmetDat,"/Volumes/External/WTC_resilience_imaging/data/behav_data/behav_data_RMET_113021.rds")
Descriptives by group
## `summarise()` has grouped output by 'trialtype', 'group'. You can override using the `.groups` argument.
Comparisons: run & trial type
Reaction time
##
## Call:
## lm(formula = meanRT ~ trialtype * run, data = rmetDat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -568.54 -166.14 0.87 158.15 919.93
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1209.95 32.18 37.598 < 2e-16 ***
## trialtypesocial 285.98 45.51 6.284 1.58e-09 ***
## runrun 2 127.80 45.51 2.808 0.0054 **
## trialtypesocial:runrun 2 -108.80 64.36 -1.690 0.0923 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 249.3 on 236 degrees of freedom
## Multiple R-squared: 0.2023, Adjusted R-squared: 0.1921
## F-statistic: 19.95 on 3 and 236 DF, p-value: 1.472e-11
##
## Call:
## lm(formula = medianRT ~ trialtype * run, data = rmetDat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -644.33 -171.66 -0.67 151.21 1214.34
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1150.74 33.46 34.387 < 2e-16 ***
## trialtypesocial 306.59 47.33 6.478 5.34e-10 ***
## runrun 2 96.50 47.33 2.039 0.0426 *
## trialtypesocial:runrun 2 -95.07 66.93 -1.420 0.1568
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 259.2 on 236 degrees of freedom
## Multiple R-squared: 0.2136, Adjusted R-squared: 0.2036
## F-statistic: 21.36 on 3 and 236 DF, p-value: 2.815e-12
Accuracy
##
## Call:
## lm(formula = accuracy ~ trialtype * run, data = rmetDat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.283333 -0.066667 0.005556 0.074074 0.250000
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.894444 0.013329 67.103 <2e-16 ***
## trialtypesocial -0.200000 0.018851 -10.610 <2e-16 ***
## runrun 2 -0.024074 0.018851 -1.277 0.203
## trialtypesocial:runrun 2 -0.009259 0.026659 -0.347 0.729
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1032 on 236 degrees of freedom
## Multiple R-squared: 0.5047, Adjusted R-squared: 0.4984
## F-statistic: 80.14 on 3 and 236 DF, p-value: < 2.2e-16
Group comparison: 2 groups
Reaction time
##
## Call:
## lm(formula = meanRT ~ group_2grp * trialtype + age, data = rmetDat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -645.55 -132.38 16.05 138.27 509.46
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 376.384 136.678 2.754 0.00636 **
## group_2grpPTSD 92.251 44.118 2.091 0.03762 *
## trialtypesocial 229.182 33.195 6.904 4.84e-11 ***
## age 16.005 2.483 6.447 6.59e-10 ***
## group_2grpPTSD:trialtypesocial 13.871 61.841 0.224 0.82273
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 215.1 on 231 degrees of freedom
## Multiple R-squared: 0.3321, Adjusted R-squared: 0.3206
## F-statistic: 28.72 on 4 and 231 DF, p-value: < 2.2e-16
##
## Call:
## lm(formula = medianRT ~ group_2grp * trialtype + age, data = rmetDat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -681.77 -133.03 8.59 137.56 630.26
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 420.727 141.366 2.976 0.00323 **
## group_2grpPTSD 81.257 45.631 1.781 0.07627 .
## trialtypesocial 257.876 34.334 7.511 1.28e-12 ***
## age 13.808 2.568 5.378 1.85e-07 ***
## group_2grpPTSD:trialtypesocial 9.782 63.962 0.153 0.87858
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 222.5 on 231 degrees of freedom
## Multiple R-squared: 0.328, Adjusted R-squared: 0.3163
## F-statistic: 28.19 on 4 and 231 DF, p-value: < 2.2e-16
Accuracy
##
## Call:
## lm(formula = accuracy ~ group_2grp * trialtype + age, data = rmetDat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30487 -0.06509 0.01435 0.06630 0.25069
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9441536 0.0650957 14.504 <2e-16 ***
## group_2grpPTSD -0.0009758 0.0210119 -0.046 0.9630
## trialtypesocial -0.1865079 0.0158099 -11.797 <2e-16 ***
## age -0.0011409 0.0011824 -0.965 0.3356
## group_2grpPTSD:trialtypesocial -0.0553221 0.0294531 -1.878 0.0616 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1025 on 231 degrees of freedom
## Multiple R-squared: 0.5074, Adjusted R-squared: 0.4988
## F-statistic: 59.48 on 4 and 231 DF, p-value: < 2.2e-16
Group comparison: 3 groups
- Resilient group is reference level.
Reaction time
##
## Call:
## lm(formula = meanRT ~ group * trialtype + age, data = rmetDat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -646.20 -130.03 20.49 140.77 508.39
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 354.3955 142.8202 2.481 0.0138
## groupLow-Exposed Control 40.5450 54.1483 0.749 0.4548
## groupSymptomatic (PTSD) 103.3639 46.7561 2.211 0.0280
## trialtypesocial 242.3736 38.7561 6.254 1.94e-09
## age 16.2143 2.5455 6.370 1.03e-09
## groupLow-Exposed Control:trialtypesocial -50.3662 75.7300 -0.665 0.5067
## groupSymptomatic (PTSD):trialtypesocial 0.6795 65.1232 0.010 0.9917
##
## (Intercept) *
## groupLow-Exposed Control
## groupSymptomatic (PTSD) *
## trialtypesocial ***
## age ***
## groupLow-Exposed Control:trialtypesocial
## groupSymptomatic (PTSD):trialtypesocial
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 215.8 on 229 degrees of freedom
## Multiple R-squared: 0.3339, Adjusted R-squared: 0.3164
## F-statistic: 19.13 on 6 and 229 DF, p-value: < 2.2e-16
##
## Call:
## lm(formula = medianRT ~ group * trialtype + age, data = rmetDat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -681.78 -132.71 8.58 132.35 639.36
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 417.251 147.879 2.822 0.0052
## groupLow-Exposed Control 12.726 56.066 0.227 0.8206
## groupSymptomatic (PTSD) 84.596 48.412 1.747 0.0819
## trialtypesocial 264.440 40.129 6.590 2.99e-10
## age 13.810 2.636 5.240 3.64e-07
## groupLow-Exposed Control:trialtypesocial -25.063 78.412 -0.320 0.7495
## groupSymptomatic (PTSD):trialtypesocial 3.218 67.430 0.048 0.9620
##
## (Intercept) **
## groupLow-Exposed Control
## groupSymptomatic (PTSD) .
## trialtypesocial ***
## age ***
## groupLow-Exposed Control:trialtypesocial
## groupSymptomatic (PTSD):trialtypesocial
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 223.4 on 229 degrees of freedom
## Multiple R-squared: 0.3283, Adjusted R-squared: 0.3107
## F-statistic: 18.65 on 6 and 229 DF, p-value: < 2.2e-16
Accuracy
##
## Call:
## lm(formula = accuracy ~ group * trialtype + age, data = rmetDat1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.29514 -0.06851 0.01423 0.06733 0.26042
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.9566450 0.0680414 14.060 <2e-16
## groupLow-Exposed Control -0.0122854 0.0257969 -0.476 0.6344
## groupSymptomatic (PTSD) -0.0045962 0.0222752 -0.206 0.8367
## trialtypesocial -0.1863799 0.0184639 -10.094 <2e-16
## age -0.0013119 0.0012127 -1.082 0.2805
## groupLow-Exposed Control:trialtypesocial -0.0004888 0.0360788 -0.014 0.9892
## groupSymptomatic (PTSD):trialtypesocial -0.0554501 0.0310255 -1.787 0.0752
##
## (Intercept) ***
## groupLow-Exposed Control
## groupSymptomatic (PTSD)
## trialtypesocial ***
## age
## groupLow-Exposed Control:trialtypesocial
## groupSymptomatic (PTSD):trialtypesocial .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1028 on 229 degrees of freedom
## Multiple R-squared: 0.5084, Adjusted R-squared: 0.4955
## F-statistic: 39.46 on 6 and 229 DF, p-value: < 2.2e-16
