rope range for linear models

[mattansb]

For logistic regression, the rope range is a function of 3/sqrt(pi), as this is the sd of the logistic distribution (on the latent scale).

However, this is the conditional sd, which is akin to sigma in a linear model.
But in a linear model, the rope range if a function of the marginal sd, no?

This seems odd to me. Have I missed something?

[DominiqueMakowski]

mmh I think the rationale we followed back then was quite simple: finding the equivalent of 0.1 * SD, and the most straightforward way was using the d -> log odds conversion formula... but indeed, there might be more appropriate ways of tackling this, what would you recommend?

[mattansb]

I thought this was based on Kruschke's book, no?

[DominiqueMakowski]

afaik he only gave recommendation for linear models and from there we derived the logistic one

[strengejacke]

afaik he only gave recommendation for linear models and from there we derived the logistic one

No, see

(from Supplement https://doi.org/10.1177/2515245918771304)

[strengejacke]

Reminder for ordinal regression: https://twitter.com/ChelseaParlett/status/1352367383243939841?s=19

[DominiqueMakowski]

#20

[mattansb]

It is sigma!

So... "fix" rope_range() with a startup warning?...

[mattansb]

It seems that he also gives the 0.1*sd(y) as another option?

How very weird - these two options give very different results!

[DominiqueMakowski]

From my empirical experience, I felt like the current default rope range for logistic was "a bit larger" than for lms (i.e., it was a bit harder to be outside the rope given similarish effects). Is the new value usually bigger or smaller than the old?

[mattansb]

1. The conditional sd will be smaller than the unconditional sd.

So for the logistic case, the default we use is expected to be rather small - a change of of 0.1*3/sqrt(pi) = 0.18 ~ OR: 1.2 is indeed rather small.

---

2. We have an inconsistency in the default rope range between models:

For linear models we use the unconditional, and for logistic we use conditional:

	Unconditional	Conditional
Logistic (any link = "logit")	--none--	3/sqrt(pi)
Gaussian (really any link = "identity")	Sy = sd(insight::get_response(model))	sigma(model)

[mattansb]

@DominiqueMakowski not 100% sure what you're simulation was doing, but you're (again 😂) using the OR to d conversation wrong - it's not for predicting group from y in a logistic regression, it's for predicting a dichotomized y from group in a logistic regression.

Also, the warnings shouldn't be for logistic models (as I noted, there is no unconditional alternative for 3/sqrt(pi)), but for linear models, where we might change from the current unconditional to the conditional.

[mattansb]

This is what I'm trying to say:

get_results <- function(n=100, d=0.5){
  df <- bayestestR::simulate_difference(n=n, d=d)
  
  m1 <- glm(I(V1>median(V1)) ~ V0, data=df, family="binomial")
  m2 <- lm(V1 ~ V0, data=df)
  
  data.frame(n=c(n, n, n),
             d=c(d, d, d), 
             type=c("LM - current (S)", "LM - sigma", "GLM (rescaled)"),
             rope=c(bayestestR::rope_range(m2)[2], 0.1 * insight::get_sigma(m2), 0.1), # rescaled, 0.1 * pi / sqrt(3) is just 0.1
             coef=c(insight::get_parameters(m2)[2, 2], insight::get_parameters(m2)[2, 2], insight::get_parameters(m1)[2, 2]))
}

library(purrr)
data <- map_dfr(seq(0, 2, length.out = 50), get_results, n = 100)


library(ggplot2)

ggplot(data, aes(x=d, y=rope, color=type)) +
  geom_line(size = 2) +
  coord_cartesian(ylim = c(0.05,0.15))

^{Created on 2021-01-26 by the reprex package (v0.3.0)}

[DominiqueMakowski]

aaaaaah I entirely misread this thread then

[mattansb]

lol

[DominiqueMakowski]

My priors made me think it was a rehash of #20, my deepest apologies 😅

[mattansb]

Classic example of bad priors :P

[mattansb]

I cleaned-up the code a bit, and put the warning in the right place.
We can finish resolving before the next CRAN release.

[DominiqueMakowski]

re-submitting, then

[DominiqueMakowski]

resubmitted ^^

[DominiqueMakowski]

Dear maintainer,

package bayestestR_0.8.2.tar.gz does not pass the incoming checks automatically, please see the following pre-tests:
Windows: https://win-builder.r-project.org/incoming_pretest/bayestestR_0.8.2_20210126_062500/Windows/00check.log
Status: OK
Debian: https://win-builder.r-project.org/incoming_pretest/bayestestR_0.8.2_20210126_062500/Debian/00check.log
Status: 1 NOTE

[mattansb]

Because of the note?? The URL works fine...

[strengejacke]

Do you have an attached text file that also contains another check log?

[DominiqueMakowski]

Lol dont worry I just pasted that here for reference but I'll look into it as soon as I reach soon :)

[DominiqueMakowski]

I'll just reply say it's prolly a false positive :)

[strengejacke]

Or remove link, merge #367, and re-submit 8-)

