This code reproduces the trial-level analyses with a 2s event duration reported in the following manuscript:

Cosme, D., Mobasser, A., & J. H. Pfeifer. If you’re happy and you know it: Neural correlates of self-evaluated psychological health and well-being

load packages

if(!require('pacman')) {
    install.packages('pacman')
}

pacman::p_load(tidyverse, gtools, GGally, sjstats, lme4, lmerTest, knitr, ggeffects, kableExtra, install = TRUE)

define aesthetics

rois = c("#006989", "#56445D", "#8EC922")
constructs = c("#FEC601", "#254E70", "#F43C13")
instructions = wesanderson::wes_palette("Darjeeling1", 2, "continuous")
valence = c("#119DA4", "#19647E")
plot_aes = theme_minimal() +
  theme(legend.position = "top",
        legend.text = element_text(size = 12),
        text = element_text(size = 16, family = "Futura Medium"),
        axis.text = element_text(color = "black"),
        axis.line = element_line(colour = "black"),
        axis.ticks.y = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        panel.background = element_blank())

define functions

model_table = function(model) {
  model %>%
    broom.mixed::tidy(., conf.int = TRUE) %>%
    filter(effect == "fixed") %>%
    rename("SE" = std.error,
           "z" = statistic,
           "p" = p.value) %>%
    mutate(term = gsub("\\(Intercept\\)", "Intercept (self-oriented well-being)", term),
           term = gsub("log\\(RT\\)", "RT", term),
           term = gsub("constructill-being", "Construct (ill-being)", term),
           term = gsub("constructsocial well-being", "Construct (social well-being)", term),
           term = gsub(":", " x ", term),
           p = ifelse(p < .001, "< .001",
                      ifelse(p == 1, "1.000", gsub("0.(.*)", ".\\1", sprintf("%.3f", p)))),
           `b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", estimate, conf.low, conf.high),
           z = abs(round(z, 2))) %>%
    select(term, `b [95% CI]`, z, p) %>%
    kable() %>%
    kableExtra::kable_styling()
}

load data

task = read.csv("../../data/task_data.csv") %>%
  mutate(responseYN = as.factor(responseYN),
         construct = recode(construct, "well-being" = "self-oriented well-being",
                            "social" = "social well-being"),
         construct = factor(construct, levels = c("self-oriented well-being", "social well-being", "ill-being")),
         valence = factor(valence, levels = c("positive", "negative")))
betas = read.csv("../../data/neuro_data_2s.csv")

exclude outliers and standardize

  • Subset self trials
  • Exclude outlier trials that are > 3 SDs from roi median across participants
trial_conditions = task %>%
  select(subjectID, trial, run, instruction)

betas_outlier = betas %>%
  left_join(., trial_conditions) %>%
  filter(instruction == "self") %>%
  group_by(roi) %>%
  mutate(median = median(meanPE, na.rm = TRUE),
         sd3 = 3*sd(meanPE, na.rm = TRUE),
         outlier = ifelse(meanPE > median + sd3 | meanPE < median - sd3, "yes", "no")) 
## Joining, by = c("subjectID", "trial", "run")
betas_outlier %>%
  group_by(outlier) %>%
  summarize(n = n()) %>%
  ungroup() %>%
  mutate(total = sum(n),
         percent = round((n / total) * 100, 2)) %>%
  kable(format = "pandoc")
outlier n total percent
no 11857 11988 98.91
yes 131 11988 1.09 
betas_ex = betas_outlier %>%
  filter(outlier == "no") %>%
  select(-c(median, sd3, outlier))

merge data and exclude participants

  • Motion exclusions: FP091
  • Technical failure: FP080, FP082
  • Non-compliance: FP021, FP049, FP085, FP121
  • Standardize within participant and ROI
data_ex = task %>%
  left_join(., betas_ex) %>%
  filter(!subjectID %in% c("FP091", "FP080", "FP082", "FP021", "FP049", "FP085", "FP121")) %>%
  filter(!is.na(meanPE)) %>%
  filter(instruction == "self") %>%
  group_by(roi, subjectID) %>%
  mutate(std = meanPE / sd(meanPE, na.rm = TRUE)) %>%
  select(-sdPE)
## Joining, by = c("trial", "subjectID", "run", "instruction")

remove missing data for MLM analysis

data_complete_mod = data_ex %>%
  select(-meanPE) %>%
  spread(roi, std) %>%
  select(subjectID, trial, run, construct, item, valence, response, responseYN, RT_original, RT, pgACC, vmPFC, VS) %>%
  na.omit()

ROI descriptives

means, SDs, and correlations

# correlations
pgACC_vmPFC = data_complete_mod %>%
   mutate(subjectID = as.factor(subjectID)) %>%
   select(subjectID, pgACC, vmPFC) %>%
   rmcorr::rmcorr(subjectID, pgACC, vmPFC, .)

pgACC_VS = data_complete_mod %>%
   mutate(subjectID = as.factor(subjectID)) %>%
   select(subjectID, pgACC, VS) %>%
   rmcorr::rmcorr(subjectID, pgACC, VS, .)

vmPFC_VS = data_complete_mod %>%
   mutate(subjectID = as.factor(subjectID)) %>%
   select(subjectID, vmPFC, VS) %>%
   rmcorr::rmcorr(subjectID, vmPFC, VS, .)

# means and SDs
mean_table = data_complete_mod %>%
  gather(roi, value, pgACC, vmPFC, VS) %>%
  group_by(roi) %>%
  summarize(M = round(mean(value, na.rm = TRUE), 2),
            SD = round(sd(value, na.rm = TRUE), 2)) %>%
  mutate(order = ifelse(roi == "pgACC", 1,
                 ifelse(roi == "vmPFC", 2, 3)),
         roi = ifelse(roi == "pgACC", "1. pgACC",
               ifelse(roi == "vmPFC", "2. vmPFC", "3. VS"))) %>%
  rename("ROI" = roi) %>%
  arrange(order) %>%
  select(-order)

# table
corr_table = data.frame(ROI = c("1. pgACC", "2. vmPFC", "3. VS"),
           `1` = c("--", paste0(round(pgACC_vmPFC[[1]],2), " [",round(pgACC_vmPFC[[4]][1],2),", ",round(pgACC_vmPFC[[4]][2],2),"]"),
                   paste0(round(pgACC_VS[[1]],2), " [",round(pgACC_VS[[4]][1],2),", ",round(pgACC_VS[[4]][2],2),"]")),
           `2` = c("", "--", paste0(round(vmPFC_VS[[1]],2), " [",round(vmPFC_VS[[4]][1],2),", ",round(vmPFC_VS[[4]][2],2),"]")),
           `3` = c("", "", "--"),
           check.names = FALSE)

mean_table %>%
  full_join(., corr_table) %>%
  knitr::kable(format = "pandoc", caption = "Repeated Measures Correlations Between ROIs (n = 3694). All correlations are statistically significant, p < .001. 95% confidence intervals are bracketed")
## Joining, by = "ROI"
Repeated Measures Correlations Between ROIs (n = 3694). All correlations are statistically significant, p < .001. 95% confidence intervals are bracketed
ROI M SD 1 2 3
1. pgACC 0.05 1.03 –
2. vmPFC 0.18 1.04 0.82 [0.81, 0.83] –
3. VS 0.19 1.03 0.41 [0.38, 0.44] 0.44 [0.41, 0.46] –

RT models

These are the primary analyses reported in the main text of the manuscript, controlling for reaction time.

run models

compare models

Determine the best fitting model.

model_0_rt = Base model including including no ROIs

model_1_rt = Main effects model including ROIs

model_2_rt = Interaction model including ROIs and their interactions with well-being construct

model_3_rt = Interaction model including ROIs and their interactions with valence rather than well-being construct

model_0_rt = glmer(responseYN ~ construct * RT + 
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model_1_rt = glmer(responseYN ~ construct * RT + pgACC + vmPFC + VS + 
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model_2_rt = glmer(responseYN ~ RT*construct + pgACC*construct + vmPFC*construct + VS*construct +
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model_3_rt = glmer(responseYN ~ RT*valence + pgACC*valence + vmPFC*valence + VS*valence +
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

# compare models
anova(model_0_rt, model_1_rt, model_2_rt, model_3_rt)

calculate R2

Compare the model R2 for the best fitting compared to the base model.

performance::r2(model_0_rt)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.516
##      Marginal R2: 0.484
performance::r2(model_2_rt)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.534
##      Marginal R2: 0.501

check variance inflation in the best fitting model

# vif
car::vif(model_2_rt)
##                      GVIF Df GVIF^(1/(2*Df))
## RT               3.990068  1        1.997516
## construct        1.499776  2        1.106641
## pgACC           11.926056  1        3.453412
## vmPFC           12.014489  1        3.466192
## VS               4.756413  1        2.180920
## RT:construct     4.126661  2        1.425278
## construct:pgACC 30.537367  2        2.350758
## construct:vmPFC 37.175639  2        2.469247
## construct:VS     6.386754  2        1.589718

model summary

summary(model_2_rt)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: responseYN ~ RT * construct + pgACC * construct + vmPFC * construct +  
##     VS * construct + (1 | subjectID)
##    Data: data_complete_mod
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   3042.3   3141.7  -1505.2   3010.3     3660 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -8.4412 -0.5133  0.1994  0.3905  2.8918 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 0.2338   0.4835  
## Number of obs: 3676, groups:  subjectID, 105
## 
## Fixed effects:
##                                  Estimate Std. Error z value
## (Intercept)                       1.91712    0.11290  16.981
## RT                               -3.20334    0.30409 -10.534
## constructsocial well-being        0.07527    0.14566   0.517
## constructill-being               -2.95940    0.12604 -23.481
## pgACC                             0.36299    0.16519   2.197
## vmPFC                            -0.45453    0.16394  -2.773
## VS                                0.13984    0.10197   1.371
## RT:constructsocial well-being     0.05326    0.40291   0.132
## RT:constructill-being             3.20795    0.36018   8.907
## constructsocial well-being:pgACC  0.00910    0.21691   0.042
## constructill-being:pgACC         -0.20030    0.19963  -1.003
## constructsocial well-being:vmPFC -0.10296    0.21833  -0.472
## constructill-being:vmPFC          0.71571    0.19862   3.603
## constructsocial well-being:VS    -0.05876    0.13957  -0.421
## constructill-being:VS            -0.29307    0.12645  -2.318
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## RT                               < 0.0000000000000002 ***
## constructsocial well-being                   0.605340    
## constructill-being               < 0.0000000000000002 ***
## pgACC                                        0.027992 *  
## vmPFC                                        0.005562 ** 
## VS                                           0.170227    
## RT:constructsocial well-being                0.894829    
## RT:constructill-being            < 0.0000000000000002 ***
## constructsocial well-being:pgACC             0.966535    
## constructill-being:pgACC                     0.315707    
## constructsocial well-being:vmPFC             0.637242    
## constructill-being:vmPFC                     0.000314 ***
## constructsocial well-being:VS                0.673751    
## constructill-being:VS                        0.020465 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation matrix not shown by default, as p = 15 > 12.
## Use print(x, correlation=TRUE)  or
##     vcov(x)        if you need it

model table

model_table(model_2_rt)
term b [95% CI] z p
Intercept (self-oriented well-being) 1.92 [1.70, 2.14] 16.98 < .001
RT -3.20 [-3.80, -2.61] 10.53 < .001
Construct (social well-being) 0.08 [-0.21, 0.36] 0.52 .605
Construct (ill-being) -2.96 [-3.21, -2.71] 23.48 < .001
pgACC 0.36 [0.04, 0.69] 2.20 .028
vmPFC -0.45 [-0.78, -0.13] 2.77 .006
VS 0.14 [-0.06, 0.34] 1.37 .170
RT x Construct (social well-being) 0.05 [-0.74, 0.84] 0.13 .895
RT x Construct (ill-being) 3.21 [2.50, 3.91] 8.91 < .001
Construct (social well-being) x pgACC 0.01 [-0.42, 0.43] 0.04 .967
Construct (ill-being) x pgACC -0.20 [-0.59, 0.19] 1.00 .316
Construct (social well-being) x vmPFC -0.10 [-0.53, 0.32] 0.47 .637
Construct (ill-being) x vmPFC 0.72 [0.33, 1.10] 3.60 < .001
Construct (social well-being) x VS -0.06 [-0.33, 0.21] 0.42 .674
Construct (ill-being) x VS -0.29 [-0.54, -0.05] 2.32 .020

marginal probabilities

all

ggeffect(model_2_rt)
## $RT
## # Predicted probabilities of responseYN
## 
##    RT | Predicted |       95% CI
## --------------------------------
## -1.50 |      0.98 | [0.97, 0.99]
## -1.00 |      0.96 | [0.94, 0.97]
## -0.50 |      0.88 | [0.86, 0.90]
##  0.00 |      0.72 | [0.69, 0.75]
##  0.50 |      0.47 | [0.42, 0.52]
##  1.00 |      0.23 | [0.18, 0.30]
## 
## $construct
## # Predicted probabilities of responseYN
## 
## construct                | Predicted |       95% CI
## ---------------------------------------------------
## self-oriented well-being |      0.90 | [0.88, 0.92]
## social well-being        |      0.91 | [0.89, 0.93]
## ill-being                |      0.27 | [0.23, 0.30]
## 
## $pgACC
## # Predicted probabilities of responseYN
## 
## pgACC | Predicted |       95% CI
## --------------------------------
##    -6 |      0.35 | [0.16, 0.59]
##    -4 |      0.49 | [0.33, 0.65]
##    -2 |      0.64 | [0.55, 0.72]
##     0 |      0.76 | [0.73, 0.79]
##     2 |      0.85 | [0.80, 0.89]
##     4 |      0.91 | [0.84, 0.95]
##     6 |      0.95 | [0.88, 0.98]
## 
## $vmPFC
## # Predicted probabilities of responseYN
## 
## vmPFC | Predicted |       95% CI
## --------------------------------
##    -6 |      0.94 | [0.84, 0.98]
##    -4 |      0.90 | [0.82, 0.95]
##    -2 |      0.85 | [0.79, 0.89]
##     0 |      0.77 | [0.74, 0.80]
##     2 |      0.67 | [0.60, 0.74]
##     4 |      0.55 | [0.40, 0.70]
## 
## $VS
## # Predicted probabilities of responseYN
## 
## VS | Predicted |       95% CI
## -----------------------------
## -6 |      0.74 | [0.59, 0.85]
## -4 |      0.75 | [0.65, 0.82]
## -2 |      0.76 | [0.70, 0.80]
##  0 |      0.76 | [0.74, 0.79]
##  2 |      0.77 | [0.73, 0.81]
##  4 |      0.78 | [0.70, 0.85]
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "model_2_rt"

pgACC

ggeffect(model_2_rt, terms = c("construct", "pgACC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.9026609 0.1193651 0.8800808 0.9213692 0 0.90 [0.88, 0.92]
social well-being 0.9060798 0.1223423 0.8835922 0.9245941 0 0.91 [0.88, 0.92]
ill-being 0.2639680 0.0841815 0.2331805 0.2972447 0 0.26 [0.23, 0.30]
self-oriented well-being 0.9302244 0.2062779 0.8989711 0.9523222 1 0.93 [0.90, 0.95]
social well-being 0.9333159 0.1930795 0.9055394 0.9533454 1 0.93 [0.91, 0.95]
ill-being 0.2967656 0.1354690 0.2444819 0.3549769 1 0.30 [0.24, 0.35]

vmPFC

ggeffect(model_2_rt, terms = c("construct", "vmPFC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.9111437 0.1264547 0.8889276 0.9292698 0 0.91 [0.89, 0.93]
social well-being 0.9157621 0.1322902 0.8934835 0.9337267 0 0.92 [0.89, 0.93]
ill-being 0.2566346 0.0864006 0.2256784 0.2902454 0 0.26 [0.23, 0.29]
self-oriented well-being 0.8668230 0.1704238 0.8233384 0.9008924 1 0.87 [0.82, 0.90]
social well-being 0.8615985 0.1487329 0.8230456 0.8928455 1 0.86 [0.82, 0.89]
ill-being 0.3095231 0.1265762 0.2591416 0.3648754 1 0.31 [0.26, 0.36]

VS

ggeffect(model_2_rt, terms = c("construct", "VS [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.9019953 0.1200613 0.8791365 0.9209197 0 0.90 [0.88, 0.92]
social well-being 0.9064356 0.1223074 0.8840293 0.9248808 0 0.91 [0.88, 0.92]
ill-being 0.2714135 0.0849806 0.2397552 0.3055717 0 0.27 [0.24, 0.31]
self-oriented well-being 0.9136814 0.1504340 0.8874138 0.9342744 1 0.91 [0.89, 0.93]
social well-being 0.9130895 0.1531348 0.8861310 0.9341399 1 0.91 [0.89, 0.93]
ill-being 0.2421944 0.1042255 0.2066948 0.2816264 1 0.24 [0.21, 0.28]

plot predicted effects

by construct

raw_values = data_complete_mod %>%
  select(subjectID, construct, RT, pgACC, vmPFC, VS) %>%
 gather(roi, x, pgACC, VS, vmPFC)

vals = seq(round(min(raw_values$x),1), round(max(raw_values$x), 1), .5)

predicted_sub = ggpredict(model_2_rt, terms = c("pgACC [vals]", "construct", "subjectID"), type = "random") %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggpredict(model_2_rt, terms = c("VS [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggpredict(model_2_rt, terms = c("vmPFC [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted = ggeffect(model_2_rt, terms = c("pgACC [vals]", "construct")) %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggeffect(model_2_rt, terms = c("VS [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggeffect(model_2_rt, terms = c("vmPFC [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .5) +
  geom_line(aes(color = group), size = 2) +
  facet_grid(~roi) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

by ROI

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = roi, group = interaction(facet, roi)), size = .25, alpha = .5) +
  geom_line(aes(color = roi), size = 2) +
  facet_grid(~group) +
  scale_color_manual(values = rois, name = "") + 
  scale_fill_manual(values = rois, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

RT

vals = seq(round(min(data_complete_mod$RT),1), round(max(data_complete_mod$RT), 1), .1)

predicted_sub = ggpredict(model_2_rt, terms = c("RT [vals]", "construct", "subjectID"), type = "random") %>%
  data.frame()

ggeffect(model_2_rt, terms = c("RT [vals]", "construct")) %>%
  data.frame() %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .5) +
  geom_line(aes(color = group), size = 1.5) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  labs(x = "\nlog-transformed grand-mean centered RT (s)", y = "probability of responding yes\n") + 
  plot_aes

no RT model

These are the original preregistered analyses not controlling for reaction time, reported in supplementary material.

run model

model_2 = glmer(responseYN ~ construct + pgACC*construct + vmPFC*construct + VS*construct +
                     (1 | subjectID), 
                family = binomial, 
                data = data_complete_mod, control = glmerControl(optimizer = "bobyqa"))

model summary

summary(model_2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: responseYN ~ construct + pgACC * construct + vmPFC * construct +  
##     VS * construct + (1 | subjectID)
##    Data: data_complete_mod
## Control: glmerControl(optimizer = "bobyqa")
## 
##      AIC      BIC   logLik deviance df.resid 
##   3313.9   3394.6  -1644.0   3287.9     3663 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -4.0888 -0.5242  0.3265  0.4124  2.9837 
## 
## Random effects:
##  Groups    Name        Variance Std.Dev.
##  subjectID (Intercept) 0.192    0.4382  
## Number of obs: 3676, groups:  subjectID, 105
## 
## Fixed effects:
##                                   Estimate Std. Error z value
## (Intercept)                       2.061415   0.103539  19.909
## constructsocial well-being       -0.067640   0.132185  -0.512
## constructill-being               -3.096894   0.118902 -26.046
## pgACC                             0.364552   0.154025   2.367
## vmPFC                            -0.401804   0.154109  -2.607
## VS                                0.106716   0.093745   1.138
## constructsocial well-being:pgACC  0.001201   0.203860   0.006
## constructill-being:pgACC         -0.207191   0.190181  -1.089
## constructsocial well-being:vmPFC -0.118854   0.206051  -0.577
## constructill-being:vmPFC          0.662684   0.190486   3.479
## constructsocial well-being:VS    -0.049635   0.127962  -0.388
## constructill-being:VS            -0.260937   0.119993  -2.175
##                                              Pr(>|z|)    
## (Intercept)                      < 0.0000000000000002 ***
## constructsocial well-being                   0.608855    
## constructill-being               < 0.0000000000000002 ***
## pgACC                                        0.017941 *  
## vmPFC                                        0.009127 ** 
## VS                                           0.254968    
## constructsocial well-being:pgACC             0.995300    
## constructill-being:pgACC                     0.275958    
## constructsocial well-being:vmPFC             0.564062    
## constructill-being:vmPFC                     0.000503 ***
## constructsocial well-being:VS                0.698097    
## constructill-being:VS                        0.029660 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) cnwll- cnstr- pgACC  vmPFC  VS     cw-:AC c-:ACC cw-:PF
## cnstrctwll- -0.622                                                        
## cnstrctll-b -0.736  0.542                                                 
## pgACC        0.258 -0.196 -0.225                                          
## vmPFC       -0.270  0.205  0.235 -0.795                                   
## VS          -0.058  0.041  0.050 -0.084 -0.191                            
## cwll-bn:ACC -0.193  0.282  0.169 -0.743  0.590  0.062                     
## cnstrc-:ACC -0.201  0.160  0.216 -0.797  0.632  0.068  0.600              
## cwll-bn:PFC  0.198 -0.371 -0.173  0.583 -0.733  0.141 -0.771 -0.471       
## cnstrc-:PFC  0.217 -0.167 -0.256  0.632 -0.796  0.152 -0.476 -0.782  0.593
## cwll-bng:VS  0.038 -0.029 -0.032  0.061  0.136 -0.720 -0.079 -0.051 -0.220
## cnstrct-:VS  0.041 -0.032 -0.087  0.066  0.146 -0.772 -0.050 -0.098 -0.110
##             c-:PFC cw-:VS
## cnstrctwll-              
## cnstrctll-b              
## pgACC                    
## vmPFC                    
## VS                       
## cwll-bn:ACC              
## cnstrc-:ACC              
## cwll-bn:PFC              
## cnstrc-:PFC              
## cwll-bng:VS -0.110       
## cnstrct-:VS -0.191  0.564

model table

model_table(model_2)
term b [95% CI] z p
Intercept (self-oriented well-being) 2.06 [1.86, 2.26] 19.91 < .001
Construct (social well-being) -0.07 [-0.33, 0.19] 0.51 .609
Construct (ill-being) -3.10 [-3.33, -2.86] 26.05 < .001
pgACC 0.36 [0.06, 0.67] 2.37 .018
vmPFC -0.40 [-0.70, -0.10] 2.61 .009
VS 0.11 [-0.08, 0.29] 1.14 .255
Construct (social well-being) x pgACC 0.00 [-0.40, 0.40] 0.01 .995
Construct (ill-being) x pgACC -0.21 [-0.58, 0.17] 1.09 .276
Construct (social well-being) x vmPFC -0.12 [-0.52, 0.28] 0.58 .564
Construct (ill-being) x vmPFC 0.66 [0.29, 1.04] 3.48 < .001
Construct (social well-being) x VS -0.05 [-0.30, 0.20] 0.39 .698
Construct (ill-being) x VS -0.26 [-0.50, -0.03] 2.17 .030

marginal probabilities

all

ggeffect(model_2)
## $construct
## # Predicted probabilities of responseYN
## 
## construct                | Predicted |       95% CI
## ---------------------------------------------------
## self-oriented well-being |      0.88 | [0.86, 0.90]
## social well-being        |      0.87 | [0.85, 0.89]
## ill-being                |      0.27 | [0.24, 0.30]
## 
## $pgACC
## # Predicted probabilities of responseYN
## 
## pgACC | Predicted |       95% CI
## --------------------------------
##    -6 |      0.31 | [0.15, 0.54]
##    -4 |      0.45 | [0.30, 0.61]
##    -2 |      0.59 | [0.51, 0.67]
##     0 |      0.73 | [0.70, 0.75]
##     2 |      0.83 | [0.77, 0.87]
##     4 |      0.90 | [0.82, 0.94]
##     6 |      0.94 | [0.86, 0.98]
## 
## $vmPFC
## # Predicted probabilities of responseYN
## 
## vmPFC | Predicted |       95% CI
## --------------------------------
##    -6 |      0.91 | [0.80, 0.97]
##    -4 |      0.87 | [0.77, 0.93]
##    -2 |      0.81 | [0.75, 0.86]
##     0 |      0.74 | [0.71, 0.76]
##     2 |      0.64 | [0.57, 0.71]
##     4 |      0.54 | [0.39, 0.68]
## 
## $VS
## # Predicted probabilities of responseYN
## 
## VS | Predicted |       95% CI
## -----------------------------
## -6 |      0.72 | [0.58, 0.83]
## -4 |      0.73 | [0.63, 0.80]
## -2 |      0.73 | [0.68, 0.77]
##  0 |      0.73 | [0.70, 0.75]
##  2 |      0.73 | [0.68, 0.77]
##  4 |      0.73 | [0.65, 0.80]
## 
## attr(,"class")
## [1] "ggalleffects" "list"        
## attr(,"model.name")
## [1] "model_2"

pgACC

ggeffect(model_2, terms = c("construct", "pgACC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.8818430 0.0992888 0.8600129 0.9006620 0 0.88 [0.86, 0.90]
social well-being 0.8711979 0.0984572 0.8479510 0.8913461 0 0.87 [0.85, 0.89]
ill-being 0.2653450 0.0801911 0.2358547 0.2970890 0 0.27 [0.24, 0.30]
self-oriented well-being 0.9148660 0.1860015 0.8818425 0.9392952 1 0.91 [0.88, 0.94]
social well-being 0.9069833 0.1737517 0.8739999 0.9320042 1 0.91 [0.87, 0.93]
ill-being 0.2971284 0.1322284 0.2459793 0.3539213 1 0.30 [0.25, 0.35]

vmPFC

ggeffect(model_2, terms = c("construct", "vmPFC [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.8910185 0.1063336 0.8690718 0.9096687 0 0.89 [0.87, 0.91]
social well-being 0.8833043 0.1088702 0.8594487 0.9035652 0 0.88 [0.86, 0.90]
ill-being 0.2579430 0.0823574 0.2282701 0.2900235 0 0.26 [0.23, 0.29]
self-oriented well-being 0.8454541 0.1529144 0.8021325 0.8807015 1 0.85 [0.80, 0.88]
social well-being 0.8180874 0.1294483 0.7772529 0.8528532 1 0.82 [0.78, 0.85]
ill-being 0.3109230 0.1237734 0.2614590 0.3651185 1 0.31 [0.26, 0.37]

VS

ggeffect(model_2, terms = c("construct", "VS [0, 1]")) %>%
  data.frame() %>%
  mutate(`b [95% CI]` = sprintf("%.2f [%0.2f, %.2f]", predicted, conf.low, conf.high)) %>%
  kable() %>%
  kableExtra::kable_styling()
x predicted std.error conf.low conf.high group b [95% CI]
self-oriented well-being 0.8817259 0.1004209 0.8596101 0.9007600 0 0.88 [0.86, 0.90]
social well-being 0.8721433 0.0990453 0.8488896 0.8922727 0 0.87 [0.85, 0.89]
ill-being 0.2727966 0.0809947 0.2424624 0.3053960 0 0.27 [0.24, 0.31]
self-oriented well-being 0.8924095 0.1289857 0.8656213 0.9143857 1 0.89 [0.87, 0.91]
social well-being 0.8783743 0.1279736 0.8489387 0.9027314 1 0.88 [0.85, 0.90]
ill-being 0.2432948 0.1008792 0.2087594 0.2815107 1 0.24 [0.21, 0.28]

plot predicted effects

by construct

raw_values = data_complete_mod %>%
  select(subjectID, construct, RT, pgACC, vmPFC, VS) %>%
 gather(roi, x, pgACC, VS, vmPFC)

vals = seq(round(min(raw_values$x),1), round(max(raw_values$x), 1), .5)

predicted_sub = ggpredict(model_2, terms = c("pgACC [vals]", "construct", "subjectID"), type = "random") %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggpredict(model_2, terms = c("VS [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggpredict(model_2, terms = c("vmPFC [vals]", "construct", "subjectID"), type = "random") %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted = ggeffect(model_2, terms = c("pgACC [vals]", "construct")) %>%
  data.frame() %>%
  mutate(roi = "pgACC") %>%
  bind_rows(ggeffect(model_2, terms = c("VS [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "VS")) %>%
  bind_rows(ggeffect(model_2, terms = c("vmPFC [vals]", "construct")) %>%
    data.frame() %>%
    mutate(roi = "vmPFC"))

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = group, group = interaction(facet, group)), size = .25, alpha = .5) +
  geom_line(aes(color = group), size = 2) +
  facet_grid(~roi) +
  scale_color_manual(values = constructs, name = "") + 
  scale_fill_manual(values = constructs, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes

by ROI

predicted %>%
  ggplot(aes(x, predicted)) +
  geom_line(data = predicted_sub, aes(color = roi, group = interaction(facet, roi)), size = .25, alpha = .5) +
  geom_line(aes(color = roi), size = 2) +
  facet_grid(~group) +
  scale_color_manual(values = rois, name = "") + 
  scale_fill_manual(values = rois, name = "") + 
  scale_x_continuous(breaks = seq(-6, 6, 2)) +
  labs(x = "\nparameter estimate (SD)", y = "probability of responding yes\n") + 
  plot_aes