This project uses behavioral and self-report data from the study “The effect of intranasal oxytocin on neural functioning in widow(er)s” (#1503744420, PI: Mary-Frances O’Connor, Co-PI: Brian Arizmendi). It was a double-blind randomized crossover study using an intranasal oxytocin manipulation and approach-avoidance task variant (grief AAT; gAAT) to investigate conflicting accounts of motivated behavior in people with complicated grief. Enrolled participants were older adults who had experienced the death of a spouse or long-term romantic partner between 6-36 months prior to their participation. Stratified sampling was used to ensure a range of grief symptom severity scores was represented in the sample.
Participants attended two identical experimental sessions, at which they received one of two intranasal sprays (oxytocin or placebo) and then took part in an fMRI scan, during which time they completed the task.
In the gAAT, five stimulus categories were presented:
Photos provided by the participant were scanned,resized, and framed with blue/yellow frames to match the other stimuli used in the task. Photos of the stranger were sex-matched to the participant’s partner.
Participants completed the task using a joystick. The stimuli were animated so that a joystick push would make the image smaller (as if they were pushing it away) and a joystick pull would make the image larger (as if they were bringing it closer).
Participants were instructed to respond during the task to the color of a photo’s frame, e.g., “push the joystick when the frame is BLUE, pull when the frame is YELLOW”). Once participants had completed the first run of the task (144 trials), they completed a second run that was identical to the first except for the instructions, which were reversed (e.g. “now push the joystick for YELLOW and pull for BLUE”). Order of instructions was counterbalanced across participants. Stimuli were presented via Inquisit 4 (2014), in a pseudorandomized order determined by genetic algorithm to optimize statistical power and psychological validity (Wager & Nichols, 2003). Each trial last 3000ms, with 2500 ms allowed for response. Intertrial interval was 500ms.
MFO:
Should it be RT to one stimulus type, RT to the contrast of one stimuli type > another, or compatibility scores? Compatibility scores have felt the easiest for me to understand and compare (and they were what at least one of our reviewers was expecting). Is it true that using compatibility scores as a DV has the advantage of reducing the number of predictors in any of our models (as they remove the need for response direction if those are summed in the DV)?
BA:
I agree with this point, and do believe that converting push/pull to one compatibility score as the DV could work well to reduce the number of predictors. This might be good considering w the power issue.
SS:
Yes, it is true that if we use bias scores we do not include “direction” as a predictor.
On the other hand, if we use the mixed effects modeling approach, we have vastly more power because we then have ~288 observations per person per session (so ~576 each subject, if using data from both OT and PL sessions) rather than 10 observations (i.e., 5 stimuli x 2 directions, using medians to compute bias) per person/session. So having “direction” as a predictor is less of an issue.
Either way, I think using the bias scores for figures is helpful in depicting approach/avoid tendencies.
Pros of “direction” model: Being able to identify whether the change in bias is driven by faster/slower push or pull is potentially useful. More power than bias scores (more observations - unless we use the trial-level bias scores, which I’m not sure that we should.)
Cons of “direction” model: May be more difficult to interpret given the additional level in the interactions (potential 4-way).
Side note: I find the term “compatibility scores” confusing (vs. “bias scores”). I know that’s what Rinck & Becker used, but I think the term “compatibility” only makes sense when you have a sense of what people’s natural response tendencies are – in their paper, they had participants who were afraid of spiders, so a negative compatibility score (faster to pull than push) indicated a compatibility or congruency between their behavior on the task and real-world behaviors. I think the thing that’s always confused me about the use of the term in this project is, compatible with what? We need a strong a priori hypothesis about which tendency the group will show in order for this term to make sense (IMO)
MFO:
MLM seems unnecessarily complicated after all. I’m generally pro-regression.
BA:
I agree that MLM appears unnecessarily complicated. My thoughts right now are that an RM-ANOVA would be the best way to go because:
It appears to be what some others (maccallum) are doing, and the creator of the AAT has published results using this method. It is likely more digestible by reviewers as well.
The RM-ANOVA and the MLM should theoretically yield very similar results for the reasons described in the previous point
Regarding regression: From my understanding, we have basically already done this in MLM. The last line of the model of the MLM is a regression with all variables included… If we did a series of regressions like Eisma did, then we cannot directly compare stimuli types. If I’m understanding the Eisma paper correctly, that group ran a series of regressions, with one regression per stimulus type, and the DV was the median RT. So those conclusions CAN say, for instance: that rumination affects RT on this stimulus, and not on this one (for example) but they CANNOT say that rumination affects RT more or less on one stimulus type than on another stimulus type, because no direct comparison is made in a series of separate regressions.
SS:
Re: BA’s last point - right. I am not a fan of running 5 different regressions. We get an inflated chance of type I error (so would need to adjust our alpha accordingly), and as Brian said, we can’t directly compare stimulus effects. And as we saw from reviewer comment in green below, they’re likely to be skeptical of that approach. Regressions also don’t account for the fact that we have multiple observations that come from the same person.
My vote is that we use the linear mixed effects approach as our main approach (“mixed effects model” is just a different term for MLM – means we include both fixed [e.g., stimulus] and random [e.g., person] effects.) I think there’s a simpler way to do it than the MLM analyses Brian and I both did previous, which were pretty complicated. I’m calling it a mixed model here to distinguish it from those previous analyses but it’s just a different term for MLM.
Advantages: Gives us way more power!!! Therefore, we are able to test more complicated models, such as the three- and four-way interactions.
We do need to account for the random effect of subject. Regression does not allow us to account for the interdependence of the data. In a mixed model, we address the interdependence by nesting observations within person.
But I also think it makes sense to run the RM-ANOVA, and see how our results differ (or don’t), both because it may be easier to understand and because it’s consistent with what Maccallum did.
I propose a single 2 (Group) x 5 (Stimulus Type) [x 2 (Direction) if RT is DV] model - in the placebo data only - to address Aims 1 and 2. This could be LME or ANOVA.
Then for Aims 3 and 4, run the full model: 2 (Group) x 2 (Treatment) x 5 (Stimulus Type) [x 2 (Direction) if RT is DV]. This could be LME or ANOVA.
If I understand correctly, the authors ran a group (2) x response (2) x stimulus (2) x OT condition (2) analysis, for five of the possible ten stimulus pairs. It would seem more appropriate to run a single omnibus 2x2x5x2 analysis with all 5 stimuli.
SS:
I agree with this reviewer.
One result was that both grief groups showed an approach pattern for idiographic spouse images. This has some methodological implications for reaction time tests and is consistent with the reward component of grief. However, the theoretical contribution of the result is minimal.
The other primary result was that reaction time was slower in the CG group overall after OT administration. This was surprising in that there was no specific effect of stimulus x group x OT (a three way interaction) that would have been very interesting. The failure to obtain a three-way interaction is itself difficult to unpack. Because of the small sample, this study is vulnerable to a type II error. Indeed, this is the principal flaw here. It is difficult to know whether the effect is not present or merely wasn’t detected.
SS:
Using the trial-level data in MLM/mixed model solves the power problem.
Note that the experimental design for the present study is not identical to the studies below.
median RT to pull grief, positive, negative and abstract, and median RT to push grief, positive, negative and abstract stimuli. Compatibility effect scores were calculated by subtracting median RT’s in the pull condition from median RT’s in the push condition (Table 2).
A 2 (Group) x 2 (Prime) x 2 (Direction) x 4 (Stimulus Type) repeated measures ANOVA indicated a main effect for Stimulus Type [F(3,49) = 39.88, p < 001, η = 0.71], a significant 2- way interaction between Direction x Stimulus Type [F(3,49) = 6.29 p< .002, η = 0.27], and a significant 3-way interaction between Group x Direction x Stimulus Type, [F(3,49) = 3.79, , p < .017, η = 0.19].
Median RT’s for all stimuli were calculated for each stimulus type. To be able to test our hypotheses, we also determined push-pull scores for each stimulus type, which are calculated by deducting the Median RT’s for pull trials from the Median RT’s for the push trials.
Our main analyses consisted of multiple regression analyses, in which rumination levels were used as a predictor of reaction times (push-pull) for each stimulus type (deceased + loss word; deceased + neutral word; stranger + loss word; stranger + neutral word; puzzle + X’s).
In both papers, bias scores were computed as (median RTpush - median RTpull) within each stimulus category.
MFO:
These were the questions I thought the field would want to know, and the reviewer seemed to like.
Rephrased as aim:
Aim 3.
To identify whether effects of intranasal oxytocin differ in high-severity [CG group] versus low-severity [Non-CG group] participants.
Analysis
Look at the condition (2) x group (2) interaction in the following model:
Using bias scores:
Using push/pull RTs:
If we see the expected interaction, then look at whether there is a within-group difference between OT and PL sessions.
See below for analysis & output…
Rephrased as aim:
Aim 4. To identify whether the response to spouse images in the CG group significantly differs under intranasal oxytocin vs. placebo.
Look at the stimulus (5) xcondition (2) x group (2) interaction in the following model (same model for Aim 3 above):
Using bias scores:
Using push/pull RTs:
If we see the expected interaction, then look at planned contrasts. We hypothesize that:
spouse,placebo,CG != spouse,oxytocin,CG
and that
spouse,placebo,NCG == spouse,oxytocin,NCG
(i.e., spouse in the placebo session is not equal to spouse in the oxytocin session, in the CG group only).
####### Aims 3/4 model (ANOVA) ######
# A model with three within-subjects factors (stimulus [5], direction [2], treatment [2]) and one between-subjects factor (group [2])
## DV: RTs on push and pull trials separately
## Predicted: Treatment x Group interaction (Aim 3), Stimulus x Direction x Treatment x Group (Aim 4)
# RM-ANOVA
m3a <- aov_ez(id = "ID", dv = "gAAT_RT", data_long, between = "group", within = c("stimulus", "direction", "treatment"), fun_aggregate = median)
Contrasts set to contr.sum for the following variables: group
m3a
Anova Table (Type 3 tests)
Response: gAAT_RT
Effect df MSE F ges p.value
1 group 1, 37 144383.07 1.75 .03 .19
2 stimulus 3.64, 134.70 2601.07 9.40 *** .01 <.0001
3 group:stimulus 3.64, 134.70 2601.07 1.65 .002 .17
4 direction 1, 37 6375.61 5.24 * .004 .03
5 group:direction 1, 37 6375.61 1.32 .001 .26
6 treatment 1, 37 27850.76 1.97 .007 .17
7 group:treatment 1, 37 27850.76 0.81 .003 .38
8 stimulus:direction 3.44, 127.37 2391.65 8.63 *** .009 <.0001
9 group:stimulus:direction 3.44, 127.37 2391.65 0.71 .0007 .57
10 stimulus:treatment 3.79, 140.20 1656.19 1.09 .0009 .36
11 group:stimulus:treatment 3.79, 140.20 1656.19 0.90 .0007 .46
12 direction:treatment 1, 37 1843.57 0.37 <.0001 .55
13 group:direction:treatment 1, 37 1843.57 7.28 * .002 .01
14 stimulus:direction:treatment 3.30, 122.16 1980.63 0.52 .0004 .69
15 group:stimulus:direction:treatment 3.30, 122.16 1980.63 0.24 .0002 .88
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1
Sphericity correction method: GG
# Stimulus is predictive in the model, but only alone or in interaction with direction.
# There is a significant three-way interaction of group x direction x treatment
# which looks like this:
emmip(m3a, treatment ~direction*group)
emmeans(m3a, ~c(treatment| direction*group), contr = "pairwise", adjust="holm")
$emmeans
direction = PUSH, group = NCG:
treatment emmean SE df lower.CL upper.CL
placebo 748 20.0 37 707 788
oxytocin 747 23.2 37 700 794
direction = PULL, group = NCG:
treatment emmean SE df lower.CL upper.CL
placebo 721 18.7 37 683 759
oxytocin 734 18.8 37 696 772
direction = PUSH, group = CG:
treatment emmean SE df lower.CL upper.CL
placebo 758 22.7 37 712 804
oxytocin 796 26.4 37 742 849
direction = PULL, group = CG:
treatment emmean SE df lower.CL upper.CL
placebo 762 21.2 37 719 805
oxytocin 779 21.4 37 736 823
Results are averaged over the levels of: stimulus
Confidence level used: 0.95
$contrasts
direction = PUSH, group = NCG:
contrast estimate SE df t.ratio p.value
placebo - oxytocin 0.377 17.5 37 0.022 0.9829
direction = PULL, group = NCG:
contrast estimate SE df t.ratio p.value
placebo - oxytocin -12.605 15.3 37 -0.825 0.4147
direction = PUSH, group = CG:
contrast estimate SE df t.ratio p.value
placebo - oxytocin -37.976 19.9 37 -1.907 0.0643
direction = PULL, group = CG:
contrast estimate SE df t.ratio p.value
placebo - oxytocin -17.500 17.4 37 -1.007 0.3206
Results are averaged over the levels of: stimulus
Pairwise comparisons indicate that the CG group showed a marginally significant difference between OT vs. PL on PUSH trials, averaged across stimuli, t = -1.91, p = .064.
####### Aims 3/4 model (LME) ######
# A model with three within-subjects factors (stimulus [5], direction [2], treatment [2]) and one between-subjects factor (group [2])
## DV: RTs on push and pull trials separately
## Predicted: Treatment x Group interaction (Aim 3), Stimulus x Direction x Treatment x Group (Aim 4)
# Linear mixed model
m3b <- lme(gAAT_RT ~ stimulus*direction*treatment*group, data = data_long, random = ~ 1|ID, method="REML")
anova(m3b)
numDF denDF F-value p-value
(Intercept) 1 21931 3371.640 <.0001
stimulus 4 21931 19.191 <.0001
direction 1 21931 27.165 <.0001
treatment 1 21931 43.261 <.0001
group 1 37 1.961 0.1698
stimulus:direction 4 21931 8.848 <.0001
stimulus:treatment 4 21931 1.448 0.2151
direction:treatment 1 21931 0.016 0.8986
stimulus:group 4 21931 1.951 0.0990
direction:group 1 21931 0.751 0.3861
treatment:group 1 21931 41.817 <.0001
stimulus:direction:treatment 4 21931 0.881 0.4741
stimulus:direction:group 4 21931 0.390 0.8163
stimulus:treatment:group 4 21931 0.338 0.8525
direction:treatment:group 1 21931 10.026 0.0015
stimulus:direction:treatment:group 4 21931 0.837 0.5017
# variance explained by the entire model
fitted <- fitted(m3b)
data_long$fitted_m3b <- as.vector(fitted)
forR2 <- lm(gAAT_RT ~ fitted_m3b, data=data_long)
summary(forR2)
Call:
lm(formula = gAAT_RT ~ fitted_m3b, data = data_long)
Residuals:
Min 1Q Median 3Q Max
-442.10 -125.60 -36.63 88.27 988.45
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.19743 11.59074 -0.535 0.593
fitted_m3b 1.00786 0.01461 68.975 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 186 on 22006 degrees of freedom
Multiple R-squared: 0.1778, Adjusted R-squared: 0.1777
F-statistic: 4758 on 1 and 22006 DF, p-value: < 2.2e-16
Like the ANOVA, the linear mixed effects model shows a significant three-way direction x treatment x group interaction. The model explains about 18% of the variance in the data (Adj. Rtreatment x group interaction (not interpreted in the presence of the higher-order interaction.)
Next, looking at the model summary tells us about some contrasts (but these are not the ones that we are really interested in so I won’t print the output table). Note that the model output includes tests of the null hypotheses that these differences are equal to zero.
summary(m3b)
# we can look at ALL pairwise comparisons using the emmeans package again, but this doesn't take into account the random effects?
emm.m3b1 <- emmeans(m3b, pairwise ~c(treatment | direction*group), adjust="holm") # emm.m3b1[["emmeans"]]@linfct gives us the contrast matrices for each pairwise comparison, this seems useful
NOTE: Results may be misleading due to involvement in interactions
emm.m3b1
$emmeans
direction = PULL, group = NCG:
treatment emmean SE df lower.CL upper.CL
placebo 759 18.3 38 722 796
oxytocin 769 18.3 38 732 806
direction = PUSH, group = NCG:
treatment emmean SE df lower.CL upper.CL
placebo 781 18.3 38 744 818
oxytocin 776 18.3 38 739 814
direction = PULL, group = CG:
treatment emmean SE df lower.CL upper.CL
placebo 791 20.8 37 749 833
oxytocin 818 20.8 37 775 860
direction = PUSH, group = CG:
treatment emmean SE df lower.CL upper.CL
placebo 793 20.8 37 751 835
oxytocin 837 20.8 37 795 879
Results are averaged over the levels of: stimulus
d.f. method: containment
Confidence level used: 0.95
$contrasts
direction = PULL, group = NCG:
contrast estimate SE df t.ratio p.value
placebo - oxytocin -9.95 4.74 21931 -2.099 0.0358
direction = PUSH, group = NCG:
contrast estimate SE df t.ratio p.value
placebo - oxytocin 4.62 4.75 21931 0.974 0.3301
direction = PULL, group = CG:
contrast estimate SE df t.ratio p.value
placebo - oxytocin -26.44 5.38 21931 -4.918 <.0001
direction = PUSH, group = CG:
contrast estimate SE df t.ratio p.value
placebo - oxytocin -44.47 5.43 21931 -8.196 <.0001
Results are averaged over the levels of: stimulus
emmip(m3b, group*direction ~treatment)
NOTE: Results may be misleading due to involvement in interactions
Pairwise comparisons show that while the CG group is significantly slower to both PUSH and PULL under OT, the NCG group also gets a bit slower to PULL (but not PUSH) under OT.
Although there is no significant interaction with stimulus, it’s interesting to visualize what it looks like:
emmip(m3b, group*direction ~treatment|stimulus) # lmem
Here we can clearly see that the CG group’s RTs are longer across stimulus types in the OT condition.
####### Aims 3/4 model (ANOVA w/bias scores) ######
# A model with two within-subjects factors (stimulus [5], treatment [2]) and one between-subjects factor (group [2])
## DV: bias scores
## Predicted: Treatment x Group interaction (Aim 3), Stimulus x Treatment x Group (Aim 4)
# RM-ANOVA
m3c <- aov_ez(id = "ID", dv = "gAAT_bias", data_bias, between = "group", within = c("stimulus", "treatment"), fun_aggregate = median)
Contrasts set to contr.sum for the following variables: group
m3c
Anova Table (Type 3 tests)
Response: gAAT_bias
Effect df MSE F ges p.value
1 group 1, 37 12751.21 1.32 .010 .26
2 stimulus 3.44, 127.37 4783.31 8.63 *** .08 <.0001
3 group:stimulus 3.44, 127.37 4783.31 0.71 .007 .57
4 treatment 1, 37 3687.15 0.37 .0008 .55
5 group:treatment 1, 37 3687.15 7.28 * .02 .01
6 stimulus:treatment 3.30, 122.16 3961.26 0.52 .004 .69
7 group:stimulus:treatment 3.30, 122.16 3961.26 0.24 .002 .88
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘+’ 0.1 ‘ ’ 1
Sphericity correction method: GG
# 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(m3c, group ~ treatment)
emmeans(m3c, ~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(m3c, ~c(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
death -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 - death 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
death - living -34.24 11.16 37 -3.067 0.0241
death - spouse -50.30 13.12 37 -3.834 0.0038
death - 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.
####### Aims 3/4 model (LME w/bias scores) ######
# A model with two within-subjects factors (stimulus [5], treatment [2]) and one between-subjects factor (group [2])
## DV: bias scores
## Predicted: Treatment x Group interaction (Aim 3), Stimulus x Treatment x Group (Aim 4)
# Linear mixed model
m3d <- lme(gAAT_bias ~ stimulus*treatment*group, data = data_bias, random = ~ 1|ID, method="REML")
anova(m3d)
numDF denDF F-value p-value
(Intercept) 1 333 6.033213 0.0145
stimulus 4 333 9.481681 <.0001
treatment 1 333 0.067814 0.7947
group 1 37 1.319279 0.2581
stimulus:treatment 4 333 0.423645 0.7916
stimulus:group 4 333 0.791315 0.5315
treatment:group 1 333 7.268341 0.0074
stimulus:treatment:group 4 333 0.216758 0.9290
# variance explained by the entire model
fitted <- fitted(m3d)
data_bias$fitted_m3d <- as.vector(fitted)
forR2 <- lm(gAAT_bias ~ fitted_m3d, data=data_bias)
summary(forR2)
Call:
lm(formula = gAAT_bias ~ fitted_m3d, data = data_bias)
Residuals:
Min 1Q Median 3Q Max
-221.580 -35.774 3.186 33.975 210.900
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.07567 3.11099 -0.989 0.323
fitted_m3d 1.21899 0.08519 14.310 <2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 56.71 on 388 degrees of freedom
Multiple R-squared: 0.3454, Adjusted R-squared: 0.3438
F-statistic: 204.8 on 1 and 388 DF, p-value: < 2.2e-16
As above, we see a main effect of stimulus and a two-way interaction of treatment x group. The model explains about 34% of the variance (adj R2 = .34).
This is what the interaction looks like:
p3 <- ggplot(data_bias, aes(fill=group, y=gAAT_bias, x=treatment)) +
geom_boxplot(aes(fill=group), outlier.alpha = 0.3)
p3 <- p3 + geom_hline(yintercept=0, linetype="solid", color="black", size=.5) + xlab("Condition") + ylab("Avoid bias Approach bias")
p4 <- ggplot(data_bias, aes(fill=treatment, y=gAAT_bias, x=group)) +
geom_boxplot(aes(fill=treatment), outlier.alpha = 0.3)
p4 <- p4 + geom_hline(yintercept=0, linetype="solid", color="black", size=.5) + xlab("Group") + ylab("Avoid bias Approach bias")
p3 # treatment on x-axis
p4 # group on x-axis
emm_m3c <- emmeans(m3c, ~ treatment|group) # compare treatment within group (ANOVA)
plot(emm_m3c, comparisons = TRUE, intervals = TRUE, horizontal = FALSE)
emm_m3d <- emmeans(m3d, ~ treatment|group) # compare treatment within group (LME)
NOTE: Results may be misleading due to involvement in interactions
plot(emm_m3d, comparisons = TRUE, intervals = TRUE, horizontal = FALSE)
https://aosmith.rbind.io/2019/04/15/custom-contrasts-emmeans/ https://cran.r-project.org/web/packages/afex/vignettes/afex_anova_example.html#running-custom-contrasts
Objects returned from emmeans can also be used to test specific contrasts. For this, we can simply create a list, where each element corresponds to one contrasts. A contrast is defined as a vector of constants on the reference grid (i.e., the object returned from emmeans, here m3). For example, we might be interested in whether there is a difference between the valid and invalid inferences in each of the two conditions.
This section in progress, no results as of yet…
m3d_ref <- emmeans(m3d, ~stimulus+group+treatment, model = "multivariate")
m3d_ref # this is the reference grid
# vector of contrast weights must line up with items as ordered in reference grid
m3a_ref <- emmeans(m3a, ~stimulus+group+treatment, model = "multivariate")
m3a_ref
m3b_ref <- emmeans(m3b, ~stimulus+group+treatment, model = "multivariate")
m3b_ref
m3c_ref <- emmeans(m3c, ~stimulus+group+treatment, model = "multivariate")
m3c_ref
spouse_pl = c(0, # neutral NCG placebo
0, # death NCG placebo
0, # living NCG placebo
.5, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
0, # death CG placebo
0, # living CG placebo
.5, # 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
0, # spouse CG oxytocin
0) # stranger CG oxytocin
death_pl = c(0, # neutral NCG placebo
.5, # death NCG placebo
0, # living NCG placebo
0, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
.5, # 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
0, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
spouse_ot = 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
0, # stranger CG placebo
0, # neutral NCG oxytocin
0, # death NCG oxytocin
0, # living NCG oxytocin
.5, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
0, # death CG oxytocin
0, # living CG oxytocin
.5, # spouse CG oxytocin
0) # stranger CG oxytocin
death_ot = 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
0, # stranger CG placebo
0, # neutral NCG oxytocin
.5, # death NCG oxytocin
0, # living NCG oxytocin
0, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
.5, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
spouse_ncg_pl = 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
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
spouse_ncg_ot = 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
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
spouse_cg_pl = 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
0, # spouse CG oxytocin
0) # stranger CG oxytocin
spouse_cg_ot = 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
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
###
living_ncg_pl = 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
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
living_ncg_ot = 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
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
living_cg_pl = 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
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
living_cg_ot = 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
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
###
death_ncg_pl = 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
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
0) # stranger CG oxytocin
death_ncg_ot = 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
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
death_cg_pl = 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
0, # death CG oxytocin
0, # living CG oxytocin
0, # spouse CG oxytocin
0) # stranger CG oxytocin
death_cg_ot = 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
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
# test the contrasts
contrast(m3d_ref, method = list("spouse NCG > death NCG (placebo)" = spouse_ncg_pl - death_ncg_pl,
"spouse CG > death CG (placebo)" = spouse_cg_pl - death_cg_pl,
"spouse NCG > death NCG (oxytocin)" = spouse_ncg_ot - death_ncg_ot,
"spouse CG > death CG (oxytocin)" = spouse_cg_ot - death_cg_ot,
"spouse > death (placebo)" = spouse_pl - death_pl,
"spouse > death (oxytocin)" = spouse_ot - death_ot), adjust="bonf")
# this is another way to do it:
ppl_vs_scenes_ot = 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
0, # stranger CG placebo
-1/4, # neutral NCG oxytocin
-1/4, # death NCG oxytocin
1/6, # living NCG oxytocin
1/6, # spouse NCG oxytocin
1/6, # stranger NCG oxytocin
-1/4, # neutral CG oxytocin
-1/4, # death CG oxytocin
1/6, # living CG oxytocin
1/6, # spouse CG oxytocin
1/6) # stranger CG oxytocin
contrast(ref, method=list(ppl_vs_scenes_ot))
#### For Aim 4 ####
# tests OT > PL for spouse stimuli within the CG group
spouse_cg_OTvPL = 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
# tests OT > PL for spouse stimuli within the NCG group
spouse_ncg_OTvPL = 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
# tests OT > PL for living stimuli within the CG group
living_cg_OTvPL = 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
# tests OT > PL for living stimuli within the NCG group
living_ncg_OTvPL = 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
# tests OT > PL for death stimuli within the CG group
death_cg_OTvPL = 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
# tests OT > PL for death stimuli within the NCG group
death_ncg_OTvPL = 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
# tests OT > PL for stranger stimuli within the CG group
stranger_cg_OTvPL = 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
# tests OT > PL for stranger stimuli within the NCG group
stranger_ncg_OTvPL = 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
contrast(m3a_ref, method = list("Effect of OT on spouse RT in CG > effect of OT on spouse RT in NCG" = spouse_cg_OTvPL - spouse_ncg_OTvPL))
contrast(m3a_ref, method = list("Effect of OT on living loved one RT in CG > effect of OT on living loved one RT in NCG" = living_cg_OTvPL - living_ncg_OTvPL))
contrast(m3a_ref, method = list("Effect of OT on generic death RT in CG > effect of OT on generic death RT in NCG" = death_cg_OTvPL - death_ncg_OTvPL))
contrast(m3a_ref, method = list("Effect of OT on stranger RT in CG > effect of OT on stranger RT in NCG" = stranger_cg_OTvPL - stranger_ncg_OTvPL))
contrast(m3b_ref, method = list("Effect of OT on spouse RT in CG > effect of OT on spouse RT in NCG" = spouse_cg_OTvPL - spouse_ncg_OTvPL))
contrast(m3b_ref, method = list("Effect of OT on living loved one RT in CG > effect of OT on living loved one RT in NCG" = living_cg_OTvPL - living_ncg_OTvPL))
contrast(m3b_ref, method = list("Effect of OT on generic death RT in CG > effect of OT on generic death RT in NCG" = death_cg_OTvPL - death_ncg_OTvPL))
contrast(m3b_ref, method = list("Effect of OT on stranger RT in CG > effect of OT on stranger RT in NCG" = stranger_cg_OTvPL - stranger_ncg_OTvPL))
contrast(m3c_ref, method = list("Effect of OT on spouse bias in CG > effect of OT on spouse bias in NCG" = spouse_cg_OTvPL - spouse_ncg_OTvPL))
contrast(m3c_ref, method = list("Effect of OT on living loved one bias in CG > effect of OT on living loved one bias in NCG" = living_cg_OTvPL - living_ncg_OTvPL))
contrast(m3c_ref, method = list("Effect of OT on generic death bias in CG > effect of OT on generic death bias in NCG" = death_cg_OTvPL - death_ncg_OTvPL))
contrast(m3c_ref, method = list("Effect of OT on stranger bias in CG > effect of OT on stranger bias in NCG" = stranger_cg_OTvPL - stranger_ncg_OTvPL))
contrast(m3d_ref, method = list("Effect of OT on spouse bias in CG > effect of OT on spouse bias in NCG" = spouse_cg_OTvPL - spouse_ncg_OTvPL))
contrast(m3d_ref, method = list("Effect of OT on living loved one bias in CG > effect of OT on living loved one bias in NCG" = living_cg_OTvPL - living_ncg_OTvPL))
contrast(m3d_ref, method = list("Effect of OT on generic death bias in CG > effect of OT on generic death bias in NCG" = death_cg_OTvPL - death_ncg_OTvPL))
contrast(m3d_ref, method = list("Effect of OT on stranger bias in CG > effect of OT on stranger bias in NCG" = stranger_cg_OTvPL - stranger_ncg_OTvPL))
c1 <- list(
spouse_death_pl = c(0, # neutral NCG placebo
-.5, # death NCG placebo
0, # living NCG placebo
.5, # spouse NCG placebo
0, # stranger NCG placebo
0, # neutral CG placebo
-.5, # death CG placebo
0, # living CG placebo
.5, # 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
0, # spouse CG oxytocin
0), # stranger CG oxytocin
spouse_death_ot = 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
0, # stranger CG placebo
0, # neutral NCG oxytocin
-.5, # death NCG oxytocin
0, # living NCG oxytocin
.5, # spouse NCG oxytocin
0, # stranger NCG oxytocin
0, # neutral CG oxytocin
-.5, # death CG oxytocin
0, # living CG oxytocin
.5, # spouse CG oxytocin
0), # stranger CG oxytocin
spouse_death_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
-.25, # death CG placebo
0, # living CG placebo
.25, # 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
-.25, # death CG oxytocin
0, # living CG oxytocin
.25, # spouse CG oxytocin
0), # stranger CG oxytocin
spouse_death_pl_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
-.5, # death CG placebo
0, # living CG placebo
.5, # 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
0, # spouse CG oxytocin
0), # stranger CG oxytocin
NCG_CG_spouse_pl = 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
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
NCG_CG_spouse_ot = 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
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
)
library(multcomp)
library(contrast)
library(sjPlot)
L3vsL2 <- contrast(m2b,
a = list(stimulus = levels(data_placebo$stimulus), direction = "PULL"),
b = list(stimulus = levels(data_placebo$stimulus), group = "CG"))
L3vsL2
?contrast
xyplot(data_placebo$gAAT_RT ~ data_placebo$direction * data_placebo$stimulus,
col=c("grey90","grey40"), las=2,
main="test")
contr <- glht(m2b, linfct=mcp(stimulus="Tukey"))
summary(contr)
coefs <- coef(m2b)
# difference between parameter estimates for stimulusdeath:directionPUSH:groupCG and stimulusspouse:directionPUSH:groupCG
unique(coefs[19] - coefs[17]) # use unique() to keep it from repeating 39 times
# difference between parameter estimates for stimulusstranger:directionPUSH:groupCG and stimulusspouse:directionPUSH:groupCG
unique(coefs[19] - coefs[20]) # use unique() to keep it from repeating 39 times
# https://aosmith.rbind.io/2019/04/15/custom-contrasts-emmeans/
# for effect coding
data_long <- data_long %>% mutate(spouse_v_death = ifelse(stimulus == "spouse", 1,
(ifelse(stimulus == "death", -1, 0))),
spouse_v_living = ifelse(stimulus == "spouse", 1,
(ifelse(stimulus == "living", -1, 0))))
data_bias <- data_bias %>% mutate(spouse_v_death = ifelse(stimulus == "spouse", 1,
(ifelse(stimulus == "death", -1, 0))),
spouse_v_living = ifelse(stimulus == "spouse", 1,
(ifelse(stimulus == "living", -1, 0))))
data_placebo <- data_placebo %>% mutate(spouse_v_death = ifelse(stimulus == "spouse", 1,
(ifelse(stimulus == "death", -1, 0))),
spouse_v_living = ifelse(stimulus == "spouse", 1,
(ifelse(stimulus == "living", -1, 0))))
data_bias_placebo <- data_bias_placebo %>% mutate(spouse_v_death = ifelse(stimulus == "spouse", 1,
(ifelse(stimulus == "death", -1, 0))),
spouse_v_living = ifelse(stimulus == "spouse", 1,
(ifelse(stimulus == "living", -1, 0))),
spouse = ifelse(stimulus == "spouse", 1, 0),
death = ifelse(stimulus == "death", 1, 0))
The scripts below generate the data files used in the analyses above.
Pre-R steps:
In Terminal, within each of the folders in turn…
find . -iname "*.iqdat" -exec bash -c 'mv "$0" "${0%\.iqdat}.tsv"' {} \;cat *.tsv > OT_gAAT_run-1.tsv (or OT_gAAT_run-2.tsv)library(tidyverse)
filter <- dplyr::filter
select <- dplyr::select
# read in the data (only importing a subset of columns)
raw_r1 <- as_data_frame(read_tsv("~/Dropbox/GLASS Lab/OT Study/data/raw-data/OT-gAAT-behavioral-data/Run1_tsv/OT_gAAT_run-1.tsv", col_types=cols_only(
date = col_integer(),
time = col_character(),
subject = col_character(),
blockcode = col_character(),
blocknum = col_character(),
trialcode = col_character(),
values.trialcode = col_character(),
values.stimulus = col_character(),
values.initialresponse = col_character(),
values.RT = col_number()
)))
raw_r2 <- as_tibble(read_tsv("~/Dropbox/GLASS Lab/OT Study/data/raw-data/OT-gAAT-behavioral-data/Run2_tsv/OT_gAAT_run-2.tsv", col_types=cols_only(
date = col_integer(),
time = col_character(),
subject = col_character(),
blockcode = col_character(),
blocknum = col_character(),
trialcode = col_character(),
values.trialcode = col_character(),
values.stimulus = col_character(),
values.initialresponse = col_character(),
values.RT = col_number()
)))
# header line repeats when I used `cat` to merge, so get those out of there
raw_r1 <- raw_r1 %>% filter(!is.na(date))
raw_r2 <- raw_r2 %>% filter(!is.na(date))
unique(raw_r1$subject)
# D118 is showing up as "test" for some reason (checked against the .iqdat file), and D130_b is missing
# both of D130's and D142's visits were entered as "D130"/"D142", so change the subject ID for the visit on 052416 to D130_b and on
raw_r1$subject[raw_r1$date == 020916] <- "D118"
raw_r1$subject[raw_r1$date == 052416] <- "D130_b"
raw_r1$subject[raw_r1$date == 111516] <- "D142_b"
# fixing some other stuff in a similar vein
raw_r1$subject[raw_r1$subject == "D101_2"] <- "D101_b"
raw_r1$subject[raw_r1$subject == "D102_B"] <- "D102_b"
raw_r1$subject[raw_r1$subject == "D117_b"] <- "D117"
raw_r1$subject[raw_r1$subject == "D117_c"] <- "D117_b"
raw_r1$subject[raw_r1$subject == "D135_c"] <- "D135_b"
length(unique(raw_r1$subject)) # 78 unique - that is correct
unique(raw_r2$subject)
# some similar issues to fix
raw_r2$subject[raw_r2$subject == "D101_2Y"] <- "D101_b"
raw_r2$subject[raw_r2$subject == "D102_B"] <- "D102_b"
raw_r2$subject[raw_r2$subject == "D107_B_2"] <- "D107_b"
raw_r2$subject[raw_r2$subject == "D117_b"] <- "D117"
raw_r2$subject[raw_r2$subject == "D117_c"] <- "D117_b"
raw_r2$subject[raw_r2$subject == "D126_B"] <- "D126_b"
raw_r2$subject[raw_r2$subject == "D135_c"] <- "D135_b"
raw_r2$subject[raw_r2$date == 111516] <- "D142_b"
length(unique(raw_r2$subject)) # 78 unique - that is correct
str(raw_r1)
str(raw_r2)
# take out ITIs
unique(raw_r1$trialcode)
mlm_r1 <- filter(raw_r1, grepl("^A",trialcode))
unique(mlm_r1$trialcode)
unique(raw_r2$trialcode)
mlm_r2 <- filter(raw_r2, grepl("^A",trialcode))
unique(mlm_r2$trialcode)
# add a "trialnum" variable, grouped by ID
mlm_r1 <- mlm_r1 %>% group_by(subject) %>% mutate(trialnum = row_number())
View(mlm_r1)
mlm_r2 <- mlm_r2 %>% group_by(subject) %>% mutate(trialnum = row_number())
View(mlm_r2)
mlm_r1 %>% group_by(subject) %>% count() # everyone has 144 trials
mlm_r2 %>% group_by(subject) %>% count() # everyone has 144 trials, except D122_b has 177
# no clue why that happened, although I do vaguely remember we had one session where the task seemed to go on and on and we force quit it?
# remove trials 145-177 for D122_b:
mlm_r2 <- filter(mlm_r2, trialnum <=144)
mlm_r2 %>% group_by(subject) %>% count()
# now everyone has 144 trials
# add a column for "run"
mlm_r1$run <- "1"
mlm_r2$run <- "2"
# add a column for "visit"
mlm_r1$visit <- as.factor(ifelse(grepl("*_b", mlm_r1$subject), "2", "1"))
mlm_r2$visit <- as.factor(ifelse(grepl("*_b", mlm_r2$subject), "2", "1"))
mlm_r1 %>% group_by(subject, visit, run) %>% count() # all looks good: everyone has 144 trials for each run at each visit
mlm_r2 %>% group_by(subject, visit, run) %>% count() # all looks good: everyone has 144 trials for each run at each visit
# now remove the "_b" from visit 2 IDs
mlm_r1 <- mlm_r1 %>% ungroup() %>%
mutate(subject = str_replace(subject, "_b", ""))
unique(mlm_r1$subject)
mlm_r2 <- mlm_r2 %>% ungroup() %>%
mutate(subject = str_replace(subject, "_b", ""))
unique(mlm_r2$subject)
# merge data from the 2 runs
bx <- bind_rows(mlm_r1, mlm_r2)
# split up by visit
bx_v1 <- bx %>% filter(visit == "1")
bx_v2 <- bx %>% filter(visit == "2")
bx_v1 %>% group_by(subject, run) %>% count() # everyone has 144 trials for each run at visit 1
bx_v2 %>% group_by(subject, run) %>% count() # everyone has 144 trials for each run at visit 2
# add columns for "cond" at visit 1 and visit 2 [condition: A/placebo or B/oxytocin]
# load in the randomization data
randomize <- readRDS("~/Dropbox/GLASS Lab/OT Study/data/cleaned-data/randomization.rds")
str(randomize)
bx_v1 <- bx_v1 %>%
rename("ID" = subject)
bx_v2 <- bx_v2 %>%
rename("ID" = subject)
# get the IDs for everyone who got treatment A at visit 1
IDs_v1_txA <- randomize %>% filter(tx_v1 == "A")
txA_list <- IDs_v1_txA$ID
txA_list
# ifelse statement: if ID = any of those listed in txA_list, make cond_v1 = "A", else make it "B"
bx_v1 <- bx_v1 %>%
mutate(cond_v1 = as.factor(ifelse(ID == "D101" | ID == "D105" | ID == "D109" | ID == "D110" | ID == "D114" | ID == "D115" | ID == "D116" | ID == "D118" | ID == "D119" | ID == "D120" | ID == "D123" | ID == "D124" | ID == "D126" | ID == "D128" | ID == "D132" | ID == "D135" | ID == "D136" | ID == "D137" | ID == "D138" | ID == "D140" | ID == "D141" | ID == "D145" | ID == "D148" | ID == "D149" | ID == "D117", "A", "B")), cond_v2 = as.factor(ifelse(ID == "D101" | ID == "D105" | ID == "D109" | ID == "D110" | ID == "D114" | ID == "D115" | ID == "D116" | ID == "D118" | ID == "D119" | ID == "D120" | ID == "D123" | ID == "D124" | ID == "D126" | ID == "D128" | ID == "D132" | ID == "D135" | ID == "D136" | ID == "D137" | ID == "D138" | ID == "D140" | ID == "D141" | ID == "D145" | ID == "D148" | ID == "D149" | ID == "D117", "B", "A")))
# ifelse statement: if ID = any of those listed in txA_list, make cond_v2 = "B", else make it "A" (people who were A at visit 1 should be B for visit 2, and vice versa)
bx_v2 <- bx_v2 %>%
mutate(cond_v2 = as.factor(ifelse(ID == "D101" | ID == "D105" | ID == "D109" | ID == "D110" | ID == "D114" | ID == "D115" | ID == "D116" | ID == "D118" | ID == "D119" | ID == "D120" | ID == "D123" | ID == "D124" | ID == "D126" | ID == "D128" | ID == "D132" | ID == "D135" | ID == "D136" | ID == "D137" | ID == "D138" | ID == "D140" | ID == "D141" | ID == "D145" | ID == "D148" | ID == "D149" | ID == "D117", "B", "A")), cond_v1 = as.factor(ifelse(ID == "D101" | ID == "D105" | ID == "D109" | ID == "D110" | ID == "D114" | ID == "D115" | ID == "D116" | ID == "D118" | ID == "D119" | ID == "D120" | ID == "D123" | ID == "D124" | ID == "D126" | ID == "D128" | ID == "D132" | ID == "D135" | ID == "D136" | ID == "D137" | ID == "D138" | ID == "D140" | ID == "D141" | ID == "D145" | ID == "D148" | ID == "D149" | ID == "D117", "A", "B")))
levels(bx_v1$cond_v1)
levels(bx_v2$cond_v2)
head(bx_v1)
txA_list
bx_v1 %>% group_by(cond_v1) %>% count(ID) # IDs and cond match txA_list *note that some IDs from the randomization dataset are listed in txA but not in the behavioral data. This is because some participants were dropped or never ended up completing their experimental session: D109 (dropped after 1st session), D116 and D136 (never enrolled), D124 (incidental finding, excluded after T1MPRAGE).
bx_v2 %>% group_by(cond_v2) %>% count(ID) # IDs and cond match txA_list
bx_long <- bind_rows(bx_v1,bx_v2)
View(bx_long)
# rename some variables
colnames(bx_long)
bx_long <- bx_long %>% rename("push_pull" = values.initialresponse,
"gAAT_RT" = values.RT)
head(bx_long)
# add a new variable for condition
bx_long <- bx_long %>% mutate(cond = as.factor(ifelse(cond_v1 == "A" & visit == "1", "A",
ifelse(cond_v2 == "A" & visit == "2", "A",
ifelse(cond_v1 == "B" & visit == "1", "B",
ifelse(cond_v2 == "B" & visit == "2", "B", NA))))))
View(bx_long) # check that it coded the "cond" column correctly
# add a new variable for stimulus category
# in values.stimulus:
# 1 = spouse, 2 = living loved one/WHOTO, 3 = stranger, 4 = nomothetic death-related, 5 = neutral images
sort(unique(bx_long$values.stimulus))
bx_long <- bx_long %>%
mutate(stim = as.factor(ifelse(values.stimulus == "1B_1.jpg"|values.stimulus == "1B_2.jpg"|values.stimulus == "1B_3.jpg"|values.stimulus == "1Y_1.jpg"|values.stimulus == "1Y_2.jpg"|values.stimulus == "1Y_3.jpg", "spouse",
ifelse(values.stimulus == "2B_1.jpg"|values.stimulus == "2B_2.jpg"|values.stimulus == "2B_3.jpg"|values.stimulus == "2Y_1.jpg"|values.stimulus == "2Y_2.jpg"|values.stimulus == "2Y_3.jpg", "living",
ifelse(values.stimulus == "3B_1.jpg"|values.stimulus == "3B_2.jpg"|values.stimulus == "3B_3.jpg"|values.stimulus == "3Y_1.jpg"|values.stimulus == "3Y_2.jpg"|values.stimulus == "3Y_3.jpg", "stranger",
ifelse(values.stimulus == "4B_1.jpg"|values.stimulus == "4B_2.jpg"|values.stimulus == "4B_3.jpg"|values.stimulus == "4Y_1.jpg"|values.stimulus == "4Y_2.jpg"|values.stimulus == "4Y_3.jpg", "death",
ifelse(values.stimulus == "5B_1.jpg"|values.stimulus == "5B_2.jpg"|values.stimulus == "5B_3.jpg"|values.stimulus == "5Y_1.jpg"|values.stimulus == "5Y_2.jpg"|values.stimulus == "5Y_3.jpg", "neutral", NA)))))))
levels(bx_long$stim)
bx_long$stim <- relevel(bx_long$stim, ref = "neutral") # make neutral the reference level
# save it
saveRDS(bx_long, "~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/bx_long.rds")
Outlier removal follows conventions from other gAAT papers, which exclude trials with RTs that are equal to or above the 99th percentile (calculated from the sample as a whole) or equal to or below 1st percentile.
# calculate percentiles separately for placebo/oxytocin conditions
bx_pl <- bx_long %>% filter(cond == "A")
bx_ot <- bx_long %>% filter(cond == "B")
p99pl <- quantile(bx_pl$gAAT_RT, 0.99) # 1716.76 ms
p01pl <- quantile(bx_pl$gAAT_RT, 0.01) # 463 ms
p99ot <- quantile(bx_ot$gAAT_RT, 0.99) # 1711.07 ms
p01ot <- quantile(bx_ot$gAAT_RT, 0.01) # 473 ms
# identify placebo condition outliers
bx_pl <- bx_pl %>%
mutate(outlier_RT = ifelse(gAAT_RT <= p01pl | gAAT_RT >= p99pl, 1, 0))
table(bx_pl$outlier_RT)
227/11005 # about 2% of trials are outliers
# identify oxytocin condition outliers
bx_ot <- bx_ot %>%
mutate(outlier_RT = ifelse(gAAT_RT <= p01ot | gAAT_RT >= p99ot, 1, 0))
table(bx_ot$outlier_RT)
229/11003 # about 2% of trials are outliers
bx_long_noOut <- bind_rows(bx_pl, bx_ot) %>% filter(outlier_RT == 0)
dim(bx_long)
dim(bx_long_noOut) # after outliers removed, there are 22008 trials total of the original 22464: 22464 - (227+229) = 22008
saveRDS(bx_long_noOut, "~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/bx_long_noOut.rds")
write_csv(bx_long_noOut, "~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/bx_long_noOut.csv")
library(psych)
# how much missing data is there?
ntrials <- bx_long_noOut %>% group_by(ID, cond) %>% count() # don't group by push_pull because peopled pulled/pushed on different numbers of trials, for ex. responses in the wrong direction
describe(ntrials$n/288) # everyone has at least 90% of the 288 trials expected at each visit (no more than 10% trials dropped because of outliers or missed responses)
describe(bx_long_noOut$gAAT_RT) # skewness is 1.17, kurtosis is 1.7
hist(bx_long_noOut$gAAT_RT) # distribution is skewed, as usual for RT data...
qqnorm(bx_long_noOut$gAAT_RT) # and skewness is reflected in QQ plot
# 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)
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 = median(neutral, na.rm = TRUE), deathApush = median(death, na.rm = TRUE), spouseApush = median(spouse, na.rm=TRUE), livingApush = median(living, na.rm=TRUE), strangerApush = median(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 = median(neutral, na.rm = TRUE), deathApull = median(death, na.rm = TRUE), spouseApull = median(spouse, na.rm=TRUE), livingApull = median(living, na.rm=TRUE), strangerApull = median(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 = median(neutral, na.rm = TRUE), deathBpush = median(death, na.rm = TRUE), spouseBpush = median(spouse, na.rm=TRUE), livingBpush = median(living, na.rm=TRUE), strangerBpush = median(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 = median(neutral, na.rm = TRUE), deathBpull = median(death, na.rm = TRUE), spouseBpull = median(spouse, na.rm=TRUE), livingBpull = median(living, na.rm=TRUE), strangerBpull = median(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)
# 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)
# which are these?
bias_long %>% filter(outlier==TRUE)
### this is another way to do it using boxplot.stats
# outlier_values <- boxplot.stats(medians_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.rds")
write_csv(bias_long, "~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/bias_long.csv")
The four extreme outliers (>3 SD) are:
To calculate trial-level bias (aka “compatibility”) scores: separate push minus pull trials, then subtract RT on that trial at run1 from RT on that trial at run2.
bx_long <- readRDS("~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/bx_long.rds")
bx_r1.1 <- bx_long %>% filter(run == "1")
bx_r2.1 <- bx_long %>% filter(run == "2")
tail(bx_r1.1)
tail(bx_r2.1)
# run1 should line up with run2: "stim" should match (i.e, spouse trials with spouse trials) but push_pull should say PUSH at one run and PULL at the other
# calculate percentiles for outlier removal (run1 and run2)
p99r1 <- quantile(bx_r1.1$gAAT_RT, 0.99)
p01r1 <- quantile(bx_r1.1$gAAT_RT, 0.01)
p99r1 # 99th percentile RTs = 1713.38ms
p01r1 # 1st percentile RTs = 472ms
p99r2 <- quantile(bx_r2.1$gAAT_RT, 0.99)
p01r2 <- quantile(bx_r2.1$gAAT_RT, 0.01)
p99r2 # 99th percentile RTs = 1712.45ms
p01r2 # 1st percentile RTs = 463ms
bx_r1.2 <- bx_r1.1 %>%
mutate(outlier_RT = ifelse(gAAT_RT <= p01r1 | gAAT_RT >= p99r1, 1, 0))
table(bx_r1.2$outlier_RT)
227/11005 # 227 trials will be dropped, 11005 will be retained...about 2% of trials.
bx_r2.2 <- bx_r2.1 %>%
mutate(outlier_RT = ifelse(gAAT_RT <= p01r2 | gAAT_RT >= p99r2, 1, 0))
table(bx_r2.2$outlier_RT)
230/11002 # 230 trials will be dropped, 11002 will be retained...about 2% of trials.
bx_r1r2 <- bind_cols(bx_r1.2, bx_r2.2)
dim(bx_r1r2) # after outliers removed, there are 11232 trials total
# make sure that trials in the two runs line up post-merge
head(bx_r1r2)
tail(bx_r1r2)
# get rid of any trials with outliers in either run1 or run2
# (need run1/run2 to be the same length so that the trials match up)
bx_r1r2.1 <- bx_r1r2 %>% filter(outlier_RT == "0" & outlier_RT1 == "0")
count(bx_r1r2.1) # 10790 obs.
10790/11232 # ~96% of trials retained = 4% data excluded
push_pull <- bx_r1r2.1 %>% filter(push_pull == "PUSH" & push_pull1 == "PULL") %>% mutate(bias = (gAAT_RT-gAAT_RT1))
pull_push <- bx_r1r2.1 %>% filter(push_pull == "PULL" & push_pull1 == "PUSH") %>% mutate(bias = (gAAT_RT1-gAAT_RT))
bx_pp <- bind_rows(push_pull,pull_push)
head(bx_pp) # can see from this that D101 (and possibly others) has only a few trials, because they failed to reverse the instructions (i.e., push_pull and push_pull1 columns have the same values instead of PUSH at one and PULL at the other)
# in essence, this removes their "error" trials
bx_pp %>% group_by(ID, cond) %>% count() # yep, this is visible in their *very low* trialcounts:
# D101 A n = 3
# D141 A n = 18
# D147 B n = 1
# D148 A n = 19
As we see from the last part of the script above, 4 participants have extremely low trial counts at one visit or another, due to failure to reverse the instructions. This will be addressed next.
# create a "drop" variable for these 4
bx_pp <- bx_pp %>% mutate(low_n_trials = ifelse(ID == "D101" & cond == "A" | ID == "D141" & cond == "A"| ID == "D148" & cond == "A" | ID == "D147" & cond == "B", 1, 0))
# check whether the rows line up
notlinedup <- bx_pp %>% filter(stim != stim1)
notlinedup
notlinedup <- bx_pp %>% filter(trialcode != trialcode1)
notlinedup
notlinedup <- bx_pp %>% filter(trialnum != trialnum1)
notlinedup
notlinedup <- bx_pp %>% filter(push_pull == push_pull1)
notlinedup
# 0 rows returned; all the rows match. Thus we conclude the rows from the two runs line up.
saveRDS(bx_pp, "~/Dropbox/GLASS Lab/OT gAAT for meeting 5-10-19/bx_long_trialbytrialBias.rds")