# Extracting Standard Errors for a combination of factorial predictors in binomial GLM

Suppose I run a binomial GLM (in R) with response variable [0,1] and 2 predictor variables that are both categorical. Let's call them a and b where a has 3 factor levels (a1,a2,a3) and b has 2 factor levels (b1,b2).

Therefore:

mod <- glm(y~a+b, family=binomial(logit),data=pretend)


Let's say the summary output of the coefficients looks like this:

Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept)        -0.62049    0.06988  -8.879  < 2e-16 ***
a a2                -2.24304    0.16111 -13.923  < 2e-16 ***
a.a3                -1.76965    0.14147 -12.509  < 2e-16 ***
s.b2                -0.79545    0.07918 -10.046  < 2e-16 ***


So I understand a1 and b1 are constrained in the intercept and the estimate values presented are log-odds ratios (I am familiar with interpreting them).

Let's say now I want to determine the estimate when a2 is true or = 1 and b1 is =1 (true), then I would just calculate it as B0 (intercept) + B1 (coeff of a2) = -0.62049 + -2.24304

Or if a3 is true and b2 then I would calculate it as B0 (intercept) + B3 (coeff of a3) + B4(coeff of b2). = -0.62049 + -1.76965 + -0.79545

But now how do I calculate the Standard errors and/or 95% Confidence intervals for the last two scenarios? Can I just add the standard errors in the same way? Or I would think I leave out the SE's of B0 (the intercept)?

So that in my first scenario where a2 is true and b1 is true, then I would just calculate it as SE of B1 (SE of a2) = 0.16111

And in my second scenario where a3 is true and b2 then I would calculate it as SE of B3 (SE of a3) + SE of B4(SE of b2). = 0.14147 + 0.07918.

Am I correct in my assumptions of interpreting the model output or am I completely wrong here?

• you're completely wrong (but this is a fine question!) You can't add standard errors, you have to add variances (and you have to take the covariances into account as well). @Jarle's answer is correct; you can also compute this directly along the lines of stats.stackexchange.com/questions/105519/… Jun 6, 2016 at 20:28

## 1 Answer

You just do something like

predict(mod,newdata=data.frame(a="a2",b="b1"),se.fit=TRUE)


See ?predict.glm for further details.

• Thank you. I did not want to use the 'predict' function because my response variable is a proportion value between 0 and 1 (fractional logit model) and I am calculating robust standard errors as suggested here: [link]stats.stackexchange.com/questions/117052/… Therefore, I am calculating my SE's after the fact, separate from the model using function 'coeftest' in package 'sandwich'. I am not sure how to bring those SE's produced together with the predict function or using the formulas supplied by @BenBolker's post.
– MiMi
Jun 7, 2016 at 3:40
• If you want the standard error of $p=logit^{-1}(\eta)$ instead of at the scale of the linear predictor $\eta$, just add the argument type="response" to predict. Jun 7, 2016 at 9:20
• No I am calculating different types of SE's - heteroskedastic robust SE's - which cannot be solved by adding 'type=response'. Thank you though.
– MiMi
Jun 7, 2016 at 9:25