1

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?

Ben Bolker
  • 34,308
  • 2
  • 93
  • 126
MiMi
  • 73
  • 1
  • 5
  • 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 http://stats.stackexchange.com/questions/105519/standard-error-of-difference-of-estimates/105524#105524 – Ben Bolker Jun 06 '16 at 20:28

1 Answers1

2

You just do something like

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

See ?predict.glm for further details.

Jarle Tufto
  • 7,989
  • 1
  • 20
  • 36
  • 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]http://stats.stackexchange.com/questions/117052/replicating-statas-robust-option-in-r/117064#117064[link] 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 07 '16 at 03: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. – Jarle Tufto Jun 07 '16 at 09: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 07 '16 at 09:25