Load libraries

Load Data

Demographic Data
##  [1] "SubjectID"     "Adate"         "Age"           "MosEmployLas" 
##  [5] "ConsentDate"   "Cohort"        "Outcome"       "DOB"          
##  [9] "SexAtBirth"    "Ethnicity"     "Race"          "MaritalStatus"
## [13] "BioChildren"   "Child1Age"     "Child2Age"     "Child3Age"    
## [17] "Child4Age"     "Child5Age"     "Child6Age"     "Dependents"   
## [21] "Education"     "EmploymentSta" "EmploymentSt0" "HoursPerWeek" 
## [25] "EmploymentFul" "EmploymentLas" "AnnualIncome"  "ExclusionSta" 
## [29] "Race_Text"     "Minority_Text"
Behavioral Data
RSA data
Create wideform summary dfs

Behavioral Data Analysis

Descriptive Tables:

Emotion Ratings

summary_allobs_emo <- describeBy(emo_resp ~ PicValence + Procedure, data=cert_mna, mat=T)

htmlTable::htmlTable(format(summary_allobs_emo, 
                            digits = 2))
item group1 group2 vars n mean sd median trimmed mad min max range skew kurtosis se
emo_resp1 1 neg reg 1 2178 3.1 1.19 3 3.2 1.5 1 5 4 -0.024 -0.92 0.026
emo_resp2 2 neut reg 1 2153 1.3 0.58 1 1.1 0.0 1 5 4 2.537 7.38 0.013
emo_resp3 3 neg watch 1 2183 3.7 1.15 4 3.8 1.5 1 5 4 -0.562 -0.59 0.025
emo_resp4 4 neut watch 1 2163 1.3 0.65 1 1.1 0.0 1 5 4 2.606 7.52 0.014

Thinking Change Ratings

summary_allobs_er <- describeBy(er_resp ~ PicValence + Procedure, data=cert_mna, mat=T)
#print(summary_allobs)
htmlTable::htmlTable(format(summary_allobs_er, 
                            digits = 2))
item group1 group2 vars n mean sd median trimmed mad min max range skew kurtosis se
er_resp1 1 neg reg 1 2180 3.1 1.15 3 3.1 1.5 1 5 4 -0.091 -0.80 0.025
er_resp2 2 neut reg 1 2168 2.7 1.13 3 2.7 1.5 1 5 4 0.159 -0.73 0.024
er_resp3 3 neg watch 1 2151 1.7 1.07 1 1.5 0.0 1 5 4 1.489 1.44 0.023
er_resp4 4 neut watch 1 2144 1.3 0.74 1 1.1 0.0 1 5 4 2.719 7.88 0.016

Visualizations:

On average, the group found negative images more negative

and they changed their thinking more during the reappraise trials

sumdf_l$Valence <- factor(sumdf_l$Valence, 
                          levels=c("neg","neut"), 
                          labels=c("Negative","Neutral"))
sumdf_l$Condition <- factor(sumdf_l$Condition, 
                            levels=c("reg","watch"), 
                            labels=c("Regulate","Watch"))
sumdf_l$Rating_Type <- factor(sumdf_l$Rating_Type, 
                            levels=c("emo_resp","er_resp"), 
                            labels=c("Emotion","Thinking Change"))


  ggplot(data = sumdf_l, aes(x = Condition, y = Rating, color = Valence)) +
  geom_boxplot(aes(fill=Valence), alpha = .5) +
  facet_wrap(~Rating_Type) +
  labs(y = "Rating", x = "Condition") +
  ggtitle("Ratings across Valence and Conditions") +
  theme_minimal() 

# emosum <- ggplot(data = sumdf, aes(x = Procedure, y = emo_resp, color = PicValence)) +
#   geom_boxplot(aes(fill=PicValence), alpha = .5) +
#   labs(y = "Emotion Rating", x = "Condition") +
#   ggtitle("Ratings across Valence and Conditions") +
#   theme(legend.position = "none") 
# 
# regsum <- ggplot(data = sumdf, aes(x = Procedure, y = er_resp, color = PicValence)) +
#   geom_boxplot(aes(fill=PicValence), alpha = .5) +
#   labs(y = "Regulation Rating", x = "Condition") +
#   ggtitle(" ") 

#grid.arrange(emosum,regsum,nrow=1)

Distribution of Emotion Ratings on Negative trials:

Purple = regulate

Blue = watch

Green = difference between reg and watch

sumdf_emo$cond22<- paste(sumdf_emo$PicValence,sumdf_emo$Procedure,sep="_")
sumdf_reg$cond22<- paste(sumdf_reg$PicValence,sumdf_reg$Procedure,sep="_")

sumdf_emo <- sumdf_emo[,-c(2,3)]
sumdf_reg<- sumdf_reg[,-c(2,3)]

sumdf_emo_w<-spread(sumdf_emo, key=cond22, value=emo_resp)
sumdf_reg_w<-spread(sumdf_reg, key=cond22, value=er_resp)

 sumdf_emo_w$diff <- sumdf_emo_w$neg_watch - sumdf_emo_w$neg_reg
 sumdf_reg_w$diff <- sumdf_reg_w$neg_watch - sumdf_reg_w$neg_reg

negemodens<-ggplot(data = sumdf_emo_w) +
  geom_density(aes(x=neg_reg), alpha = .4, fill = "plum4")+
    geom_density(aes(x=neg_watch), alpha = .4, fill = "cadetblue4") +
      geom_density(aes(x=diff), alpha = .4, fill = "palegreen4")
negemodens + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))

Distribution of Thinking Change Ratings on Negative trials:

