Logistic regression: Strange standard errors from glm() in R

To my surprise I found that standard errors and thus Wald confidence intervals became smaller when I removed the intercept from a simple logistic regression model, using glm() and R.

# load an object named "my.df" to the global enviroment
# fit a model with intercept to data
my.fit <- glm(deprived.of.education ~ religion, data = my.df, family = binomial("logit"))
# fit a model without any intercept to data
my.fit.without.intercept <- glm(deprived.of.education ~ 0 + religion, data = my.df, family = binomial("logit"))

# inspect the first fit
summary(my.fit)$coefficients # Estimate Std. Error z value Pr(>|z|) # (Intercept) -2.8718056 0.03175130 -90.44687 0.000000e+00 # religionChristianity 0.4934891 0.03234887 15.25522 1.519805e-52 # religionHinduism 0.5257316 0.03376535 15.57015 1.161317e-54 # religionIslam 1.5734832 0.03231692 48.68914 0.000000e+00 # religionNonreligious 1.5975456 0.03555164 44.93592 0.000000e+00 # inspect the second fit summary(my.fit.without.intercept)$coefficients
#                       Estimate  Std. Error    z value Pr(>|z|)
# religionBuddhism     -2.871806 0.031751299  -90.44687        0
# religionChristianity -2.378317 0.006189045 -384.27842        0
# religionHinduism     -2.346074 0.011487113 -204.23530        0
# religionIslam        -1.298322 0.006019850 -215.67354        0
# religionNonreligious -1.274260 0.015992939  -79.67642        0


I understand why the z values are different, because the null hypotheses in the two cases are different. In the first case, with the intercept, the null is "same as the reference category", while without the intercept, the null becomes "zero".
But I do not understand the large difference in standard errors between the two models.

Without the intercept, the standard errors seem to vary with n of each level, i.e. there are many cases of "Christianity" and "Islam", and they have small standard errors, but with the intercept, there is essentially no variation in the standard errors.

Could someone please explain the reason for the differences in the magnitude of the standard errors between the two models?

I would like to calculate probabilities and confidence intervals around them, and I have done so using the estimates from the first model. If I would do that with the estimates from the second model, the confidence intervals would be much smaller, but would they be reliable?

Your coefficients, even when they share common names, are not the same, i.e. their interpretation is different.

In the first model, the effect of religionChristianity is a variation in the outcome wrt the baseline (religionBuddhism), a relative variation. In the second model the effect of religionChristianity is an absolute variation.

The effects are numerically equal, $-2.8718056+0.4934891=-2.378317+5e-07$, but in the first case the effect is a sum of two effects, i.e. you should compare the joint significance of Intercept and religionChristianity in the first model with the significance of religionChristianity in the second one. You should compare a joint confidence interval (first model) with a simple one (second model).

The simple CI for religionChristianity is:

> confint(my.fit.without.intercept)
Waiting for profiling to be done...
2.5 %    97.5 %
...
religionChristianity -2.390467 -2.366207


There are several ways to compute joint intervals. Using arm:

> library(arm)
> n.sims <- 1000
> sim.i <- sim(my.fit, n.sims)
> intercept.plus.christianity <- sim.i@coef[,1] + sim.i@coef[,2]
> quantile(intercept.plus.christianity, c(0.025, 0.975))
2.5%     97.5%
-2.390826 -2.366828


Can you see any significant (relevant) difference?

• No, I don't. Thanks for the code, it gives a practical way to get to my goal, and it proves that both sets of standard errors are correct. However, I would like to have the formula too. For the second, non-intercept model, getting eg. the lower bound for "Christianity" is simple enough: $-2.390448 = -2.378317 - 1.96 \times 0.006189045$. My naive idea was to create the "combined" interval for the first model by $-2.8718056 + 0.4934891 - 1.96 * 0.03234887$, but that gave a much larger confidence interval. – Hans Ekbrand May 22 '14 at 20:49
• Your answer is useful, and I am thankful for it, but to close the question I would like to see a formula that arrives at the lower bound for "Christianity" from the estimates in the intercept model. – Hans Ekbrand May 22 '14 at 20:54
• Sorry, but... why? The models are numerically equivalent (this is what I wanted to highlight), but statistically different, they address different scientific questions. If you are interested in "relative" effects you have to a) choose a baseline (Buddhism or something else, but you shoud explain why one or another) and b) run the first model and not bother to run the second one. If you are interested in "absolute" effects (no baseline) you should run the second model and not bother to run the first one. IMHO you can and should compare several models only if they address the same question. – Sergio May 22 '14 at 21:22
• Because I want to learn. What brought me into this was the desire to understand when I could do without an intercept term. Do you imply that there is no formula for me to understand in this case? Is sampling the only way? – Hans Ekbrand May 22 '14 at 21:43
• Joint confidence "intervals" are actually confidence regions (see, e.g., here and here), this is why sampling is so convenient. – Sergio May 22 '14 at 22:04