Question on betabinomial family

[jsocolar]

posterior_predict never yields (just) the linear predictor including random effects, it yields that plus a simulated residual.

Do you have any specific questions about the documentation for the transform argument of posterior_linpred (see ?brms::posterior_linpred)?

[jd_c]

@ssalimi

@jd_c yes, exactly. Different observation have different size and incident and proportions should be different. But, posterior linpred or dividing on soze doesn’t provide the proportion! I appreciate any hints to get real proportion. Thanks

@ssalimi

repeated measures have different size and incident so the proportions shouldn’t be the same.

@ssalimi

I am still confused why posterior_linpred for intercept with various incident and trials in different time-point provides a fixed number.

Hi @ssalimi . Based on these comments, I think you might be a little confused regarding what you expect the model to give you in terms of predictions. @jsocolar has outlined why the predictions are the same. Here is some code below that might help clear up the confusion.

library(brms)
data("cbpp", package = "lme4")
cbpp$obs <- 1:length(cbpp$size)

#intercept only model
m1 <- brm(incidence | trials(size) ~ 1, data=cbpp, family=beta_binomial, cores=4)

#varying intercept for herd
m2 <- brm(incidence | trials(size) ~ 1 + (1|herd), data=cbpp, family=beta_binomial, cores=4)

#varying intercept for every observation
m3 <- brm(incidence | trials(size) ~ 1 + (1|obs), data=cbpp, family=beta_binomial, cores=4)

#predicted proportions from intercept only model - all proportions are the same for all rows
p1 <- posterior_linpred(m1, transform=TRUE)
p1_table<-as.data.frame(posterior_summary(p1))
head(p1_table)

#predicted proportions from the model with a varying intercept for herd - all proportions within each herd are the same for all rows within each herd
p2 <- posterior_linpred(m2, transform=TRUE)
p2_table<-as.data.frame(posterior_summary(p2))
head(p2_table)

#predicted proportions from the model with a varying intercept for every observation in the data - a different proportion for each row of data
p3 <- posterior_linpred(m3, transform=TRUE)
p3_table<-as.data.frame(posterior_summary(p3))
head(p3_table)

#the observed proportions in the data
cbpp$prop <- cbpp$incidence/cbpp$size
head(cbpp$prop)

As you can see when you run the code, if you include only an intercept, then all the predicted proportions will be the same across all observations. When you include the varying intercept for herd, the proportions will be the same across observations within the same herd but vary between herds. The only way to get a different predicted proportion for every row of data (which is what I think you were expecting the model to do, based on your comments) is to include a varying intercept for every row of data. However, this would likely be a much too flexible model (and as you can see from all of the warnings, it does not run well).

[ssalimi]

@jd_c Thank you for the nice summary. I agree with you, and I also ran m3. I ran it using more channels and iterations and in my data, it worked well using loo and pp_check.
The point is with posterior_predict and other models than beta_binomial I could get different values for each individual. lin_pred (summary=TRUE) provides a linear relationship I assume. I am concerned that it is enough information for a model. I am considering using the random intercept model (m3) and getting the fitt value and then comparing two models with loo_compare. Thanks for all the hints. Or figure out how I can apply posterior_predict to get random intercept proportions for each individual. Thanks for all.

[jsocolar]

@ssalimi

The point is with posterior_predict and other models than beta_binomial I could get different values for each individual.

posterior_predict always gives values from the posterior predictive distribution–the linear predictor transformed via the inverse link plus a simulated residual. It works just the same for the beta binomial family.

@ssalimi

I am considering using the random intercept model (m3) and getting the fitt value and then comparing two models with loo_compare.

There is no way that loo_compare will be able to handle m3; unmarginalized observation-level random effects are just too flexible for PSIS-LOO. Note also that the beta-binomial already contains an observation-level random effect; it’s just that this random effect is:

1. beta-distributed on the identity scale, rather than Gaussian on the link scale [technically I guess we’d say that its distributed as a shifted beta on the identity scale, centered around its own mean]
2. marginalized out.

If you want to recover the posterior distribution for the identity-scale linear predictor, inclusive of this beta-distributed random effect, the procedure for doing so is

@jsocolar

This calculation refers to the proportion #2 mentioned above. We can generate a posterior distribution for this proportion ppp by sampling from the following probability density function iteration-wise over aaa and bbb which are fitted parameters in the model:

\frac{\beta(p | a, b)\times B(y | p, t)}{\int_0^1\beta(p | a, b)\times B(y | p, t)dp}β(p|a,b)×B(y|p,t)∫10β(p|a,b)×B(y|p,t)dp\frac{\beta(p | a, b)\times B(y | p, t)}{\int_0^1\beta(p | a, b)\times B(y | p, t)dp}, where aaa and bbb are the fitted parameters of the beta distribution, \betaβ\beta is the beta PDF, yyy is the observed count, ttt is the number of trials in the binomial distribution, and BBB is the binomial PMF.

Alternatively, you could switch to a model that is binomial with an observation-level random effect that is Gaussian on the link scale (i.e. switch the family to binomial and include the observation level random effect in the formula). Then posterior_linpred(transform = TRUE) (not posterior_predict) would

@ssalimi

get random intercept proportions for each individual.

based on the Gaussian random intercept.

Unless you have a very good reason for believing that neither the beta-binomial nor the binomial-with-Gaussian-random-effect is sufficiently flexible by itself, I think it is very unlikely that the additional flexibility of m3, which combines the beta-binomial with the Gaussian random effect, would be useful in practice.

[jsocolar]

In summary:

• If you want predictions of potential new data that shares covariates with the actual data, use posterior_predict. You can convert these to proportions by dividing through by the number of trials. This will give predictions of hypothetical observed sample proportions from new data, but not estimates of the true underlying population proportions (well, not good estimates anyway).
• If you want estimates of the true underlying population proportions, then:
  i. Do not use posterior_predict.
  ii. Begin by formalizing your ideas about how these proportions vary. For observations with identical covariates, are the modeled proportions (the underlying population-level proportions, not the observed sample proportions) supposed to be the same?
  iii. Since you do not expect or want these proportions to be the same, you need a model for how the proportions vary, even when covariates are identical. In this thread, we have four candidate models for how the proportions could vary. 1) The observation-level sampling distribution could be binomial, and the between-observation variation in the proportion could be beta. This is the beta-binomial model. 2) The observation-level sampling distribution could be binomial, and the between-observation variation in the proportion could be Gaussian on the logit scale. This is the binomial model with a Gaussian random effect. 3) The observation-level sampling distribution could be beta-binomial, and the between-observation variation in the proportion could be Gaussian on the logit scale. This is m3, treating the beta part as observation-level variation. 4) The observation-level sampling variation could be binomial, and the the between-observation variation in the proportion could be lGaussian-beta (where lGaussian refers to a variable whose logit is Gaussian). This is m3 treating the beta part as between-observation variation.

Each of these assumptions requires a different technique for fitting the model and/or extracting posterior estimates of the population proportions that underlie the observed observation-level samples. And some of them, particularly the latter two (based on m3) are so flexible that they might be difficult to fit.

[ssalimi]

@jsocolar @jd_c Thank you both for patiently discussing all the points. Now I am confident that m3 model with posterior_linpred(transform=TRUE) provides what I am looking for. I appreciate all hints and points. It has been quite interesting.