easystats
Revise the way we calculate % in ROPE (SGPV) for frequentist models?
This idea came up after a chat with Isabella on Twitter.
See following example at https://easystats.github.io/parameters/reference/equivalence_test.lm.html
library(parameters)
data(qol_cancer)
model <- lm(QoL ~ time + age + education, data = qol_cancer)
equivalence_test(model)
#> # TOST-test for Practical Equivalence
#>
#> ROPE: [-1.99 1.99]
#>
#> Parameter | 90% CI | % in ROPE | H0 | p
#> ------------------------------------------------------------------
#> (Intercept) | [59.33, 68.41] | 0% | Rejected | > .999
#> time | [-0.76, 2.53] | 83.52% | Undecided | 0.137
#> age | [-0.26, 0.32] | 100% | Accepted | < .001
#> education [mid] | [ 5.13, 12.39] | 0% | Rejected | 0.999
#> education [high] | [10.14, 18.57] | 0% | Rejected | > .999
If we look at the parameter time, we see a CI of (-0.76, 2.53) and 83.52% inside ROPE. The latter is equivalent to the SGPV. Both SGPV and % in ROPE are calculated the same way (see GitHub repo from Daniel: https://github.com/Lakens/TOST_vs_SGPV/tree/master/functions, related paper https://open.lnu.se/index.php/metapsychology/article/view/933).
However, although we actually want the overlap of the interval (probability mass) with the ROPE, it's calculated (for simplicity) based on the overlap of the range.
What we are actually doing is:
# hypothetical min/max
est.hi <- 2.53
est.lo <- -0.76
x <- bayestestR::distribution_uniform(
n = 500,
min = est.lo,
max = est.hi
)
# actual (simulated) min/max
min(x)
#> [1] -0.75671
max(x)
#> [1] 2.52671
bayestestR::equivalence_test(x, c(-2, 2), ci = 1)
#> # Test for Practical Equivalence
#>
#> ROPE: [-2.00 2.00]
#>
#> H0 | inside ROPE | 100% HDI
#> --------------------------------------
#> Undecided | 83.80 % | [-0.76 2.53]
bayestestR::rope(x, c(-2, 2), ci = 1) |> plot()
But what we actually want, is:
# hypothetical min/max
est.hi <- 2.53
est.lo <- -0.76
x <- bayestestR::distribution_normal(
n = 500,
mean = est.hi - ((est.hi - est.lo) / 2),
sd = (est.hi - est.lo) / 6.1
)
# actual (simulated) min/max
min(x)
#> [1] -0.7816991
max(x)
#> [1] 2.551699
bayestestR::equivalence_test(x, c(-2, 2), ci = 1)
#> # Test for Practical Equivalence
#>
#> ROPE: [-2.00 2.00]
#>
#> H0 | inside ROPE | 100% HDI
#> -------------------------------------
#> Accepted | 98.00 % | [-0.78 2.55]
bayestestR::rope(x, c(-2, 2), ci = 1) |> plot()
So my question is, to compute the % in ROPE for frequentist models, should we switch to simulating the CI and calculate the % in ROPE based on the actual probability mass of that interval, than just the range? Wouldn't this be more precise? Or am I missing something here?
@easystats/core-team
Tagging @vincentarelbundock too, maybe this is of interest for you, too?
However, although we actually want the overlap of the interval (probability mass) with the ROPE, it's calculated (for simplicity) based on the overlap of the range.
We've discussed this elsewhere I think, and concluded that giving "posterior" probabilities estimates from MLE models was a bad idea? Why can't we just have the p-value from the equivalence test?
Not sure we discussed it in this context, though. But we may decide on dropping the % in rope completely, but I also see that it might be useful for frequentist case.
If it's so useful, might as well switch to a framework that can provide with posterior probabilities 😉
I'm not following after a quick read but will circle back.
…perennial Brenton refrain about confidence distributions and frequentist posterior distributions 😜…
I don't have a strong view.
…perennial Brenton refrain about confidence distributions and frequentist posterior distributions 😜…
I read that book that you recommend and liked it a lot. Thanks for that!
Although, this kind of interpretation feels very attractive, I will say that my impressionistic sense is that it seems a bit "exotic" in the stats community. For software defaults and labels, I'd lean toward more traditional/conservative language, or at least use labels which make it clear that the suggested interpretation departs somewhat from canon.
I think labeling is more a side topic here. ;-) Currently, we report the % in ROPE, which is equivalent to the SGPV (which, in turn, is proposed by strong proponents of the Neyman-Pearson philosophy). However, the % in ROPE usually refers to the probability mass of the interval distribution that is inside the ROPE. Currently, we just divide the ranges:
ci <- c(-0.76, 2.53)
# rope <- c(-2, 2)
overlap <- c(-0.76, 2)
100 * diff(overlap) / diff(ci)
#> [1] 83.89058
That is, we consider the CI as a uniform distribution with equal density for all values inside this interval. Usually, you would consider the CI as roughly normal distribution, so the % in ROPE of that interval cannot be calculated as if it was uniformly distributed. Thus, my suggestion to simulate a normal distribution (based on the lower/upper limit) and then calculate the proportion of that probability mass that is inside the ROPE range. This would no longer be in line with the proposed calculation of the SGPV. But it would give a more realistic picture of the "% in ROPE".
We could, by default, just report the p-value of the equivalence test and then add columns like SGPV and % in ROPE when requested by the user, if these measures are too avantgarde for some of you to be shown by default :-P
if these measures are too avantgarde for some of you to be shown by default :-P
Haha! again, don't listen to me much here. I don't really use this function and don't have a strong view.
I think this can work... I have implemented a % in ROPE based on a simulated normal distribution of the CI. The "new" results from % inside ROPE are very close to what we would get from a bootstrapped model or a model where we simulate estimates. We now have the SGPV, calculated in the "classical" way as proposed by Lakens, and the % inside ROPE, which assumes a normally distributed confidence interval and where the proportion of the probability mass of this distribution inside the ROPE is shown. See the close results from simulated model and normal model.
library(parameters)
data(qol_cancer)
model <- lm(QoL ~ time + age + education + phq4 + hospital, data = qol_cancer)
equivalence_test(model, metric = "all")
#> # TOST-test for Practical Equivalence
#>
#> ROPE: [-1.99 1.99]
#>
#> Parameter | 90% CI | SGPV | % in ROPE | Equivalence | p
#> -----------------------------------------------------------------------------
#> (Intercept) | [56.21, 69.05] | < .001 | 0% | Rejected | > .999
#> time | [-0.18, 2.52] | 0.802 | 85.68% | Undecided | 0.160
#> age | [-0.40, 0.07] | > .999 | 100% | Accepted | < .001
#> education [mid] | [ 3.12, 9.12] | < .001 | 0% | Rejected | 0.988
#> education [high] | [ 3.92, 11.06] | < .001 | 0% | Rejected | 0.994
#> phq4 | [-5.75, -4.70] | < .001 | 0% | Rejected | > .999
#> hospital | [-1.52, 8.95] | 0.335 | 27.05% | Undecided | 0.743
x <- simulate_model(model)
equivalence_test(x, ci = 0.9)
#> # Test for Practical Equivalence
#>
#> ROPE: [-1.99 1.99]
#>
#> Parameter | H0 | inside ROPE | 90% HDI
#> -------------------------------------------------------
#> (Intercept) | Rejected | 0.00 % | [55.86 69.00]
#> time | Undecided | 88.22 % | [-0.19 2.65]
#> age | Accepted | 100.00 % | [-0.38 0.09]
#> educationmid | Rejected | 0.00 % | [ 3.00 9.08]
#> educationhigh | Rejected | 0.00 % | [ 3.96 11.05]
#> phq4 | Rejected | 0.00 % | [-5.72 -4.74]
#> hospital | Undecided | 25.78 % | [-1.37 9.19]
model <- glm(vs ~ wt + mpg + cyl, data = mtcars, family = "binomial")
equivalence_test(model, metric = "all")
#> # TOST-test for Practical Equivalence
#>
#> ROPE: [-0.18 0.18]
#>
#> Parameter | 90% CI | SGPV | % in ROPE | Equivalence | p
#> ------------------------------------------------------------------------
#> (Intercept) | [-14.74, 30.58] | 0.008 | 0.95% | Undecided | 0.991
#> wt | [ -0.59, 5.11] | 0.064 | 3.79% | Undecided | 0.964
#> mpg | [ -0.48, 0.61] | 0.332 | 42.74% | Undecided | 0.593
#> cyl | [ -5.25, -0.34] | < .001 | 1.68% | Rejected | 0.983
x <- simulate_model(model)
equivalence_test(x, ci = 0.9)
#> # Test for Practical Equivalence
#>
#> ROPE: [-0.18 0.18]
#>
#> Parameter | H0 | inside ROPE | 90% HDI
#> ------------------------------------------------------
#> (Intercept) | Undecided | 0.89 % | [-15.52 30.46]
#> wt | Undecided | 4.22 % | [ -0.33 5.15]
#> mpg | Undecided | 44.22 % | [ -0.48 0.63]
#> cyl | Rejected | 0.00 % | [ -5.46 -0.48]
# for predictions
model <- lm(QoL ~ time + age + education + phq4 + hospital, data = qol_cancer)
ggeffects::ggpredict(model, "education") |>
equivalence_test(metric = "all")
#> # TOST-test for Practical Equivalence
#>
#> ROPE: [-1.99 1.99]
#>
#> education | Contrast | 95% CI | SGPV | % in ROPE | Equivalence | p
#> ---------------------------------------------------------------------------------
#> low-mid | -6.12 | [ -9.69, -2.56] | < .001 | 0.42% | Rejected | 0.988
#> low-high | -7.49 | [-11.74, -3.24] | < .001 | 0% | Rejected | 0.994
#> mid-high | -1.37 | [ -4.59, 1.86] | 0.596 | 61.05% | Undecided | 0.353