[DominiqueMakowski]

On CRAN

[bbolker]

help ... it's great that this warning exists (I guess), but it's super-opaque to figure out what to do about it if it pops up (apparently) out of nowhere ...

I'm getting it from running dotwhisker::dwplot() on a stanreg object. dw_tidy is calling broomExtra::tidy_parameters which is eventually making its way down the chain to rope. The message says it's about binomial models but is actually about linear models ...?

reprex:

library(rstanarm)
library(dotwhisker)
fit <- stan_lm(mpg ~ wt + qsec + cyl, data = mtcars, prior = R2(0.75),
               refresh=0, seed=101)
dwplot(fit)

Session info has 111 packages (if you want I can post the whole thing) but this is dotwhisker-0.6.0, rstanarm 2.21.1, bayestestR 0.8.2

Should this be posted as an issue somewhere upstream ... ?

Warning in doTryCatch(return(expr), name, parentenv, handler): Note that the default rope range for binomial models might change in future versions (see #364 set it explicitly to preserve current results.

where 1: rope(x, range = rope_range, ci = rope_ci, ...)
where 2: .describe_posterior(posteriors, centrality = centrality, dispersion = dispersion, 
    ci = ci, ci_method = ci_method, test = test, rope_range = rope_range, 
    rope_ci = rope_ci, bf_prior = bf_prior, BF = BF, effects = effects, 
    component = component, parameters = parameters, ...)
where 3: describe_posterior.stanreg(model, centrality = centrality, dispersion = dispersion, 
    ci = ci, ci_method = ci_method, test = test, rope_range = rope_range, 
    rope_ci = rope_ci, bf_prior = bf_prior, diagnostic = diagnostic, 
    priors = priors, ...)
where 4: bayestestR::describe_posterior(model, centrality = centrality, 
    dispersion = dispersion, ci = ci, ci_method = ci_method, 
    test = test, rope_range = rope_range, rope_ci = rope_ci, 
    bf_prior = bf_prior, diagnostic = diagnostic, priors = priors, 
    ...)
where 5: .extract_parameters_bayesian(model, centrality = centrality, 
    dispersion = dispersion, ci = ci, ci_method = ci_method, 
    test = test, rope_range = rope_range, rope_ci = rope_ci, 
    bf_prior = bf_prior, diagnostic = diagnostic, priors = priors, 
    effects = effects, standardize = standardize, ...)
where 6: model_parameters.stanreg(x, ...)
where 7: parameters::model_parameters(x, ...)
where 8: standardize_names(parameters::model_parameters(x, ...), style = "broom")
where 9: doTryCatch(return(expr), name, parentenv, handler)
where 10: tryCatchOne(expr, names, parentenv, handlers[[1L]])
where 11: tryCatchList(expr, classes, parentenv, handlers)
where 12: tryCatch(expr = standardize_names(parameters::model_parameters(x, 
    ...), style = "broom"), error = function(e) NULL)
where 13: broomExtra::tidy_parameters(x, conf.int = TRUE, ...)
where 14: dw_tidy(x, ci, by_2sd, margins, ...)
where 15: dwplot(lms)

[strengejacke]

The message says it's about binomial models but is actually about linear models ...?

Yes, this was a misnomer fixed in this commit, but not yet in the CRAN version:

e70ec3e#diff-8a9dfe10a23f01c6067d481572ad8f67365347367637e87548e6fa48f1119d54

You can use verbose argument, or better, @IndrajeetPatil do you pass down verbose to rope_range?

[mattansb]

I think it would perhaps be better to have this in the startup message, what do you think?

Also, @bbolker while you're here... Care to chime in? 😅 Should the default rope be based on sd(y) or sigma?

[strengejacke]

I think it would perhaps be better to have this in the startup message, what do you think?

Or we resolve this issue and can then remove the message ;-)

[strengejacke]

The question is whether the conditional "rope" makes more sense than the unconditional one... The idea behind the ROPE is to have a relationship/association/effectsize between predictors and outcome that can be considered as "large enough" to be relevant. From my understanding, "large enough" refers to a significant change in the outcome - i.e. the unconditional rope range (+/- 0.1 * SD(y)).

[DominiqueMakowski]

From my understanding, "large enough" refers to a significant change in the outcome - i.e. the unconditional rope range (+/- 0.1 * SD(y))

I spontaneously tend to agree

[strengejacke]

One could argue that, due to multiple regression, the "raw" effect of a predictor (i.e. unconditional rope range) is not of interest, since coefficients are "adjusted" for each other. That could mean we should use the conditional rope range. However, in real worlds, we always have confounders and no "isolated" associations, but we still want to have a relevant change in our outcome - even if conditioned on / adjusted for other covariates. Then we would still look for +/- 0.1 * SD(y) - unless there's a good reason that the unconditional rope range is biased, and the conditional rope range is preferred.

[mattansb]

I think there are reasons for each of the options (and as you know I think that defaults should be avoided anyway)... But,

