Obtaining stratified odds ratios of logistic regression with two interaction terms

[Martin Mueller]

Dear Statalist Forum,

I have a dataset of a case/control study (coded: 1 for cases), with effect modifiers sex (coded: 1 for male) and age_gp (coded: 15, 26, 40, 65) and a confounder (coded: 1,2,3,4,5); see below a random sample of the dataset with dataex command.

 * Example generated by -dataex-. To install: ssc install dataex clear input byte(case bin_outcome) float(sex_m1 age_gp) byte confounder 1 1 0 26 3 0 0 0 65 2 0 0 1 40 3 0 0 0 65 1 0 0 1 15 3 0 1 0 40 2 0 0 1 65 1 0 0 0 15 3 0 0 0 65 2 0 0 1 26 3 end

I ran the following logistic regression:

logistic bin_outcome case##age_gp case##sex_m1 i.confounder, base

I was wondering if the way and my interpretation how I got the stratified odds ratios is correct:

// Odds ratio for cases vs. controls in age-group 15-25 for females adjusted for triage for bin_outcome
lincom 1.case+ 1.case#15.age_gp + 1.case#0.sex

// Odds ratio for cases vs. controls in in age-group 65-120 for males adjusted for triage
lincom 1.case+ 1.case#65.age_gp + 1.case#1.sex

Thank you for your expertise in advance.
Kind regards
Martin

[Martin Mueller]

Sorry.

* Example generated by -dataex-. To install: ssc install dataex
clear
input byte(case bin_outcome sex_m1) float(age_gp) byte triage
0 1 0 26 3
0 0 0 65 2
0 0 1 40 3
0 0 0 65 1
0 0 1 15 3
0 1 0 40 2
0 0 1 65 1
0 0 0 15 3
0 0 0 65 2
0 0 1 26 3
end

[Martin Mueller]

Just to close this post. Here is the correct answer and a beautiful derivation of the correct lincom command (Thank you Clyde Schechter)!

Run the following logisitic regression:

logistic bin_outcome case##age_gp##sex_m1 i.confounder, base

An extra explicit age_gp#sex_m1 term in the regression is not needed because having used the ## interaction operator in the regression, all of the lower-level interactions among case, age_gp and sex_m1, including age_gp#sex_m1, are automatically generated. If you look at the Stata output from your logistic regression command, you will see that it's there.

To calculate the case:non-case odds ratio we can do the numerator and denominator separately.

For a non-case, E(log odds outcome) = _cons + 15.age_gp + 0.sex + 15.age_gp#0.sex + *. For a case, E(log odds outcome) = _cons + 15.age_gp + 0.sex + 1.case + 1.case#15.age_gp + 1.case#0.sex + 15.age_gp#0.sex + 1.case#15.age_gp#0.sex. So the log odds ratio will be the difference between these. Many of the terms are common to both and cancel. So we are left with:

log(odds ratio) = 1.case + 1.case#15.age_gp + 1.case#0.sex + 1.case#15.age_gp#0.sex

If your original command was -logistic- (as opposed to -logit-), -lincom- will automatically exponentiate this for you, so you can get the odds ratio with

lincom 1.case + 1.case#15.age_gp + 1.case#0.sex + 1.case#15.age_gp#0.sex.

PS: * ignored the terms involving the other confounders because, unless they, too, have an interaction with case in the model (which you do not suggest is the case), they will appear in both the numerator and denominator, and thus drop out.