MLM ADA #1 v2

About

Overview

This ADA uses data from the project “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, designed to address conflicting accounts of motivated behavior in complicated grief using a novel approach-avoidance task variant (grief AAT; gAAT) that allows for comparison between participant-provided photos and nomothetic grief stimuli. Participants were 40 bereaved adults ages 55-80, who experienced the death of a spouse or long-term romantic partner between six months and three years prior to their participation.

Task

Participants attended two identical experimental sessions, at which they received one of two intranasal sprays (treatment A/B) and then took part in an fMRI scan, which involved a structural sequence, the gAAT, and a resting state sequence.

For the gAAT, five stimulus categories were presented: (1) deceased spouse, (2) living loved one [another attachment figure], (3) stranger, (4) generic grief-related scenes, and (5) neutral scenes. Photos provided by the participant were scanned and resized to a standard 750-pixel height and width, and framed in colored frames to match the stimuli used for non-idiographic categories (i.e., stranger, grief-related, neutral). Photos of the stranger were sex-matched to the participant’s partner.

During the task, participants were instructed to follow rules based on an irrelevant characteristic of the stimuli, namely, the color of the frame around the photo (e.g., “push the joystick when the frame is blue, pull when the frame is yellow”). The task was completed in two runs of approximately seven minutes each. In the second run, the instructions were reversed (e.g. “now pull the joystick for blue framed photos”). 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.

Participants completed the task using an MRI-compatible joystick. The stimuli were animated so that a joystick push would make the image smaller (as if they were pushing it away from them) and a joystick pull would make the image larger (as if they were bringing it nearer to them).

Variables

The primary dependent variable of interest for the present analysis is reaction time bias, computed as RT_push-RT_pull. Negative bias scores (faster to push than pull) indicate relative avoidance bias, while positive scores (faster to pull than to push) indicate relative approach bias. Because the two runs were identical in every way except for the instruction provided, each trial had both a push RT and a pull RT (e.g., run1 trial1 was push, run2 trial1 was pull, so subtract RT on run2 trial 1 from run1 trial1 to get the difference - this is further elaborated in the data cleaning script below), unless they failed to reverse the instructions.

Main predictors of interest included trial number (in order to establish whether peoples’ response bias changed over the 7-minute run), stimulus (neutral, spouse, stranger, living loved one, death), and condition (treatment A/B - oxytocin or placebo [I’m still blinded]).

Participants also completed a number of self-report measures (state and trait) and provided biological samples for OXTR genotyping and cortisol assay, some of which were used as predictors.

Equations

Level 1:

We first performed MLM to estimate the effect of trial number (i.e., the random effect of time) on response bias (push RT-pull RT) on a given trial. The level-1 model estimated, for each participant j (j = 1–35), a regression line that predicted each participant’s response bias in milliseconds on each trial i from the trial number. This model was represented as follows:

[bias_i]_j = b_0j + b_1j [trialnumber_i] + r_ij

where b_0j, the intercept, is interpreted as participant j’s mean response bias for a given trial; b_1j, the slope, represents the effect of trial number on response bias for each participant j (positive coefficients represent the predictive magnitude of trial number on response bias for a given participant); and r_ij is the residual error term.

Level 2:

The level-2 models estimated the average effects for the entire sample. The models were represented as follows:

b_0j = γ₀₀ + μ_0j

b_1j = γ₁₀ + μ_1j

where the intercept, γ₀₀, is interpreted as the average response bias for the entire sample; γ₁₀ is the average effect of represents the effect of trial number on response bias for each participant j (positive coefficients represent the predictive magnitude of trial number on response bias for the sample as a whole), and thus, the primary estimate of interest; and μ_0j and μ_1j are the residual error terms.

Adding the fixed effect of stimulus type resulted in the following model:

b_0j = γ₀₀ + γ₀₁(stimulus_i) + μ_0j

b_1j = γ₁₀ + γ₁₁(stimulus_i) + μ_1j

Data cleaning

Note: Went back to the raw data due to some inconsistencies in the OT_ADAA_Omnibus.xlsx file from BA that I found.

Pre-R steps:

1. Removed all .iqdat files under 25KB or so in order to filter out pilot data and false starts.
   
2. Removed non-study data (i.e., tests and OT undergrad files).
   
