# load packages
library(afex)
Loading required package: lme4
Loading required package: Matrix
************
Welcome to afex. For support visit: http://afex.singmann.science/
- Functions for ANOVAs: aov_car(), aov_ez(), and aov_4()
- Methods for calculating p-values with mixed(): 'KR', 'S', 'LRT', and 'PB'
- 'afex_aov' and 'mixed' objects can be passed to emmeans() for follow-up tests
- NEWS: library('emmeans') now needs to be called explicitly!
- Get and set global package options with: afex_options()
- Set orthogonal sum-to-zero contrasts globally: set_sum_contrasts()
- For example analyses see: browseVignettes("afex")
************
Attaching package: ‘afex’
The following object is masked from ‘package:lme4’:
lmer
library(emmeans)
library(tidyverse)
Registered S3 method overwritten by 'rvest':
method from
read_xml.response xml2
[30m── [1mAttaching packages[22m ──────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──[39m
[30m[32m✔[30m [34mggplot2[30m 3.2.0 [32m✔[30m [34mpurrr [30m 0.3.2
[32m✔[30m [34mtibble [30m 2.1.3 [32m✔[30m [34mdplyr [30m 0.8.1
[32m✔[30m [34mtidyr [30m 0.8.3 [32m✔[30m [34mstringr[30m 1.4.0
[32m✔[30m [34mreadr [30m 1.3.1 [32m✔[30m [34mforcats[30m 0.4.0[39m
[30m── [1mConflicts[22m ─────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
[31m✖[30m [34mtidyr[30m::[32mexpand()[30m masks [34mMatrix[30m::expand()
[31m✖[30m [34m.GlobalEnv[30m::[32mfilter()[30m masks [34mdplyr[30m::filter(), [34mstats[30m::filter()
[31m✖[30m [34mdplyr[30m::[32mlag()[30m masks [34mstats[30m::lag()[39m
library(viridis)
Loading required package: viridisLite
library(psych)
Attaching package: ‘psych’
The following objects are masked from ‘package:ggplot2’:
%+%, alpha
library(ggbeeswarm)
library(nlme)
Attaching package: ‘nlme’
The following object is masked from ‘package:dplyr’:
collapse
The following object is masked from ‘package:lme4’:
lmList
library(jtools)
afex_options(emmeans_model = "multivariate")
filter <- dplyr::filter
select <- dplyr::select
# load data
# NOTE that this is the dataset that calculates bias by subtracting the median pull reaction time in each category from the median push reaction time in that category (NOT the trial-by-trial bias subtracting run1 from run2)
data_bias <- readRDS("~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/bias_long.rds")
varstoadd <- readRDS("~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/varstoadd.rds")
# mean-center tot_icg and tot_bdi for use later
meancent <- function(x) { x - mean(x, na.rm=TRUE) } #simple worker function to mean center a variable
varstoadd <- varstoadd %>% mutate_at(vars("tot_icg", "tot_bdi"), funs(cm=meancent))
funs() is soft deprecated as of dplyr 0.8.0
please use list() instead
# Before:
funs(name = f(.))
# After:
list(name = ~ f(.))
[90mThis warning is displayed once per session.[39m
data_bias <- left_join(data_bias, varstoadd, by="ID") %>%
rename(treatment = cond,
gAAT_bias = bias) %>%
mutate(treatment = recode(treatment,
A = "placebo",
B = "oxytocin"),
stimulus = recode_factor(stimulus,
death = "grief"),
stimulus = fct_relevel(stimulus, "neutral")) # make neutral the reference level
Detecting old grouped_df format, replacing `vars` attribute by `groups`
# subset the sessions
data_bias_placebo <- data_bias %>%
filter(treatment == "placebo")
data_bias_ot <- data_bias %>%
filter(treatment == "oxytocin")
What are the comparisons we care most about?
To address the question, What experimental methods should we be using with this group?, we look at the contrasts of spouse > stranger
, death > neutral
, and (spouse > stranger) > (death > neutral)
. Results confirm what we see in the plots: there is a significant difference between RT bias in spouse and stranger conditions, but not death and neutral conditions. Further, the difference between those two contrasts is itself (just) statistically significant:
q2_ref <- emmeans(q2, ~stimulus+group, model = "multivariate")
q2_ref # this is the reference grid
stimulus group emmean SE df lower.CL upper.CL
neutral NCG 17.07 10.8 37 -4.88 39.02
grief NCG 1.45 19.1 37 -37.16 40.07
living NCG 40.73 15.4 37 9.53 71.92
spouse NCG 52.73 15.1 37 22.10 83.35
stranger NCG 19.57 11.4 37 -3.43 42.57
neutral CG -15.68 12.3 37 -40.65 9.29
grief CG -27.00 21.7 37 -70.93 16.93
living CG 23.32 17.5 37 -12.17 58.81
spouse CG 21.79 17.2 37 -13.05 56.64
stranger CG -20.76 12.9 37 -46.93 5.40
Confidence level used: 0.95
# vector of contrast weights must line up with items as ordered in reference grid
spouse_v_str_placebo = c(
0, # neutral NCG
0, # death NCG
0, # living NCG
.5, # spouse NCG
-.5, # stranger NCG
0, # neutral CG
0, # death CG
0, # living CG
.5, # spouse CG
-.5) # stranger CG
death_v_neu_placebo = c(
-.5, # neutral NCG
.5, # death NCG
0, # living NCG
0, # spouse NCG
0, # stranger NCG
-.5, # neutral CG
.5, # death CG
0, # living CG
0, # spouse CG
0) # stranger CG
# test the contrasts
q2_contrasts <- emmeans::contrast(q2_ref, method=list("spouse > stranger (placebo)" = spouse_v_str_placebo,
"death > neutral (placebo)" = death_v_neu_placebo,
"[spouse > stranger] > [death > neutral]" = spouse_v_str_placebo-death_v_neu_placebo))
q2_contrasts
contrast estimate SE df t.ratio p.value
spouse > stranger (placebo) 37.9 10.8 37 3.497 0.0012
death > neutral (placebo) -13.5 16.8 37 -0.803 0.4271
[spouse > stranger] > [death > neutral] 51.3 21.7 37 2.361 0.0236
# for plotting
q2_singlecontrasts <- emmeans::contrast(q2_ref, method=list("spouse > stranger \n(placebo)" = spouse_v_str_placebo,
"death > neutral \n(placebo)" = death_v_neu_placebo))
# vertical
plot(q2_singlecontrasts, comparisons = T, horizontal=F, ylab="Avoidance bias Approach bias") +
geom_hline(yintercept=0, linetype="solid", color="black", size=.2)
# horizontal
plot(q2_singlecontrasts, comparisons = T, horizontal=T, xlab="Avoidance bias Approach bias") +
geom_vline(xintercept=0, linetype="solid", color="black", size=.2)
#### Q1 (model without group)
q1_ref <- emmeans(q1, ~stimulus, model = "multivariate")
q1_ref # this is the reference grid
stimulus emmean SE df lower.CL upper.CL
neutral 2.79 8.45 38 -14.31 19.9
grief -10.95 14.31 38 -39.92 18.0
living 33.14 11.50 38 9.87 56.4
spouse 39.24 11.48 38 16.01 62.5
stranger 1.99 9.02 38 -16.27 20.2
Confidence level used: 0.95
# vector of contrast weights must line up with items as ordered in reference grid
death_v_neu_q1 = c(
-1, # neutral
1, # death
0, # living
0, # spouse
0) # stranger
spouse_v_str_q1 = c(
0, # neutral
0, # death
0, # living
1, # spouse
-1) # stranger
# test the contrasts
q1_contrasts <- emmeans::contrast(q1_ref, method=list("spouse > stranger (placebo)" = spouse_v_str_q1,
"death > neutral (placebo)" = death_v_neu_q1,
"[spouse > stranger] > [death > neutral]" = spouse_v_str_q1-death_v_neu_q1))
q1_contrasts
contrast estimate SE df t.ratio p.value
spouse > stranger (placebo) 37.3 10.6 38 3.508 0.0012
death > neutral (placebo) -13.7 16.4 38 -0.837 0.4078
[spouse > stranger] > [death > neutral] 51.0 21.3 38 2.397 0.0216
# for plotting
q1_singlecontrasts <- emmeans::contrast(q2_ref, method=list("spouse > stranger \n(placebo)" = spouse_v_str_placebo,
"death > neutral \n(placebo)" = death_v_neu_placebo))
# vertical
plot(q2_singlecontrasts, comparisons = T, horizontal=F, ylab="Avoidance bias Approach bias") +
geom_hline(yintercept=0, linetype="solid", color="black", size=.2)
# horizontal
plot(q2_singlecontrasts, comparisons = T, horizontal=T, xlab="Avoidance bias Approach bias") +
geom_vline(xintercept=0, linetype="solid", color="black", size=.2)
This isn’t the most attractive plot ever, but it does allow us to visualize the contrast results:
The blue bars are confidence intervals for the EMMs, and the red arrows are for the comparisons among them. If an arrow from one mean overlaps an arrow from another group, the difference is not “significant,” based on the correction for multiple comparisons and the value of alpha (which defaults to 0.05).
Aim 3.
To identify whether effects of intranasal oxytocin differ in high-severity [CG group] versus low-severity [Non-CG group] participants.
# A model in the full dataset (both sessions) only with two within-subjects factors (stimulus [5], treatment [2])
# and one between-subjects factor (group [2])
## DV: RT bias (push - pull)
## Predicted: Treatment x Group interaction (Aim 3), Stimulus x Treatment x Group (Aim 4)
# RM-ANOVA
q3 <- aov_ez(id = "ID", dv = "gAAT_bias", data_bias, between = c("group"), within = c("stimulus", "treatment"), factorize=T)
Contrasts set to contr.sum for the following variables: group
summary(q3)
Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
Sum Sq num Df Error SS den Df F value Pr(>F)
(Intercept) 66795 1 471795 37 5.2383 0.02790 *
group 16822 1 471795 37 1.3193 0.25809
stimulus 142185 4 609249 148 8.6350 2.723e-06 ***
group:stimulus 11688 4 609249 148 0.7098 0.58646
treatment 1347 1 136424 37 0.3652 0.54931
group:treatment 26838 1 136424 37 7.2789 0.01045 *
stimulus:treatment 6832 4 483926 148 0.5224 0.71944
group:stimulus:treatment 3202 4 483926 148 0.2448 0.91244
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Mauchly Tests for Sphericity
Test statistic p-value
stimulus 0.74057 0.301957
group:stimulus 0.74057 0.301957
stimulus:treatment 0.61959 0.049689
group:stimulus:treatment 0.61959 0.049689
Greenhouse-Geisser and Huynh-Feldt Corrections
for Departure from Sphericity
GG eps Pr(>F[GG])
stimulus 0.86061 1.101e-05 ***
group:stimulus 0.86061 0.5664
stimulus:treatment 0.82544 0.6850
group:stimulus:treatment 0.82544 0.8818
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
HF eps Pr(>F[HF])
stimulus 0.9596356 4.078504e-06
group:stimulus 0.9596356 5.809212e-01
stimulus:treatment 0.9159706 7.037024e-01
group:stimulus:treatment 0.9159706 8.990019e-01
# Stimulus is predictive in the model, but only alone.
# There is a significant two-way interaction of group x treatment
# which looks like this:
emmip(q3, group ~ treatment)
emmeans(q3, ~ c(treatment|group), contr = "pairwise", adjust="holm") #pairwise comparisons of treatment within group
$emmeans
group = NCG:
treatment emmean SE df lower.CL upper.CL
placebo 26.31 7.88 37 10.35 42.3
oxytocin 13.33 9.35 37 -5.61 32.3
group = CG:
treatment emmean SE df lower.CL upper.CL
placebo -3.66 8.96 37 -21.82 14.5
oxytocin 16.81 10.63 37 -4.73 38.4
Results are averaged over the levels of: stimulus
Confidence level used: 0.95
$contrasts
group = NCG:
contrast estimate SE df t.ratio p.value
placebo - oxytocin 13.0 8.19 37 1.586 0.1214
group = CG:
contrast estimate SE df t.ratio p.value
placebo - oxytocin -20.5 9.31 37 -2.198 0.0343
Results are averaged over the levels of: stimulus
There is a significant group
x treatment
interaction, such that the CG group became more approach-biased when given oxytocin. This is consistent with BA’s previous results.
Pairwise comparisons (treatment
within group
) confirm that the CG group is significantly more approach-biased in the OT condition.
We can also look at the pairwise comparisons for the main effect of stimulus, if that is useful/interesting in any way…
emmeans(q3, ~stimulus, contr = "pairwise", adjust="holm") # pairwise comparisons of stimulus
$emmeans
stimulus emmean SE df lower.CL upper.CL
neutral -0.63 6.96 37 -14.73 13.5
grief -6.58 9.48 37 -25.77 12.6
living 27.67 9.45 37 8.51 46.8
spouse 43.72 9.85 37 23.75 63.7
stranger 1.79 7.49 37 -13.39 17.0
Results are averaged over the levels of: group, treatment
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
neutral - grief 5.94 10.28 37 0.578 1.0000
neutral - living -28.30 8.94 37 -3.165 0.0217
neutral - spouse -44.35 9.57 37 -4.636 0.0004
neutral - stranger -2.42 8.42 37 -0.288 1.0000
grief - living -34.24 11.16 37 -3.067 0.0241
grief - spouse -50.30 13.12 37 -3.834 0.0038
grief - stranger -8.37 10.96 37 -0.763 1.0000
living - spouse -16.05 11.07 37 -1.449 0.6226
living - stranger 25.88 9.95 37 2.600 0.0666
spouse - stranger 41.93 9.30 37 4.508 0.0006
Results are averaged over the levels of: group, treatment
P value adjustment: holm method for 10 tests
We see that averaged across levels of group
and treatment
, the sample as a whole does show more approach bias for spouse compared to death, neutral, and stranger stimuli - but not living loved one.
As shown in the plot below, OT has differential effects in the two groups that generally occur across the different stimuli.
emmip(q3, stimulus ~ treatment|group) # non-significant stimulus x group x treatment interaction, shows how overall OT has differential effects in the two groups across stimuli
# create a bar plot with error bars (95% CI)
p4 <- ggplot(data_bias, aes(group, gAAT_bias, fill=treatment)) +
stat_summary(geom = "bar", fun.y = mean, position = "dodge", color="black") +
stat_summary(geom = "errorbar", fun.data = mean_cl_normal, position=position_dodge(.9), width = 0.25, color="black") +
geom_hline(yintercept=0, linetype="solid", color="black", size=.5) +
xlab("Group") + ylab("Avoid bias Approach bias") + theme_bw(base_size=14) + theme(legend.title = element_blank(), legend.text = element_text(size=14), axis.text.x=element_text(size=14, face = "bold"))
p4
ggsave("~/Desktop/p4.png", plot = p4, dpi=600, width = 8, height = 6)
# show as a boxplot:
p4b <- ggplot(data_bias, aes(group, gAAT_bias, fill=treatment)) +
geom_hline(yintercept=0, linetype="solid", color="black", size=.5) +
geom_boxplot(aes(fill=treatment), outlier.alpha = 0.3) +
xlab("Group") + ylab("Avoid bias Approach bias")
p4b
p5 <- ggplot(data_bias, aes(fill=treatment, y=gAAT_bias, x=group)) +
stat_summary(geom = "bar", fun.y = mean, position = "dodge", color="black") +
stat_summary(geom = "errorbar", fun.data = mean_cl_normal, position=position_dodge(.9), width = 0.25, color="black") +
facet_grid(. ~ stim_ordered) + theme(strip.text.x = element_text(size=12, color="red",
face="bold.italic"),
strip.text.y = element_text(size=12, color="red",
face="bold.italic")) + geom_hline(yintercept=0, linetype="solid", color="black", size=.5) + xlab("Group") + ylab("Avoid bias Approach bias") + theme(legend.title = element_blank(), legend.text = element_text(size=14), axis.text.x=element_text(size=14, face = "bold"))
p5
p6 <- ggplot(data_bias, aes(fill=treatment, y=gAAT_bias, x=stim_ordered)) +
stat_summary(geom = "bar", fun.y = mean, position = "dodge", color="black") +
stat_summary(geom = "errorbar", fun.data = mean_cl_normal, position=position_dodge(.9), width = 0.25, color="black") +
facet_grid(. ~ group) + theme(strip.text.x = element_text(size=12, color="red",
face="bold.italic"),
strip.text.y = element_text(size=12, color="red",
face="bold.italic")) + geom_hline(yintercept=0, linetype="solid", color="black", size=.5) + xlab("Group") + ylab("Avoid bias Approach bias")
p6
q3lme <- lme(gAAT_bias ~ stimulus*tot_icg_cm*treatment, random = ~ 1|ID, data_bias, method="REML")
summary(q3lme)
Linear mixed-effects model fit by REML
Data: data_bias
AIC BIC logLik
4299.202 4385.299 -2127.601
Random effects:
Formula: ~1 | ID
(Intercept) Residual
StdDev: 30.15073 60.44021
Fixed effects: gAAT_bias ~ stimulus * tot_icg_cm * treatment
Value Std.Error DF t-value p-value
(Intercept) 2.83597 10.815612 333 0.2622105 0.7933
stimulusgrief -13.72181 13.687066 333 -1.0025386 0.3168
stimulusliving 30.33701 13.687066 333 2.2164732 0.0273
stimulusspouse 36.43912 13.687066 333 2.6623030 0.0081
stimulusstranger -0.79288 13.687066 333 -0.0579295 0.9538
tot_icg_cm -1.18720 0.867564 37 -1.3684357 0.1794
treatmentoxytocin -3.73823 13.687066 333 -0.2731211 0.7849
stimulusgrief:tot_icg_cm -0.62914 1.097894 333 -0.5730399 0.5670
stimulusliving:tot_icg_cm 0.26402 1.097894 333 0.2404817 0.8101
stimulusspouse:tot_icg_cm 0.27738 1.097894 333 0.2526480 0.8007
stimulusstranger:tot_icg_cm -0.42778 1.097894 333 -0.3896384 0.6971
stimulusgrief:treatmentoxytocin 13.00832 19.356434 333 0.6720411 0.5020
stimulusliving:treatmentoxytocin -7.28588 19.356434 333 -0.3764061 0.7069
stimulusspouse:treatmentoxytocin 13.81421 19.356434 333 0.7136752 0.4759
stimulusstranger:treatmentoxytocin 6.90453 19.356434 333 0.3567044 0.7215
tot_icg_cm:treatmentoxytocin 0.95617 1.097894 333 0.8709140 0.3844
stimulusgrief:tot_icg_cm:treatmentoxytocin 1.24112 1.552657 333 0.7993540 0.4247
stimulusliving:tot_icg_cm:treatmentoxytocin 1.22172 1.552657 333 0.7868604 0.4319
stimulusspouse:tot_icg_cm:treatmentoxytocin 0.55253 1.552657 333 0.3558597 0.7222
stimulusstranger:tot_icg_cm:treatmentoxytocin -0.20482 1.552657 333 -0.1319158 0.8951
Correlation:
(Intr) stmlsg stmlsl stmlssp stmlsst tt_cg_ trtmnt stmlsg:__ stmlsl:__ stmlssp:__
stimulusgrief -0.633
stimulusliving -0.633 0.500
stimulusspouse -0.633 0.500 0.500
stimulusstranger -0.633 0.500 0.500 0.500
tot_icg_cm -0.003 0.002 0.002 0.002 0.002
treatmentoxytocin -0.633 0.500 0.500 0.500 0.500 0.002
stimulusgrief:tot_icg_cm 0.002 -0.003 -0.001 -0.001 -0.001 -0.633 -0.001
stimulusliving:tot_icg_cm 0.002 -0.001 -0.003 -0.001 -0.001 -0.633 -0.001 0.500
stimulusspouse:tot_icg_cm 0.002 -0.001 -0.001 -0.003 -0.001 -0.633 -0.001 0.500 0.500
stimulusstranger:tot_icg_cm 0.002 -0.001 -0.001 -0.001 -0.003 -0.633 -0.001 0.500 0.500 0.500
stimulusgrief:treatmentoxytocin 0.447 -0.707 -0.354 -0.354 -0.354 -0.001 -0.707 0.002 0.001 0.001
stimulusliving:treatmentoxytocin 0.447 -0.354 -0.707 -0.354 -0.354 -0.001 -0.707 0.001 0.002 0.001
stimulusspouse:treatmentoxytocin 0.447 -0.354 -0.354 -0.707 -0.354 -0.001 -0.707 0.001 0.001 0.002
stimulusstranger:treatmentoxytocin 0.447 -0.354 -0.354 -0.354 -0.707 -0.001 -0.707 0.001 0.001 0.001
tot_icg_cm:treatmentoxytocin 0.002 -0.001 -0.001 -0.001 -0.001 -0.633 -0.003 0.500 0.500 0.500
stimulusgrief:tot_icg_cm:treatmentoxytocin -0.001 0.002 0.001 0.001 0.001 0.447 0.002 -0.707 -0.354 -0.354
stimulusliving:tot_icg_cm:treatmentoxytocin -0.001 0.001 0.002 0.001 0.001 0.447 0.002 -0.354 -0.707 -0.354
stimulusspouse:tot_icg_cm:treatmentoxytocin -0.001 0.001 0.001 0.002 0.001 0.447 0.002 -0.354 -0.354 -0.707
stimulusstranger:tot_icg_cm:treatmentoxytocin -0.001 0.001 0.001 0.001 0.002 0.447 0.002 -0.354 -0.354 -0.354
stmlsst:__ stmlsg: stmlsl: stmlssp: stmlsst: tt_c_: stmlsg:__: stmlsl:__: stmlssp:__:
stimulusgrief
stimulusliving
stimulusspouse
stimulusstranger
tot_icg_cm
treatmentoxytocin
stimulusgrief:tot_icg_cm
stimulusliving:tot_icg_cm
stimulusspouse:tot_icg_cm
stimulusstranger:tot_icg_cm
stimulusgrief:treatmentoxytocin 0.001
stimulusliving:treatmentoxytocin 0.001 0.500
stimulusspouse:treatmentoxytocin 0.001 0.500 0.500
stimulusstranger:treatmentoxytocin 0.002 0.500 0.500 0.500
tot_icg_cm:treatmentoxytocin 0.500 0.002 0.002 0.002 0.002
stimulusgrief:tot_icg_cm:treatmentoxytocin -0.354 -0.003 -0.001 -0.001 -0.001 -0.707
stimulusliving:tot_icg_cm:treatmentoxytocin -0.354 -0.001 -0.003 -0.001 -0.001 -0.707 0.500
stimulusspouse:tot_icg_cm:treatmentoxytocin -0.354 -0.001 -0.001 -0.003 -0.001 -0.707 0.500 0.500
stimulusstranger:tot_icg_cm:treatmentoxytocin -0.707 -0.001 -0.001 -0.001 -0.003 -0.707 0.500 0.500 0.500
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.85308052 -0.59976610 -0.02327937 0.56193906 3.63273721
Number of Observations: 390
Number of Groups: 39
anova(q3lme)
numDF denDF F-value p-value
(Intercept) 1 333 6.036778 0.0145
stimulus 4 333 9.584126 <.0001
tot_icg_cm 1 37 1.341922 0.2541
treatment 1 333 0.068546 0.7936
stimulus:tot_icg_cm 4 333 0.991220 0.4124
stimulus:treatment 4 333 0.428223 0.7882
tot_icg_cm:treatment 1 333 9.562118 0.0022
stimulus:tot_icg_cm:treatment 4 333 0.373410 0.8276
We see a similar result when using the continuous tot_icg
variable instead of group
- a significant two-way interaction of treatment and grief severity, shown in the plot below:
library(effects)
Loading required package: carData
lattice theme set by effectsTheme()
See ?effectsTheme for details.
sd(data_bias_placebo$tot_icg)
[1] 12.49869
ef <- effect("tot_icg_cm*treatment", q3lme)
NOTE: tot_icg_cm:treatment is not a high-order term in the model
summary(ef)
tot_icg_cm*treatment effect
treatment
tot_icg_cm placebo oxytocin
-20 39.094412 10.27877
-8 23.610717 13.01446
4 8.127023 15.75016
20 -12.517903 19.39776
30 -25.420982 21.67751
Lower 95 Percent Confidence Limits
treatment
tot_icg_cm placebo oxytocin
-20 14.960970 -13.854675
-8 8.441897 -2.154356
4 -5.252641 2.370497
20 -36.591257 -4.675596
30 -58.614498 -11.516010
Upper 95 Percent Confidence Limits
treatment
tot_icg_cm placebo oxytocin
-20 63.227854 34.41221
-8 38.779538 28.18328
4 21.506687 29.12983
20 11.555450 43.47111
30 7.772534 54.87102
ef <- as.data.frame(ef)
ggplot(ef, aes(tot_icg_cm, color=treatment, fit)) + geom_point(size=2) + geom_errorbar(aes(ymin=fit-se, ymax=fit+se), width=0.6) + geom_line() +
xlab("Grief severity (mean-centered ICG)") + ylab("RT bias (fit)") + theme_bw(base_size=14) + theme(legend.title = element_blank(), legend.text = element_text(size=14), axis.text.x=element_text(size=14)) +
geom_vline(xintercept=0, linetype="dotted", color="black", size=.2)
# make a data frame with columns for the fixed-effect parameter estimates and fixed-effect parameter standard errors
df<-data.frame(cbind(summary(q3lme)$tTable[,1], summary(q3lme)$tTable[,2]))
names(df)<-c("Estimate","SE")
df
Estimate SE
(Intercept) 2.8359674 10.8156118
stimulusgrief -13.7218119 13.6870658
stimulusliving 30.3370146 13.6870658
stimulusspouse 36.4391163 13.6870658
stimulusstranger -0.7928845 13.6870658
tot_icg_cm -1.1872049 0.8675636
treatmentoxytocin -3.7382265 13.6870658
stimulusgrief:tot_icg_cm -0.6291373 1.0978944
stimulusliving:tot_icg_cm 0.2640236 1.0978944
stimulusspouse:tot_icg_cm 0.2773808 1.0978944
stimulusstranger:tot_icg_cm -0.4277819 1.0978944
stimulusgrief:treatmentoxytocin 13.0083201 19.3564341
stimulusliving:treatmentoxytocin -7.2858802 19.3564341
stimulusspouse:treatmentoxytocin 13.8142074 19.3564341
stimulusstranger:treatmentoxytocin 6.9045258 19.3564341
tot_icg_cm:treatmentoxytocin 0.9561716 1.0978944
stimulusgrief:tot_icg_cm:treatmentoxytocin 1.2411226 1.5526572
stimulusliving:tot_icg_cm:treatmentoxytocin 1.2217245 1.5526572
stimulusspouse:tot_icg_cm:treatmentoxytocin 0.5525280 1.5526572
stimulusstranger:tot_icg_cm:treatmentoxytocin -0.2048200 1.5526572
ef <- effect("tot_icg_cm:treatment", q3lme, xlevels=list(tot_icg_cm=c((-12.54-12.54), -12.54, 0, 12.54, (12.54+12.54))))
NOTE: tot_icg_cm:treatment is not a high-order term in the model
summary(ef)
tot_icg_cm*treatment effect
treatment
tot_icg_cm placebo oxytocin
-25.08 45.649176 9.120655
-12.54 29.468715 11.979459
0 13.288254 14.838262
12.54 -2.892207 17.697066
25.08 -19.072667 20.555870
Lower 95 Percent Confidence Limits
treatment
tot_icg_cm placebo oxytocin
-25.08 16.9727932 -19.5557276
-12.54 11.3589273 -6.1303290
0 0.5380056 2.0880137
12.54 -20.9517806 -0.3625079
25.08 -47.6856457 -8.0571086
Upper 95 Percent Confidence Limits
treatment
tot_icg_cm placebo oxytocin
-25.08 74.325559 37.79704
-12.54 47.578503 30.08925
0 26.038503 27.58851
12.54 15.167367 35.75664
25.08 9.540311 49.16885
# plot(predictorEffects(q3lme, "tot_icg_cm", xlevels=5))
# variance explained by the entire model
fitted_q3lme <- fitted(q3lme)
data_bias$fitted_q3lme <- as.vector(fitted_q3lme)
forR2 <- lm(gAAT_bias ~ fitted_q3lme, data=data_bias)
summary(forR2)
Call:
lm(formula = gAAT_bias ~ fitted_q3lme, data = data_bias)
Residuals:
Min 1Q Median 3Q Max
-226.853 -35.677 0.759 34.325 219.691
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -2.95372 3.08714 -0.957 0.339
fitted_q3lme 1.21031 0.08332 14.526 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 56.42 on 388 degrees of freedom
Multiple R-squared: 0.3523, Adjusted R-squared: 0.3506
F-statistic: 211 on 1 and 388 DF, p-value: < 2.2e-16
# ICC (variance explained by between-person variance)
# function from the `reghelper` package
round(reghelper::ICC(q3lme),3)
[1] 0.199
Aim 4. To identify whether the response to spouse images in the CG group significantly differs under intranasal oxytocin vs. placebo (i.e., spouse in the placebo session is significantly different than spouse in the oxytocin session, in the CG group only).
As shown about, there was no three-way interaction of stimulus
x group
x treatment
.
In this question, we are specifically looking at effects of group
x treatment
within stimulus:spouse
.
q4_ref <- emmeans(q3, ~stimulus+group+treatment, model = "multivariate")
q4_ref # this is the reference grid
stimulus group treatment emmean SE df lower.CL upper.CL
neutral NCG placebo 17.07 10.8 37 -4.88 39.02
grief NCG placebo 1.45 19.1 37 -37.16 40.07
living NCG placebo 40.73 15.4 37 9.53 71.92
spouse NCG placebo 52.73 15.1 37 22.10 83.35
stranger NCG placebo 19.57 11.4 37 -3.43 42.57
neutral CG placebo -15.68 12.3 37 -40.65 9.29
grief CG placebo -27.00 21.7 37 -70.93 16.93
living CG placebo 23.32 17.5 37 -12.17 58.81
spouse CG placebo 21.79 17.2 37 -13.05 56.64
stranger CG placebo -20.76 12.9 37 -46.93 5.40
neutral NCG oxytocin 6.20 12.5 37 -19.04 31.45
grief NCG oxytocin -9.93 12.1 37 -34.54 14.67
living NCG oxytocin 14.57 16.6 37 -19.06 48.20
spouse NCG oxytocin 43.86 15.5 37 12.50 75.23
stranger NCG oxytocin 11.93 14.0 37 -16.53 40.39
neutral CG oxytocin -10.12 14.2 37 -38.84 18.60
grief CG oxytocin 9.18 13.8 37 -18.81 37.17
living CG oxytocin 32.06 18.9 37 -6.20 70.32
spouse CG oxytocin 56.50 17.6 37 20.82 92.18
stranger CG oxytocin -3.56 16.0 37 -35.93 28.82
Confidence level used: 0.95
# vector of contrast weights must line up with items as ordered in reference grid
spouseOT_v_spousePL_CG = c(
0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
0, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
-1, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
0, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
1, # spouse CG oxytocin
0) # stranger CG oxytocin
spouseOT_v_spousePL_NCG = c(
0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
-1, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
0, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
1, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
strangerOT_v_strangerPL_CG = c(
0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
0, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
0, # spouse CG placebo
-1, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
0, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
1) # stranger CG oxytocin
strangerOT_v_strangerPL_NCG = c(
0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
0, # spouse NCG placebo
-1, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
0, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
0, # spouse NCG oxytocin
1, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
deathOT_v_deathPL_CG = c(
0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
0, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
-1, # death CG placebo
0, # living CG placebo
0, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
0, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
1, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
deathOT_v_deathPL_NCG = c(
0, # neutral NCG placebo
-1, # death NCG placebo
0, # living NCG placebo
0, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
0, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
1, # death NCG oxytocin
0, # living NCG oxytocin
0, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
livingOT_v_livingPL_CG = c(
0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
0, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
-1, # living CG placebo
0, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
0, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
1, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
livingOT_v_livingPL_NCG = c(
0, # neutral NCG placebo
0, # death NCG placebo
-1, # living NCG placebo
0, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
0, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
1, # living NCG oxytocin
0, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
# test the contrasts
q4_contrasts <- emmeans::contrast(q4_ref, method=list("CG: spouse OT > spouse PL" = spouseOT_v_spousePL_CG,
"NCG: spouse OT > spouse PL" = spouseOT_v_spousePL_NCG))
q4_contrasts
contrast estimate SE df t.ratio p.value
CG: spouse OT > spouse PL 34.71 18.3 37 1.896 0.0657
NCG: spouse OT > spouse PL -8.86 16.1 37 -0.551 0.5850
There is a non-significant effect of oxytocin on spouse-related bias in the CG group, in which oxytocin increases approach bias relative to placebo, t = 1.90, p = .066. There is no effect of oxytocin on spouse-related bias in the NCG group.
# the vector of contrasts below is equivalent to
# q4con <- contrast(q4_ref, method=list("Q4 contrast" = spouseOT_v_spousePL_CG - spouseOT_v_spousePL_NCG))
CGspouseOT_v_spousePL__v__NCGspouseOT_v_spousePL = c(
0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
1, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
-1, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
-1, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
1, # spouse CG oxytocin
0) # stranger CG oxytocin
q4con <- emmeans::contrast(q4_ref, method=list("[CG spouseOT > spousePL] > [NCG spouseOT > spousePL]" = CGspouseOT_v_spousePL__v__NCGspouseOT_v_spousePL))
q4con
contrast estimate SE df t.ratio p.value
[CG spouseOT > spousePL] > [NCG spouseOT > spousePL] 43.6 24.4 37 1.788 0.0819
# for plotting
q4_singlecontrasts <- emmeans::contrast(q4_ref, method=list("CG: spouse OT > spouse PL" = spouseOT_v_spousePL_CG,
"NCG: spouse OT > spouse PL" = spouseOT_v_spousePL_NCG))
q4_singlecontrasts
contrast estimate SE df t.ratio p.value
CG: spouse OT > spouse PL 34.71 18.3 37 1.896 0.0657
NCG: spouse OT > spouse PL -8.86 16.1 37 -0.551 0.5850
# vertical
plot(q4_singlecontrasts, comparisons = T, horizontal=F, ylab="Avoidance bias Approach bias") +
geom_hline(yintercept=0, linetype="solid", color="black", size=.2)
# horizontal
plot(q4_singlecontrasts, comparisons = T, horizontal=T, xlab="Avoidance bias Approach bias") +
geom_vline(xintercept=0, linetype="solid", color="black", size=.2)
The plots above show the group comparison of oxytocin effects on spouse-related bias. The difference between CG and NCG was not statistically significant, t = 1.79, p = .082.
# load in behavioral data (long dataset containing trial-level reaction times)
data_long <- readRDS("~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/bx_long_noOut.rds")
# read in subset variables from master dataset to add to behavioral dataset
varstoadd <- readRDS("~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/varstoadd.rds")
# The model requires a stacked/long dataset containing the following, at minimum:
## one variable `ID` with participant IDs
## one numeric variable `gAAT_RT` with TRIAL-LEVEL reaction times
## one factor `stimulus` with 5 levels (spouse, stranger, loved one, grief scenes, neutral scenes)
## one factor `treatment` with 2 levels (oxytocin, placebo)
## one factor `direction` with 2 levels (push, pull)
## one factor `group` with 2 factors (CG, NCG) *OR* one numeric variable `tot_icg` with individual grief severity scores (time-invariant)
meancent <- function(x) { x - mean(x, na.rm=TRUE) } #simple worker function to mean center a variable
varstoadd <- varstoadd %>% mutate_at(vars("tot_icg", "tot_bdi"), funs(cm=meancent))
data_long <- left_join(data_long, varstoadd, ID="ID") %>% rename(treatment = cond,
stimulus = stim,
direction = push_pull) %>%
mutate(treatment = recode(treatment,
A = "placebo",
B = "oxytocin"),
stimulus = fct_relevel(stimulus, "neutral"))
Joining, by = "ID"
data_long <- data_long %>% mutate(spouse = as.factor(ifelse(stimulus == "spouse", 1, 0))) # dummy-coded variable for spouse
####### Aims 3/4 model (LME) ######
# A model with three within-subjects factors (spouse [2], direction [2], treatment [2]) and one between-subjects factor (group [2])
## DV: RTs on push and pull trials separately
# Linear mixed model
q4b <- lme(gAAT_RT ~ direction*spouse*group*treatment, data = data_long, random = ~ 1|ID, method="REML")
anova(q4b)
numDF denDF F-value p-value
(Intercept) 1 21955 3375.408 <.0001
direction 1 21955 26.625 <.0001
spouse 1 21955 36.814 <.0001
group 1 37 1.957 0.1702
treatment 1 21955 43.210 <.0001
direction:spouse 1 21955 19.721 <.0001
direction:group 1 21955 0.843 0.3586
spouse:group 1 21955 3.526 0.0604
direction:treatment 1 21955 0.010 0.9199
spouse:treatment 1 21955 2.348 0.1254
group:treatment 1 21955 41.697 <.0001
direction:spouse:group 1 21955 0.253 0.6153
direction:spouse:treatment 1 21955 0.001 0.9729
direction:group:treatment 1 21955 10.097 0.0015
spouse:group:treatment 1 21955 0.051 0.8222
direction:spouse:group:treatment 1 21955 0.992 0.3192
summary(q4b)
Linear mixed-effects model fit by REML
Data: data_long
AIC BIC logLik
292704.2 292848.2 -146334.1
Random effects:
Formula: ~1 | ID
(Intercept) Residual
StdDev: 84.36568 186.5109
Fixed effects: gAAT_RT ~ direction * spouse * group * treatment
Value Std.Error DF t-value p-value
(Intercept) 761.2938 18.371921 21955 41.43790 0.0000
directionPUSH 15.0130 5.280863 21955 2.84290 0.0045
spouse1 -9.7174 8.439500 21955 -1.15142 0.2496
groupCG 27.7401 27.817548 37 0.99721 0.3251
treatmentoxytocin 6.2350 5.282687 21955 1.18027 0.2379
directionPUSH:spouse1 36.8890 11.971335 21955 3.08144 0.0021
directionPUSH:groupCG -17.0591 7.987169 21955 -2.13581 0.0327
spouse1:groupCG 23.1558 12.786235 21955 1.81100 0.0702
directionPUSH:treatmentoxytocin -12.1292 7.461599 21955 -1.62555 0.1041
spouse1:treatmentoxytocin 16.6552 11.885122 21955 1.40135 0.1611
groupCG:treatmentoxytocin 19.7527 7.967166 21955 2.47926 0.0132
directionPUSH:spouse1:groupCG -19.2106 18.130061 21955 -1.05960 0.2893
directionPUSH:spouse1:treatmentoxytocin -11.4352 16.903268 21955 -0.67651 0.4987
directionPUSH:groupCG:treatmentoxytocin 27.2741 11.298759 21955 2.41391 0.0158
spouse1:groupCG:treatmentoxytocin -15.5399 18.054943 21955 -0.86070 0.3894
directionPUSH:spouse1:groupCG:treatmentoxytocin 25.5584 25.658271 21955 0.99611 0.3192
Correlation:
(Intr) drPUSH spous1 gropCG trtmnt drPUSH:1 drPUSH:CG sp1:CG drPUSH: sps1:t grpCG:
directionPUSH -0.144
spouse1 -0.090 0.314
groupCG -0.660 0.095 0.060
treatmentoxytocin -0.144 0.502 0.314 0.095
directionPUSH:spouse1 0.064 -0.441 -0.705 -0.042 -0.221
directionPUSH:groupCG 0.095 -0.661 -0.208 -0.142 -0.332 0.292
spouse1:groupCG 0.060 -0.207 -0.660 -0.089 -0.207 0.465 0.309
directionPUSH:treatmentoxytocin 0.102 -0.708 -0.222 -0.067 -0.708 0.312 0.468 0.147
spouse1:treatmentoxytocin 0.064 -0.223 -0.710 -0.042 -0.445 0.501 0.147 0.469 0.315
groupCG:treatmentoxytocin 0.096 -0.333 -0.208 -0.143 -0.663 0.147 0.497 0.310 0.469 0.295
directionPUSH:spouse1:groupCG -0.042 0.291 0.466 0.063 0.146 -0.660 -0.441 -0.705 -0.206 -0.331 -0.219
directionPUSH:spouse1:treatmentoxytocin -0.045 0.312 0.499 0.030 0.313 -0.708 -0.207 -0.330 -0.441 -0.703 -0.207
directionPUSH:groupCG:treatmentoxytocin -0.067 0.467 0.147 0.101 0.468 -0.206 -0.707 -0.219 -0.660 -0.208 -0.705
spouse1:groupCG:treatmentoxytocin -0.042 0.147 0.467 0.063 0.293 -0.330 -0.219 -0.708 -0.207 -0.658 -0.441
directionPUSH:spouse1:groupCG:treatmentoxytocin 0.030 -0.206 -0.329 -0.044 -0.206 0.467 0.311 0.498 0.291 0.463 0.311
drPUSH:1:CG drPUSH:1: dPUSH:CG: s1:CG:
directionPUSH
spouse1
groupCG
treatmentoxytocin
directionPUSH:spouse1
directionPUSH:groupCG
spouse1:groupCG
directionPUSH:treatmentoxytocin
spouse1:treatmentoxytocin
groupCG:treatmentoxytocin
directionPUSH:spouse1:groupCG
directionPUSH:spouse1:treatmentoxytocin 0.468
directionPUSH:groupCG:treatmentoxytocin 0.311 0.292
spouse1:groupCG:treatmentoxytocin 0.499 0.463 0.311
directionPUSH:spouse1:groupCG:treatmentoxytocin -0.707 -0.659 -0.440 -0.704
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.3275393 -0.6751297 -0.1941000 0.4703443 5.2964704
Number of Observations: 22008
Number of Groups: 39
# how to manually specify contrasts:
# https://stats.stackexchange.com/questions/348641/specify-contrasts-for-lme-with-interactions
### CONTRASTS ###
q4_ref_lme <- emmeans(q4b, ~spouse+group+treatment+direction, model = "multivariate")
q4_ref_lme # this is the reference grid
spouse group treatment direction emmean SE df lower.CL upper.CL
0 NCG placebo PULL 761 18.4 38 724 798
1 NCG placebo PULL 752 19.5 38 712 791
0 CG placebo PULL 789 20.9 37 747 831
1 CG placebo PULL 802 22.2 37 757 847
0 NCG oxytocin PULL 768 18.4 38 730 805
1 NCG oxytocin PULL 774 19.5 38 735 814
0 CG oxytocin PULL 815 20.9 37 773 857
1 CG oxytocin PULL 830 22.2 37 785 875
0 NCG placebo PUSH 776 18.4 38 739 813
1 NCG placebo PUSH 803 19.5 38 764 843
0 CG placebo PUSH 787 20.9 37 745 829
1 CG placebo PUSH 818 22.2 37 773 863
0 NCG oxytocin PUSH 770 18.4 38 733 808
1 NCG oxytocin PUSH 803 19.5 38 763 842
0 CG oxytocin PUSH 828 20.9 37 786 870
1 CG oxytocin PUSH 874 22.3 37 829 920
d.f. method: containment
Confidence level used: 0.95
# vector of contrast weights must line up with items as ordered in reference grid
spouseOT_v_spousePL__NCG_PULL = c(
0, # 0 NCG placebo PULL 761 18.4 38 724 798
-1, # 1 NCG placebo PULL 752 19.5 38 712 791
0, # 0 CG placebo PULL 789 20.9 37 747 831
0, # 1 CG placebo PULL 802 22.2 37 757 847
0, # 0 NCG oxytocin PULL 768 18.4 38 730 805
1, # 1 NCG oxytocin PULL 774 19.5 38 735 814
0, # 0 CG oxytocin PULL 815 20.9 37 773 857
0, # 1 CG oxytocin PULL 830 22.2 37 785 875
0, # 0 NCG placebo PUSH 776 18.4 38 739 813
0, # 1 NCG placebo PUSH 803 19.5 38 764 843
0, # 0 CG placebo PUSH 787 20.9 37 745 829
0, # 1 CG placebo PUSH 818 22.2 37 773 863
0, # 0 NCG oxytocin PUSH 770 18.4 38 733 808
0, # 1 NCG oxytocin PUSH 803 19.5 38 763 842
0, # 0 CG oxytocin PUSH 828 20.9 37 786 870
0) # 1 CG oxytocin PUSH 874 22.3 37 829 920
spouseOT_v_spousePL__NCG_PUSH = c(
0, # 0 NCG placebo PULL 761 18.4 38 724 798
0, # 1 NCG placebo PULL 752 19.5 38 712 791
0, # 0 CG placebo PULL 789 20.9 37 747 831
0, # 1 CG placebo PULL 802 22.2 37 757 847
0, # 0 NCG oxytocin PULL 768 18.4 38 730 805
0, # 1 NCG oxytocin PULL 774 19.5 38 735 814
0, # 0 CG oxytocin PULL 815 20.9 37 773 857
0, # 1 CG oxytocin PULL 830 22.2 37 785 875
0, # 0 NCG placebo PUSH 776 18.4 38 739 813
-1, # 1 NCG placebo PUSH 803 19.5 38 764 843
0, # 0 CG placebo PUSH 787 20.9 37 745 829
0, # 1 CG placebo PUSH 818 22.2 37 773 863
0, # 0 NCG oxytocin PUSH 770 18.4 38 733 808
1, # 1 NCG oxytocin PUSH 803 19.5 38 763 842
0, # 0 CG oxytocin PUSH 828 20.9 37 786 870
0) # 1 CG oxytocin PUSH 874 22.3 37 829 920
spouseOT_v_spousePL__CG_PULL = c(
0, # 0 NCG placebo PULL 761 18.4 38 724 798
0, # 1 NCG placebo PULL 752 19.5 38 712 791
0, # 0 CG placebo PULL 789 20.9 37 747 831
-1, # 1 CG placebo PULL 802 22.2 37 757 847
0, # 0 NCG oxytocin PULL 768 18.4 38 730 805
0, # 1 NCG oxytocin PULL 774 19.5 38 735 814
0, # 0 CG oxytocin PULL 815 20.9 37 773 857
1, # 1 CG oxytocin PULL 830 22.2 37 785 875
0, # 0 NCG placebo PUSH 776 18.4 38 739 813
0, # 1 NCG placebo PUSH 803 19.5 38 764 843
0, # 0 CG placebo PUSH 787 20.9 37 745 829
0, # 1 CG placebo PUSH 818 22.2 37 773 863
0, # 0 NCG oxytocin PUSH 770 18.4 38 733 808
0, # 1 NCG oxytocin PUSH 803 19.5 38 763 842
0, # 0 CG oxytocin PUSH 828 20.9 37 786 870
0) # 1 CG oxytocin PUSH 874 22.3 37 829 920
spouseOT_v_spousePL__CG_PUSH = c(
0, # 0 NCG placebo PULL 761 18.4 38 724 798
0, # 1 NCG placebo PULL 752 19.5 38 712 791
0, # 0 CG placebo PULL 789 20.9 37 747 831
0, # 1 CG placebo PULL 802 22.2 37 757 847
0, # 0 NCG oxytocin PULL 768 18.4 38 730 805
0, # 1 NCG oxytocin PULL 774 19.5 38 735 814
0, # 0 CG oxytocin PULL 815 20.9 37 773 857
0, # 1 CG oxytocin PULL 830 22.2 37 785 875
0, # 0 NCG placebo PUSH 776 18.4 38 739 813
0, # 1 NCG placebo PUSH 803 19.5 38 764 843
0, # 0 CG placebo PUSH 787 20.9 37 745 829
-1, # 1 CG placebo PUSH 818 22.2 37 773 863
0, # 0 NCG oxytocin PUSH 770 18.4 38 733 808
0, # 1 NCG oxytocin PUSH 803 19.5 38 763 842
0, # 0 CG oxytocin PUSH 828 20.9 37 786 870
1) # 1 CG oxytocin PUSH 874 22.3 37 829 920
# test the contrasts
q4_contrasts_PULL <- emmeans::contrast(q4_ref_lme, method=list("OT > PL (spouse/NCG/PULL)" = spouseOT_v_spousePL__NCG_PULL,
"OT > PL (spouse/CG/PULL)" = spouseOT_v_spousePL__CG_PULL,
"[OT > PL (spouse/CG/PULL)] > [OT > PL (spouse/NCG/PULL)]" = spouseOT_v_spousePL__CG_PULL-spouseOT_v_spousePL__NCG_PULL))
q4_contrasts_PUSH <- emmeans::contrast(q4_ref_lme, method=list("OT > PL (spouse/NCG/PUSH)" = spouseOT_v_spousePL__NCG_PUSH,
"OT > PL (spouse/CG/PUSH)" = spouseOT_v_spousePL__CG_PUSH,
"[OT > PL (spouse/CG/PUSH)] > [OT > PL (spouse/NCG/PUSH)]" = spouseOT_v_spousePL__CG_PUSH-spouseOT_v_spousePL__NCG_PUSH))
q4_contrasts_PULL
contrast estimate SE df t.ratio p.value
OT > PL (spouse/NCG/PULL) 22.89 10.6 21955 2.150 0.0316
OT > PL (spouse/CG/PULL) 27.10 12.2 21955 2.219 0.0265
[OT > PL (spouse/CG/PULL)] > [OT > PL (spouse/NCG/PULL)] 4.21 16.2 21955 0.260 0.7949
q4_contrasts_PUSH
contrast estimate SE df t.ratio p.value
OT > PL (spouse/NCG/PUSH) -0.674 10.8 21955 -0.062 0.9502
OT > PL (spouse/CG/PUSH) 56.371 12.3 21955 4.581 <.0001
[OT > PL (spouse/CG/PUSH)] > [OT > PL (spouse/NCG/PUSH)] 57.045 16.4 21955 3.484 0.0005
# contrast push-pull
q4_contrasts_PUSH_PULL <- emmeans::contrast(q4_ref_lme, method=list(
"( [OT>PL (spouse/CG/PUSH)] - [OT>PL (spouse/CG/PULL)] ) > ( [OT>PL (spouse/NCG/PUSH)] - [OT>PL (spouse/NCG/PULL)] )" = (spouseOT_v_spousePL__CG_PUSH-spouseOT_v_spousePL__CG_PULL) - (spouseOT_v_spousePL__NCG_PUSH-spouseOT_v_spousePL__NCG_PULL)))
q4_contrasts_PUSH_PULL
contrast estimate SE
( [OT>PL (spouse/CG/PUSH)] - [OT>PL (spouse/CG/PULL)] ) > ( [OT>PL (spouse/NCG/PUSH)] - [OT>PL (spouse/NCG/PULL)] ) 52.8 23
df t.ratio p.value
21955 2.293 0.0218
# contrasts (push-pull by group) for plotting:
# these are just for plotting, don't care about whether the push-pull difference is significant
q4_contrasts_PUSH_PULL_plot <- emmeans::contrast(q4_ref_lme, method=list(
"( [OT>PL (spouse/CG/PUSH)] > \n [OT>PL (spouse/CG/PULL)] )" = spouseOT_v_spousePL__CG_PUSH-spouseOT_v_spousePL__CG_PULL,
"( [OT>PL (spouse/NCG/PUSH)] > \n [OT>PL (spouse/NCG/PULL)] )" = spouseOT_v_spousePL__NCG_PUSH-spouseOT_v_spousePL__NCG_PULL))
q4_contrasts_PUSH_PULL_plot
contrast estimate SE df t.ratio p.value
( [OT>PL (spouse/CG/PUSH)] > \n [OT>PL (spouse/CG/PULL)] ) 29.3 17.3 21955 1.688 0.0914
( [OT>PL (spouse/NCG/PUSH)] > \n [OT>PL (spouse/NCG/PULL)] ) -23.6 15.2 21955 -1.554 0.1203
# horizontal
plot(q4_contrasts_PUSH_PULL_plot, comparisons = T, horizontal=T, xlab="Avoidance bias Approach bias") +
geom_vline(xintercept=0, linetype="solid", color="black", size=.2)
data_long_spouse <- data_long %>% filter(spouse == 1) # get spouse data only to visualize
# create a bar plot with error bars (95% CI)
p7 <- ggplot(data_long_spouse, aes(fill=treatment, y=gAAT_RT, x=group)) +
stat_summary(geom = "bar", fun.y = mean, position = "dodge", color="black") +
stat_summary(geom = "errorbar", fun.data = mean_cl_normal, position=position_dodge(.9), width = 0.25, color="black") +
facet_grid(. ~ direction) + theme(strip.text.x = element_text(size=12, color="black",
face="bold.italic")) +
geom_hline(yintercept=0, linetype="solid", color="black", size=.5) + xlab("Group") + ylab("Mean RT to spouse photos (ms)")
p7
FYI, we get the same exact results w/dummy-coded spouse
as using stimulus
and doing contrasts only on the spouse category.
# these are the contrast vectors I'd use w/stimulus in the model instead of dummy-coded spouse
# Linear mixed model
q4b <- lme(gAAT_RT ~ direction*stimulus*group*treatment, data = data_long, random = ~ 1|ID, method="REML")
anova(q4b)
numDF denDF F-value p-value
(Intercept) 1 21931 3371.640 <.0001
direction 1 21931 26.683 <.0001
stimulus 4 21931 19.311 <.0001
group 1 37 1.956 0.1702
treatment 1 21931 43.266 <.0001
direction:stimulus 4 21931 8.848 <.0001
direction:group 1 21931 0.783 0.3764
stimulus:group 4 21931 1.937 0.1013
direction:treatment 1 21931 0.015 0.9021
stimulus:treatment 4 21931 1.455 0.2129
group:treatment 1 21931 41.817 <.0001
direction:stimulus:group 4 21931 0.392 0.8147
direction:stimulus:treatment 4 21931 0.879 0.4755
direction:group:treatment 1 21931 10.006 0.0016
stimulus:group:treatment 4 21931 0.343 0.8490
direction:stimulus:group:treatment 4 21931 0.837 0.5017
q4_ref_lme <- emmeans(q4b, ~stimulus+group+treatment+direction, model = "multivariate")
q4_ref_lme # this is the reference grid
stimulus group treatment direction emmean SE df lower.CL upper.CL
neutral NCG placebo PULL 740 19.6 38 700 779
death NCG placebo PULL 761 19.7 38 721 801
living NCG placebo PULL 757 19.4 38 718 796
spouse NCG placebo PULL 752 19.5 38 712 791
stranger NCG placebo PULL 784 19.3 38 745 823
neutral CG placebo PULL 773 22.3 37 728 818
death CG placebo PULL 794 22.3 37 749 839
living CG placebo PULL 774 22.0 37 729 819
spouse CG placebo PULL 802 22.2 37 757 847
stranger CG placebo PULL 812 21.9 37 768 856
neutral NCG oxytocin PULL 762 19.6 38 723 802
death NCG oxytocin PULL 763 19.6 38 724 803
living NCG oxytocin PULL 764 19.4 38 725 803
spouse NCG oxytocin PULL 774 19.5 38 735 814
stranger NCG oxytocin PULL 779 19.4 38 740 818
neutral CG oxytocin PULL 811 22.3 37 766 856
death CG oxytocin PULL 815 22.4 37 769 860
living CG oxytocin PULL 793 22.1 37 748 838
spouse CG oxytocin PULL 830 22.2 37 785 875
stranger CG oxytocin PULL 840 22.0 37 795 884
neutral NCG placebo PUSH 763 19.6 38 723 802
death NCG placebo PUSH 769 19.6 38 729 809
living NCG placebo PUSH 791 19.4 38 752 830
spouse NCG placebo PUSH 803 19.5 38 764 843
stranger NCG placebo PUSH 779 19.4 38 740 818
neutral CG placebo PUSH 781 22.3 37 736 826
death CG placebo PUSH 776 22.4 37 731 821
living CG placebo PUSH 793 22.1 37 748 837
spouse CG placebo PUSH 818 22.2 37 773 863
stranger CG placebo PUSH 796 22.1 37 751 840
neutral NCG oxytocin PUSH 766 19.6 38 726 806
death NCG oxytocin PUSH 764 19.7 38 725 804
living NCG oxytocin PUSH 770 19.4 38 731 810
spouse NCG oxytocin PUSH 803 19.6 38 763 842
stranger NCG oxytocin PUSH 779 19.3 38 740 818
neutral CG oxytocin PUSH 814 22.3 37 769 859
death CG oxytocin PUSH 834 22.4 37 789 879
living CG oxytocin PUSH 829 22.1 37 784 874
spouse CG oxytocin PUSH 874 22.3 37 829 920
stranger CG oxytocin PUSH 835 22.0 37 790 880
d.f. method: containment
Confidence level used: 0.95
# vector of contrast weights must line up with items as ordered in reference grid
spouseOT_v_spousePL__NCG_PULL = c(
0 , # neutral NCG placebo PULL
0 , # death NCG placebo PULL
0 , # living NCG placebo PULL
-1 , # spouse NCG placebo PULL
0 , # stranger NCG placebo PULL
0 , # neutral CG placebo PULL
0 , # death CG placebo PULL
0 , # living CG placebo PULL
0 , # spouse CG placebo PULL
0 , # stranger CG placebo PULL
0 , # neutral NCG oxytocin PULL
0 , # death NCG oxytocin PULL
0 , # living NCG oxytocin PULL
1 , # spouse NCG oxytocin PULL
0 , # stranger NCG oxytocin PULL
0 , # neutral CG oxytocin PULL
0 , # death CG oxytocin PULL
0 , # living CG oxytocin PULL
0 , # spouse CG oxytocin PULL
0 , # stranger CG oxytocin PULL
0 , # neutral NCG placebo PUSH
0 , # death NCG placebo PUSH
0 , # living NCG placebo PUSH
0 , # spouse NCG placebo PUSH
0 , # stranger NCG placebo PUSH
0 , # neutral CG placebo PUSH
0 , # death CG placebo PUSH
0 , # living CG placebo PUSH
0 , # spouse CG placebo PUSH
0 , # stranger CG placebo PUSH
0 , # neutral NCG oxytocin PUSH
0 , # death NCG oxytocin PUSH
0 , # living NCG oxytocin PUSH
0 , # spouse NCG oxytocin PUSH
0 , # stranger NCG oxytocin PUSH
0 , # neutral CG oxytocin PUSH
0 , # death CG oxytocin PUSH
0 , # living CG oxytocin PUSH
0 , # spouse CG oxytocin PUSH
0 ) # stranger CG oxytocin PUSH
spouseOT_v_spousePL__CG_PULL = c(
0 , # neutral NCG placebo PULL
0 , # death NCG placebo PULL
0 , # living NCG placebo PULL
0 , # spouse NCG placebo PULL
0 , # stranger NCG placebo PULL
0 , # neutral CG placebo PULL
0 , # death CG placebo PULL
0 , # living CG placebo PULL
-1 , # spouse CG placebo PULL
0 , # stranger CG placebo PULL
0 , # neutral NCG oxytocin PULL
0 , # death NCG oxytocin PULL
0 , # living NCG oxytocin PULL
0 , # spouse NCG oxytocin PULL
0 , # stranger NCG oxytocin PULL
0 , # neutral CG oxytocin PULL
0 , # death CG oxytocin PULL
0 , # living CG oxytocin PULL
1 , # spouse CG oxytocin PULL
0 , # stranger CG oxytocin PULL
0 , # neutral NCG placebo PUSH
0 , # death NCG placebo PUSH
0 , # living NCG placebo PUSH
0 , # spouse NCG placebo PUSH
0 , # stranger NCG placebo PUSH
0 , # neutral CG placebo PUSH
0 , # death CG placebo PUSH
0 , # living CG placebo PUSH
0 , # spouse CG placebo PUSH
0 , # stranger CG placebo PUSH
0 , # neutral NCG oxytocin PUSH
0 , # death NCG oxytocin PUSH
0 , # living NCG oxytocin PUSH
0 , # spouse NCG oxytocin PUSH
0 , # stranger NCG oxytocin PUSH
0 , # neutral CG oxytocin PUSH
0 , # death CG oxytocin PUSH
0 , # living CG oxytocin PUSH
0 , # spouse CG oxytocin PUSH
0 ) # stranger CG oxytocin PUSH
spouseOT_v_spousePL__NCG_PUSH = c(
0 , # neutral NCG placebo PULL
0 , # death NCG placebo PULL
0 , # living NCG placebo PULL
0 , # spouse NCG placebo PULL
0 , # stranger NCG placebo PULL
0 , # neutral CG placebo PULL
0 , # death CG placebo PULL
0 , # living CG placebo PULL
0 , # spouse CG placebo PULL
0 , # stranger CG placebo PULL
0 , # neutral NCG oxytocin PULL
0 , # death NCG oxytocin PULL
0 , # living NCG oxytocin PULL
0 , # spouse NCG oxytocin PULL
0 , # stranger NCG oxytocin PULL
0 , # neutral CG oxytocin PULL
0 , # death CG oxytocin PULL
0 , # living CG oxytocin PULL
0 , # spouse CG oxytocin PULL
0 , # stranger CG oxytocin PULL
0 , # neutral NCG placebo PUSH
0 , # death NCG placebo PUSH
0 , # living NCG placebo PUSH
-1 , # spouse NCG placebo PUSH
0 , # stranger NCG placebo PUSH
0 , # neutral CG placebo PUSH
0 , # death CG placebo PUSH
0 , # living CG placebo PUSH
0 , # spouse CG placebo PUSH
0 , # stranger CG placebo PUSH
0 , # neutral NCG oxytocin PUSH
0 , # death NCG oxytocin PUSH
0 , # living NCG oxytocin PUSH
1 , # spouse NCG oxytocin PUSH
0 , # stranger NCG oxytocin PUSH
0 , # neutral CG oxytocin PUSH
0 , # death CG oxytocin PUSH
0 , # living CG oxytocin PUSH
0 , # spouse CG oxytocin PUSH
0 ) # stranger CG oxytocin PUSH
spouseOT_v_spousePL__CG_PUSH = c(
0 , # neutral NCG placebo PULL
0 , # death NCG placebo PULL
0 , # living NCG placebo PULL
0 , # spouse NCG placebo PULL
0 , # stranger NCG placebo PULL
0 , # neutral CG placebo PULL
0 , # death CG placebo PULL
0 , # living CG placebo PULL
0 , # spouse CG placebo PULL
0 , # stranger CG placebo PULL
0 , # neutral NCG oxytocin PULL
0 , # death NCG oxytocin PULL
0 , # living NCG oxytocin PULL
0 , # spouse NCG oxytocin PULL
0 , # stranger NCG oxytocin PULL
0 , # neutral CG oxytocin PULL
0 , # death CG oxytocin PULL
0 , # living CG oxytocin PULL
0 , # spouse CG oxytocin PULL
0 , # stranger CG oxytocin PULL
0 , # neutral NCG placebo PUSH
0 , # death NCG placebo PUSH
0 , # living NCG placebo PUSH
0 , # spouse NCG placebo PUSH
0 , # stranger NCG placebo PUSH
0 , # neutral CG placebo PUSH
0 , # death CG placebo PUSH
0 , # living CG placebo PUSH
-1 , # spouse CG placebo PUSH
0 , # stranger CG placebo PUSH
0 , # neutral NCG oxytocin PUSH
0 , # death NCG oxytocin PUSH
0 , # living NCG oxytocin PUSH
0 , # spouse NCG oxytocin PUSH
0 , # stranger NCG oxytocin PUSH
0 , # neutral CG oxytocin PUSH
0 , # death CG oxytocin PUSH
0 , # living CG oxytocin PUSH
1 , # spouse CG oxytocin PUSH
0 ) # stranger CG oxytocin PUSH
# test the contrasts
q4_contrasts_PULL <- emmeans::contrast(q4_ref_lme, method=list("OT > PL (spouse/NCG/PULL)" = spouseOT_v_spousePL__NCG_PULL,
"OT > PL (spouse/CG/PULL)" = spouseOT_v_spousePL__CG_PULL,
"[OT > PL (spouse/CG/PULL)] > [OT > PL (spouse/NCG/PULL)]" = spouseOT_v_spousePL__CG_PULL-spouseOT_v_spousePL__NCG_PULL))
q4_contrasts_PUSH <- emmeans::contrast(q4_ref_lme, method=list("OT > PL (spouse/NCG/PUSH)" = spouseOT_v_spousePL__NCG_PUSH,
"OT > PL (spouse/CG/PUSH)" = spouseOT_v_spousePL__CG_PUSH,
"[OT > PL (spouse/CG/PUSH)] > [OT > PL (spouse/NCG/PUSH)]" = spouseOT_v_spousePL__CG_PUSH-spouseOT_v_spousePL__NCG_PUSH))
q4_contrasts_PULL
contrast estimate SE df t.ratio p.value
OT > PL (spouse/NCG/PULL) 22.89 10.6 21931 2.153 0.0313
OT > PL (spouse/CG/PULL) 27.10 12.2 21931 2.221 0.0263
[OT > PL (spouse/CG/PULL)] > [OT > PL (spouse/NCG/PULL)] 4.21 16.2 21931 0.260 0.7949
q4_contrasts_PUSH
contrast estimate SE df t.ratio p.value
OT > PL (spouse/NCG/PUSH) -0.68 10.8 21931 -0.063 0.9497
OT > PL (spouse/CG/PUSH) 56.37 12.3 21931 4.586 <.0001
[OT > PL (spouse/CG/PUSH)] > [OT > PL (spouse/NCG/PUSH)] 57.05 16.4 21931 3.488 0.0005
# contrast push-pull
q4_contrasts_PUSH_PULL <- emmeans::contrast(q4_ref_lme, method=list(
"(CG group: Effect of oxytocin vs. placebo on the difference between spouse-PUSH and spouse-PULL trials) vs. (NCG group: Effect of oxytocin vs. placebo on the difference between spouse-PUSH and spouse-PULL trials)" = (spouseOT_v_spousePL__CG_PUSH-spouseOT_v_spousePL__CG_PULL) - (spouseOT_v_spousePL__NCG_PUSH-spouseOT_v_spousePL__NCG_PULL)))
q4_contrasts_PUSH_PULL
contrast
(CG group: Effect of oxytocin vs. placebo on the difference between spouse-PUSH and spouse-PULL trials) vs. (NCG group: Effect of oxytocin vs. placebo on the difference between spouse-PUSH and spouse-PULL trials)
estimate SE df t.ratio p.value
52.8 23 21931 2.297 0.0217
# contrasts (push-pull by group) for plotting:
# these are just for plotting, don't care about whether the push-pull difference is significant
q4_contrasts_PUSH_PULL_plot <- emmeans::contrast(q4_ref_lme, method=list(
"NCG" = spouseOT_v_spousePL__NCG_PUSH-spouseOT_v_spousePL__NCG_PULL,
"CG" = spouseOT_v_spousePL__CG_PUSH-spouseOT_v_spousePL__CG_PULL))
q4_contrasts_PUSH_PULL_plot
contrast estimate SE df t.ratio p.value
NCG -23.6 15.2 21931 -1.556 0.1197
CG 29.3 17.3 21931 1.690 0.0910
# contrast A = CG group: \n Effect of oxytocin vs. placebo on the \n difference between spouse-PUSH \n and spouse-PULL trials
# contrast B = NCG group: \n Effect of oxytocin vs. placebo on the \n difference between spouse-PUSH \n and spouse-PULL trials
# horizontal
plot(q4_contrasts_PUSH_PULL_plot, comparisons = T, horizontal=F, xlab = "Estimated marginal means of the effect of oxytocin (vs. placebo) \n on the difference between spouse-PUSH and spouse-PULL trials", ylab="Avoidance bias Approach bias", intervals = FALSE, ylab()) +
geom_hline(yintercept=0, linetype="solid", color="black", size=.2) +
theme_bw() +
theme(axis.title.x = element_text(size = rel(1.2), margin = margin(t = 20, r = 0, b = 0, l = 0)),
axis.title.y = element_text(size = rel(1.)),
axis.text.x = element_text(size = rel(2)))
For the first section of the results/Table 1, showing the group comparison on a number of demographic and other variables.
Redid these after removing D149, who completed the gAAT outside of the scanner (was in the parent study but not gAAT study - recruited specifically for resting state after BA had already completed dissertation).
varstoadd <- varstoadd %>% filter(ID !="D149")
varstoadd %>% count(group) # 22 NCG, 17 CG in the whole sample
data_bias_placebo %>% filter(ID !="D149") %>% group_by(group) %>% count()/5 # double-checked and yep, that matches 22 NCG, 17 NCG with behavioral data
‘/’ not meaningful for factors
group n
1 NA 22
2 NA 17
#we had excluded D149 who completed the task outside of the scanner (NCG male recruited specifically for resting state)
# descriptives in the full sample (not broken down by group)
describe(varstoadd$age_yrs)
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 39 69.34 6.53 70.25 69.5 7.29 57.11 79.73 22.62 -0.27 -1.01 1.05
describe(varstoadd$timesincedeath/30.417) # in months
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 39 15.41 8.27 13.51 14.54 6.53 4.96 36 31.04 1.02 -0.02 1.32
describe(varstoadd$tot_icg)
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 39 23.38 12.63 22 22.79 10.38 4 51 47 0.36 -0.76 2.02
# descriptives by group w/t-tests
# age
describeBy(varstoadd$age_yrs, group=varstoadd$group)
Descriptive statistics by group
group: NCG
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 22 68.96 6.54 69.74 69.06 5.81 57.39 79.73 22.34 -0.18 -1.05 1.39
---------------------------------------------------------------------------------------------------
group: CG
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 17 69.83 6.68 71.33 70.02 7.29 57.11 79.72 22.61 -0.37 -1.16 1.62
t.test(varstoadd$age_yrs ~ varstoadd$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
Welch Two Sample t-test
data: varstoadd$age_yrs by varstoadd$group
t = -0.40662, df = 34.181, p-value = 0.6868
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-5.212653 3.474216
sample estimates:
mean in group NCG mean in group CG
68.96301 69.83223
# relationship length
describeBy(varstoadd$yrs_together, group=varstoadd$group)
Descriptive statistics by group
group: NCG
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 22 36.89 11.32 35 36.36 13.34 21 59 38 0.32 -1.11 2.41
---------------------------------------------------------------------------------------------------
group: CG
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 17 38.85 13.86 45 39.1 14.83 17 57 40 -0.28 -1.55 3.36
t.test(varstoadd$yrs_together ~ varstoadd$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
Welch Two Sample t-test
data: varstoadd$yrs_together by varstoadd$group
t = -0.47528, df = 30.569, p-value = 0.638
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-10.410269 6.477114
sample estimates:
mean in group NCG mean in group CG
36.88636 38.85294
# timesincedeath
describeBy((varstoadd$timesincedeath)/30.417, group=varstoadd$group)
Descriptive statistics by group
group: NCG
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 22 16.7 8.7 14.99 15.69 8.46 7.33 36 28.67 0.93 -0.48 1.86
---------------------------------------------------------------------------------------------------
group: CG
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 17 13.74 7.61 9.99 13.03 4.92 4.96 33.17 28.21 1.01 0.07 1.84
t.test(((varstoadd$timesincedeath)/30.417) ~ varstoadd$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
Welch Two Sample t-test
data: ((varstoadd$timesincedeath)/30.417) by varstoadd$group
t = 1.132, df = 36.377, p-value = 0.265
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-2.342493 8.265862
sample estimates:
mean in group NCG mean in group CG
16.69820 13.73651
# ICG
describeBy(varstoadd$tot_icg, group=varstoadd$group)
Descriptive statistics by group
group: NCG
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 22 14.27 6.36 16 14.33 8.15 4 24 20 -0.19 -1.47 1.36
---------------------------------------------------------------------------------------------------
group: CG
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 17 35.18 7.99 34 34.73 8.9 26 51 25 0.59 -1.09 1.94
t.test(varstoadd$tot_icg ~ varstoadd$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
Welch Two Sample t-test
data: varstoadd$tot_icg by varstoadd$group
t = -8.8371, df = 29.996, p-value = 7.508e-10
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-25.73467 -16.07282
sample estimates:
mean in group NCG mean in group CG
14.27273 35.17647
# BDI
describeBy(varstoadd$tot_bdi, group=varstoadd$group)
Descriptive statistics by group
group: NCG
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 22 6 4.86 5 5.44 2.97 0 19 19 1.01 0.38 1.04
---------------------------------------------------------------------------------------------------
group: CG
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 17 16.41 6.7 16 16.6 5.93 2 28 26 -0.1 -0.61 1.62
t.test(varstoadd$tot_bdi ~ varstoadd$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
Welch Two Sample t-test
data: varstoadd$tot_bdi by varstoadd$group
t = -5.4028, df = 28.116, p-value = 9.113e-06
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-14.358525 -6.465004
sample estimates:
mean in group NCG mean in group CG
6.00000 16.41176
# prescription meds
describeBy(varstoadd$meds_total, group=varstoadd$group)
Descriptive statistics by group
group: NCG
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 22 2.95 2.24 3 2.83 1.48 0 7 7 0.47 -0.79 0.48
---------------------------------------------------------------------------------------------------
group: CG
vars n mean sd median trimmed mad min max range skew kurtosis se
X1 1 17 2.76 2.49 2 2.67 2.97 0 7 7 0.41 -1.36 0.6
t.test(varstoadd$meds_total ~ varstoadd$group,
alternative = c("two.sided"),
mu = 0, paired = FALSE, var.equal = FALSE,
conf.level = 0.95)
Welch Two Sample t-test
data: varstoadd$meds_total by varstoadd$group
t = 0.24687, df = 32.54, p-value = 0.8066
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.375544 1.755223
sample estimates:
mean in group NCG mean in group CG
2.954545 2.764706
# chi-square tests
## sex
chisq.test(varstoadd$group, y = varstoadd$sex_m, correct = FALSE, simulate.p.value = TRUE)
Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)
data: varstoadd$group and varstoadd$sex_m
X-squared = 5.2901, df = NA, p-value = 0.02199
# use the simulated p value because if you don't, with small cell sizes, many of the expected values will be very small
# and therefore the approximations of p may not be correct (R will warn you about this)
varstoadd %>% group_by(group) %>%
count(sex_m) %>%
mutate(`(\\%)` = prop.table(n)*100)
## race
chisq.test(varstoadd$group, y = varstoadd$race, correct = FALSE, simulate.p.value = TRUE)
Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)
data: varstoadd$group and varstoadd$race
X-squared = 1.3282, df = NA, p-value = 0.4378
varstoadd %>% group_by(group) %>%
count(race) %>%
mutate(`(\\%)` = prop.table(n)*100)
## ethnicity
chisq.test(varstoadd$group, y = varstoadd$ethnicity_hisp, correct = FALSE, simulate.p.value = TRUE)
Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)
data: varstoadd$group and varstoadd$ethnicity_hisp
X-squared = 2.7281, df = NA, p-value = 0.1939
varstoadd %>% group_by(group) %>%
count(ethnicity_hisp) %>%
mutate(`(\\%)` = prop.table(n)*100)
## psychoactive meds
chisq.test(varstoadd$group, y = varstoadd$meds_psychoactive, correct = FALSE, simulate.p.value = FALSE)
Pearson's Chi-squared test
data: varstoadd$group and varstoadd$meds_psychoactive
X-squared = 0.28966, df = 1, p-value = 0.5904
varstoadd %>% group_by(group) %>%
count(meds_psychoactive) %>%
mutate(`(\\%)` = prop.table(n)*100)
# make variables for retired/non-retired; college grad. vs. < college grad.
varstoadd <- varstoadd %>% mutate(employment_retired = as.factor(ifelse(employment == "retired", 1, 0)))
varstoadd <- varstoadd %>% mutate(education_collegegrad = recode_factor(education, "less than high school" = "no", "high school grad" = "no", "some college" = "no", "college grad" = "college grad", "1 year grad school" = "college grad", "2 years grad school" = "college grad", "3 years grad school" = "college grad", "4+ years grad school" = "college grad", .ordered = TRUE))
table(varstoadd$education_collegegrad)
no college grad
16 23
## education (college grad vs. < college grad)
chisq.test(varstoadd$group, y = varstoadd$education_collegegrad, correct = FALSE, simulate.p.value = FALSE)
Pearson's Chi-squared test
data: varstoadd$group and varstoadd$education_collegegrad
X-squared = 0.00028336, df = 1, p-value = 0.9866
varstoadd %>% group_by(group) %>%
count(education_collegegrad) %>%
mutate(`(\\%)` = prop.table(n)*100)
## employment (retired vs. not)
chisq.test(varstoadd$group, y = varstoadd$employment_retired, correct = FALSE, simulate.p.value = TRUE)
Pearson's Chi-squared test with simulated p-value (based on 2000 replicates)
data: varstoadd$group and varstoadd$employment_retired
X-squared = 0.16819, df = NA, p-value = 0.7206
varstoadd %>% group_by(group) %>%
count(employment_retired) %>%
mutate(`(\\%)` = prop.table(n)*100)
I wondered whether the LME results differed from the ANOVA results due to better power or because the ANOVA uses medians (so estimated marginal means are means of the medians) whereas the estimated marginal means from the LME are means of means.
# load the data (trial-level data with outliers at 1st and 99th percentile RTs removed)
bx_long_noOut <- readRDS("~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/bx_long_noOut.rds")
# creating a wide dataset with the bias variables, which will then be converted back to a long dataset
# The arguments to spread():
# - data: Data object
# - key: Name of column containing the new column names
# - value: Name of column containing values
data_wide <- spread(bx_long_noOut, key=stim, value=gAAT_RT)
colnames(data_wide)
[1] "date" "time" "ID" "blockcode" "blocknum" "trialcode"
[7] "values.trialcode" "values.stimulus" "push_pull" "trialnum" "run" "visit"
[13] "cond_v1" "cond_v2" "cond" "outlier_RT" "neutral" "death"
[19] "living" "spouse" "stranger"
View(data_wide)
# subset by condition and direction
data_pushA <- data_wide %>% select(ID, cond, push_pull, visit, neutral, death, living, spouse, stranger) %>% filter(cond == "A" & push_pull == "PUSH") %>% group_by(ID) %>% mutate(neutralApush = mean(neutral, na.rm = TRUE), deathApush = mean(death, na.rm = TRUE), spouseApush = mean(spouse, na.rm=TRUE), livingApush = mean(living, na.rm=TRUE), strangerApush = mean(stranger, na.rm=TRUE))
data_pullA <- data_wide %>% select(ID, cond, push_pull, visit, neutral, death, living, spouse, stranger) %>% filter(cond == "A" & push_pull == "PULL") %>% group_by(ID) %>% mutate(neutralApull = mean(neutral, na.rm = TRUE), deathApull = mean(death, na.rm = TRUE), spouseApull = mean(spouse, na.rm=TRUE), livingApull = mean(living, na.rm=TRUE), strangerApull = mean(stranger, na.rm=TRUE))
data_pushB <- data_wide %>% select(ID, cond, push_pull, visit, neutral, death, living, spouse, stranger) %>% filter(cond == "B" & push_pull == "PUSH") %>% group_by(ID) %>% mutate(neutralBpush = mean(neutral, na.rm = TRUE), deathBpush = mean(death, na.rm = TRUE), spouseBpush = mean(spouse, na.rm=TRUE), livingBpush = mean(living, na.rm=TRUE), strangerBpush = mean(stranger, na.rm=TRUE))
data_pullB <- data_wide %>% select(ID, cond, push_pull, visit, neutral, death, living, spouse, stranger) %>% filter(cond == "B" & push_pull == "PULL") %>% group_by(ID) %>% mutate(neutralBpull = mean(neutral, na.rm = TRUE), deathBpull = mean(death, na.rm = TRUE), spouseBpull = mean(spouse, na.rm=TRUE), livingBpull = mean(living, na.rm=TRUE), strangerBpull = mean(stranger, na.rm=TRUE))
# summarize by ID
data_pushA1 <- data_pushA[!duplicated(data_pushA$ID), ] %>% select(-c(neutral, death, living, spouse, stranger, visit, cond, push_pull))
data_pushB1 <- data_pushB[!duplicated(data_pushB$ID), ] %>% select(-c(neutral, death, living, spouse, stranger, visit, cond, push_pull))
data_pullA1 <- data_pullA[!duplicated(data_pullA$ID), ] %>% select(-c(neutral, death, living, spouse, stranger, visit, cond, push_pull))
data_pullB1 <- data_pullB[!duplicated(data_pullB$ID), ] %>% select(-c(neutral, death, living, spouse, stranger, visit, cond, push_pull))
# merge the 4 subsets together
join1 <- left_join(data_pushA1, data_pushB1, by="ID")
join2 <- left_join(join1, data_pullA1, by = "ID")
bias <- left_join(join2, data_pullB1, by = "ID")
head(bias)
bias <- bias %>% mutate(bias_neu_A = neutralApush - neutralApull,
bias_dea_A = deathApush - deathApull,
bias_spo_A = spouseApush - spouseApull,
bias_str_A = strangerApush - strangerApull,
bias_liv_A = livingApush - livingApull,
bias_neu_B = neutralBpush - neutralBpull,
bias_dea_B = deathBpush - deathBpull,
bias_spo_B = spouseBpush - spouseBpull,
bias_str_B = strangerBpush - strangerBpull,
bias_liv_B = livingBpush - livingBpull
)
# drop the condition/stim specific variables
bias <- bias %>% select(-c(neutralApush, deathApush, spouseApush, livingApush, strangerApush, neutralBpush, deathBpush, spouseBpush, livingBpush, strangerBpush, neutralApull, deathApull, spouseApull, livingApull, strangerApull, neutralBpull, deathBpull, spouseBpull, livingBpull, strangerBpull))
colnames(bias)
[1] "ID" "bias_neu_A" "bias_dea_A" "bias_spo_A" "bias_str_A" "bias_liv_A" "bias_neu_B" "bias_dea_B" "bias_spo_B"
[10] "bias_str_B" "bias_liv_B"
# turn it into a long dataset
bias_long <- gather(bias,
key = "stim",
value = "bias", -c(ID))
head(bias_long)
# create new columns for stimulus category and condition
bias_long$cond <- as.factor(ifelse(grepl("*_A", bias_long$stim), "A", "B"))
bias_long$stimulus <- as.factor(ifelse(grepl("*neu*", bias_long$stim), "neutral",
ifelse(grepl("*dea*", bias_long$stim), "death",
ifelse(grepl("*spo*", bias_long$stim), "spouse",
ifelse(grepl("*str*", bias_long$stim), "stranger",
ifelse(grepl("*liv*", bias_long$stim), "living", NA))))))
bias_long <- bias_long %>% select(-stim) # no longer need the "stim" variable
#describe(bias_long$bias) # skewness is -.01, kurtosis is 1.34, M = 14.04, SD = 70.01
hist(bias_long$bias) # normal distribution
qqnorm(bias_long$bias) # definitely some outliers
# check out the outliers
## function from https://stackoverflow.com/questions/12866189/calculating-the-outliers-in-r
outfun <- function(x) {
abs(x-mean(x,na.rm=TRUE)) > 3*sd(x,na.rm=TRUE)
}
# add a variable for outlier = T/F
bias_long$outlier <- outfun(bias_long$bias)
table(bias_long$outlier) # 4 observations where outlier = TRUE (>3SD)
FALSE TRUE
389 1
# which are these?
bias_long %>% filter(outlier==TRUE)
### this is another way to do it using boxplot.stats
# outlier_values <- boxplot.stats(means_long$bias)$out
# boxplot(data_long$bias, main="Outliers", boxwex=0.3)
# mtext(paste("Outlying values: ", paste(round(outlier_values, digits=2), collapse=", ")), cex=0.9, side=1)
saveRDS(bias_long, "~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/bias_long_means.rds")
write_csv(bias_long, "~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/bias_long_means.csv")
data_bias <- readRDS("~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/bias_long_means.rds")
varstoadd <- readRDS("~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/varstoadd.rds")
# mean-center tot_icg and tot_bdi for use later
meancent <- function(x) { x - mean(x, na.rm=TRUE) } #simple worker function to mean center a variable
varstoadd <- varstoadd %>% mutate_at(vars("tot_icg", "tot_bdi"), funs(cm=meancent))
data_bias <- left_join(data_bias, varstoadd, by="ID") %>%
rename(treatment = cond,
gAAT_bias = bias) %>%
mutate(treatment = recode(treatment,
A = "placebo",
B = "oxytocin"),
stimulus = fct_relevel(stimulus, "neutral")) # make neutral the reference level
# A model in the full dataset (both sessions) only with two within-subjects factors (stimulus [5], treatment [2])
# and one between-subjects factor (group [2])
## DV: RT bias (push - pull)
## Predicted: Treatment x Group interaction (Aim 3), Stimulus x Treatment x Group (Aim 4)
# RM-ANOVA
q3 <- aov_ez(id = "ID", dv = "gAAT_bias", data_bias, between = c("group"), within = c("stimulus", "treatment"), factorize=T)
Contrasts set to contr.sum for the following variables: group
summary(q3)
HF eps > 1 treated as 1
Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
Sum Sq num Df Error SS den Df F value Pr(>F)
(Intercept) 63128 1 381107 37 6.1289 0.018001 *
group 2079 1 381107 37 0.2019 0.655844
stimulus 82982 4 493956 148 6.2158 0.000119 ***
group:stimulus 3210 4 493956 148 0.2404 0.915054
treatment 468 1 131120 37 0.1320 0.718442
group:treatment 24669 1 131120 37 6.9612 0.012114 *
stimulus:treatment 9030 4 368504 148 0.9066 0.461817
group:stimulus:treatment 9771 4 368504 148 0.9810 0.419843
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Mauchly Tests for Sphericity
Test statistic p-value
stimulus 0.75078 0.33886
group:stimulus 0.75078 0.33886
stimulus:treatment 0.77231 0.42403
group:stimulus:treatment 0.77231 0.42403
Greenhouse-Geisser and Huynh-Feldt Corrections
for Departure from Sphericity
GG eps Pr(>F[GG])
stimulus 0.88237 0.0002559 ***
group:stimulus 0.88237 0.8958272
stimulus:treatment 0.89372 0.4535478
group:stimulus:treatment 0.89372 0.4140424
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
HF eps Pr(>F[HF])
stimulus 0.9868407 0.0001295954
group:stimulus 0.9868407 0.9131148508
stimulus:treatment 1.0010882 0.4618169043
group:stimulus:treatment 1.0010882 0.4198428526
# Stimulus is predictive in the model, but only alone.
# There is a significant two-way interaction of group x treatment
# which looks like this:
emmip(q3, group ~ treatment)
emmeans(q3, ~ c(treatment|group), contr = "pairwise") #pairwise comparisons of treatment within group
$emmeans
group = NCG:
treatment emmean SE df lower.CL upper.CL
placebo 22.07 7.29 37 7.2921 36.9
oxytocin 8.24 8.52 37 -9.0282 25.5
group = CG:
treatment emmean SE df lower.CL upper.CL
placebo 1.38 8.30 37 -15.4367 18.2
oxytocin 19.62 9.70 37 -0.0219 39.3
Results are averaged over the levels of: stimulus
Confidence level used: 0.95
$contrasts
group = NCG:
contrast estimate SE df t.ratio p.value
placebo - oxytocin 13.8 8.03 37 1.723 0.0932
group = CG:
contrast estimate SE df t.ratio p.value
placebo - oxytocin -18.2 9.13 37 -1.998 0.0531
Results are averaged over the levels of: stimulus
# using ICG (mean-centered)
q3lme <- lme(gAAT_bias ~ stimulus*tot_icg_cm*treatment, random = ~ 1|ID, data_bias, method="REML")
summary(q3lme)
Linear mixed-effects model fit by REML
Data: data_bias
AIC BIC logLik
4220.377 4306.474 -2088.189
Random effects:
Formula: ~1 | ID
(Intercept) Residual
StdDev: 27.14454 54.32694
Fixed effects: gAAT_bias ~ stimulus * tot_icg_cm * treatment
Value Std.Error DF t-value p-value
(Intercept) 16.443918 9.724764 333 1.6909324 0.0918
stimulusdeath -19.933433 12.302676 333 -1.6202517 0.1061
stimulusliving 9.871666 12.302676 333 0.8023998 0.4229
stimulusspouse 18.328365 12.302676 333 1.4897868 0.1372
stimulusstranger -25.083510 12.302676 333 -2.0388661 0.0423
tot_icg_cm -0.298530 0.780062 37 -0.3827007 0.7041
treatmentoxytocin -13.611336 12.302676 333 -1.1063719 0.2694
stimulusdeath:tot_icg_cm -0.959373 0.986847 333 -0.9721601 0.3317
stimulusliving:tot_icg_cm -0.302443 0.986847 333 -0.3064743 0.7594
stimulusspouse:tot_icg_cm -0.869538 0.986847 333 -0.8811277 0.3789
stimulusstranger:tot_icg_cm -0.651229 0.986847 333 -0.6599091 0.5098
stimulusdeath:treatmentoxytocin 25.540577 17.398612 333 1.4679663 0.1431
stimulusliving:treatmentoxytocin 6.087726 17.398612 333 0.3498972 0.7266
stimulusspouse:treatmentoxytocin 16.202373 17.398612 333 0.9312452 0.3524
stimulusstranger:treatmentoxytocin 20.743229 17.398612 333 1.1922347 0.2340
tot_icg_cm:treatmentoxytocin 0.708320 0.986847 333 0.7177604 0.4734
stimulusdeath:tot_icg_cm:treatmentoxytocin 0.718906 1.395612 333 0.5151185 0.6068
stimulusliving:tot_icg_cm:treatmentoxytocin 1.705593 1.395612 333 1.2221107 0.2225
stimulusspouse:tot_icg_cm:treatmentoxytocin 1.073324 1.395612 333 0.7690700 0.4424
stimulusstranger:tot_icg_cm:treatmentoxytocin 0.014154 1.395612 333 0.0101420 0.9919
Correlation:
(Intr) stmlsd stmlsl stmlssp stmlsst tt_cg_ trtmnt stmlsd:__ stmlsl:__ stmlssp:__
stimulusdeath -0.633
stimulusliving -0.633 0.500
stimulusspouse -0.633 0.500 0.500
stimulusstranger -0.633 0.500 0.500 0.500
tot_icg_cm -0.003 0.002 0.002 0.002 0.002
treatmentoxytocin -0.633 0.500 0.500 0.500 0.500 0.002
stimulusdeath:tot_icg_cm 0.002 -0.003 -0.001 -0.001 -0.001 -0.633 -0.001
stimulusliving:tot_icg_cm 0.002 -0.001 -0.003 -0.001 -0.001 -0.633 -0.001 0.500
stimulusspouse:tot_icg_cm 0.002 -0.001 -0.001 -0.003 -0.001 -0.633 -0.001 0.500 0.500
stimulusstranger:tot_icg_cm 0.002 -0.001 -0.001 -0.001 -0.003 -0.633 -0.001 0.500 0.500 0.500
stimulusdeath:treatmentoxytocin 0.447 -0.707 -0.354 -0.354 -0.354 -0.001 -0.707 0.002 0.001 0.001
stimulusliving:treatmentoxytocin 0.447 -0.354 -0.707 -0.354 -0.354 -0.001 -0.707 0.001 0.002 0.001
stimulusspouse:treatmentoxytocin 0.447 -0.354 -0.354 -0.707 -0.354 -0.001 -0.707 0.001 0.001 0.002
stimulusstranger:treatmentoxytocin 0.447 -0.354 -0.354 -0.354 -0.707 -0.001 -0.707 0.001 0.001 0.001
tot_icg_cm:treatmentoxytocin 0.002 -0.001 -0.001 -0.001 -0.001 -0.633 -0.003 0.500 0.500 0.500
stimulusdeath:tot_icg_cm:treatmentoxytocin -0.001 0.002 0.001 0.001 0.001 0.447 0.002 -0.707 -0.354 -0.354
stimulusliving:tot_icg_cm:treatmentoxytocin -0.001 0.001 0.002 0.001 0.001 0.447 0.002 -0.354 -0.707 -0.354
stimulusspouse:tot_icg_cm:treatmentoxytocin -0.001 0.001 0.001 0.002 0.001 0.447 0.002 -0.354 -0.354 -0.707
stimulusstranger:tot_icg_cm:treatmentoxytocin -0.001 0.001 0.001 0.001 0.002 0.447 0.002 -0.354 -0.354 -0.354
stmlsst:__ stmlsd: stmlsl: stmlssp: stmlsst: tt_c_: stmlsd:__: stmlsl:__: stmlssp:__:
stimulusdeath
stimulusliving
stimulusspouse
stimulusstranger
tot_icg_cm
treatmentoxytocin
stimulusdeath:tot_icg_cm
stimulusliving:tot_icg_cm
stimulusspouse:tot_icg_cm
stimulusstranger:tot_icg_cm
stimulusdeath:treatmentoxytocin 0.001
stimulusliving:treatmentoxytocin 0.001 0.500
stimulusspouse:treatmentoxytocin 0.001 0.500 0.500
stimulusstranger:treatmentoxytocin 0.002 0.500 0.500 0.500
tot_icg_cm:treatmentoxytocin 0.500 0.002 0.002 0.002 0.002
stimulusdeath:tot_icg_cm:treatmentoxytocin -0.354 -0.003 -0.001 -0.001 -0.001 -0.707
stimulusliving:tot_icg_cm:treatmentoxytocin -0.354 -0.001 -0.003 -0.001 -0.001 -0.707 0.500
stimulusspouse:tot_icg_cm:treatmentoxytocin -0.354 -0.001 -0.001 -0.003 -0.001 -0.707 0.500 0.500
stimulusstranger:tot_icg_cm:treatmentoxytocin -0.707 -0.001 -0.001 -0.001 -0.003 -0.707 0.500 0.500 0.500
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-4.12394988 -0.63794412 -0.05353412 0.65825302 2.86524706
Number of Observations: 390
Number of Groups: 39
anova(q3lme) # yup, essentially same result as with medians dataset
numDF denDF F-value p-value
(Intercept) 1 333 6.512304 0.0112
stimulus 4 333 7.090156 <.0001
tot_icg_cm 1 37 0.131607 0.7188
treatment 1 333 0.000766 0.9779
stimulus:tot_icg_cm 4 333 1.003989 0.4055
stimulus:treatment 4 333 0.729390 0.5724
tot_icg_cm:treatment 1 333 10.217601 0.0015
stimulus:tot_icg_cm:treatment 4 333 0.541989 0.7050
# create a bar plot with error bars (95% CI)
p4 <- ggplot(data_bias, aes(group, gAAT_bias, fill=treatment)) +
stat_summary(geom = "bar", fun.y = mean, position = "dodge", color="black") +
stat_summary(geom = "errorbar", fun.data = mean_cl_normal, position=position_dodge(.9), width = 0.25, color="black") +
geom_hline(yintercept=0, linetype="solid", color="black", size=.5) +
xlab("Group") + ylab("Avoid bias Approach bias")
p4
# show as a boxplot:
p4b <- ggplot(data_bias, aes(group, gAAT_bias, fill=treatment)) +
geom_hline(yintercept=0, linetype="solid", color="black", size=.5) +
geom_boxplot(aes(fill=treatment), outlier.alpha = 0.3) +
xlab("Group") + ylab("Avoid bias Approach bias")
p4b
q4_ref <- emmeans(q3, ~stimulus+group+treatment, model = "multivariate")
q4_ref # this is the reference grid
stimulus group treatment emmean SE df lower.CL upper.CL
neutral NCG placebo 22.492 11.5 37 -0.795 45.8
death NCG placebo 8.212 15.0 37 -22.102 38.5
living NCG placebo 33.547 12.5 37 8.283 58.8
spouse NCG placebo 50.738 15.2 37 19.893 81.6
stranger NCG placebo -4.630 10.9 37 -26.742 17.5
neutral CG placebo 8.593 13.1 37 -17.898 35.1
death CG placebo -18.733 17.0 37 -53.218 15.8
living CG placebo 16.909 14.2 37 -11.831 45.6
spouse CG placebo 14.017 17.3 37 -21.072 49.1
stranger CG placebo -13.904 12.4 37 -39.058 11.3
neutral NCG oxytocin 3.709 12.0 37 -20.560 28.0
death NCG oxytocin 1.218 12.6 37 -24.291 26.7
living NCG oxytocin 5.863 13.8 37 -22.138 33.9
spouse NCG oxytocin 29.622 12.6 37 4.133 55.1
stranger NCG oxytocin 0.796 13.2 37 -26.037 27.6
neutral CG oxytocin 1.731 13.6 37 -25.877 29.3
death CG oxytocin 17.799 14.3 37 -11.219 46.8
living CG oxytocin 35.667 15.7 37 3.813 67.5
spouse CG oxytocin 47.430 14.3 37 18.434 76.4
stranger CG oxytocin -4.507 15.1 37 -35.032 26.0
Confidence level used: 0.95
# vector of contrast weights must line up with items as ordered in reference grid
spouseOT_v_spousePL_CG = c(
0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
0, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
-1, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
0, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
1, # spouse CG oxytocin
0) # stranger CG oxytocin
spouseOT_v_spousePL_NCG = c(
0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
-1, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
0, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
1, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
strangerOT_v_strangerPL_CG = c(
0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
0, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
0, # spouse CG placebo
-1, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
0, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
1) # stranger CG oxytocin
strangerOT_v_strangerPL_NCG = c(
0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
0, # spouse NCG placebo
-1, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
0, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
0, # spouse NCG oxytocin
1, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
deathOT_v_deathPL_CG = c(
0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
0, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
-1, # death CG placebo
0, # living CG placebo
0, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
0, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
1, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
deathOT_v_deathPL_NCG = c(
0, # neutral NCG placebo
-1, # death NCG placebo
0, # living NCG placebo
0, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
0, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
1, # death NCG oxytocin
0, # living NCG oxytocin
0, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
livingOT_v_livingPL_CG = c(
0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
0, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
-1, # living CG placebo
0, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
0, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
1, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
livingOT_v_livingPL_NCG = c(
0, # neutral NCG placebo
0, # death NCG placebo
-1, # living NCG placebo
0, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
0, # spouse CG placebo
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
1, # living NCG oxytocin
0, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
# test the contrasts
q4_contrasts <- emmeans::contrast(q4_ref, method=list("CG: spouse OT > spouse PL" = spouseOT_v_spousePL_CG,
"NCG: spouse OT > spouse PL" = spouseOT_v_spousePL_NCG))
q4_contrasts
contrast estimate SE df t.ratio p.value
CG: spouse OT > spouse PL 33.4 16.2 37 2.065 0.0460
NCG: spouse OT > spouse PL -21.1 14.2 37 -1.484 0.1462
library(tidyverse)
# removed the 1 participant who did not complete the gAAT in the scanner, per MFO (not in the "gAAT study", although he was in the parent study) so for the parent study, add 1 to the completers and all boxes flowing into that
# create an empty 100x100 grid
data <- tibble(x= 0:100, y= 0:100)
data %>%
ggplot(aes(x, y)) +
scale_x_continuous(minor_breaks = seq(10, 100, 10)) +
scale_y_continuous(minor_breaks = seq(10, 100, 10)) +
theme_linedraw() ->
p
p
# for everything that follows, you will use x/y, xmin/ymin, and xmax/ymax values to resize and arrange elements as you wish, using the grid coordinates as reference
# also, because R is weird about how it displays graphics (e.g., making the image bigger does NOT just scale everything up, including enlarging any text), you'll notice that some of the text overflows its box when displayed in RStudio
# that's because I went back and optimized the text size for when I save it out as high-res PNG file at the end
# add the first box and text
p +
geom_rect(xmin = 30, xmax=64, ymin=94, ymax=100, color='black',
fill='white', size=0.25) +
annotate('text', x=48, y=97,label= 'Prospective n = 163', size=3) ->
p
# then add some more boxes and text: those not screened and those excluded
p +
geom_rect(xmin = 30, xmax=64, ymin=73, ymax=80, color='black',
fill='white', size=0.25) +
annotate('text', x = 48, y=77,label= 'Screened n = 142', size=3) +
geom_rect(xmin = 70, xmax=100, ymin=85, ymax=100, color='black',
fill = 'white', size=0.25) +
annotate('text', x= 72, y=93, hjust=0, label= '21 not screened: \n 16 Unable to reach \n 5 No response after initial contact', size=3) +
geom_rect(xmin = 70, xmax=100, ymin=39, ymax=80, color='black',
fill='white', size=0.25) +
annotate('text', x= 72, y=60, hjust=0, label= '98 did not meet criteria:\n 24 MRI safety \n 11 Not interested \n 11 ICG score \n 11 Loss <6 or >36 months \n 11 Age <55 or >80 years \n 8 Claustrophobia \n 6 Medications \n 4 Medical condition \n 5 Other (e.g., out of town) \n 7 Unspecified (but did not\n complete full screen)', size=3) ->
p2
p2
# add a third box: n participants enrolled
p2 + geom_rect(xmin = 30, xmax=64, ymin=50, ymax=60, color='black',
fill='white', size=0.25) +
annotate('text', x= 48, y=55,label= 'Enrolled n = 44', size=3) -> p3
p3
# and some more boxes with n completers/n drops
p3 +
geom_rect(xmin = 30, xmax=64, ymin=20, ymax=30, color='black',
fill='white', size=0.25) +
annotate('text', x= 48, y=25,label= 'Completed n = 39', size=3) +
geom_rect(xmin = 70, xmax=100, ymin=10, ymax=30, color='black',
fill='white', size=0.25) +
annotate('text', x= 72, y=20, hjust=0, label= '5 dropped: \n 2 Undisclosed metal \n 1 Incidental finding \n 1 Back pain in MRI \n 1 Nausea', size=3) -> p4
p4
# finally, the last boxes: those who completed follow-up and those who did not
# p4 +
# geom_rect(xmin = 30, xmax=64, ymin=-2, ymax=6, color='black',
# fill='white', size=0.25) +
# annotate('text', x= 48, y=2,label= '33 Completed follow-up', size=3) +
# geom_rect(xmin = 70, xmax=100, ymin=-2, ymax=6, color='black',
# fill='white', size=0.25) +
# annotate('text', x= 72, y=2, hjust=0, label= '6 Lost to follow up \n1 Not assessed', size=3) -> p5
#
# p5
# if not including follow-up:
p5 <- p4
# now time to add the arrows
# in the following, x and y = the point on the plot where the arrow starts
# the point where the arrow finishes is defined by xend and yend
p5+
geom_segment(
# first vertical arrow
x=50, xend=50, y=94, yend=80.5,
size=0.15, linejoin = "mitre", lineend = "butt",
arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# first horizonal arrow
geom_segment(
x=50, xend=69.7, y=89, yend=89,
size=0.15, linejoin = "mitre", lineend = "butt",
arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# second vertical arrow
geom_segment(
x=50, xend=50, y=73, yend=60.5,
size=0.15, linejoin = "mitre", lineend = "butt",
arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# second horizonal arrow
geom_segment(
x=50, xend=69.7, y=67, yend=67,
size=0.15, linejoin = "mitre", lineend = "butt",
arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# third vertical arrow
geom_segment(
x=50, xend=50, y=50, yend=30.5,
size=0.15, linejoin = "mitre", lineend = "butt",
arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# diagonal arrow
geom_segment(
x=64, xend=69.7, y=40, yend=25,
size=0.15, linejoin = "mitre", lineend = "butt",
arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# short horizonal line, no arrow
geom_segment(
x=50, xend=64, y=40, yend=40,
size=0.15, linejoin = "mitre", lineend = "butt") +
# # fourth vertical arrow
# geom_segment(
# x=50, xend=50, y=20, yend=6.5,
# size=0.15, linejoin = "mitre", lineend = "butt",
# arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# 2nd diagonal arrow
# geom_segment(
# x=64, xend=69.7, y=12, yend=1,
# size=0.15, linejoin = "mitre", lineend = "butt",
# arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# # 2nd short horizonal line, no arrow
# geom_segment(
# x=50, xend=64, y=12, yend=12,
# size=0.15, linejoin = "mitre", lineend = "butt") +
# then get rid of the grid...
theme_void() +
# ...and add a title (centered)
ggtitle("CONSORT Diagram") +
theme(plot.title = element_text(hjust=.5)) -> p6
p6
# finally...save it as a high-resolution (dpi=600) PNG file
# width and height are in inches
# you'll have to play around with those (as well as the text size, potentially) to see which make your plot look best
ggsave("~/Desktop/consort.png", plot = p6, dpi=600, width = 8, height = 6)
icg <- readRDS("~/Dropbox/GLASS Lab/OT Study/data/master-dataset/ot-fmri_master-dataset_020719.rds")
icg <- subset(icg, select=c(ID, icg_1:icg_19)) %>%
filter(!ID %in% c("D149", "D147", "D142"))
icg <- subset(icg, select=c(icg_1:icg_19))
alpha(icg)
Reliability analysis
Call: alpha(x = icg)
raw_alpha std.alpha G6(smc) average_r S/N ase mean sd median_r
0.93 0.92 0.96 0.39 12 0.016 1.2 0.67 0.37
lower alpha upper 95% confidence boundaries
0.89 0.93 0.96
Reliability if an item is dropped:
raw_alpha std.alpha G6(smc) average_r S/N alpha se var.r med.r
icg_1 0.92 0.91 0.96 0.37 11 0.018 0.035 0.36
icg_2 0.92 0.92 0.96 0.39 11 0.017 0.038 0.37
icg_3 0.92 0.91 0.96 0.37 11 0.018 0.035 0.37
icg_4 0.92 0.92 0.96 0.38 11 0.017 0.040 0.37
icg_5 0.93 0.92 0.96 0.40 12 0.016 0.038 0.39
icg_6 0.92 0.92 0.96 0.38 11 0.017 0.037 0.36
icg_7 0.92 0.92 0.96 0.38 11 0.018 0.035 0.37
icg_8 0.92 0.92 0.96 0.38 11 0.017 0.037 0.36
icg_9 0.92 0.92 0.96 0.38 11 0.017 0.038 0.37
icg_10 0.92 0.92 0.96 0.38 11 0.017 0.040 0.37
icg_11 0.93 0.93 0.96 0.42 13 0.016 0.032 0.39
icg_12 0.93 0.93 0.96 0.41 12 0.016 0.037 0.39
icg_13 0.92 0.91 0.96 0.37 11 0.018 0.038 0.36
icg_14 0.93 0.92 0.96 0.41 12 0.016 0.036 0.39
icg_15 0.93 0.92 0.96 0.39 12 0.016 0.041 0.39
icg_16 0.92 0.92 0.96 0.38 11 0.017 0.040 0.36
icg_17 0.92 0.92 0.96 0.38 11 0.017 0.037 0.37
icg_18 0.92 0.92 0.96 0.39 12 0.016 0.041 0.38
icg_19 0.92 0.92 0.95 0.38 11 0.017 0.039 0.36
Item statistics
n raw.r std.r r.cor r.drop mean sd
icg_1 37 0.85 0.84 0.84 0.82 1.41 0.90
icg_2 37 0.65 0.63 0.62 0.60 1.76 0.95
icg_3 37 0.85 0.83 0.83 0.81 1.19 1.22
icg_4 37 0.71 0.71 0.70 0.66 2.49 1.15
icg_5 37 0.43 0.43 0.39 0.37 1.86 0.95
icg_6 37 0.75 0.73 0.73 0.70 1.46 1.22
icg_7 37 0.79 0.78 0.78 0.76 1.57 1.14
icg_8 37 0.75 0.75 0.75 0.72 1.27 1.04
icg_9 37 0.74 0.75 0.75 0.71 0.68 0.75
icg_10 37 0.68 0.70 0.69 0.64 0.86 0.95
icg_11 37 0.23 0.28 0.24 0.19 0.24 0.55
icg_12 37 0.36 0.38 0.33 0.29 0.70 0.94
icg_13 37 0.82 0.81 0.80 0.79 1.95 1.29
icg_14 37 0.36 0.39 0.37 0.31 0.59 0.80
icg_15 37 0.53 0.56 0.54 0.49 0.24 0.60
icg_16 37 0.72 0.71 0.70 0.67 0.78 1.25
icg_17 37 0.68 0.67 0.66 0.63 0.81 1.10
icg_18 37 0.60 0.59 0.58 0.53 0.95 1.10
icg_19 37 0.75 0.74 0.74 0.70 2.00 1.11
Non missing response frequency for each item
0 1 2 3 4 miss
icg_1 0.16 0.38 0.35 0.11 0.00 0
icg_2 0.08 0.32 0.38 0.19 0.03 0
icg_3 0.41 0.19 0.27 0.08 0.05 0
icg_4 0.08 0.08 0.30 0.35 0.19 0
icg_5 0.05 0.30 0.43 0.16 0.05 0
icg_6 0.27 0.27 0.24 0.16 0.05 0
icg_7 0.19 0.32 0.27 0.16 0.05 0
icg_8 0.27 0.32 0.30 0.08 0.03 0
icg_9 0.49 0.35 0.16 0.00 0.00 0
icg_10 0.43 0.35 0.14 0.08 0.00 0
icg_11 0.81 0.14 0.05 0.00 0.00 0
icg_12 0.57 0.22 0.16 0.05 0.00 0
icg_13 0.14 0.27 0.27 0.16 0.16 0
icg_14 0.59 0.22 0.19 0.00 0.00 0
icg_15 0.84 0.08 0.08 0.00 0.00 0
icg_16 0.65 0.11 0.11 0.08 0.05 0
icg_17 0.57 0.16 0.19 0.05 0.03 0
icg_18 0.49 0.22 0.16 0.14 0.00 0
icg_19 0.08 0.27 0.30 0.27 0.08 0
# does RT differ by sex?
t.test(gAAT_bias ~ sex_m, data=data_bias_placebo)
Welch Two Sample t-test
data: gAAT_bias by sex_m
t = 2.4111, df = 84.075, p-value = 0.01808
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
5.167684 53.808290
sample estimates:
mean in group female mean in group male
21.560714 -7.927273
t.test(gAAT_bias ~ sex_m, data=data_bias_ot)
Welch Two Sample t-test
data: gAAT_bias by sex_m
t = 2.7309, df = 122.76, p-value = 0.007246
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
7.345009 46.041355
sample estimates:
mean in group female mean in group male
22.375000 -4.318182
# yes, RT bias differs by sex
# RM-ANOVA
q3b <- aov_ez(id = "ID", dv = "gAAT_bias", data_bias, between = "group", within = c("stimulus", "treatment"), covariate = c("sex_m", "tot_bdi_cm"), factorize= FALSE, print.formula = TRUE)
Formula send to aov_car: gAAT_bias ~ group+sex_m+tot_bdi_cm + Error(ID/(stimulus * treatment))
Numerical variables NOT centered on 0 (i.e., likely bogus results): tot_bdi_cmContrasts set to contr.sum for the following variables: group, sex_m
anova(q3b)
Anova Table (Type 3 tests)
Response: gAAT_bias
num Df den Df MSE F ges Pr(>F)
group 1.0000 35.00 9748.1 0.0446 0.000343 0.83403
sex_m 1.0000 35.00 9748.1 4.0711 0.030400 0.05134 .
tot_bdi_cm 1.0000 35.00 9748.1 0.0015 0.000012 0.96934
stimulus 3.5865 125.53 3516.0 3.6099 0.034715 0.01044 *
group:stimulus 3.5865 125.53 3516.0 0.2876 0.002857 0.86709
sex_m:stimulus 3.5865 125.53 3516.0 2.9653 0.028694 0.02662 *
tot_bdi_cm:stimulus 3.5865 125.53 3516.0 1.1142 0.010979 0.35056
treatment 1.0000 35.00 3646.1 0.0000 0.000000 0.99943
group:treatment 1.0000 35.00 3646.1 2.0144 0.005769 0.16466
sex_m:treatment 1.0000 35.00 3646.1 0.2407 0.000693 0.62677
tot_bdi_cm:treatment 1.0000 35.00 3646.1 0.6706 0.001928 0.41839
stimulus:treatment 3.5646 124.76 2850.3 0.9712 0.007735 0.41917
group:stimulus:treatment 3.5646 124.76 2850.3 0.4737 0.003788 0.73356
sex_m:stimulus:treatment 3.5646 124.76 2850.3 0.9381 0.007474 0.43648
tot_bdi_cm:stimulus:treatment 3.5646 124.76 2850.3 0.3122 0.002500 0.84904
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# this looks at all interaction effects of covariates, which we don't want
# so run using lme? see below
# LME
# using continuous ICG
q3blme_c <- lme(gAAT_bias ~ stimulus*tot_icg_cm*treatment + sex_m + tot_bdi_cm, random = ~ 1|ID, data_bias, method="REML")
# using groups
q3blme <- lme(gAAT_bias ~ stimulus*group*treatment + sex_m + tot_bdi_cm, random = ~ 1|ID, data_bias, method="REML")
anova(q3blme)
numDF denDF F-value p-value
(Intercept) 1 333 6.894147 0.0090
stimulus 4 333 7.013386 <.0001
group 1 35 0.213292 0.6471
treatment 1 333 0.000758 0.9781
sex_m 1 35 4.094047 0.0507
tot_bdi_cm 1 35 0.001498 0.9693
stimulus:group 4 333 0.268928 0.8979
stimulus:treatment 4 333 0.721492 0.5777
group:treatment 1 333 8.267916 0.0043
stimulus:group:treatment 4 333 0.818662 0.5140
summary(q3blme)
Linear mixed-effects model fit by REML
Data: data_bias
AIC BIC logLik
4151.676 4245.47 -2051.838
Random effects:
Formula: ~1 | ID
(Intercept) Residual
StdDev: 26.00842 54.62346
Fixed effects: gAAT_bias ~ stimulus * group * treatment + sex_m + tot_bdi_cm
Value Std.Error DF t-value p-value
(Intercept) 13.85979 14.22509 333 0.9743203 0.3306
stimulusdeath -14.28030 16.46959 333 -0.8670703 0.3865
stimulusliving 11.05495 16.46959 333 0.6712339 0.5025
stimulusspouse 28.24617 16.46959 333 1.7150493 0.0873
stimulusstranger -27.12252 16.46959 333 -1.6468236 0.1005
groupCG -6.18870 22.12154 35 -0.2797588 0.7813
treatmentoxytocin -18.78342 16.46959 333 -1.1404906 0.2549
sex_m.L -17.07842 8.46430 35 -2.0176999 0.0513
tot_bdi_cm 0.03474 0.89756 35 0.0387085 0.9693
stimulusdeath:groupCG -13.04515 24.94542 333 -0.5229476 0.6014
stimulusliving:groupCG -2.73861 24.94542 333 -0.1097842 0.9126
stimulusspouse:groupCG -22.82165 24.94542 333 -0.9148637 0.3609
stimulusstranger:groupCG 4.62601 24.94542 333 0.1854452 0.8530
stimulusdeath:treatmentoxytocin 11.78932 23.29152 333 0.5061634 0.6131
stimulusliving:treatmentoxytocin -8.90076 23.29152 333 -0.3821457 0.7026
stimulusspouse:treatmentoxytocin -2.33273 23.29152 333 -0.1001535 0.9203
stimulusstranger:treatmentoxytocin 24.20961 23.29152 333 1.0394174 0.2994
groupCG:treatmentoxytocin 11.92161 24.94542 333 0.4779079 0.6330
stimulusdeath:groupCG:treatmentoxytocin 31.60410 35.27815 333 0.8958548 0.3710
stimulusliving:groupCG:treatmentoxytocin 34.52079 35.27815 333 0.9785318 0.3285
stimulusspouse:groupCG:treatmentoxytocin 42.60694 35.27815 333 1.2077431 0.2280
stimulusstranger:groupCG:treatmentoxytocin -7.95117 35.27815 333 -0.2253852 0.8218
Correlation:
(Intr) stmlsd stmlsl stmlssp stmlsst gropCG trtmnt sx_m.L tt_bd_ stmlsd:CG stmlsl:CG
stimulusdeath -0.579
stimulusliving -0.579 0.500
stimulusspouse -0.579 0.500 0.500
stimulusstranger -0.579 0.500 0.500 0.500
groupCG -0.710 0.372 0.372 0.372 0.372
treatmentoxytocin -0.579 0.500 0.500 0.500 0.500 0.372
sex_m.L 0.322 0.000 0.000 0.000 0.000 -0.205 0.000
tot_bdi_cm 0.291 0.000 0.000 0.000 0.000 -0.433 0.000 0.058
stimulusdeath:groupCG 0.382 -0.660 -0.330 -0.330 -0.330 -0.564 -0.330 0.000 0.000
stimulusliving:groupCG 0.382 -0.330 -0.660 -0.330 -0.330 -0.564 -0.330 0.000 0.000 0.500
stimulusspouse:groupCG 0.382 -0.330 -0.330 -0.660 -0.330 -0.564 -0.330 0.000 0.000 0.500 0.500
stimulusstranger:groupCG 0.382 -0.330 -0.330 -0.330 -0.660 -0.564 -0.330 0.000 0.000 0.500 0.500
stimulusdeath:treatmentoxytocin 0.409 -0.707 -0.354 -0.354 -0.354 -0.263 -0.707 0.000 0.000 0.467 0.233
stimulusliving:treatmentoxytocin 0.409 -0.354 -0.707 -0.354 -0.354 -0.263 -0.707 0.000 0.000 0.233 0.467
stimulusspouse:treatmentoxytocin 0.409 -0.354 -0.354 -0.707 -0.354 -0.263 -0.707 0.000 0.000 0.233 0.233
stimulusstranger:treatmentoxytocin 0.409 -0.354 -0.354 -0.354 -0.707 -0.263 -0.707 0.000 0.000 0.233 0.233
groupCG:treatmentoxytocin 0.382 -0.330 -0.330 -0.330 -0.330 -0.564 -0.660 0.000 0.000 0.500 0.500
stimulusdeath:groupCG:treatmentoxytocin -0.270 0.467 0.233 0.233 0.233 0.399 0.467 0.000 0.000 -0.707 -0.354
stimulusliving:groupCG:treatmentoxytocin -0.270 0.233 0.467 0.233 0.233 0.399 0.467 0.000 0.000 -0.354 -0.707
stimulusspouse:groupCG:treatmentoxytocin -0.270 0.233 0.233 0.467 0.233 0.399 0.467 0.000 0.000 -0.354 -0.354
stimulusstranger:groupCG:treatmentoxytocin -0.270 0.233 0.233 0.233 0.467 0.399 0.467 0.000 0.000 -0.354 -0.354
stmlssp:CG stmlsst:CG stmlsd: stmlsl: stmlssp: stmlsst: grpCG: stmlsd:CG: stmlsl:CG:
stimulusdeath
stimulusliving
stimulusspouse
stimulusstranger
groupCG
treatmentoxytocin
sex_m.L
tot_bdi_cm
stimulusdeath:groupCG
stimulusliving:groupCG
stimulusspouse:groupCG
stimulusstranger:groupCG 0.500
stimulusdeath:treatmentoxytocin 0.233 0.233
stimulusliving:treatmentoxytocin 0.233 0.233 0.500
stimulusspouse:treatmentoxytocin 0.467 0.233 0.500 0.500
stimulusstranger:treatmentoxytocin 0.233 0.467 0.500 0.500 0.500
groupCG:treatmentoxytocin 0.500 0.500 0.467 0.467 0.467 0.467
stimulusdeath:groupCG:treatmentoxytocin -0.354 -0.354 -0.660 -0.330 -0.330 -0.330 -0.707
stimulusliving:groupCG:treatmentoxytocin -0.354 -0.354 -0.330 -0.660 -0.330 -0.330 -0.707 0.500
stimulusspouse:groupCG:treatmentoxytocin -0.707 -0.354 -0.330 -0.330 -0.660 -0.330 -0.707 0.500 0.500
stimulusstranger:groupCG:treatmentoxytocin -0.354 -0.707 -0.330 -0.330 -0.330 -0.660 -0.707 0.500 0.500
stmlssp:CG:
stimulusdeath
stimulusliving
stimulusspouse
stimulusstranger
groupCG
treatmentoxytocin
sex_m.L
tot_bdi_cm
stimulusdeath:groupCG
stimulusliving:groupCG
stimulusspouse:groupCG
stimulusstranger:groupCG
stimulusdeath:treatmentoxytocin
stimulusliving:treatmentoxytocin
stimulusspouse:treatmentoxytocin
stimulusstranger:treatmentoxytocin
groupCG:treatmentoxytocin
stimulusdeath:groupCG:treatmentoxytocin
stimulusliving:groupCG:treatmentoxytocin
stimulusspouse:groupCG:treatmentoxytocin
stimulusstranger:groupCG:treatmentoxytocin 0.500
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-3.93650355 -0.66336641 -0.03033002 0.62162957 2.79026574
Number of Observations: 390
Number of Groups: 39
anova(q3blme_c)
numDF denDF F-value p-value
(Intercept) 1 333 6.930476 0.0089
stimulus 4 333 7.090157 <.0001
tot_icg_cm 1 35 0.140058 0.7105
treatment 1 333 0.000766 0.9779
sex_m 1 35 4.354869 0.0443
tot_bdi_cm 1 35 0.021043 0.8855
stimulus:tot_icg_cm 4 333 1.003989 0.4055
stimulus:treatment 4 333 0.729390 0.5724
tot_icg_cm:treatment 1 333 10.217602 0.0015
stimulus:tot_icg_cm:treatment 4 333 0.541989 0.7050
# Linear mixed model
q4c <- lme(gAAT_RT ~ spouse*group*treatment*direction + sex_m + tot_bdi_cm, data = data_long, random = ~ 1|ID, method="REML")
anova(q4c)
numDF denDF F-value p-value
(Intercept) 1 21955 3297.907 <.0001
spouse 1 21955 36.621 <.0001
group 1 35 1.905 0.1762
treatment 1 21955 43.259 <.0001
direction 1 21955 26.776 <.0001
sex_m 1 35 0.019 0.8904
tot_bdi_cm 1 35 1.143 0.2924
spouse:group 1 21955 3.556 0.0594
spouse:treatment 1 21955 2.291 0.1302
group:treatment 1 21955 41.572 <.0001
spouse:direction 1 21955 19.866 <.0001
group:direction 1 21955 0.850 0.3565
treatment:direction 1 21955 0.004 0.9492
spouse:group:treatment 1 21955 0.053 0.8188
spouse:group:direction 1 21955 0.254 0.6146
spouse:treatment:direction 1 21955 0.001 0.9731
group:treatment:direction 1 21955 10.092 0.0015
spouse:group:treatment:direction 1 21955 0.992 0.3193
summary(q4c)
Linear mixed-effects model fit by REML
Data: data_long
AIC BIC logLik
292695.3 292855.3 -146327.7
Random effects:
Formula: ~1 | ID
(Intercept) Residual
StdDev: 85.35971 186.5109
Fixed effects: gAAT_RT ~ spouse * group * treatment * direction + sex_m + tot_bdi_cm
Value Std.Error DF t-value p-value
(Intercept) 747.5275 24.82848 21955 30.107661 0.0000
spouse1 -9.7150 8.43950 21955 -1.151135 0.2497
groupCG 57.2625 40.03833 35 1.430191 0.1615
treatmentoxytocin 6.2354 5.28269 21955 1.180345 0.2379
directionPUSH 15.0140 5.28086 21955 2.843093 0.0045
sex_m.L -4.7215 23.23893 35 -0.203173 0.8402
tot_bdi_cm -2.6213 2.46427 35 -1.063728 0.2947
spouse1:groupCG 23.1553 12.78624 21955 1.810958 0.0702
spouse1:treatmentoxytocin 16.6537 11.88512 21955 1.401222 0.1612
groupCG:treatmentoxytocin 19.7543 7.96717 21955 2.479459 0.0132
spouse1:directionPUSH 36.8846 11.97134 21955 3.081075 0.0021
groupCG:directionPUSH -17.0575 7.98717 21955 -2.135609 0.0327
treatmentoxytocin:directionPUSH -12.1284 7.46160 21955 -1.625445 0.1041
spouse1:groupCG:treatmentoxytocin -15.5399 18.05494 21955 -0.860702 0.3894
spouse1:groupCG:directionPUSH -19.2068 18.13006 21955 -1.059392 0.2894
spouse1:treatmentoxytocin:directionPUSH -11.4344 16.90327 21955 -0.676464 0.4988
groupCG:treatmentoxytocin:directionPUSH 27.2716 11.29876 21955 2.413681 0.0158
spouse1:groupCG:treatmentoxytocin:directionPUSH 25.5533 25.65827 21955 0.995910 0.3193
Correlation:
(Intr) spous1 gropCG trtmnt drPUSH sx_m.L tt_bd_ sp1:CG sps1:t grpCG: s1:PUS gCG:PU
spouse1 -0.067
groupCG -0.779 0.042
treatmentoxytocin -0.107 0.314 0.066
directionPUSH -0.107 0.314 0.066 0.502
sex_m.L 0.506 0.000 -0.312 0.000 0.000
tot_bdi_cm 0.457 0.000 -0.657 0.000 0.000 0.058
spouse1:groupCG 0.044 -0.660 -0.062 -0.207 -0.207 0.000 0.000
spouse1:treatmentoxytocin 0.048 -0.710 -0.029 -0.445 -0.223 0.000 0.000 0.469
groupCG:treatmentoxytocin 0.071 -0.208 -0.099 -0.663 -0.333 0.000 0.000 0.310 0.295
spouse1:directionPUSH 0.047 -0.705 -0.029 -0.221 -0.441 0.000 0.000 0.465 0.501 0.147
groupCG:directionPUSH 0.071 -0.208 -0.099 -0.332 -0.661 0.000 0.000 0.309 0.147 0.497 0.292
treatmentoxytocin:directionPUSH 0.076 -0.222 -0.047 -0.708 -0.708 0.000 0.000 0.147 0.315 0.469 0.312 0.468
spouse1:groupCG:treatmentoxytocin -0.031 0.467 0.044 0.293 0.147 0.000 0.000 -0.708 -0.658 -0.441 -0.330 -0.219
spouse1:groupCG:directionPUSH -0.031 0.466 0.044 0.146 0.291 0.000 0.000 -0.705 -0.331 -0.219 -0.660 -0.441
spouse1:treatmentoxytocin:directionPUSH -0.033 0.499 0.021 0.313 0.312 0.000 0.000 -0.330 -0.703 -0.207 -0.708 -0.207
groupCG:treatmentoxytocin:directionPUSH -0.050 0.147 0.070 0.468 0.467 0.000 0.000 -0.219 -0.208 -0.705 -0.206 -0.707
spouse1:groupCG:treatmentoxytocin:directionPUSH 0.022 -0.329 -0.031 -0.206 -0.206 0.000 0.000 0.498 0.463 0.311 0.467 0.311
t:PUSH sp1:CG: s1:CG:P s1::PU gCG::P
spouse1
groupCG
treatmentoxytocin
directionPUSH
sex_m.L
tot_bdi_cm
spouse1:groupCG
spouse1:treatmentoxytocin
groupCG:treatmentoxytocin
spouse1:directionPUSH
groupCG:directionPUSH
treatmentoxytocin:directionPUSH
spouse1:groupCG:treatmentoxytocin -0.207
spouse1:groupCG:directionPUSH -0.206 0.499
spouse1:treatmentoxytocin:directionPUSH -0.441 0.463 0.468
groupCG:treatmentoxytocin:directionPUSH -0.660 0.311 0.311 0.292
spouse1:groupCG:treatmentoxytocin:directionPUSH 0.291 -0.704 -0.707 -0.659 -0.440
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-2.3278425 -0.6750483 -0.1943970 0.4703152 5.2966698
Number of Observations: 22008
Number of Groups: 39
q4c_ref <- emmeans(q4c, ~spouse+group+treatment+direction, model = "multivariate")
q4c_ref
spouse group treatment direction emmean SE df lower.CL upper.CL
0 NCG placebo PULL 747 25.1 35 696 798
1 NCG placebo PULL 737 25.9 35 685 790
0 CG placebo PULL 804 25.6 35 752 856
1 CG placebo PULL 818 26.7 35 763 872
0 NCG oxytocin PULL 753 25.0 35 702 804
1 NCG oxytocin PULL 760 25.9 35 708 813
0 CG oxytocin PULL 830 25.6 35 778 882
1 CG oxytocin PULL 845 26.7 35 791 899
0 NCG placebo PUSH 762 25.0 35 711 813
1 NCG placebo PUSH 789 25.9 35 737 842
0 CG placebo PUSH 802 25.6 35 750 854
1 CG placebo PUSH 833 26.7 35 779 888
0 NCG oxytocin PUSH 756 25.0 35 705 807
1 NCG oxytocin PUSH 789 25.9 35 736 841
0 CG oxytocin PUSH 843 25.6 35 791 895
1 CG oxytocin PUSH 890 26.7 35 835 944
Results are averaged over the levels of: sex_m
d.f. method: containment
Confidence level used: 0.95
spouseOT_v_spousePL__NCG_PULL = c(
0, # 0 NCG placebo PULL 761 18.4 38 724 798
-1, # 1 NCG placebo PULL 752 19.5 38 712 791
0, # 0 CG placebo PULL 789 20.9 37 747 831
0, # 1 CG placebo PULL 802 22.2 37 757 847
0, # 0 NCG oxytocin PULL 768 18.4 38 730 805
1, # 1 NCG oxytocin PULL 774 19.5 38 735 814
0, # 0 CG oxytocin PULL 815 20.9 37 773 857
0, # 1 CG oxytocin PULL 830 22.2 37 785 875
0, # 0 NCG placebo PUSH 776 18.4 38 739 813
0, # 1 NCG placebo PUSH 803 19.5 38 764 843
0, # 0 CG placebo PUSH 787 20.9 37 745 829
0, # 1 CG placebo PUSH 818 22.2 37 773 863
0, # 0 NCG oxytocin PUSH 770 18.4 38 733 808
0, # 1 NCG oxytocin PUSH 803 19.5 38 763 842
0, # 0 CG oxytocin PUSH 828 20.9 37 786 870
0) # 1 CG oxytocin PUSH 874 22.3 37 829 920
spouseOT_v_spousePL__NCG_PUSH = c(
0, # 0 NCG placebo PULL 761 18.4 38 724 798
0, # 1 NCG placebo PULL 752 19.5 38 712 791
0, # 0 CG placebo PULL 789 20.9 37 747 831
0, # 1 CG placebo PULL 802 22.2 37 757 847
0, # 0 NCG oxytocin PULL 768 18.4 38 730 805
0, # 1 NCG oxytocin PULL 774 19.5 38 735 814
0, # 0 CG oxytocin PULL 815 20.9 37 773 857
0, # 1 CG oxytocin PULL 830 22.2 37 785 875
0, # 0 NCG placebo PUSH 776 18.4 38 739 813
-1, # 1 NCG placebo PUSH 803 19.5 38 764 843
0, # 0 CG placebo PUSH 787 20.9 37 745 829
0, # 1 CG placebo PUSH 818 22.2 37 773 863
0, # 0 NCG oxytocin PUSH 770 18.4 38 733 808
1, # 1 NCG oxytocin PUSH 803 19.5 38 763 842
0, # 0 CG oxytocin PUSH 828 20.9 37 786 870
0) # 1 CG oxytocin PUSH 874 22.3 37 829 920
spouseOT_v_spousePL__CG_PULL = c(
0, # 0 NCG placebo PULL 761 18.4 38 724 798
0, # 1 NCG placebo PULL 752 19.5 38 712 791
0, # 0 CG placebo PULL 789 20.9 37 747 831
-1, # 1 CG placebo PULL 802 22.2 37 757 847
0, # 0 NCG oxytocin PULL 768 18.4 38 730 805
0, # 1 NCG oxytocin PULL 774 19.5 38 735 814
0, # 0 CG oxytocin PULL 815 20.9 37 773 857
1, # 1 CG oxytocin PULL 830 22.2 37 785 875
0, # 0 NCG placebo PUSH 776 18.4 38 739 813
0, # 1 NCG placebo PUSH 803 19.5 38 764 843
0, # 0 CG placebo PUSH 787 20.9 37 745 829
0, # 1 CG placebo PUSH 818 22.2 37 773 863
0, # 0 NCG oxytocin PUSH 770 18.4 38 733 808
0, # 1 NCG oxytocin PUSH 803 19.5 38 763 842
0, # 0 CG oxytocin PUSH 828 20.9 37 786 870
0) # 1 CG oxytocin PUSH 874 22.3 37 829 920
spouseOT_v_spousePL__CG_PUSH = c(
0, # 0 NCG placebo PULL 761 18.4 38 724 798
0, # 1 NCG placebo PULL 752 19.5 38 712 791
0, # 0 CG placebo PULL 789 20.9 37 747 831
0, # 1 CG placebo PULL 802 22.2 37 757 847
0, # 0 NCG oxytocin PULL 768 18.4 38 730 805
0, # 1 NCG oxytocin PULL 774 19.5 38 735 814
0, # 0 CG oxytocin PULL 815 20.9 37 773 857
0, # 1 CG oxytocin PULL 830 22.2 37 785 875
0, # 0 NCG placebo PUSH 776 18.4 38 739 813
0, # 1 NCG placebo PUSH 803 19.5 38 764 843
0, # 0 CG placebo PUSH 787 20.9 37 745 829
-1, # 1 CG placebo PUSH 818 22.2 37 773 863
0, # 0 NCG oxytocin PUSH 770 18.4 38 733 808
0, # 1 NCG oxytocin PUSH 803 19.5 38 763 842
0, # 0 CG oxytocin PUSH 828 20.9 37 786 870
1) # 1 CG oxytocin PUSH 874 22.3 37 829 920
# test the contrasts
q4c_contrasts_PULL <- emmeans::contrast(q4c_ref, method=list("OT > PL (spouse/NCG/PULL)" = spouseOT_v_spousePL__NCG_PULL,
"OT > PL (spouse/CG/PULL)" = spouseOT_v_spousePL__CG_PULL,
"[OT > PL (spouse/CG/PULL)] > [OT > PL (spouse/NCG/PULL)]" = spouseOT_v_spousePL__CG_PULL-spouseOT_v_spousePL__NCG_PULL))
q4c_contrasts_PUSH <- emmeans::contrast(q4c_ref, method=list("OT > PL (spouse/NCG/PUSH)" = spouseOT_v_spousePL__NCG_PUSH,
"OT > PL (spouse/CG/PUSH)" = spouseOT_v_spousePL__CG_PUSH,
"[OT > PL (spouse/CG/PUSH)] > [OT > PL (spouse/NCG/PUSH)]" = spouseOT_v_spousePL__CG_PUSH-spouseOT_v_spousePL__NCG_PUSH))
q4c_contrasts_PULL
contrast estimate SE df t.ratio p.value
OT > PL (spouse/NCG/PULL) 22.89 10.6 21955 2.150 0.0316
OT > PL (spouse/CG/PULL) 27.10 12.2 21955 2.219 0.0265
[OT > PL (spouse/CG/PULL)] > [OT > PL (spouse/NCG/PULL)] 4.21 16.2 21955 0.260 0.7948
Results are averaged over the levels of: sex_m
q4c_contrasts_PUSH
contrast estimate SE df t.ratio p.value
OT > PL (spouse/NCG/PUSH) -0.674 10.8 21955 -0.062 0.9503
OT > PL (spouse/CG/PUSH) 56.365 12.3 21955 4.580 <.0001
[OT > PL (spouse/CG/PUSH)] > [OT > PL (spouse/NCG/PUSH)] 57.039 16.4 21955 3.483 0.0005
Results are averaged over the levels of: sex_m
# contrast push-pull
q4c_contrasts_PUSH_PULL <- emmeans::contrast(q4c_ref, method=list(
"( [OT>PL (spouse/CG/PUSH)] - [OT>PL (spouse/CG/PULL)] ) > ( [OT>PL (spouse/NCG/PUSH)] - [OT>PL (spouse/NCG/PULL)] )" = (spouseOT_v_spousePL__CG_PUSH-spouseOT_v_spousePL__CG_PULL) - (spouseOT_v_spousePL__NCG_PUSH-spouseOT_v_spousePL__NCG_PULL)))
q4_contrasts_PUSH_PULL
contrast
(CG group: Effect of oxytocin vs. placebo on the difference between spouse-PUSH and spouse-PULL trials) vs. (NCG group: Effect of oxytocin vs. placebo on the difference between spouse-PUSH and spouse-PULL trials)
estimate SE df t.ratio p.value
52.8 23 21931 2.297 0.0217
# contrasts (push-pull by group) for plotting:
# these are just for plotting, don't care about whether the push-pull difference is significant
q4_contrasts_PUSH_PULL_plot <- emmeans::contrast(q4c_ref, method=list(
"( [OT>PL (spouse/CG/PUSH)] > \n [OT>PL (spouse/CG/PULL)] )" = spouseOT_v_spousePL__CG_PUSH-spouseOT_v_spousePL__CG_PULL,
"( [OT>PL (spouse/NCG/PUSH)] > \n [OT>PL (spouse/NCG/PULL)] )" = spouseOT_v_spousePL__NCG_PUSH-spouseOT_v_spousePL__NCG_PULL))
q4_contrasts_PUSH_PULL_plot
contrast estimate SE df t.ratio p.value
( [OT>PL (spouse/CG/PUSH)] > \n [OT>PL (spouse/CG/PULL)] ) 29.3 17.3 21955 1.688 0.0915
( [OT>PL (spouse/NCG/PUSH)] > \n [OT>PL (spouse/NCG/PULL)] ) -23.6 15.2 21955 -1.554 0.1203
Results are averaged over the levels of: sex_m
m_ref <- emmeans(q3blme, ~stimulus+group+treatment, model = "multivariate")
m_ref # this is the reference grid
stimulus group treatment emmean SE df lower.CL upper.CL
neutral NCG placebo 13.867 14.3 35 -15.13 42.9
death NCG placebo -0.413 14.3 35 -29.41 28.6
living NCG placebo 24.922 14.3 35 -4.07 53.9
spouse NCG placebo 42.113 14.3 35 13.12 71.1
stranger NCG placebo -13.255 14.3 35 -42.25 15.7
neutral CG placebo 7.679 15.6 35 -23.97 39.3
death CG placebo -19.647 15.6 35 -51.29 12.0
living CG placebo 15.995 15.6 35 -15.65 47.6
spouse CG placebo 13.103 15.6 35 -18.54 44.7
stranger CG placebo -14.818 15.6 35 -46.46 16.8
neutral NCG oxytocin -4.916 14.3 35 -33.91 24.1
death NCG oxytocin -7.407 14.3 35 -36.40 21.6
living NCG oxytocin -2.762 14.3 35 -31.76 26.2
spouse NCG oxytocin 20.997 14.3 35 -8.00 50.0
stranger NCG oxytocin -7.829 14.3 35 -36.82 21.2
neutral CG oxytocin 0.817 15.6 35 -30.83 32.5
death CG oxytocin 16.885 15.6 35 -14.76 48.5
living CG oxytocin 34.753 15.6 35 3.11 66.4
spouse CG oxytocin 46.515 15.6 35 14.87 78.2
stranger CG oxytocin -5.421 15.6 35 -37.07 26.2
Results are averaged over the levels of: sex_m
d.f. method: containment
Confidence level used: 0.95
# vector of contrast weights must line up with items as ordered in reference grid
# ref grid is the same, so can use same contrast vectors as for Q4
# test the contrasts
m_contrasts <- emmeans::contrast(m_ref, method=list("CG: spouse OT > spouse PL" = spouseOT_v_spousePL_CG,
"NCG: spouse OT > spouse PL" = spouseOT_v_spousePL_NCG))
m_contrasts
contrast estimate SE df t.ratio p.value
CG: spouse OT > spouse PL 33.4 18.7 333 1.783 0.0754
NCG: spouse OT > spouse PL -21.1 16.5 333 -1.282 0.2007
Results are averaged over the levels of: sex_m
q4c_ref_lme <- emmeans(q4c, ~spouse+group+treatment+direction, model = "multivariate")
q4c_ref_lme # this is the reference grid
spouse group treatment direction emmean SE df lower.CL upper.CL
0 NCG placebo PULL 747 25.1 35 696 798
1 NCG placebo PULL 737 25.9 35 685 790
0 CG placebo PULL 804 25.6 35 752 856
1 CG placebo PULL 818 26.7 35 763 872
0 NCG oxytocin PULL 753 25.0 35 702 804
1 NCG oxytocin PULL 760 25.9 35 708 813
0 CG oxytocin PULL 830 25.6 35 778 882
1 CG oxytocin PULL 845 26.7 35 791 899
0 NCG placebo PUSH 762 25.0 35 711 813
1 NCG placebo PUSH 789 25.9 35 737 842
0 CG placebo PUSH 802 25.6 35 750 854
1 CG placebo PUSH 833 26.7 35 779 888
0 NCG oxytocin PUSH 756 25.0 35 705 807
1 NCG oxytocin PUSH 789 25.9 35 736 841
0 CG oxytocin PUSH 843 25.6 35 791 895
1 CG oxytocin PUSH 890 26.7 35 835 944
Results are averaged over the levels of: sex_m
d.f. method: containment
Confidence level used: 0.95
# vector of contrast weights must line up with items as ordered in reference grid
spouseOT_v_spousePL__NCG_PULL = c(
0, # 0 NCG placebo PULL 761 18.4 38 724 798
-1, # 1 NCG placebo PULL 752 19.5 38 712 791
0, # 0 CG placebo PULL 789 20.9 37 747 831
0, # 1 CG placebo PULL 802 22.2 37 757 847
0, # 0 NCG oxytocin PULL 768 18.4 38 730 805
1, # 1 NCG oxytocin PULL 774 19.5 38 735 814
0, # 0 CG oxytocin PULL 815 20.9 37 773 857
0, # 1 CG oxytocin PULL 830 22.2 37 785 875
0, # 0 NCG placebo PUSH 776 18.4 38 739 813
0, # 1 NCG placebo PUSH 803 19.5 38 764 843
0, # 0 CG placebo PUSH 787 20.9 37 745 829
0, # 1 CG placebo PUSH 818 22.2 37 773 863
0, # 0 NCG oxytocin PUSH 770 18.4 38 733 808
0, # 1 NCG oxytocin PUSH 803 19.5 38 763 842
0, # 0 CG oxytocin PUSH 828 20.9 37 786 870
0) # 1 CG oxytocin PUSH 874 22.3 37 829 920
spouseOT_v_spousePL__NCG_PUSH = c(
0, # 0 NCG placebo PULL 761 18.4 38 724 798
0, # 1 NCG placebo PULL 752 19.5 38 712 791
0, # 0 CG placebo PULL 789 20.9 37 747 831
0, # 1 CG placebo PULL 802 22.2 37 757 847
0, # 0 NCG oxytocin PULL 768 18.4 38 730 805
0, # 1 NCG oxytocin PULL 774 19.5 38 735 814
0, # 0 CG oxytocin PULL 815 20.9 37 773 857
0, # 1 CG oxytocin PULL 830 22.2 37 785 875
0, # 0 NCG placebo PUSH 776 18.4 38 739 813
-1, # 1 NCG placebo PUSH 803 19.5 38 764 843
0, # 0 CG placebo PUSH 787 20.9 37 745 829
0, # 1 CG placebo PUSH 818 22.2 37 773 863
0, # 0 NCG oxytocin PUSH 770 18.4 38 733 808
1, # 1 NCG oxytocin PUSH 803 19.5 38 763 842
0, # 0 CG oxytocin PUSH 828 20.9 37 786 870
0) # 1 CG oxytocin PUSH 874 22.3 37 829 920
spouseOT_v_spousePL__CG_PULL = c(
0, # 0 NCG placebo PULL 761 18.4 38 724 798
0, # 1 NCG placebo PULL 752 19.5 38 712 791
0, # 0 CG placebo PULL 789 20.9 37 747 831
-1, # 1 CG placebo PULL 802 22.2 37 757 847
0, # 0 NCG oxytocin PULL 768 18.4 38 730 805
0, # 1 NCG oxytocin PULL 774 19.5 38 735 814
0, # 0 CG oxytocin PULL 815 20.9 37 773 857
1, # 1 CG oxytocin PULL 830 22.2 37 785 875
0, # 0 NCG placebo PUSH 776 18.4 38 739 813
0, # 1 NCG placebo PUSH 803 19.5 38 764 843
0, # 0 CG placebo PUSH 787 20.9 37 745 829
0, # 1 CG placebo PUSH 818 22.2 37 773 863
0, # 0 NCG oxytocin PUSH 770 18.4 38 733 808
0, # 1 NCG oxytocin PUSH 803 19.5 38 763 842
0, # 0 CG oxytocin PUSH 828 20.9 37 786 870
0) # 1 CG oxytocin PUSH 874 22.3 37 829 920
spouseOT_v_spousePL__CG_PUSH = c(
0, # 0 NCG placebo PULL 761 18.4 38 724 798
0, # 1 NCG placebo PULL 752 19.5 38 712 791
0, # 0 CG placebo PULL 789 20.9 37 747 831
0, # 1 CG placebo PULL 802 22.2 37 757 847
0, # 0 NCG oxytocin PULL 768 18.4 38 730 805
0, # 1 NCG oxytocin PULL 774 19.5 38 735 814
0, # 0 CG oxytocin PULL 815 20.9 37 773 857
0, # 1 CG oxytocin PULL 830 22.2 37 785 875
0, # 0 NCG placebo PUSH 776 18.4 38 739 813
0, # 1 NCG placebo PUSH 803 19.5 38 764 843
0, # 0 CG placebo PUSH 787 20.9 37 745 829
-1, # 1 CG placebo PUSH 818 22.2 37 773 863
0, # 0 NCG oxytocin PUSH 770 18.4 38 733 808
0, # 1 NCG oxytocin PUSH 803 19.5 38 763 842
0, # 0 CG oxytocin PUSH 828 20.9 37 786 870
1) # 1 CG oxytocin PUSH 874 22.3 37 829 920
# test the contrasts
q4c_contrasts_PULL <- emmeans::contrast(q4c_ref_lme, method=list("OT > PL (spouse/NCG/PULL)" = spouseOT_v_spousePL__NCG_PULL,
"OT > PL (spouse/CG/PULL)" = spouseOT_v_spousePL__CG_PULL,
"[OT > PL (spouse/CG/PULL)] > [OT > PL (spouse/NCG/PULL)]" = spouseOT_v_spousePL__CG_PULL-spouseOT_v_spousePL__NCG_PULL))
q4c_contrasts_PUSH <- emmeans::contrast(q4c_ref_lme, method=list("OT > PL (spouse/NCG/PUSH)" = spouseOT_v_spousePL__NCG_PUSH,
"OT > PL (spouse/CG/PUSH)" = spouseOT_v_spousePL__CG_PUSH,
"[OT > PL (spouse/CG/PUSH)] > [OT > PL (spouse/NCG/PUSH)]" = spouseOT_v_spousePL__CG_PUSH-spouseOT_v_spousePL__NCG_PUSH))
q4c_contrasts_PULL
contrast estimate SE df t.ratio p.value
OT > PL (spouse/NCG/PULL) 22.89 10.6 21955 2.150 0.0316
OT > PL (spouse/CG/PULL) 27.10 12.2 21955 2.219 0.0265
[OT > PL (spouse/CG/PULL)] > [OT > PL (spouse/NCG/PULL)] 4.21 16.2 21955 0.260 0.7948
Results are averaged over the levels of: sex_m
q4c_contrasts_PUSH
contrast estimate SE df t.ratio p.value
OT > PL (spouse/NCG/PUSH) -0.674 10.8 21955 -0.062 0.9503
OT > PL (spouse/CG/PUSH) 56.365 12.3 21955 4.580 <.0001
[OT > PL (spouse/CG/PUSH)] > [OT > PL (spouse/NCG/PUSH)] 57.039 16.4 21955 3.483 0.0005
Results are averaged over the levels of: sex_m
# contrast push-pull
q4c_contrasts_PUSH_PULL <- emmeans::contrast(q4c_ref_lme, method=list(
"( [OT>PL (spouse/CG/PUSH)] - [OT>PL (spouse/CG/PULL)] ) > ( [OT>PL (spouse/NCG/PUSH)] - [OT>PL (spouse/NCG/PULL)] )" = (spouseOT_v_spousePL__CG_PUSH-spouseOT_v_spousePL__CG_PULL) - (spouseOT_v_spousePL__NCG_PUSH-spouseOT_v_spousePL__NCG_PULL)))
q4c_contrasts_PUSH_PULL
contrast estimate SE
( [OT>PL (spouse/CG/PUSH)] - [OT>PL (spouse/CG/PULL)] ) > ( [OT>PL (spouse/NCG/PUSH)] - [OT>PL (spouse/NCG/PULL)] ) 52.8 23
df t.ratio p.value
21955 2.293 0.0219
Results are averaged over the levels of: sex_m
# contrasts (push-pull by group) for plotting:
# these are just for plotting, don't care about whether the push-pull difference is significant
q4c_contrasts_PUSH_PULL_plot <- emmeans::contrast(q4c_ref_lme, method=list(
"( [OT>PL (spouse/CG/PUSH)] > \n [OT>PL (spouse/CG/PULL)] )" = spouseOT_v_spousePL__CG_PUSH-spouseOT_v_spousePL__CG_PULL,
"( [OT>PL (spouse/NCG/PUSH)] > \n [OT>PL (spouse/NCG/PULL)] )" = spouseOT_v_spousePL__NCG_PUSH-spouseOT_v_spousePL__NCG_PULL))
q4c_contrasts_PUSH_PULL_plot
contrast estimate SE df t.ratio p.value
( [OT>PL (spouse/CG/PUSH)] > \n [OT>PL (spouse/CG/PULL)] ) 29.3 17.3 21955 1.688 0.0915
( [OT>PL (spouse/NCG/PUSH)] > \n [OT>PL (spouse/NCG/PULL)] ) -23.6 15.2 21955 -1.554 0.1203
Results are averaged over the levels of: sex_m
# horizontal
plot(q4c_contrasts_PUSH_PULL_plot, comparisons = T, horizontal=T, xlab="Avoidance bias Approach bias") +
geom_vline(xintercept=0, linetype="solid", color="black", size=.2)
library(tidyverse)
# create an empty 100x100 grid
dat <- tibble(x= 0:400, y= 0:400)
dat %>%
ggplot(aes(x, y)) +
scale_x_continuous(minor_breaks = seq(10, 400, 10)) +
scale_y_continuous(minor_breaks = seq(10, 400, 10)) +
theme_linedraw() +
geom_vline(xintercept=10, linetype="solid", color="black", size=.5) -> p
# for everything that follows, you will use x/y, xmin/ymin, and xmax/ymax values to resize and arrange elements as you wish, using the grid coordinates as reference
# also, because RStudio is weird about how it displays graphics (e.g., making the image bigger does NOT just scale everything up, including enlarging any text), you'll notice that some of the text overflows its box when displayed in RStudio
# that's because I went back and optimized the text size for when I save it out as high-res PNG file at the end
# add the first box and text
p <- p +
annotate('text', x=52, y=375,label= 'CG group', size=3) +
geom_rect(xmin = 12, xmax=51, ymin=383, ymax=400, color='black',
fill='white', size=0.25) +
annotate('text', x=32, y=393,label= 'PUSH SPOUSE trials', size=3) +
geom_rect(xmin = 58, xmax=98, ymin=383, ymax=400, color='black',
fill='white', size=0.25) +
annotate('text', x=78, y=393,label= 'PULL SPOUSE trials', size=3) +
annotate('text', x=152, y=375,label= 'NCG group', size=3) +
geom_rect(xmin = 112, xmax=151, ymin=383, ymax=400, color='black',
fill='white', size=0.25) +
annotate('text', x=132, y=393,label= 'PUSH SPOUSE trials', size=3) +
geom_rect(xmin = 158, xmax=198, ymin=383, ymax=400, color='black',
fill='white', size=0.25) +
annotate('text', x=178, y=393,label= 'PULL SPOUSE trials', size=3)
p
# then add some more boxes and text: those not screened and those excluded
p +
geom_rect(xmin = 30, xmax=64, ymin=73, ymax=80, color='black',
fill='white', size=0.25) +
annotate('text', x = 48, y=77,label= 'Screened n = 143', size=3) +
geom_rect(xmin = 70, xmax=100, ymin=85, ymax=100, color='black',
fill = 'white', size=0.25) +
annotate('text', x= 72, y=93, hjust=0, label= '21 not screened: \n 16 Unable to reach \n 5 No response after initial contact', size=3) +
geom_rect(xmin = 70, xmax=100, ymin=39, ymax=80, color='black',
fill='white', size=0.25) +
annotate('text', x= 72, y=60, hjust=0, label= '98 did not meet criteria:\n 24 MRI safety \n 11 Not interested \n 11 ICG score \n 11 Loss <6 or >36 months \n 11 Age <55 or >80 years \n 8 Claustrophobia \n 6 Medications \n 4 Medical condition \n 5 Other (e.g., out of town) \n 7 Unspecified (but did not\n complete full screen)', size=3) ->
p2
p2
# add a third box: n participants enrolled
p2 + geom_rect(xmin = 30, xmax=64, ymin=50, ymax=60, color='black',
fill='white', size=0.25) +
annotate('text', x= 48, y=55,label= 'Enrolled n = 45', size=3) -> p3
p3
# and some more boxes with n completers/n drops
p3 +
geom_rect(xmin = 30, xmax=64, ymin=20, ymax=30, color='black',
fill='white', size=0.25) +
annotate('text', x= 48, y=25,label= 'Completed n = 40', size=3) +
geom_rect(xmin = 70, xmax=100, ymin=10, ymax=30, color='black',
fill='white', size=0.25) +
annotate('text', x= 72, y=20, hjust=0, label= '5 dropped: \n 2 Undisclosed metal \n 1 Incidental finding \n 1 Back pain in MRI \n 1 Nausea', size=3) -> p4
p4
# finally, the last boxes: those who completed follow-up and those who did not
# p4 +
# geom_rect(xmin = 30, xmax=64, ymin=-2, ymax=6, color='black',
# fill='white', size=0.25) +
# annotate('text', x= 48, y=2,label= '33 Completed follow-up', size=3) +
# geom_rect(xmin = 70, xmax=100, ymin=-2, ymax=6, color='black',
# fill='white', size=0.25) +
# annotate('text', x= 72, y=2, hjust=0, label= '6 Lost to follow up \n1 Not assessed', size=3) -> p5
#
# p5
# if not including follow-up:
p5 <- p4
# now time to add the arrows
# in the following, x and y = the point on the plot where the arrow starts
# the point where the arrow finishes is defined by xend and yend
p5+
geom_segment(
# first vertical arrow
x=50, xend=50, y=94, yend=80.5,
size=0.15, linejoin = "mitre", lineend = "butt",
arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# first horizonal arrow
geom_segment(
x=50, xend=69.7, y=89, yend=89,
size=0.15, linejoin = "mitre", lineend = "butt",
arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# second vertical arrow
geom_segment(
x=50, xend=50, y=73, yend=60.5,
size=0.15, linejoin = "mitre", lineend = "butt",
arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# second horizonal arrow
geom_segment(
x=50, xend=69.7, y=67, yend=67,
size=0.15, linejoin = "mitre", lineend = "butt",
arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# third vertical arrow
geom_segment(
x=50, xend=50, y=50, yend=30.5,
size=0.15, linejoin = "mitre", lineend = "butt",
arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# diagonal arrow
geom_segment(
x=64, xend=69.7, y=40, yend=25,
size=0.15, linejoin = "mitre", lineend = "butt",
arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# short horizonal line, no arrow
geom_segment(
x=50, xend=64, y=40, yend=40,
size=0.15, linejoin = "mitre", lineend = "butt") +
# # fourth vertical arrow
# geom_segment(
# x=50, xend=50, y=20, yend=6.5,
# size=0.15, linejoin = "mitre", lineend = "butt",
# arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# 2nd diagonal arrow
# geom_segment(
# x=64, xend=69.7, y=12, yend=1,
# size=0.15, linejoin = "mitre", lineend = "butt",
# arrow = arrow(length = unit(1, "mm"), type= "closed")) +
# # 2nd short horizonal line, no arrow
# geom_segment(
# x=50, xend=64, y=12, yend=12,
# size=0.15, linejoin = "mitre", lineend = "butt") +
# then get rid of the grid...
theme_void() +
# ...and add a title (centered)
ggtitle("CONSORT Diagram") +
theme(plot.title = element_text(hjust=.5)) -> p6
p6
# finally...save it as a high-resolution (dpi=600) PNG file
# width and height are in inches
# you'll have to play around with those (as well as the text size, potentially) to see which make your plot look best
ggsave("~/Desktop/consort.png", plot = p6, dpi=600, width = 8, height = 6)