negemoregdens<-ggplot(data = sumdf_reg_w) +
  geom_density(aes(x=neg_reg), alpha = .4, fill = "plum1")+
    geom_density(aes(x=neg_watch), alpha = .4, fill = "cadetblue1") +
      geom_density(aes(x=diff), alpha = .4, fill = "palegreen1")
negemoregdens + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))

#Decompose within and between person effects for model
cert_mna$emo_win <- calc.mcent(emo_resp, subject, data=cert_mna)
cert_mna$emo_bw <- calc.mean(emo_resp, subject, data=cert_mna, expand=TRUE)

cert_mna$reg_win <- calc.mcent(er_resp, subject, data=cert_mna)
cert_mna$reg_bw <- calc.mean(er_resp, subject, data=cert_mna, expand=TRUE)

Individuals also rate the negative images more negative during Watch compared to Regulate

cert_mna_neg <- subset(cert_mna, cert_mna$PicValence=="neg")
sumdf_neg <- subset(sumdf, sumdf$PicValence=="neg")

# win <- ggplot(data = cert_mna_neg, aes(x = Procedure, y= emo_win)) + 
#         geom_smooth(aes(group = subject, color = subject), alpha = .4, size = .5, method = "lm", se = F) +
#         geom_smooth(aes(group = 1), color = "black",size=1.5, method = "lm", se = F) +
#         scale_color_gradient(low = "indianred4", high = "indianred2") +
#         ylim(1,3.25) +
#         ylab("Within-person Centered Emotion Ratings") +
#         xlab("Condition") +
#         theme_minimal()
#   
# win <- win + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
# panel.background = element_blank(), axis.line = element_line(colour = "black"))

bw <- ggplot(data = sumdf_neg, aes(x = Procedure, y= emo_resp, group=subject,color=subject)) +
       geom_line(alpha = .5, size = .5) +
        geom_smooth(aes(group = 1), color = "black",size=1.5, method = "lm", se = F) +
        scale_color_gradient(low = "deeppink4", high = "deeppink2") +
        #scale_fill_gradient(low = "indianred4", high = "indianred2") 
        xlab("Condition") +
        ylab("Person-mean Emotion Ratings on Negative Trials")
bw <- bw + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))

bw

Individuals also more thinking change during Regulate compared to Watch

bw2 <- ggplot(data = sumdf_neg, aes(x = Procedure, y= er_resp, group=subject,color=subject)) +
       geom_line(alpha = .5, size = .5) +
        geom_smooth(aes(group = 1), color = "black",size=1.5, method = "lm", se = F) +
        scale_color_gradient(low = "darkblue", high = "skyblue3") +
        #scale_fill_gradient(low = "indianred4", high = "indianred2") 
        xlab("Condition") +
        ylab("Person-mean Thinking Change Ratings on Negative Trials")
bw2 <- bw2 + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))

bw2

Confirm visual effects with HLM

longest_emo <- subset(longest,longest$Rating_Type=="emo_resp")
longest_reg <- subset(longest,longest$Rating_Type=="er_resp")

longest_emo$emo_win <- calc.mcent(Rating, Subject, data=longest_emo)
longest_emo$emo_bw <- calc.mean(Rating, Subject, data=longest_emo, expand=TRUE)

longest_reg$reg_win <- calc.mcent(Rating, Subject, data=longest_reg)
longest_reg$reg_bw <- calc.mean(Rating, Subject, data=longest_reg, expand=TRUE)

model1 <- lmer(Rating ~ Valence*Condition + (1|Subject), data=longest_emo)
#summary(model1)
#htmlTable::htmlTable(format(model1, digits = 2))
#sjPlot::tab_model(model1, p.val = "kr", show.df = TRUE)

gtsummary::tbl_regression(model1)
Characteristic Beta 95% CI1 p-value
Valence


    neg — —
    neut -1.9 -1.9, -1.8 <0.001
Condition


    reg — —
    watch 0.56 0.51, 0.61 <0.001
Valence * Condition


    neut * watch -0.54 -0.61, -0.47 <0.001
1 CI = Confidence Interval

How correlated are Emotion ratings and Thinking Change ratings for Negative Regulate trials?

Visually, there is a negative relationship (higher thinking change ~ lower emotions)
sumdf_neg_reg <- subset(sumdf_neg,sumdf_neg$Procedure=="reg")
emo_er <- ggplot(data = sumdf_neg_reg, aes(x = er_resp, y= emo_resp, group=subject,color=subject)) +
       geom_point(alpha = .5, size = 1.5) +
        geom_smooth(aes(group = 1), color = "black",size=1.5, method = "lm", se = F) +
        scale_color_gradient(low = "orchid4", high = "orchid2") +
        #scale_fill_gradient(low = "indianred4", high = "indianred2") 
        xlab("Emotion Rating") +
        ylab("Thinking Change")
emo_er <- emo_er + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))

emo_er
## `geom_smooth()` using formula = 'y ~ x'

Statisically, it is not significant

cor.test(sumdf_neg_reg$emo_resp,sumdf_neg_reg$er_resp)
## 
##  Pearson's product-moment correlation
## 
## data:  sumdf_neg_reg$emo_resp and sumdf_neg_reg$er_resp
## t = -1.4378, df = 111, p-value = 0.1533
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.3121514  0.0507794
## sample estimates:
##        cor 
## -0.1352188

RSA Data Analysis:

ROIs from Neurosynth
ROIs from Neurosynth
Across all people and all runs: RDMs for each ROI – NeutralW, Neutral R, NegW, Neg R
Across all people and all runs: RDMs for each ROI – NeutralW, Neutral R, NegW, Neg R

Fisher Transform the correlations (don’t convert back to a distance)

Untransformed values…

