I use betareg betaregresion to model degradation percentage over time: perc_degr~time. To calculate half-life (assuming logit link fucntion), I do as follows: th = -coef(model)[1]/coef(model)[2]. I also used the mvrnorm command from the MASS package to do Monte Carlo simulations to get an idea of the distribution of the values and later calculate confidence intervals:

sims_coef = mvrnorm(n=10000, mu=coefs, Sigma=vcovmat, empirical=T)
th_sims = (-sims_coef[,1])/sims_coef[,2]

Now I understood that betareg gives me the precision paraeter phi. How would I incorporate this value to simulate the distribution/confidence intervals for the half-life?
Thanks a lot!

  • $\begingroup$ Could you elaborate on what you mean by "incorporate"? $\endgroup$
    – whuber
    Sep 1, 2020 at 13:39
  • $\begingroup$ @whuber I'd like to use phi to make the predictions more realistic/more appropriate for betaregression.... So far I'm assuming a normal distribution around the coefficients, and I though maybe phi can be used to calculate the CIs of the half-life... $\endgroup$
    – OnLeRo
    Sep 1, 2020 at 14:25
  • $\begingroup$ When you calculate coef(betareg_model) this gives you the full coefficient vector, i.e., including the coefficients for both $\mu$ and $\phi$. And the vcov(betareg_model) gives you the matching covariance matrix. Then you can simulate from the combined coefficient distribution and just compute those predictions you are interested in. So for your th you simply discard the simulated coefficient(s) for $\phi$. $\endgroup$ Sep 1, 2020 at 16:17
  • $\begingroup$ @AchimZeileis OK, thanks. So in my model perc ~ time I'm only interested in coef(betareg_model)[1:2] and vcov(betareg_model)[1:2,1:2] which I pass to mvrnorn to get the distribution, right? $\endgroup$
    – OnLeRo
    Sep 2, 2020 at 8:08

1 Answer 1


The strategy is to simulate the coefficients from the joint normal distribution of all parameters, including the $\phi$ parameter. The reason is that the coefficients for $\mu$ and $\phi$ are not orthogonal, i.e., depend on each other. After having obtained the full simulated parameter vector, you can than compute any quanity you want from the relevant subset of parameters.

I would recommend though to use a log-link for the $\phi$ in this case because the normal approximation will typically be much better on the log-scale because the parameters are unbounded. The easiest way to achieve that is to specify a two-part formula: perc_degr ~ time | 1 in your case, declaring that $\mu$ depends on time but $\phi$ is constant. Moreover, I would encourage you to consider a model with varying $\phi$ though: perc_degr ~ time | time where both parameters are allowed to change over time.

As a worked example, consider the following illustration based on a (too) simple model for the ReadingSkills data available in the package:

data("ReadingSkills", package = "betareg")
m <- betareg(accuracy ~ iq | 1, data = ReadingSkills)

Then we can simulate 10,000 draws from the full parameter vector (using mvtnorm rather than MASS):

cf <- rmvnorm(n = 10000, mean = coef(m), sigma = vcov(m))

And then we can compute any quantity we want. For example, we can compute the negative ratio of coefficients for $\mu$ and compare it with the quantity at the estimated coefficients:

hist(-cf[,1]/cf[,2], main = "")
abline(v = -coef(m)[1]/coef(m)[2], col = 2, lwd = 2)

Histogram ratio

In this particular model this plot looks rather bad because we have to extrapolate quite a bit. Hence I'm also showing another quantity which looks a bit more reasonable on this data set.

Let's also compute the predicted $\mu$ at regressor iq set to -1, 0, and 1. So we set up a suitable 10,000 x 3 matrix and compute the predicted $\mu$ with the inverse link from the fitted model object:

x <- matrix(rep(c(-1, 0, 1), each = nrow(cf)), ncol = 3)
mu <- m$link$mean$linkinv(cf[,1] + cf[,2] * x)
mu0 <- predict(m, newdata = data.frame(iq = c(-1, 0, 1)), type = "response")

The fitted quantity from the original model is also added for comparison:

for(i in 1:3) {
  hist(mu[, i], main = "", xlab = paste("x =", x[1,i]), xlim = c(0.5, 1))
  abline(v = mu0[i], col = 2, lwd = 2)

enter image description here

  • $\begingroup$ Thanks! This helps a lot. In my actual model I already set phi to be variable with time... Also, what's the difference between rmvnorm and mvrnorm? I get very similar values... $\endgroup$
    – OnLeRo
    Sep 4, 2020 at 9:41
  • $\begingroup$ The mvtnorm package has more functionality and is more robust in certain cases. Hence, I would generally prefer that over the older/simpler implementation in MASS. $\endgroup$ Sep 4, 2020 at 10:51
  • $\begingroup$ OK, thanks for the explanation, will do so $\endgroup$
    – OnLeRo
    Sep 4, 2020 at 11:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.