Beta regression: Monte Carlo simulations for coefficients

Question

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!

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:

library("betareg")
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):

library("mvtnorm")
set.seed(1)
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)

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)
}