# Make df ever longer by gathering the difference pairwise condition comparisons and letting there be just 1 dissim variable
cert_rsa_l <- gather(cert_rsa, key = "comparison", value = "dissimilarity", 
                  NtW_by_NtW:NgR_by_NgR, factor_key = TRUE,
                  -subject, -ROI, -session)

# for the wide format data...
cert_rsa$NtR_by_NtW_similarity_fz <- fisherz(1 - cert_rsa$NtR_by_NtW)
cert_rsa$NgW_by_NtW_similarity_fz <- fisherz(1 - cert_rsa$NgW_by_NtW)
cert_rsa$NgW_by_NtR_similarity_fz <- fisherz(1 - cert_rsa$NgW_by_NtR)
cert_rsa$NgR_by_NtW_similarity_fz <- fisherz(1 - cert_rsa$NgR_by_NtW)
cert_rsa$NgR_by_NgW_similarity_fz <- fisherz(1 - cert_rsa$NgR_by_NgW)
cert_rsa$NtR_by_NtR_similarity_fz <- fisherz(1 - cert_rsa$NgR_by_NtR)

# IF we want fisher transformed values:
cert_rsa_l$dissimilarity_corr <- (1 - cert_rsa_l$dissimilarity)
  cert_rsa_l$dissimilarity_fz <- (1 - fisherz(cert_rsa_l$dissimilarity_corr))
cert_rsa_l$similarity_fz <- fisherz(cert_rsa_l$dissimilarity_corr)

cert_rsa_l_red <- subset(cert_rsa_l,
                         !(cert_rsa_l$comparison %in% c("NtW_by_NtW", "NtR_by_NtR", "NgW_by_NgW", "NgR_by_NgR")))
plot(cert_rsa_l_red$dissimilarity)

hist(cert_rsa_l_red$dissimilarity)

F Transformed (representing similarity now):

#plot(cert_rsa_l_red$dissimilarity_fz)
#hist(cert_rsa_l_red$dissimilarity_fz)

plot(cert_rsa_l_red$similarity_fz)

hist(cert_rsa_l_red$similarity_fz)

Descriptive Table

(commented out for ease of interpretation)

summary_rsa <- describeBy(similarity_fz ~ comparison + ROI, data=cert_rsa_l_red, mat=T)
#print(summary_allobs)
#htmlTable::htmlTable(format(summary_rsa, 
#                            digits = 2))

Visualizations:

How do the conditions compare to each other across ROIs?

 cert_rsa_l_red$comparison <- as.factor(cert_rsa_l_red$comparison)
 rois1234<-subset(cert_rsa_l_red,cert_rsa_l_red$ROI<5)
 rois5678<-subset(cert_rsa_l_red,cert_rsa_l_red$ROI>4 & cert_rsa_l_red$ROI <9)
 rois91011<-subset(cert_rsa_l_red,cert_rsa_l_red$ROI>8)

 allROIs<-ggplot(data = rois1234, aes(x = comparison, y = similarity_fz, color = comparison)) +
  geom_violin(aes(fill=comparison), alpha = .5) +
  stat_summary(fun=mean, color="black", geom="point", size = 1) + 
  facet_wrap(~ROI, scales = "free", ncol = 2) +
  labs(y = "", x = "Comparison") +
  #ggtitle("Ratings across Valence and Conditions") +
  theme_minimal() 

allROIs + theme(
  axis.text.x = element_text(angle = 45, hjust = 1, size=8),  # rotate x-axis labels if needed
  strip.text = element_text(size = 8),  # adjust the size of facet labels
  strip.background = element_blank())

 allROIs<-ggplot(data = rois5678, aes(x = comparison, y = similarity_fz, color = comparison)) +
  geom_violin(aes(fill=comparison), alpha = .5) +
  stat_summary(fun=mean, color="black", geom="point", size = 1) + 
  facet_wrap(~ROI, scales = "free", ncol = 2) +
  labs(y = "Dissimilarity (lower is more alike)", x = "Comparison") +
  #ggtitle("Ratings across Valence and Conditions") +
  theme_minimal() 

allROIs + theme(
  axis.text.x = element_text(angle = 45, hjust = 1, size=8),  # rotate x-axis labels if needed
  strip.text = element_text(size = 8),  # adjust the size of facet labels
  strip.background = element_blank())

 allROIs<-ggplot(data = rois91011, aes(x = comparison, y = similarity_fz, color = comparison)) +
  geom_violin(aes(fill=comparison), alpha = .5) +
  stat_summary(fun=mean, color="black", geom="point", size = 1) + 
  facet_wrap(~ROI, scales = "free", ncol = 2) +
  labs(y = "Dissimilarity (lower is more alike)", x = "Comparison") +
  #ggtitle("Ratings across Valence and Conditions") +
  theme_minimal() 

allROIs + theme(
  axis.text.x = element_text(angle = 45, hjust = 1, size=8),  # rotate x-axis labels if needed
  strip.text = element_text(size = 8),  # adjust the size of facet labels
  strip.background = element_blank())

These full plots are a little hard to take in though

Another way visualize: Plot the density of

Negative_Regulate vs Negative Watch (purple)

Negative_Regulate vs Neutral Watch (gold)

(notice Neg_Reg is more similar to Neutral_Watch than Negative_Watch)

# further reduce conditions to 3 from 6
cert_rsa_l_redcons <- subset(cert_rsa_l_red,
                         (cert_rsa_l_red$comparison %in% c("NgR_by_NgW", "NtR_by_NtW", "NgR_by_NtW")))

# using the wide form dataset, plot the density distribution for the conditions (not transformed btw)
 
