0
$\begingroup$

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?

$\endgroup$
3
  • 1
    $\begingroup$ You should be able to run emmeans::emmeans(model, ~ tree, type = "response") or something similar to get upper and lower bounds. $\endgroup$
    – Mark White
    Oct 18, 2018 at 3:01
  • $\begingroup$ @MarkWhite Thanks! That's exactly what I ended up doing. $\endgroup$ Oct 19, 2018 at 1:05
  • $\begingroup$ @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)$? $\endgroup$ Oct 19, 2018 at 1:38

2 Answers 2

3
$\begingroup$

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.

$\endgroup$
2
$\begingroup$

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\}$

$\endgroup$
6
  • $\begingroup$ 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? $\endgroup$
    – Mark White
    Oct 18, 2018 at 2:58
  • 1
    $\begingroup$ 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. $\endgroup$
    – user158565
    Oct 18, 2018 at 3:05
  • 1
    $\begingroup$ 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. $\endgroup$
    – user158565
    Oct 18, 2018 at 3:33
  • $\begingroup$ Yup, I did a little simulation and found this to be the case: gist.github.com/markhwhiteii/cb5635d75589dcbedc48859b2a1f3dd4 $\endgroup$
    – Mark White
    Oct 18, 2018 at 14:55
  • $\begingroup$ 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? $\endgroup$
    – Mark White
    Oct 18, 2018 at 15:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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