Closed Dallak closed 2 years ago
I'm thinking about this. But meanwhile, please note the issue I just posted for bayestestR: https://github.com/easystats/bayestestR/issues/559. I think there are serious problems with using rope
on a model, including the example you show.
Also, I question whether ROPE is useful for your example with summary(r3)
, because it is usually not that interesting to test estimated means against zero. Doing so with summary(pairs(r3))
would seem to make more sense.
Thanks @rvlenth for this input. All variables are z-scored in my data. The aim of why I'm using it is that I want to test how robust each contrast is. I know the example I provided with summary(r3, point.est = mean)
is not idea (my apologies), but in my real analysis I'm using your approach, and I was able to indirectly calculate the ROPE. For my analysis, I guess the ROPE is helpful in quantifying the robustness of each contrast. Many thanks again for your time and the helpful discussion on the other post.
Contrast Estimate Q2.5 Q97.5 Rope
1 fin - ini -0.4130980 -0.685896 -0.1401137 0.0000000
2 fin - med -0.0202940 -0.342101 0.3002460 0.4921786
3 ini - med 0.3934325 0.146396 0.6463180 0.0000000
Well, maybe you closed the issue, but I still intend to think about providing equivalence assessment in the summary of a Bayesian model. It probably will not be ROPE exactly. But Ido provide equivalence testing for frequentist models, and it would be good to not leave the Bayesian out.
Russ
Sent from my iPad
On Aug 17, 2022, at 7:06 AM, Abdulrahman Dallak @.***> wrote:
Thanks @rvlenthhttps://github.com/rvlenth for this input. All variables are z-scored in my data. The aim of why I'm using it is that I want to test how robust each contrast is. I know the example I provided with summary(r3, point.est = mean) is not idea, but in my real analysis I'm using your approach, and I was able to indirectly calculate the ROPE. For my analysis, I guess the ROPE is helpful in quantifying the robustness of each contrast. Many thanks again for your time and the helpful discussion on the other post.
Contrast Estimate Q2.5 Q97.5 Rope 1 fin - ini -0.4130980 -0.685896 -0.1401137 0.0000000 2 fin - med -0.0202940 -0.342101 0.3002460 0.4921786 3 ini - med 0.3934325 0.146396 0.6463180 0.0000000
— Reply to this email directly, view it on GitHubhttps://github.com/rvlenth/emmeans/issues/370#issuecomment-1217922068, or unsubscribehttps://github.com/notifications/unsubscribe-auth/AGMJPLYZAM7DYXRZ6BSLFO3VZTIUTANCNFSM56SOKWAQ. You are receiving this because you were mentioned.Message ID: @.***>
FWIW, it seems pretty straightforward to compute ROPE areas with emmGrid
objects, so long as you specify the range
manually. Just use coda::as.mcmc
and (ifnecessary) pick one matrix or rbind them. Here's an example, starting with an equivalence-testing example from a vignette:
library(emmeans)
### Vignette example:
### https://cran.r-project.org/web/packages/emmeans/vignettes/confidence-intervals.html
### (Near the end)
# Frequentist model
pigs.lm = lm(log(conc) ~ source + factor(percent), data = pigs)
EMMF = emmeans(pigs.lm, "source")
test(pairs(EMMF), delta = log(1.25), adjust = "none")
## contrast estimate SE df t.ratio p.value
## fish - soy -0.273 0.0529 23 0.937 0.8209
## fish - skim -0.402 0.0542 23 3.308 0.9985
## soy - skim -0.130 0.0530 23 -1.765 0.0454
##
## Results are averaged over the levels of: percent
## Results are given on the log (not the response) scale.
## Statistics are tests of equivalence with a threshold of 0.22314
## P values are left-tailed
# Comparable Bayesian model
set.seed(123)
library(rstanarm)
## ... messages deleted
pigs.stan = stan_glm(log(conc) ~ source + factor(percent), data = pigs)
## ... messages deleted
EMMB = emmeans(pigs.stan, "source")
post = coda::as.mcmc(pairs(EMMB))
X = do.call(rbind, post) # combine into one matrix
apply(X, 2, bayestestR::rope, range = log(1.25) * c(-1, 1))
## $`contrast fish - soy`
## # Proportion of samples inside the ROPE [-0.22, 0.22]:
##
## inside ROPE
## -----------
## 16.45 %
##
##
## $`contrast fish - skim`
## # Proportion of samples inside the ROPE [-0.22, 0.22]:
##
## inside ROPE
## -----------
## 0.00 %
##
##
## $`contrast soy - skim`
## # Proportion of samples inside the ROPE [-0.22, 0.22]:
##
## inside ROPE
## -----------
## 97.42 %
# Straight calculation of posterior areas
apply(X, 2, \(x) mean(abs(x) < log(1.25)))
## contrast fish - soy contrast fish - skim contrast soy - skim
## 0.18125 0.00150 0.95050
Created on 2022-08-17 by the reprex package (v2.0.1)
Those last results are close to what you get from rope()
, except for some kind of scaling.
What I am inclined to do is allow the delta
argument (that you see in the frequentist example) for Bayesian models, and if non-zero, you get the posterior probability of $Pr(|X| < \delta|$ that we see at the end of the output. That preserves a consistent interface between frequentist and Bayesian models, and gives you something darn close to what you want. Unlike bayestestR::rope()
, you will need to specify the range explicitly -- but in my mind that is not a bad thing anyway. If you want the range that rope()
provides, specify delta = SD/10
where SD is the SD of the response.
This was pretty easy to add. For the example above, I get
> test(pairs(EMMB), adjust = "none", delta=log(1.25))
contrast estimate lower.HPD upper.HPD p.equiv
fish - soy -0.273 -0.384 -0.1626 0.1812
fish - skim -0.403 -0.513 -0.2830 0.0015
soy - skim -0.129 -0.240 -0.0134 0.9505
Results are averaged over the levels of: percent
Point estimate displayed: median
'p.equiv' based on posterior P(|lin. pred.| < 0.2231)
Results are given on the log (not the response) scale.
HPD interval probability: 0.95
> test(pairs(EMMB), adjust = "none", delta=log(1.25), type = "response")
contrast ratio lower.HPD upper.HPD p.equiv
fish / soy 0.761 0.679 0.847 0.1812
fish / skim 0.668 0.593 0.747 0.0015
soy / skim 0.879 0.777 0.977 0.9505
Results are averaged over the levels of: percent
Point estimate displayed: median
'p.equiv' based on posterior P(|lin. pred.| < 0.2231)
Results are back-transformed from the log scale
HPD interval probability: 0.95
Many thanks Russell for this thorough answer! It means a lot to me. I tried your approach, and it seems working and producing what I'm looking for. Much obliged for your time and assistance.
Dear all,
I'm using this great package for calculating the marginal means for a brms model. However, I'm wondering if it possible to calculate also the rope for emmean between (-0.1, 0.1) for all predictors.
Here is an example for what I mean:
This is a model summary for a brms model where rope (last column) was calculated using
rope(model)
from the packagebayestestR
. I would like to use the same thing and report the rope for the following table. Is it currently supported?Thanks in advance!