cert_rsa_rois1234<-subset(cert_rsa,cert_rsa$ROI<5)
 cert_rsa_rois5678<-subset(cert_rsa,cert_rsa$ROI>4 & cert_rsa$ROI <9)
 cert_rsa_rois91011<-subset(cert_rsa,cert_rsa$ROI>8)

 density1<-ggplot(data = cert_rsa_rois1234) +
  geom_density(aes(x=NgR_by_NgW_similarity_fz), fill="slateblue2", alpha = .5) +
  geom_density(aes(x=NgR_by_NtW_similarity_fz), fill="goldenrod", alpha = .5) +
  facet_wrap(~ROI, scales = "free", ncol = 2) +
  theme_minimal() 

density1 + theme(
  axis.text.x = element_text(angle = 45, hjust = 1, size=8),  # rotate x-axis labels if needed
  strip.text = element_text(size = 8),  # adjust the size of facet labels
  strip.background = element_blank())

density1<-ggplot(data = cert_rsa_rois5678) +
  geom_density(aes(x=NgR_by_NgW_similarity_fz), fill="slateblue2", alpha = .5) +
  geom_density(aes(x=NgR_by_NtW_similarity_fz), fill="goldenrod", alpha = .5) +
  facet_wrap(~ROI, scales = "free", ncol = 2) +
  theme_minimal() 

density1 + theme(
  axis.text.x = element_text(angle = 45, hjust = 1, size=8),  # rotate x-axis labels if needed
  strip.text = element_text(size = 8),  # adjust the size of facet labels
  strip.background = element_blank())

density1<-ggplot(data = cert_rsa_rois91011) +
  geom_density(aes(x=NgR_by_NgW_similarity_fz), fill="slateblue2", alpha = .5) +
  geom_density(aes(x=NgR_by_NtW_similarity_fz), fill="goldenrod", alpha = .5) +
  facet_wrap(~ROI, scales = "free", ncol = 2) +
  theme_minimal() 

density1 + theme(
  axis.text.x = element_text(angle = 45, hjust = 1, size=8),  # rotate x-axis labels if needed
  strip.text = element_text(size = 8),  # adjust the size of facet labels
  strip.background = element_blank())

Is this difference statistically significant?

Yes, all comparisons are less similar than reference (NtR_by_NtW)

There is a lot of overlap between the distributions at the group level, what are the correlations like between conditions?

First 4 ROIs plotted as examples: significant correlations but also unaccounted for variance as well

for (roi in 1:4){
  tmpdata <- subset(cert_rsa, cert_rsa$ROI == roi)
  #create correlation matrix for fz similarity vals of interest
  cor1_tab <- rcorr(as.matrix(tmpdata[,14:19]),type="pearson")
  
  r1 <- as.data.frame(cor1_tab$r)
  #print(r1)
  p1 <- as.data.frame(cor1_tab$P)
  p.mat1 <- cor_pmat(tmpdata[,14:19])

  # FDR correct the lower triangle values
  p.mat.tri  <- p.mat1[lower.tri(p.mat1)]
  p.mat.tri.fdr <- p.adjust(p.mat.tri, method="fdr")
  p.mat1[lower.tri(p.mat1)] <- p.mat.tri.fdr


print(paste("plot of correlations for ROI #", roi,sep=""))
 plot_roi <- ggcorrplot(r1, method = "square", ggtheme = ggplot2::theme_minimal, title = " ", show.legend = TRUE, legend.title = "Corr", show.diag = FALSE,
  colors = c("turquoise4","white", "violetred4"), outline.color = "black",
  hc.order = F, hc.method = "pairwise", lab = T,
   p.mat = p.mat1, sig.level = 0.05,
  insig = c("pch"), pch = 4, pch.col = "black",
  pch.cex = 5, tl.cex = 10, tl.col = "black", tl.srt = 45,
  digits = 2)
 print(plot_roi)

}
## [1] "plot of correlations for ROI #1"

## [1] "plot of correlations for ROI #2"

## [1] "plot of correlations for ROI #3"

## [1] "plot of correlations for ROI #4"

RSA & Behavioral Data Analysis:

Merge the RSA and behavioral data together:

## Prep behavioral data for merge: 
#calc mean for each person each condition:
personmn_emo <- describeBy(emo_resp ~ PicValence + Procedure + subject, data=cert_mna, mat=T)
personmn_emo <- personmn_emo[,c(2:4,6,7)]
personmn_emo$condition <- paste(personmn_emo$group1, personmn_emo$group2, sep = "_")
personmn_n <- personmn_emo[,c(3,4,6)]
personmn_emo <- personmn_emo[,c(3,5,6)]

personmn_think <- describeBy(er_resp ~ PicValence + Procedure + subject, data=cert_mna, mat=T)
personmn_think <- personmn_think[,c(2:4,6,7)]
personmn_think$condition <- paste(personmn_think$group1, personmn_think$group2, sep = "_")
personmn_think <- personmn_think[,c(3,5,6)]

#spread to wide form 
colnames(personmn_n) <- c("subject","number","ntrials")
personmn_n_wide = personmn_n %>% 
  spread(ntrials, number, sep = "_") #ntrials actually = conditions, but labelled weird for naming purposes 

colnames(personmn_emo) <- c("subject","value","AvgEmo")
personmn_emo_wide = personmn_emo %>% 
  spread(AvgEmo, value, sep = "_") #ntrials actually = conditions, but labelled weird for naming purposes 

colnames(personmn_think) <- c("subject","value","AvgThink")
personmn_think_wide = personmn_think %>% 
  spread(AvgThink, value, sep = "_") #ntrials actually = conditions, but labelled weird for naming purposes 