3. Organized raw .iqdat files for each visit by run 1/run 2 (two separate folders) based on date/time stamp.

   In Terminal, within each of the folders in turn…

4. Changed file extensions to .tsv (tab-separated) from .iqdat: find . -iname "*.iqdat" -exec bash -c 'mv "$0" "${0%\.iqdat}.tsv"' {} \;
   
5. Merged the files: cat *.tsv > OT_gAAT_run-1.tsv (or OT_gAAT_run-2.tsv) 5. Merged the files: mlr --tsvlite cat *.tsv > OT_gAAT_run-1.tsv (or OT_gAAT_run-2.tsv)
   5a. mlr --tsvlite drops the header from all but the first file, so column names are not repeated throughout the dataset.
   5b. This requires Miller, which can be installed via Homebrew on OSX (http://johnkerl.org/miller/doc/build.html) tried to make this work but no luck

Notes:

• D117: Had 3 visits due to scanner technical issues. Use “_b" and “_c" as visit 1 and visit 2 respectively.
• D114: Has multiple files from visit 2. Use the last two.
• D135: Use files from 6/06/16 and 6/17/18 (scanner issues on 06/13/16). Use “_b" and “_c" as visit 1 and visit 2 respectively.
• D147: Visit 1 run 2 is partial data for some reason (I do have notes that she switched back to PUSH-YELLOW during run 2, but did not re-run because running way over time.)

Creating the MLM dataset

Import data

I do recall one scan session where Inquisit had some kind of hiccup and the task kept going beyond the end of the fMRI sequence so we had to force-quit it. D122 has complete data from the first 144 trials, so we can just drop the extra ones.

Bias scores

Separate out push minus pull, then subtract run1 from run2 to generate a score quantifying relative approach or avoidance bias. The top and bottom 1% of reaction times (RTs) will be removed as part of conventional data cleaning for this task.

Transformations

Because bias is quite non-normal, looking at the plots above, analyses will use the transformed bias score (bias_t) instead.

Checking the data

Add time-invariant variables

Next step is to add demographics and other self-reported and genetic data variables that we might potentially want to include in the model but that do not change over time. Looking at the codebook for the master dataset, these include:

Demographics:

• age
• age_yrs
• sex_m
• ethnicity_hisp
• race
• timesincedeath

Genetic variables:

• rs2254298
• rs2268498
• rs53576
• rs2254298_A.carrier
• rs2268498_C.carrier
• rs53576_A.carrier

Self-report measures:

• tot_icg [grief]
• tot_bdi [depression]
• tot_bgq [brief grief questionnaire with impairment item]
• tot_ecrrs_global_anx [anxious attachment - global]
• tot_ecrrs_global_avoid [avoidant attachment - global]
• tot_ecrrs_spouse_anx [anxious attachment - spouse]
• tot_ecrrs_spouse_avoid [avoidant attachment - spouse]
• tot_ysl [yearning for the deceased]

Exploratory analysis and checking assumptions

First step is to overview the entire dataset using summarytools.

None of the variables have any missing data, apart from rs2254298 and rs2254298_A.carrier, and I’m aware from the data cleaning that I did for the master dataset that one participant’s genotype for that SNP was unable to be identified.

Descriptives

Frequencies

Checking distribution shapes and normality

Got a preview of this above when we used summarytools, but here’s a closer look at some histograms and QQ plots.

The most problematic of these appear to be anxious/avoidant attachment style for the spouse (tot_ecrrs_spouse_anx and tot_ecrrs_spouse_avoid), global anxious attachment style (tot_ecrrs_global_anx), depressive symptoms (tot_bdi), and time since spouse’s death (timesincedeath - in days).

Centering L2 predictors

Evaluate missing data

One person is missing rs2254298 genotype (rs2254298 and rs2254298_A.carrier), which we already knew about. No one else is missing data on any variable.

Scale reliabilities

This dataset does not include the item-level data, but we can get that from the master dataset we loaded in earlier.

All reliabilities look good or excellent (a >=.86) except for the Brief Grief Questionnaire, where a = .72. It’s a short screening measure with only 7 items, so makes sense that its reliability would be lower.

Drop the 4 people with low trial counts

Recall that these participants had bias scores on <14% of trials, likely due to failing to reverse the instructions on run 2.

OLS estimates for exploratory purposes

Peoples’ intercepts do look fairly different from each other…

Stem plot for the intercepts based on the OLS estimates

Stem plot for the slopes based on the OLS estimates

Note that these are all extremely close to zero - probably because of the lack of time effect.

R-square for OLS models

The correlation between intercepts and slopes is -0.475, so people who show less avoidance bias overall change less over the course of the task.

Graphing multilevel data for exploratory purposes

There does not seem to be either a linear or curvilinear effect of trial number (i.e., time) on RT bias in the task.

Analyses (NHST)

Remember, the effect of time refers to trialnumber (this is NOT the time variable; time refers to time of day when task was run.)

The dependent variable, bias_t, is the measure of relative behavioral approach/avoidance bias, based on reaction time to push or pull a joystick in response to various images. We transformed it because it was not normally distributed, and use the transformed variable for all analyses.

Unconditional means model

The overall mean bias score (transformed) is about -0.002, with a standard deviation of 0.023.

ICC is about 1.6% of the variance. This means that 1.6% of the variance in bias is explained by people being different from one another, while the remainder is explained by…something else.

Add fixed and random intercept and slope of time

The average bias at trial number = 0 is -0.016 and gets increasingly approach-biased by a nonsignificant fraction of a percent with each unit of time. There is no significant effect of trial number (time) on RT bias score, t = 0.74, p = 0.457.

99.99% of variance in bias is due to the between-person variance, after controlling for time. This is consistent with the lack of a significant effect of time.

Adding stimulus as a fixed effect

There were 5 different types of stimuli, presented in a randomized order (consistent across participants/sessions) determined by genetic algorithm. Everyone saw the same stimuli in the same order, so it is a fixed effect. The different categories of images were: Death (generic), Deceased (personal photo), Other living loved one (personal photo), Stranger (sex-matched to spouse), and Neutral (objects). The factor stim was created such that Neutral was the reference level.

NO effect of time. There is a significant fixed effect of spouse, t = 3.55, p = .0004. There is also a non-significant effect of stranger, t = -1.73, p = .084. This indicates that people tended to show more approach bias to photos of their spouse (across conditions). The parameter estimate is 0.12 with an SE of 0.03. They also tended to show greater more avoidance bias to photos of the sex-matched stranger.

95% CIs don’t cross zero for spouse, but do for all other stimulus types. This supports the significant effect of spouse stimuli on reaction time bias.

The ICC is identical to that in model 2.

Comparing the model with the fixed effect of stimulus (m3) to the model without (m2), m3 is better (p < .0001) as evidenced by lower AIC, BIC, and lL values.

Take out fixed and random effect of time

Because there’s no significant effect of time, makes sense to drop it from the model. Now we’re basically doing regression but accounting for whatever interdependence exists. Note that these results will NOT match the results of the regressions, because the regressions are based on bias score medians as the DV, rather than means.

Fixed effects are still more or less the same - significant for spouse, marginal for stranger - although the parameter estimates and test statistics are slightly different.

As above, spouse is the only stimulus category whose 95 CIs don’t cross zero, supporting a significant effect of spouse on bias - in this case, viewing spouse images was associated with a relative approach bias.

Is m4 any better than m3?

Nope. Models do not significantly differ (p = 0.851), though the AIC and BIC values are a little smaller for m4 than m3. So let’s forget time and go with m4…

Adding condition as a fixed effect

Participants completed the task at two different sessions, where they were randomized to receive treatment A at one and treatment B at the other (order counterbalanced across participants). The two treatments are oxytocin vs. placebo.

The fixed effect of condition B is not significant, t = 0.15, p = 0.884.

The two models are not significantly different, but the one with the fixed effect of condition (m5) has slightly higher AIC/BIC. Not compelling enough to select that model over m4, though.

Adding the stimulus*condition interaction

We might expect that peoples’ tendency to approach or avoid might differ depending on which treatment they got, but only for certain stimuli.

Interesting. With the interaction in the model, there is a significant effect of death and stranger but not spouse. 95% CIs indicate that viewing the generic death and stranger stimuli are both associated with greater relative avoidance bias.

The model with the interaction in it is not significantly better than either the one with just the fixed effect of condition or with no effect of condition at all, so makes sense to go with the simplest model which only includes stimulus as a fixed effect.

Adding other plausible predictors

There are a few other variables that we might expect to impact approach/avoidance bias, namely grief severity, OXTR genotype, time since death, and attachment style (avoidant, anxious).

None of the added fixed effects are significant in the new model, and the AIC/BIC values are higher, so take ’em out.

Final model

A case of significant, but not meaningful: The R² for the entire model is .02, indicating that 2% of the variance in bias is explained by the model. So this model is not very good at predicting approach/avoidance bias from the fixed effect of stimulus within person. Oh well.

Check final model assumptions

Checking normality of the residuals

The QQ plot looks pretty okay, albeit a few outliers. D113 is the person who most falls off the line, and does so at both ends. The residuals are also pretty evenly distributed across time, though many fall outside of the +/- 2 range.

Getting an error here for the random slopes Error in [.data.frame(ranef(m4), 72) : undefined columns selected. I think this may be because I don’t have any random slopes, since I took time out of the model. Indeed, when I go into m4, it doesn’t look like a second element for the slopes exists. Let’s test this…

My suspicion was correct: the error message above is because there are no random slopes for my final model.

Checking homoscedasticity

Note that again there are a lot of residuals falling outside of the +/-2 range.

Also, there’s something weird about the L2 plot. Not sure why it shows only 3 of 5 stimulus types.

Analyses (Bayesian)

Unconditional means model

First step is to look at the unconditional means model. We know from the exploratory stuff in previous sections that bias_t is normally distributed, so we can tell brms to use normal likelihood (Gaussian).

library(brms)
# b1N <- brm(bias_t ~ 1 + (1|ID), data=data, family=gaussian()) 
# saveRDS(b1N, "mlm_ada1_b1N.rds")
summary(b1N)

 Family: gaussian 
  Links: mu = identity; sigma = identity 
Formula: bias_t ~ 1 + (1 | ID) 
   Data: data (Number of observations: 9289) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~ID (Number of levels: 35) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     0.13      0.02     0.10     0.17       1184 1.00

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept    -0.00      0.02    -0.05     0.05       1712 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     0.99      0.01     0.98     1.01       4000 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

Rhat = 1 across the board, so that’s good.

plot(b1N) # check trace plots

I think these plots look okay? b_Intercept and SD_ID__Intercept plots are not exactly on top of each other, but do have a lot of overlap in their distributions.

pp_check(b1N, nsamples=20)  # posterior predictive check using normal likelihood

pp_check(b1N, type="error_hist", nsamples=20, binwidth=1) # error histogram (histogram of prediction residuals)

The observed distribution is a very good fit to the predicted data (based on the normal distribution). The error histogram also shows that residuals are normally distributed around 0. See this page for some helpful documentation.

# check ICC to get sense of how variance is distributed between- and within- people
iSd <- as.numeric(summary(b1N)$random$ID[1,1]) # finds 1st row, 1st column of random effects inside the summary of b1N
iVar <- iSd^2
residSd <- as.numeric(summary(b1N)$spec_pars[[1,1]])
residVar <- residSd^2
icc1 <- iVar/(iVar+residVar) 
icc1 # this is the unconditional growth model - just happens in this dataset that it's a nuisance model

[1] 0.01749842

For the unconditional growth model, the ICC is 0.017, meaning that only about 2% of the variance in bias is because people on average perform differently from each other on the task.

Time as fixed and random effect

## check whether we need time as a control variable
# b2 <- brm(bias_t ~ 1 + trialnum + (1 + trialnum|ID), data=data, family=skew_normal(), control=list(adapt_delta = .95))
# saveRDS(b2, "mlm_ada1_b2.rds")
summary(b2) 

 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: bias_t ~ 1 + trialnum + (1 + trialnum | ID) 
   Data: data (Number of observations: 9289) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~ID (Number of levels: 35) 
                        Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)               0.12      0.03     0.07     0.17       1651 1.00
sd(trialnum)                0.00      0.00     0.00     0.00        510 1.01
cor(Intercept,trialnum)     0.22      0.47    -0.75     0.95       2352 1.00

Population-Level Effects: 
          Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept    -0.02      0.03    -0.07     0.04       1819 1.00
trialnum      0.00      0.00    -0.00     0.00       4000 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     0.99      0.01     0.98     1.01       4000 1.00
alpha    -0.02      0.35    -0.64     0.60       4000 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

Note: Estimate = mean of posterior distribution for each of the parameters.

Group-level effects: When trialnum is 0, the standard deviation of the intercepts is 0.12. There is no random effect of time: the slopes are not at all variable (people are not different from each other). Zero is a credible value for the correlation, because the 95% CIs (credible intervals) cross zero.

Note that Rhat for sd(trialnum) in the random effects is 1.01. Per the brms documentation, this does not seem to indicate a convergence issue due to lack of error messages and the fact that the Rhat for the group-level effect in their example is 1.01: “If ‘Rhat’ is considerably greater than 1, the algorithm has not yet converged and it is necessary to run more iterations and / or set stronger priors.” However, it does indicate that I would probably want to run it again with a higher number of iterations or stronger priors.

Population-level effects: There is no fixed effect of time, either, given that the upper and lower CI limits are -0.0, 0.0.

plot(b2)

SD_ID__trialnum and cor_ID__Intercept__trialnum are totally skewed.

pp_check(b2, nsamples=20)

pp_check(b2, type="error_hist", nsamples=20, binwidth=1) 

Again, the trace plot shows that the observed distribution is a very good fit to the predicted data based on the normal distribution. The error histogram also shows that residuals are normally distributed around 0.

SKIPPING THE CHUNK BELOW… I can’t get sjstats to install (issue with the specific TMB package version- freaking dependencies are the bane of my existence with all these different packages), and no point in running cross validation on a model that we’re going to drop anyway due to lack of both fixed/random effects of time.

# install.packages("sjstats")
# library(sjstats)
# equi_test(b2) # is there a fixed effect?
# waic(b1N, b2) # does WAIC think we should keep time? If so, given there is no fixed effect, this is due to the random effect accounting for variance.
# b1 <- add_ic(b1, ic="R2")
# b2 <- add_ic(b2, ic="R2") 
# mean(b1$R2)
# mean(b2$R2) # does R2 thinks we should keep time? 
# kfold(b1SN, b2, K=10) # use CV to check for over-fitting
# this takes too long, so reduce n(folds) or only perform for final model

Fixed effect of stimulus

# b3 <- brm(bias_t ~ 1 + stim + (1 |ID), data=data, family=skew_normal(), control=list(adapt_delta = .95))
# saveRDS(b3, "mlm_ada1_b3.rds")
summary(b3) 

 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: bias_t ~ 1 + stim + (1 | ID) 
   Data: data (Number of observations: 9289) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~ID (Number of levels: 35) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     0.13      0.02     0.10     0.17        981 1.00

Population-Level Effects: 
             Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept       -0.01      0.03    -0.07     0.05       1156 1.00
stimdeath       -0.05      0.03    -0.12     0.01       2955 1.00
stimliving       0.04      0.03    -0.02     0.11       2423 1.00
stimspouse       0.12      0.03     0.05     0.18       2704 1.00
stimstranger    -0.06      0.03    -0.12     0.01       2560 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     0.99      0.01     0.98     1.00       4000 1.00
alpha    -0.04      0.35    -0.65     0.58       3285 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

Similar to the results of the NHST analysis, there is a notable fixed effect of spouse (estimate = 0.12, est. error = 0.03, 95% CI[0.05,0.18]) on bias, such that people show more approach bias on trials where the stimulus is their spouse.

plot(b3)

marginal_effects(b3)

The high density intervals illustrate the positive effect of spouse on bias, which we also see in the marginal effects plot. In the analysis above, the 95% CI for death and stranger stimuli barely crossed zero, which corresponds with what we see in the marginal effects plot, where the upper limit of the error bar is very close to zero (or below). This also corresponds with the effects we saw in NHST, where there were significant or marginally significant effects of death and stranger in addition to spouse.

pp_check(b3, nsamples=20)

pp_check(b3, type="error_hist", nsamples=20, binwidth=1) # looks like we have improved the residuals; if we rerun this random sampling at least some of the time we no longer see the lower end outliers

# equi_test(b3) # TRY THIS CODE WHEN I GET SJTOOLS TO INSTALL (what does ROPE tells us about the size of effect?)
# waic(b2, b3) # TRY THIS CODE WHEN I GET SJTOOLS TO INSTALL

b2 <- add_ic(b2, ic="R2") 
b3 <- add_ic(b3, ic="R2") 
round(mean(b2$R2), digits=3)

[1] 0.017

round(mean(b3$R2), digits=3)

[1] 0.02

In both models, we have explained about 2 percent of the variance (1.7% vs. 2.0%).

Adding the stimulus*condition interaction (no time)

One of the original study hypotheses was that there would be an interaction (actually, a 3-way interaction with group, but not going to go there for ADA #1; also we may be underpowered for that) between stimulus and condition, such that people might display different response biases for different stimulus categories in condition A vs. B, or vice versa.

# b4 <- brm(bias_t ~ 1 + stim*cond + (1|ID), data=data, family=skew_normal(), control=list(adapt_delta = .95))
# saveRDS(b4, "mlm_ada1_b4.rds")
summary(b4) 

 Family: skew_normal 
  Links: mu = identity; sigma = identity; alpha = identity 
Formula: bias_t ~ 1 + stim * cond + (1 | ID) 
   Data: data (Number of observations: 9289) 
Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup samples = 4000

Group-Level Effects: 
~ID (Number of levels: 35) 
              Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sd(Intercept)     0.13      0.02     0.10     0.17       1088 1.00

Population-Level Effects: 
                   Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
Intercept              0.02      0.04    -0.06     0.10       1357 1.00
stimdeath             -0.11      0.05    -0.20    -0.02       1815 1.00
stimliving             0.04      0.05    -0.05     0.13       1791 1.00
stimspouse             0.06      0.05    -0.03     0.16       2030 1.00
stimstranger          -0.09      0.05    -0.18    -0.00       2068 1.00
condB                 -0.05      0.05    -0.15     0.04       1546 1.00
stimdeath:condB        0.11      0.07    -0.02     0.24       1760 1.00
stimliving:condB       0.01      0.06    -0.12     0.13       1959 1.00
stimspouse:condB       0.11      0.07    -0.02     0.24       2012 1.00
stimstranger:condB     0.07      0.06    -0.06     0.20       2009 1.00

Family Specific Parameters: 
      Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
sigma     0.99      0.01     0.98     1.00       4000 1.00
alpha    -0.05      0.35    -0.64     0.58       3415 1.00

Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
is a crude measure of effective sample size, and Rhat is the potential 
scale reduction factor on split chains (at convergence, Rhat = 1).

# b4 <- readRDS("mlm_ada1_b4.rds")

Got some warnings while running:

Rejecting initial value: Gradient evaluated at the initial value is not finite. Stan can't start sampling from this initial value. Rejecting initial value: Log probability evaluates to log(0), i.e. negative infinity. Stan can't start sampling from this initial value.

Apparently setting priors may help with this.

Here, we see notable effects of death and stranger stimuli, but not condition nor its interaction with stimulus type.

plot(b4)

marginal_effects(b4) 

pp_check(b4, nsamples=20)

pp_check(b4, type="error_hist", nsamples=20, binwidth=1) # residuals seem about the same

# equi_test(b4) 

Nothing terribly new here, except for the marginal effects. Looking at the plot for stim, the distributions for spouse and stranger and spouse and death don’t seem to overlap, and living loved one and death just barely overlap. Unsurprisingly, people tend to show more approach bias for spouse and living loved one photos compared to death and stranger photos. The third marginal effects plot illustrates the lack of condition or stimulus*condition interaction effects.

The variance accounted for shows a negligable increase with the stimulus*condition interaction in the model (goes from 0.020 to 0.022). The negligable effect is shown in the plot - note no real differences

Using priors

# First find out what default priors are. 
get_prior(bias_t ~ 1 + stim*cond + (1|ID), data=data, family=skew_normal())

brms has set some weakly informative default priors on some of the special parameters. The default for the random coefficients is a half student-t with 3 degrees of freedom.

Let’s assume some weak prior knowledge from previous research. It would actually be somewhat difficult to get solid priors from previous research on this task, given that no previous work has used this specific variant in this population with this specific manipulation (i.e., oxytocin vs. placebo) so I’m making these numbers up.

Based on the error message here (The following priors do not correspond to any model parameter…), I clearly have no idea what I’m doing and need to go back to that lesson (get_prior was not helpful to me, but thanks anyway, brms).