Interpreting coefficients in multiple logistic regression

Question

So I thought I had a handle on interpreting exponentiated coefficients from logistic regression in terms of odds and odds ratios, but I'm having more difficulty than I expected.

I'm playing with some example data I generated in R, in about the simplest possible situation. I made 12 observations of two binary predictors and a binary response:

 df <- data.frame(
   x1 = factor(c("High","High","Low","Low","Low","Low","High","High","High","Low", "Low", "High"), 
   levels = c("Low", "High")), 
   x2 = factor(c("No","No","No","No","No","Yes","Yes","Yes","Yes","Yes", "Yes", "No")),
   y = c(0,1,1,0,1,0,1,0,1,1,1,0))

If I use just one predictor, things work like I expect.

 xtabs(~ y + x1, df)

		x1
		Low	High
y	0	2	3
	1	4	3

 mod1 <- glm(data = df, y ~ x1, "binomial")
 exp(mod1$coefficients)

(Intercept)	x1High
2.0	0.5

I can interpret the exponentiated intercept of 2 as the odds of y = 1 in the reference level of x1 (4/2 = 2), and the exponentiated coefficient 0.5 of x1High as the odds ratio (3/3)/(4/2) = 0.5. So far so good, and I can get the same sort of intelligible results if I only use x2 as a predictor, though I omit those results in the interest of space.

But when I put both x1 and x2 in the model it all goes to hell, and I realize that I must not know what's going on:

xtabs(~ y + x1 + x2, df)

x2 = No
		x1
		Low	High
y	0	1	2
	1	2	1

x2 = Yes
		x1
		Low	High
y	0	1	1
	1	2	2

mod3 <- glm(data = df, y ~ x1 + x2, "binomial")
exp(mod3$coefficients)

(Intercept)	x1High	x2Yes
1.4294	0.4894	2.0433

Now I'm clearly off on the wrong foot, because I naively suppose that the exponentiated intercept should correspond to the odds of y=1 when the predictors are both at their reference levels -- when x1 = Low and x2 = No. But looking at the relevant cells for that condition in the crosstabs, it would appear the odds in that condition are 2/1 = 2, not 1.4294. Changing data the response values of data points that aren't in the reference level of either category also changes the fitted intercept, so that whole interpretation seems straight out the window. And I also can't figure out any way of working with odds ratios that yield either of the coefficients on x1High or x2Yes.

I'd appreciate any pointers about how I can relate any odds or odds ratios that could be manually computed from the above cross-tab to any of these exponentiated logistic regression coefficients. I started from a more complex applied problem and reduced it as far as I could to what's shown here, and I'm temporarily stuck at this point, so I thought I'd ask for some help getting unstuck.

Thanks!

P.S. I did check that my contrasts option for unordered factors is set to the default contr.treatment, like I expected.

Answer 1

It seems you are expecting your stratified cross-tabs and the covariate-adjusted logistic regression model to yield the same implied odds and odds ratios, but they do not contain the same information or make the same assumptions about the data. To reproduce the results of a stratified analysis using logistic regression, you need to include an interaction between x1 and x2 in the model. This makes the model saturated, meaning it will perfectly reproduce the results of the cross-tab. See below:

mod4 <- glm(data = df, y ~ x1*x2, "binomial")
exp(mod4$coefficients)

# (Intercept)       x1High        x2Yes x1High:x2Yes 
#        2.00         0.25         1.00         4.00 

The odds are indeed 2 for being in the reference category of each variable. See my answer here for essentially the same question.