prsnmn_all <- merge(personmn_n_wide,personmn_emo_wide, by = "subject")
prsnmn_all <- merge(prsnmn_all,personmn_think_wide, by = "subject")
#rsa and behavior
  #long
  cert_rsa_l_red_beh <- merge(cert_rsa_l_red, prsnmn_all, by="subject", all.x = T)
  #wide
  cert_rsa_beh_w <- merge(cert_rsa, prsnmn_all, by="subject", all.x = T)
  
#now add demogrpahic data
  #long
  cert_rsa_beh_dem <- merge(cert_rsa_l_red_beh, MNA, by.x ="subject", by.y = "SubjectID", all.x = T)
  #wide
  cert_rsa_beh_dem_w <- merge(cert_rsa_beh_w, MNA, by.x ="subject", by.y = "SubjectID", all.x = T)
#created an averaged value across sessions:
  #long 
  cert_rsa_beh_dem_NegRxNegW <- subset(cert_rsa_beh_dem, cert_rsa_beh_dem$comparison == "NgR_by_NgW")
  cert_rsa_beh_dem_NegRxNegW$ROI <- as.factor(cert_rsa_beh_dem_NegRxNegW$ROI)
  
  avgsim2sess <- describeBy(similarity_fz ~ ROI + subject, data=cert_rsa_beh_dem_NegRxNegW, mat=T)
  avgsim2sess <- avgsim2sess[,c(2,3,6)]
  colnames(avgsim2sess) <- c("ROI","subject","avgsimilarity_fz")
  
  cert_rsa_beh_dem_NegRxNegW_avg <- merge(cert_rsa_beh_dem_NegRxNegW,avgsim2sess, by=c("subject", "ROI"), all.x = T)
  cert_rsa_beh_dem_NegRxNegW_avg <- subset(cert_rsa_beh_dem_NegRxNegW_avg,cert_rsa_beh_dem_NegRxNegW_avg$session==1)

  #wide
  avgsim2sess <- describeBy(NgR_by_NgW_similarity_fz ~ ROI + subject, data=cert_rsa_beh_dem_w, mat=T)
  avgsim2sess <- avgsim2sess[,c(2,3,6)]
  colnames(avgsim2sess) <- c("ROI","subject","NgR_by_NgW_avgsimilarity_fz")
  cert_rsa_beh_dem_w_avg <- merge(cert_rsa_beh_dem_w,avgsim2sess, by=c("subject", "ROI"), all.x = T)
  avgsim2sess <- describeBy(NtR_by_NtW_similarity_fz ~ ROI + subject, data=cert_rsa_beh_dem_w, mat=T)
  avgsim2sess <- avgsim2sess[,c(2,3,6)]
  colnames(avgsim2sess) <- c("ROI","subject","NtR_by_NtW_avgsimilarity_fz")
  cert_rsa_beh_dem_w_avg <- merge(cert_rsa_beh_dem_w_avg,avgsim2sess, by=c("subject", "ROI"), all.x = T)
  cert_rsa_beh_dem_w_avg <- subset(cert_rsa_beh_dem_w_avg,cert_rsa_beh_dem_w_avg$session==1)

Build primary model (and repeat for each ROI)

Thinking Change on NegReg trials ~ NegR by NegW Similarity + EmotionRating on NegReg trials + covariates(cohort + age + sexAOB)

(Averaged over sessions/runs)

No mixed effects models, models with multiple obs per person (ROIs, sessions, etc.) were not converging…
# Primary model: 
# Dissimilarity between Neg watch and Neg regulate ~  thinking change on Neg regulate + behavioral response rate + age + sex
### Thinking Change on NegR trials ~ ROI sim between Neg watch and Neg regulate + behavioral response rate + age + sex
# also think about adding emotion rating as another variable (it is not significantly corr with thinking change)

R1 <- lm(AvgThink_neg_reg ~ NgR_by_NgW_avgsimilarity_fz + Cohort + AvgEmo_neg_reg + Age + SexAtBirth, data = subset(cert_rsa_beh_dem_w_avg, cert_rsa_beh_dem_w_avg$ROI==1))
R2<- lm(AvgThink_neg_reg ~ NgR_by_NgW_avgsimilarity_fz + Cohort + AvgEmo_neg_reg + Age + SexAtBirth, data = subset(cert_rsa_beh_dem_w_avg, cert_rsa_beh_dem_w_avg$ROI==2))
R3<- lm(AvgThink_neg_reg ~ NgR_by_NgW_avgsimilarity_fz + Cohort + AvgEmo_neg_reg + Age + SexAtBirth, data = subset(cert_rsa_beh_dem_w_avg, cert_rsa_beh_dem_w_avg$ROI==3))
R4<- lm(AvgThink_neg_reg ~ NgR_by_NgW_avgsimilarity_fz + Cohort + AvgEmo_neg_reg + Age + SexAtBirth, data = subset(cert_rsa_beh_dem_w_avg, cert_rsa_beh_dem_w_avg$ROI==4))
R5<- lm(AvgThink_neg_reg ~ NgR_by_NgW_avgsimilarity_fz + Cohort + AvgEmo_neg_reg + Age + SexAtBirth, data = subset(cert_rsa_beh_dem_w_avg, cert_rsa_beh_dem_w_avg$ROI==5))
R6<- lm(AvgThink_neg_reg ~ NgR_by_NgW_avgsimilarity_fz + Cohort + AvgEmo_neg_reg + Age + SexAtBirth, data = subset(cert_rsa_beh_dem_w_avg, cert_rsa_beh_dem_w_avg$ROI==6))
R7<- lm(AvgThink_neg_reg ~ NgR_by_NgW_avgsimilarity_fz + Cohort + AvgEmo_neg_reg + Age + SexAtBirth, data = subset(cert_rsa_beh_dem_w_avg, cert_rsa_beh_dem_w_avg$ROI==7))
R8<- lm(AvgThink_neg_reg ~ NgR_by_NgW_avgsimilarity_fz + Cohort + AvgEmo_neg_reg + Age + SexAtBirth, data = subset(cert_rsa_beh_dem_w_avg, cert_rsa_beh_dem_w_avg$ROI==8))
R9<- lm(AvgThink_neg_reg ~ NgR_by_NgW_avgsimilarity_fz + Cohort + AvgEmo_neg_reg + Age + SexAtBirth, data = subset(cert_rsa_beh_dem_w_avg, cert_rsa_beh_dem_w_avg$ROI==9))
R10<- lm(AvgThink_neg_reg ~ NgR_by_NgW_avgsimilarity_fz + Cohort + AvgEmo_neg_reg + Age + SexAtBirth, data = subset(cert_rsa_beh_dem_w_avg, cert_rsa_beh_dem_w_avg$ROI==10))
R11<- lm(AvgThink_neg_reg ~ NgR_by_NgW_avgsimilarity_fz + Cohort + AvgEmo_neg_reg + Age + SexAtBirth, data = subset(cert_rsa_beh_dem_w_avg, cert_rsa_beh_dem_w_avg$ROI==11))

