Extracting draws from posterior after using emmeans and hpd.summary

Question

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:

• a way to convert the summary_emm object into an emmGrid object (or any other object that gather_emmeans_draws accepts), OR,
• another function that allows me to draw from the posterior of an emmeans object in the way gather_emmeans_draws does. The function posterior_samples from brms unfortunately does not work in this specific case, since I do not have the contrast of interest in my summary model output, OR,
• another function that allows me to specify pairwise comparisons in the way emmeans does, AND a function that allows me to extract posterior draws of this.

Any ideas are highly appreciated!

Answer 1

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.