saudiwin / ordbetareg_pack

Repository for R package ordbetareg, used to fit continuous data with lower and upper bounds.
Other
17 stars 3 forks source link

Trouble Specifying No Random Intercept #8

Closed AdamChapnik1 closed 1 year ago

AdamChapnik1 commented 1 year ago

I have been trying to specify a model formula with ordbetareg() that does not include an intercept by using binary dummies explicitly in the formula. I have found a workaround by running through brm() directly, where family = ord_beta_reg but this is not ideal.

Here is an example (I am aware that this distribution isn't actually ordered beta, but it is faster to do it this way and it shouldn't cause any errors). I want, ideally, a formula like bf(Y ~ 0 + X1 + X2).

Setup:

library(ordbetareg) library(brms)

Data:

Y = c(rbeta(900, 1.5, 5), rbinom(100, 1, 0.15)) X1 = rbinom(1000, 1, 0.4) X2 = 1 - X1 df = cbind(Y, X1, X2) %>% as.data.frame()

This is what I would like to do:

fit1 <- ordbetareg( bf(Y ~ 0 + X1 + X2), data = df, iter = 500, warmup = 100, chains = 1, cores = 1 )

Note that the output of fit1 includes a coefficient for the intercept. This runs, but if you include a prior like extra_prior = set_prior("normal(0.1, 0.1)", class = "b") (which I would like to do) then there is an error about specifying duplicate priors.

I am aware that the documentation for ordbetareg() says "Please avoid using 0 or Intercept in the formula definition." I figured this might just be a feature of the function, so I tried just including the dummy in the intercept and specifying a prior for the intercept.

fit2 <- ordbetareg( bf(Y ~ 1 + X1), data = df, extra_prior = set_prior("normal(0.1, 0.1)", class = "Intercept"), iter = 500, warmup = 100, chains = 1, cores = 1 )

This returns an error: Error in nlist(prior, class, coef, group, resp, dpar, nlpar, lb, ub, check) : object 'Intercept' not found. If you run the same line without the extra_prior = argument, there is no issue, and the intercept coefficient is included in the output.

This is my workaround, which works without any issues:

fit3 <- brm( bf(Y ~ 0 + X1 + X2), data = df, family = ord_beta_reg, prior = set_prior("normal(0.1, 0.1)", class = "b"), iter = 1000, warmup = 500, chains = 1, cores = 1, stanvar = stanvar )

saudiwin commented 1 year ago

So the issue is that you're trying to zero out the Intercept?

AdamChapnik commented 1 year ago

Yes. It seems like the default intercept for ordbetareg() is the issue, because even fit2 in my example shows the formula actually being used has a redundant term for the intercept being added.

saudiwin commented 1 year ago

Hi - I have just pushed code that should allow you to put any prior on the intercept:

library(ordbetareg)

Y = c(rbeta(900, 1.5, 5), rbinom(100, 1, 0.15))
X1 = rbinom(1000, 1, 0.4)
X2 = 1 - X1
df = cbind(Y, X1, X2) %>% as.data.frame()

fit1 <- ordbetareg(
  bf(Y ~ 0 + X1 + X2),intercept_prior_mean =0,
  intercept_prior_SD=0.01,
  data = df,
  iter = 500, warmup = 100,
  chains = 1, cores = 1
)
summary(fit1)
AdamChapnik commented 1 year ago

Great, thank you!

On Fri, Feb 3, 2023 at 8:24 AM Robert Kubinec @.***> wrote:

Hi - I have just pushed code that should allow you to put any prior on the intercept:

library(ordbetareg)

Y = c(rbeta(900, 1.5, 5), rbinom(100, 1, 0.15)) X1 = rbinom(1000, 1, 0.4) X2 = 1 - X1 df = cbind(Y, X1, X2) %>% as.data.frame()

fit1 <- ordbetareg( bf(Y ~ 0 + X1 + X2),intercept_prior_mean =0, intercept_prior_SD=0.01, data = df, iter = 500, warmup = 100, chains = 1, cores = 1 ) summary(fit1)

— Reply to this email directly, view it on GitHub https://github.com/saudiwin/ordbetareg_pack/issues/8#issuecomment-1415867739, or unsubscribe https://github.com/notifications/unsubscribe-auth/ARKDHL42OXYWYXFR3JMMG3LWVUBKTANCNFSM6AAAAAAQRFNYHY . You are receiving this because you commented.Message ID: @.***>