#summary(R1)
gtsummary::tbl_regression(R1)
Characteristic Beta 95% CI1 p-value
NgR_by_NgW_avgsimilarity_fz 0.15 -0.16, 0.47 0.3
Cohort


    MDD — —
    HC 0.04 -0.29, 0.38 0.8
AvgEmo_neg_reg -0.10 -0.27, 0.08 0.3
Age 0.03 -0.01, 0.08 0.2
SexAtBirth 0.24 -0.06, 0.55 0.12
1 CI = Confidence Interval 
gtsummary::tbl_regression(R2)
Characteristic Beta 95% CI1 p-value
NgR_by_NgW_avgsimilarity_fz 0.16 -0.12, 0.43 0.3
Cohort


    MDD — —
    HC 0.08 -0.24, 0.40 0.6
AvgEmo_neg_reg -0.09 -0.27, 0.08 0.3
Age 0.04 -0.01, 0.08 0.2
SexAtBirth 0.25 -0.05, 0.56 0.10
1 CI = Confidence Interval 
gtsummary::tbl_regression(R3)
Characteristic Beta 95% CI1 p-value
NgR_by_NgW_avgsimilarity_fz -0.13 -0.40, 0.14 0.4
Cohort


    MDD — —
    HC 0.10 -0.22, 0.42 0.5
AvgEmo_neg_reg -0.12 -0.30, 0.06 0.2
Age 0.03 -0.01, 0.08 0.2
SexAtBirth 0.26 -0.04, 0.56 0.091
1 CI = Confidence Interval 
gtsummary::tbl_regression(R4)
Characteristic Beta 95% CI1 p-value
NgR_by_NgW_avgsimilarity_fz 0.14 -0.24, 0.52 0.5
Cohort


    MDD — —
    HC 0.10 -0.22, 0.42 0.5
AvgEmo_neg_reg -0.12 -0.29, 0.06 0.2
Age 0.03 -0.02, 0.08 0.2
SexAtBirth 0.28 -0.03, 0.59 0.071
1 CI = Confidence Interval 
gtsummary::tbl_regression(R5)
Characteristic Beta 95% CI1 p-value
NgR_by_NgW_avgsimilarity_fz -0.17 -0.49, 0.16 0.3
Cohort


    MDD — —
    HC 0.14 -0.19, 0.48 0.4
AvgEmo_neg_reg -0.11 -0.29, 0.06 0.2
Age 0.03 -0.02, 0.08 0.2
SexAtBirth 0.27 -0.03, 0.57 0.078
1 CI = Confidence Interval 
gtsummary::tbl_regression(R6)
Characteristic Beta 95% CI1 p-value
NgR_by_NgW_avgsimilarity_fz -0.02 -0.31, 0.26 0.9
Cohort


    MDD — —
    HC 0.10 -0.23, 0.43 0.6
AvgEmo_neg_reg -0.10 -0.28, 0.07 0.2
Age 0.04 -0.01, 0.09 0.15
SexAtBirth 0.27 -0.04, 0.57 0.087
1 CI = Confidence Interval 
gtsummary::tbl_regression(R7)
Characteristic Beta 95% CI1 p-value
NgR_by_NgW_avgsimilarity_fz 0.09 -0.19, 0.37 0.5
Cohort


    MDD — —
    HC 0.10 -0.22, 0.43 0.5
AvgEmo_neg_reg -0.09 -0.27, 0.08 0.3
Age 0.03 -0.02, 0.08 0.2
SexAtBirth 0.27 -0.03, 0.58 0.079
1 CI = Confidence Interval 
gtsummary::tbl_regression(R8)
Characteristic Beta 95% CI1 p-value
NgR_by_NgW_avgsimilarity_fz -0.17 -0.48, 0.15 0.3
Cohort


    MDD — —
    HC 0.11 -0.22, 0.43 0.5
AvgEmo_neg_reg -0.08 -0.26, 0.10 0.4
Age 0.03 -0.02, 0.08 0.2
SexAtBirth 0.25 -0.05, 0.56 0.10
1 CI = Confidence Interval 
gtsummary::tbl_regression(R9)
Characteristic Beta 95% CI1 p-value
NgR_by_NgW_avgsimilarity_fz 0.10 -0.17, 0.37 0.5
Cohort


    MDD — —
    HC 0.05 -0.29, 0.39 0.8
