# Beta regression estimates and confidence intervals on response scale

I'm using GLMMadmb for my beta regression. I'm having a bit of trouble. I understand beta regression uses the logit link function and I know how to get from logits to probability. Here's selected output from the regression:

Call:
glmmadmb(formula = pctrans ~ tree + (1 | tree:trayid), data = BT.data,
family = "beta")

AIC: -1149.8

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.16290    0.15136   -7.68  1.6e-14 ***
tree1       -0.01288    0.21310   -0.06     0.95
tree2       -0.26647    0.21345   -1.25     0.21
tree3       -0.03386    0.21318   -0.16     0.87
tree4        0.00934    0.21285    0.04     0.96
tree5        0.82941    0.21220    3.91  9.3e-05 ***


Converting the estimates for any variable, aside from the intercept, requires adding the intercept estimate to its own estimate before converting from logits to probability.

When I use confint() on the model, I get:

                 2.5 %     97.5 %
(Intercept) -1.4595573 -0.8662370
tree1       -0.4305435  0.4047932
tree2       -0.6848254  0.1518833
tree3       -0.4516802  0.3839700
tree4       -0.4078356  0.4265210
tree5        0.4135036  1.2453123


I'm assuming a similar procedure must be done here? e.g.,

exp(-1.459 + (-0.431)) / (1 + exp(-1.459 + (-0.431)))


After doing so, all of my confidence intervals are lopsided. Is that to be expected? Am I missing something critical?

• You should be able to run emmeans::emmeans(model, ~ tree, type = "response") or something similar to get upper and lower bounds. Oct 18, 2018 at 3:01
• @MarkWhite Thanks! That's exactly what I ended up doing. Oct 19, 2018 at 1:05
• @MarkWhite It seems that for beta and GLMM regression, emmeans uses asymptotic confidence intervals. Any thoughts on that? Obviously I'd rather have more accuracy, but I also think in this case the confidence intervals might not be critically important to the take home messages of my work. Also, how do I find $Cov(\hat \alpha, \hat \beta_1)$? Oct 19, 2018 at 1:38

Indeed as a_statistician said, assuming that there is a vcov() and a coef() method available for objects produced by glmmadmb(), you could do something like this:

X <- rbind(c(1, 1, 0, 0, 0, 0))
V <- vcov(fm)
logit <- c(X %*% coef(fm))
se_logit <- sqrt(diag(X %*% V %*% t(X)))

probs <- plogis(logit)
lower <- plogis(logit - 1.96 * se_logit)
upper <- plogis(logit + 1.96 * se_logit)


where fm is the name of your fitted model.

No, CI can not be manipulated in that way. See Confidence intervals of linear model with several factors. I answered this very similar question today.

In your situation, let $$X_i=$$ tree$$i$$, your model is $$g(\mu) = \alpha + \sum_{i=1}^5\beta_iX_i$$ To get the CI for $$\mu|X_1=1$$, need to get CI for $$\alpha + \beta_1$$ first, which is $$\hat \alpha +\hat\beta_1 \pm \sqrt{Var(\hat \alpha +\hat\beta_1)}$$, where $$Var(\hat \alpha +\hat\beta_1) = Var(\hat \alpha) +Var(\hat\beta_1)+2Cov(\hat \alpha, \hat\beta_1) = 0.15136^2 + 0.21310^2+2Cov(\hat \alpha, \hat\beta_1)$$.

You did not provide $$Cov(\hat \alpha, \hat\beta_1)$$ in your question.

Then convert this CI into response acale by $$\left\{\frac{exp(L)}{1+exp(L)}, \frac{exp(H)}{1+exp(H)}\right\}$$

• Why couldn't OP do inv_logit <- function(x) exp(x) / (1 + exp(x)); lb <- inv_logit(-1.16290 - 0.4305435); ub <- inv_logit(-1.16290 + 0.4047932) where lb is the lower bound and ub is the upper for $y | X_1 = 1$? That is, always take the point estimate of the intercept and use the lower and upper bound for $\beta_1$ to increment or decrement it? Oct 18, 2018 at 2:58
• Let $z=f(x,y)$ be function of parameters $x,y$. You can use $\hat z = f(\hat x, \hat y)$ as the point estimate of $z$, although it is biased if $f$ is not linear. (but asymptotically unbiased). But for CI, $z_l=f(x_l,y_l)$ and $z_u=f(x_u,y_u)$ are wrong. Oct 18, 2018 at 3:05
• I create an example: Let $X\sim N(10,4)$ and $Y\sim N(10,4)$ and independent with each other. We have 95% of $X$ will be in (10-2*1.96, 10+2*1.96), smae as $Y$. Let $Z=X+Y$, can we use (20-4*1.96, 20+4*1.96) cover 95% of $Z$? We know that $Z\sim N(20,8)$, (20-1.96*$\sqrt 8$, 20+1.96*$\sqrt 8$) covers 95% of $Z$. Obviously, (20-4*1.96, 20+4*1.96) is too wide. Oct 18, 2018 at 3:33
• Yup, I did a little simulation and found this to be the case: gist.github.com/markhwhiteii/cb5635d75589dcbedc48859b2a1f3dd4 Oct 18, 2018 at 14:55
• Actually, is this more due to confint() using the profile likelihood method for CIs instead of the typical confint.default() method that uses the standard 1.96 * the standard error? Oct 18, 2018 at 15:45