I have a dataset from participants that provided liking ratings (on a scale from 0-100) of stimuli associated with rewards of different magnitudes (factor pval, with levels small/medium/large) and delay (factor time, with levels delayed/immediate). A subset of the data looks like this:
structure(list(ppnr = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 2L,
2L, 2L, 2L, 2L, 2L), .Label = c("7", "8"), class = "factor"),
time = structure(c(2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L, 2L, 1L,
2L, 1L), .Label = c("del", "imm"), class = "factor"), pval = structure(c(1L,
1L, 2L, 2L, 3L, 3L, 1L, 1L, 2L, 2L, 3L, 3L), .Label = c("pval_L",
"pval_M", "pval_S"), class = "factor"), rating = c(66, 55,
81, 11, 30, 11, 18, 28, 85, 61, 6, 76), stimJPG = structure(c(1L,
2L, 3L, 5L, 4L, 6L, 5L, 1L, 3L, 2L, 6L, 4L), .Label = c("pStim01.jpg",
"pStim02.jpg", "pStim03.jpg", "pStim04.jpg", "pStim05.jpg",
"pStim06.jpg"), class = "factor")), row.names = 283:294, class = "data.frame")
ppnr time pval rating stimJPG
283 7 imm pval_L 66 pStim01.jpg
284 7 del pval_L 55 pStim02.jpg
285 7 imm pval_M 81 pStim03.jpg
286 7 del pval_M 11 pStim05.jpg
287 7 imm pval_S 30 pStim04.jpg
288 7 del pval_S 11 pStim06.jpg
289 8 imm pval_L 18 pStim05.jpg
290 8 del pval_L 28 pStim01.jpg
291 8 imm pval_M 85 pStim03.jpg
292 8 del pval_M 61 pStim02.jpg
293 8 imm pval_S 6 pStim06.jpg
294 8 del pval_S 76 pStim04.jpg
To investigate whether the ratings were influenced by the time and magnitude of the reward associated with the cues, I ran the following model in brms:
n_chains <- 4
n_threads <- 2
options(contrasts = c("contr.sum", "contr.poly"))
model <- brm(rating ~ time*pval + (1 + time + pval | ppnr) + (1 + time * pval | stimJPG), data = data, backend = "cmdstanr", chains = n_chains, cores = n_chains, threads = threading(n_threads), iter = 4000, warmup = 2000, control = list(adapt_delta = 0.9999, max_treedepth = 15))
Next, I wanted to draw samples from the posterior distribution for two specific contrasts (i.e., two pairwise comparisons). First, I obtained the estimates for those contrasts using emmeans. I could in principle use the function gather_emmeans_draws (from the tidybayes package) to draw samples from the posterior of these contrasts without problem. However, going a step back, emmeans uses the median as point estimate for Bayesian models, whereas I would like to use the mean. Obtaining the mean is possible by using hpd.summary on the emmeans object. However, this converts the emmGrid object that is created by emmeans into a summary_emm object. Unfortunately, gather_emmeans_draws() does not accept summery_emm objects, but only seems to accept emmGrid objects (or S4 objects in general). See:
emm_withmedian <- emmeans(model, pairwise ~ pval * time)$contrasts
emm_withmean <- hpd.summary(emm_withmedian, point.est = mean)
#This results in all pairwise comparisons, but I am only interested in 2 of these:
contrast estimate lower.HPD upper.HPD
pval_L del - pval_M del 14.79 3.87 26.17
pval_L del - pval_S del 26.67 11.69 42.55
pval_L del - pval_L imm -5.85 -17.98 7.67
pval_L del - pval_M imm 9.51 -3.17 22.61
pval_L del - pval_S imm 17.70 4.23 31.75
pval_M del - pval_S del 11.89 -1.45 26.84
pval_M del - pval_L imm -20.64 -33.83 -7.43
pval_M del - pval_M imm -5.28 -18.19 6.96
pval_M del - pval_S imm 2.91 -9.40 16.33
pval_S del - pval_L imm -32.53 -47.46 -18.05
pval_S del - pval_M imm -17.16 -29.95 -3.68
pval_S del - pval_S imm -8.98 -22.43 5.10
pval_L imm - pval_M imm 15.36 4.28 26.58
pval_L imm - pval_S imm 23.55 9.58 39.50
pval_M imm - pval_S imm 8.19 -4.94 22.43
Point estimate displayed: mean
HPD interval probability: 0.95
#I then want to draw from the posterior, but that's where it goes wrong:
posteriorsamples <- gather_emmeans_draws(emm_withmean)
Error in as_tibble(object@grid) :
trying to get slot "grid" from an object (class "summary_emm") that is not an S4 object
#Just for comparison's sake, if I would do the following, it would be no problem, because it uses an emmGrid object as input:
posteriorsamples2 <- gather_emmeans_draws(emm_withmedian)
Thus, it seems that I can only draw from the posterior if I work directly from the emmGrid object (emm_withmedian), forcing me to use the median instead of the mean.
I already tried converting the summary_emm object into an emmGrid object using as.emmGrid(), but that does not work, and gives me the following error: Error in nrow(V) : argument "V" is missing, with no default.
I already looked into both error messages, but haven't found a way to work around them. I also made sure to update all packages used, but that also didn't help.
Thus, I am looking for:
Any ideas are highly appreciated!
Regarding the first question: As is true of most summary methods, the returned object is just a summary, and it doesn't contain the information to convert it back to an object like the one that was summarized. However, the original emmGrid object does have all the needed content.
The other barrier is trying to work from the contrasts you don't want rather than getting the ones you do want. It is usually best to do the means and contrasts in two separate steps. It is quite simple to do:
EMM <- emmeans(model, ~ pval * time)
CON <- contrast(EMM, list(c1 = c(<desired coefs>), c2 = c(<desired coefs>)))
draws <- coda::as.mcmc(CON)
The latter does the same computation internally, making the linfct slot the identity matrix and the post.beta slot equal to the estimates.
Either way, the fact that the summary you see shows the median has no importance -- draws contains the posterior draws of the desired contrasts, and if you want a point estimate of those draws, you may use the mean, median, or whatever else you want.