AvgEmo_neg_reg -0.10 -0.27, 0.07 0.2
Age 0.04 -0.01, 0.08 0.15
SexAtBirth 0.26 -0.04, 0.57 0.088
1 CI = Confidence Interval 
gtsummary::tbl_regression(R10)
Characteristic Beta 95% CI1 p-value
NgR_by_NgW_avgsimilarity_fz 0.15 -0.08, 0.37 0.2
Cohort


    MDD — —
    HC 0.08 -0.23, 0.40 0.6
AvgEmo_neg_reg -0.10 -0.28, 0.07 0.2
Age 0.03 -0.01, 0.08 0.2
SexAtBirth 0.25 -0.05, 0.55 0.10
1 CI = Confidence Interval 
gtsummary::tbl_regression(R11)
Characteristic Beta 95% CI1 p-value
NgR_by_NgW_avgsimilarity_fz -0.15 -0.38, 0.08 0.2
Cohort


    MDD — —
    HC 0.13 -0.19, 0.46 0.4
AvgEmo_neg_reg -0.11 -0.28, 0.06 0.2
Age 0.03 -0.02, 0.08 0.2
SexAtBirth 0.28 -0.02, 0.58 0.068
1 CI = Confidence Interval

Visually, there are not strong negative relationships as we hypthothesized

sim_thinking <- ggplot(data = subset(cert_rsa_beh_dem_w_avg,cert_rsa_beh_dem_w_avg$ROI<5), aes(x = NgR_by_NgW_avgsimilarity_fz, y = AvgThink_neg_reg)) +
          geom_point(aes(color=subject), alpha = .5, size = 1.5) +
          geom_smooth(aes(group = 1), color = "black",size=1.5, method = "lm", se = F) +
          #scale_color_gradient(low = "orchid4", high = "orchid2") +
         facet_wrap(~ROI, ncol = 2)+
        xlab("Negative Reg vs Watch similarity") +
        ylab("Thinking Change")
sim_thinking <- sim_thinking + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))

sim_thinking
## `geom_smooth()` using formula = 'y ~ x'

sim_thinking <- ggplot(data = subset(cert_rsa_beh_dem_w_avg,cert_rsa_beh_dem_w_avg$ROI>4 & cert_rsa_beh_dem_w_avg$ROI<9), aes(x = NgR_by_NgW_avgsimilarity_fz, y = AvgThink_neg_reg)) +
          geom_point(aes(color=subject), alpha = .5, size = 1.5) +
          geom_smooth(aes(group = 1), color = "black",size=1.5, method = "lm", se = F) +
          #scale_color_gradient(low = "orchid4", high = "orchid2") +
         facet_wrap(~ROI, ncol = 2)+
        xlab("Negative Reg vs Watch similarity") +
        ylab("Thinking Change")
sim_thinking <- sim_thinking + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))

sim_thinking
## `geom_smooth()` using formula = 'y ~ x'

sim_thinking <- ggplot(data = subset(cert_rsa_beh_dem_w_avg,cert_rsa_beh_dem_w_avg$ROI>8), aes(x = NgR_by_NgW_avgsimilarity_fz, y = AvgThink_neg_reg)) +
          geom_point(aes(color=subject), alpha = .5, size = 1.5) +
          geom_smooth(aes(group = 1), color = "black",size=1.5, method = "lm", se = F) +
          #scale_color_gradient(low = "orchid4", high = "orchid2") +
         facet_wrap(~ROI, ncol = 2)+
        xlab("Negative Reg vs Watch similarity") +
        ylab("Thinking Change")
sim_thinking <- sim_thinking + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
panel.background = element_blank(), axis.line = element_line(colour = "black"))

sim_thinking
## `geom_smooth()` using formula = 'y ~ x'

example of convergence issue …