My main concern (which is how this whole thing started #364 (comment)) is that for GLMs, we return the "conditional"-based rope, and so there is an issue of consistency if we return the unconditional-rope only for gaussian models.

[strengejacke]

We could keep this inconsistency top secret, delete this issue and, to be sure, the whole internet as well... 😬

[DominiqueMakowski]

should we have it as an option to toggle? Fixed vs. Conditional? In the first case, it would be 0.1 * SD(y) and the equivalent based on the log-odds -> transformation for logit, and in the second we base it on sigma?

[mattansb]

So:

• rope(..., range = "default") == rope(..., range = "conditional")
• rope(..., range = "marginal") = sigma based - would only be for identity-link quantitative models (so basically just linear models).

(There is no equivalent for the marginal option for GLMs... 😢)

[DominiqueMakowski]

The above seems good!

[DominiqueMakowski]

So that one?

Model Type	Marginal	Conditional (default)
Logistic (any link = "logit")	NA	3/sqrt(pi)
Gaussian (any link = "identity")	sigma(model)	Sy = sd(insight::get_response(model))

[mattansb]

Yes. So we should add a note in the startup message that the default rope for linear models is changed.

[DominiqueMakowski]

but currently the default is sd(insight::get_response(model)) so nothing will change?

[mattansb]

Ah, you switched my table around!

Should be this:

Model Type	Marginal	Conditional (default)
Logistic (any link = "logit")	NA	3/sqrt(pi)
Gaussian (any link = "identity")	Sy = sd(insight::get_response(model))	sigma(model)

[DominiqueMakowski]

Oops my bad. Something that will need to be explained and clarified that might bug people's intuition is that for lms, the better the model the lower sigma and the smaller the rope, but for glms the rope doesn't change depending on the model quality, it's always the same value

[mattansb]

That's true - for the same data.

However, for different samples, over different ranges of Xs (say you sample everyone and I only sample children) sigma is (or should be) the same.

I can add a note to the relevant vignette, that you can reference in the docs and startup message. Sounds good?

[DominiqueMakowski]

I can add a note to the relevant vignette, that you can reference in the docs and startup message. Sounds good?

Yeah

But wouldn't 3/sqrt(pi) also apply for the marginal ROPE (for logit)? I'm trying to wrap my head around the implications

[mattansb]

Nope... GLMs don't have a straight forward way of getting the marginal variance of the latent variable 🤷‍♂️

[DominiqueMakowski]

We could keep this inconsistency top secret, delete this issue and, to be sure, the whole internet as well... 😬

I'm starting to think that this is actually a not-so-bad alternative 😁

[DominiqueMakowski]

To be honest, while I understand the crux of this issue (though I still have trouble understanding why 3/sqrt(pi) is conditional and why is it fundamentally different from the marginal linear version), the current behaviour (i.e., "marginal" for linear and "conditional" for logit) is really straightforward to explain and transparent to implement

We discuss(ed) it clearly in the docs, posts and papers, so it's really weird that nobody else had seemed to notice the discrepancy (especially if it leads to some problems)... I mean the Bayesian framework have been under so much scrutiny lately, and ROPEs are even used outside of that context (btw I wonder what would be Lacken's perspective on that)...

This relative absence of agitation surrounding our previous default suggests that it's not a timely issue: we can probablyt run a bit more with it

So I'd suggest adding a method = "legacy" argument, which corresponds to our original (inconsistent) default. We can explain in the docs a bit about the issue and say that additional methods might be implemented in the future (and we link to that discussion). Then, we _{"delete this issue and, to be sure, the whole internet as well"} let this issue sit around for some more time until Mat' has a eureka moment in his bath, realizing that he's been wrong all along 😁 Or until we gather some more thoughts on this

[mattansb]

Yes, I also find it weird that no one has pointed this out. However, the might be because (1) _serious people don't use the default values, or (2) they tend towards BFs more than ROPEs?

I agree that we can let this steam some more - no need to add the method argument at all right now, I think.

---

though I still have trouble understanding why 3/sqrt(pi) pi / sqrt(3) is conditional and why is it fundamentally different from the marginal linear version

This is the definition of a logistic regression:

P(y==1) = 1/(exp(-η)), where η is a latent variable (we only see it's dichotomized binary form) that is conditionally normally distributed ::: η ~ N(bX, pi / sqrt(3)) (where bX is the linear part we usually care about).

This is also why log(OR) * sqrt(3) / pi is "equivalent" to Cohen's d:

• log(OR) is the distance between two groups, standardized by the conditional variance pi / sqrt(3) (on the latent scale).
• M1-M2 is the distance between two groups, standardized by the conditional variance sigma.

---

Here's some more complexity:

In a (G)LMM, the default rope range should probably be conditioned on the level of the effect.

1. Level 1 effects sigma - a within-person effect should be compared to the variance within people, no?
2. Level 2 effects Var(intercept) - a between-person effect should be compared to the variance between people, no?
3. ...

---

And remember kids - don't use defaults!

[DominiqueMakowski]

thanks for the explanations, (and your patience), it's much clearer now ☺️