summary(lmer(AvgThink_neg_reg ~ similarity_fz*ROI + Cohort + Age + SexAtBirth + (1|subject), data = cert_rsa_beh_dem_NegRxNegW))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 2.78609 (tol = 0.002, component 1)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, : Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: AvgThink_neg_reg ~ similarity_fz * ROI + Cohort + Age + SexAtBirth +  
##     (1 | subject)
##    Data: cert_rsa_beh_dem_NegRxNegW
## 
## REML criterion at convergence: -48965.3
## 
## Scaled residuals: 
##        Min         1Q     Median         3Q        Max 
## -3.035e-06 -4.212e-07  1.940e-08  6.414e-07  3.683e-06 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  subject  (Intercept) 2.506e-02 1.583e-01
##  Residual             3.679e-13 6.065e-07
## Number of obs: 2024, groups:  subject, 92
## 
## Fixed effects:
##                       Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)         -4.299e-01  1.319e-01  5.487e-02  -3.260    0.834    
## similarity_fz        5.451e-13  7.016e-08  1.535e+01   0.000    1.000    
## ROI2                 7.687e-13  1.184e-07  1.535e+01   0.000    1.000    
## ROI3                 8.539e-13  1.218e-07  1.535e+01   0.000    1.000    
## ROI4                 8.358e-13  1.178e-07  1.535e+01   0.000    1.000    
## ROI5                 8.581e-13  1.179e-07  1.535e+01   0.000    1.000    
## ROI6                 8.562e-13  1.165e-07  1.535e+01   0.000    1.000    
## ROI7                 8.206e-13  1.130e-07  1.535e+01   0.000    1.000    
## ROI8                 8.784e-13  1.114e-07  1.535e+01   0.000    1.000    
## ROI9                 8.117e-13  1.131e-07  1.535e+01   0.000    1.000    
## ROI10                8.073e-13  1.111e-07  1.535e+01   0.000    1.000    
## ROI11                8.862e-13  1.163e-07  1.535e+01   0.000    1.000    
## CohortHC             2.331e-01  3.704e-02  1.495e+01   6.293 1.46e-05 ***
## Age                  1.428e-01  5.652e-03  5.896e-02  25.265    0.731    
## SexAtBirth           4.458e-01  3.529e-02  6.164e+00  12.633 1.23e-05 ***
## similarity_fz:ROI2  -4.700e-13  9.646e-08  1.535e+01   0.000    1.000    
## similarity_fz:ROI3  -5.586e-13  9.634e-08  1.535e+01   0.000    1.000    
## similarity_fz:ROI4  -5.398e-13  1.072e-07  1.535e+01   0.000    1.000    
## similarity_fz:ROI5  -5.665e-13  1.007e-07  1.535e+01   0.000    1.000    
## similarity_fz:ROI6  -5.626e-13  9.444e-08  1.535e+01   0.000    1.000    
## similarity_fz:ROI7  -5.205e-13  9.670e-08  1.535e+01   0.000    1.000    
## similarity_fz:ROI8  -6.055e-13  1.031e-07  1.535e+01   0.000    1.000    
## similarity_fz:ROI9  -5.170e-13  8.728e-08  1.535e+01   0.000    1.000    
## similarity_fz:ROI10 -5.056e-13  9.066e-08  1.535e+01   0.000    1.000    
## similarity_fz:ROI11 -5.904e-13  9.053e-08  1.535e+01   0.000    1.000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 25 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it
## optimizer (nloptwrap) convergence code: 0 (OK)
## Model failed to converge with max|grad| = 2.78609 (tol = 0.002, component 1)
## Model is nearly unidentifiable: very large eigenvalue
##  - Rescale variables?
summary(lmer(AvgThink_neg_reg ~ avgsimilarity_fz*ROI + Cohort + Age + SexAtBirth + (1|subject), data = cert_rsa_beh_dem_NegRxNegW_avg))
## Warning in optwrap(optimizer, devfun, getStart(start, rho$pp), lower =
## rho$lower, : convergence code -4 from nloptwrap: NLOPT_ROUNDOFF_LIMITED:
## Roundoff errors led to a breakdown of the optimization algorithm. In this case,
## the returned minimum may still be useful. (e.g. this error occurs in NEWUOA if
## one tries to achieve a tolerance too close to machine precision.)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: 
## AvgThink_neg_reg ~ avgsimilarity_fz * ROI + Cohort + Age + SexAtBirth +  
##     (1 | subject)
##    Data: cert_rsa_beh_dem_NegRxNegW_avg
## 
## REML criterion at convergence: -19101.1
## 
## Scaled residuals: 
##        Min         1Q     Median         3Q        Max 
## -2.345e-05 -4.524e-06 -4.002e-07  5.236e-06  1.988e-05 
## 
## Random effects:
##  Groups   Name        Variance  Std.Dev. 
##  subject  (Intercept) 4.123e-02 2.031e-01
##  Residual             2.529e-11 5.029e-06
## Number of obs: 1012, groups:  subject, 92
## 
## Fixed effects:
##                          Estimate Std. Error         df t value Pr(>|t|)    
## (Intercept)             2.144e+00  1.713e-01  4.826e-01  12.515   0.1903    
## avgsimilarity_fz        9.713e-12  1.151e-06  8.313e+01   0.000   1.0000    
## ROI2                    6.747e-13  1.766e-06  8.313e+01   0.000   1.0000    
## ROI3                    1.611e-11  1.803e-06  8.313e+01   0.000   1.0000    
## ROI4                    6.558e-12  1.829e-06  8.313e+01   0.000   1.0000    
## ROI5                    1.746e-11  1.771e-06  8.313e+01   0.000   1.0000    
## ROI6                    1.089e-11  1.757e-06  8.313e+01   0.000   1.0000    
## ROI7                    5.066e-12  1.689e-06  8.313e+01   0.000   1.0000    
## ROI8                    1.833e-11  1.654e-06  8.313e+01   0.000   1.0000    
## ROI9                    5.240e-12  1.742e-06  8.313e+01   0.000   1.0000    
## ROI10                   3.775e-12  1.657e-06  8.313e+01   0.000   1.0000    
## ROI11                   1.866e-11  1.714e-06  8.313e+01   0.000   1.0000    
## CohortHC                1.228e-01  4.756e-02  8.469e+01   2.581   0.0116 *  
## Age                     3.499e-02  7.335e-03  4.814e-01   4.771   0.3030    
## SexAtBirth              2.498e-01  4.544e-02  5.450e+01   5.498 1.05e-06 ***
## avgsimilarity_fz:ROI2   1.075e-12  1.541e-06  8.313e+01   0.000   1.0000    
## avgsimilarity_fz:ROI3  -1.472e-11  1.517e-06  8.313e+01   0.000   1.0000    
## avgsimilarity_fz:ROI4  -4.001e-12  1.825e-06  8.313e+01   0.000   1.0000    
## avgsimilarity_fz:ROI5  -1.740e-11  1.624e-06  8.313e+01   0.000   1.0000    
## avgsimilarity_fz:ROI6  -9.720e-12  1.539e-06  8.313e+01   0.000   1.0000    
## avgsimilarity_fz:ROI7  -2.235e-12  1.558e-06  8.313e+01   0.000   1.0000    
## avgsimilarity_fz:ROI8  -2.134e-11  1.639e-06  8.313e+01   0.000   1.0000    
## avgsimilarity_fz:ROI9  -4.069e-12  1.476e-06  8.313e+01   0.000   1.0000    
## avgsimilarity_fz:ROI10 -1.063e-12  1.460e-06  8.313e+01   0.000   1.0000    
## avgsimilarity_fz:ROI11 -1.731e-11  1.427e-06  8.313e+01   0.000   1.0000    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 25 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it