Appendix 3C: OR and MEM Comparison

// Health and Retirement Study 2006 //
use "https://tdmize.github.io/data/data/hrs06", clear

(hrs06_dmv.dta | Health & Retirement Study 2006 | DMV Workshop)

logit diabetes c.age i.black c.bmi i.woman i.ed4cat, or cformat(%6.4f)

Iteration 0:  Log likelihood = -8471.9299  
Iteration 1:  Log likelihood = -7829.5253  
Iteration 2:  Log likelihood = -7801.8913  
Iteration 3:  Log likelihood = -7801.8369  
Iteration 4:  Log likelihood = -7801.8369  

Logistic regression                                    Number of obs =  17,012
                                                       LR chi2(7)    = 1340.19
                                                       Prob > chi2   =  0.0000
Log likelihood = -7801.8369                            Pseudo R2     =  0.0791

------------------------------------------------------------------------------
    diabetes | Odds ratio   Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         age |     1.0221     0.0020   10.972   0.000       1.0181      1.0261
             |
       black |
      Black  |     1.6465     0.0856    9.593   0.000       1.4870      1.8231
         bmi |     1.1060     0.0038   29.064   0.000       1.0985      1.1135
             |
       woman |
      Women  |     0.7363     0.0300   -7.508   0.000       0.6797      0.7975
             |
      ed4cat |
   12 years  |     0.7591     0.0390   -5.364   0.000       0.6864      0.8395
13-15 years  |     0.6194     0.0375   -7.905   0.000       0.5500      0.6975
  16+ years  |     0.5426     0.0341   -9.715   0.000       0.4796      0.6138
             |
       _cons |     0.0045     0.0009  -26.745   0.000       0.0030      0.0067
------------------------------------------------------------------------------
Note: _cons estimates baseline odds.

est store diabmod

library(haven)           # import Stata data and labels
library(dplyr)           # data manipulation
library(marginaleffects) # predictions and marginal effects
library(questionr)       # odds ratios

# Health and Retirement Study 2006
hrs <- read_dta("https://tdmize.github.io/data/data/hrs06.dta")

# Treat binary/nominal variables as factors using their value labels
hrs <- hrs |>
  mutate(black = as_factor(black), woman = as_factor(woman), ed4cat = as_factor(ed4cat))

diabmod <- glm(diabetes ~ age + black + bmi + woman + ed4cat,
  family = binomial(link = "logit"), data = hrs)
odds.ratio(diabmod, level = .95, digits = 4)

                         OR     2.5 % 97.5 %         p    
(Intercept)       0.0044862 0.0030136 0.0067 < 2.2e-16 ***
age               1.0220590 1.0180880 1.0261 < 2.2e-16 ***
blackBlack        1.6464795 1.4864500 1.8224 < 2.2e-16 ***
bmi               1.1059647 1.0985064 1.1135 < 2.2e-16 ***
womanWomen        0.7362919 0.6797383 0.7975 5.987e-14 ***
ed4cat12 years    0.7591201 0.6864286 0.8396 8.126e-08 ***
ed4cat13-15 years 0.6193869 0.5498380 0.6973 2.686e-15 ***
ed4cat16+ years   0.5425925 0.4794171 0.6136 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Keep the OR for Black (1.6465) and the OR for BMI (1.1060) in mind for the following sections

Compute MEM and OR for effect of race

* Compute OR for effect of race
est restore diabmod

(results diabmod are active now)

mtable, at(black=(0 1)) decimal(7) atmeans post

Expression: Pr(diabetes), predict()

           |    black     Pr(y)
 ----------+-------------------
         1 |        0  0.1668597
         2 |        1  0.2479810

Specified values of covariates

           |                            1.          2.        3.        4.
           |      age       bmi     woman      ed4cat    ed4cat    ed4cat
 ----------+-------------------------------------------------------------
   Current |   68.798   27.8424   .583059    .3410534    .21185   .217023

mlincom 2 - 1 // MEM

             |   lincom    pvalue        ll        ul 
-------------+---------------------------------------
           1 |    0.081     0.000     0.063     0.099

di 0.2479810 / (1 - 0.2479810) // Odds when black = 1

.32975364

di 0.1668597 / (1 - 0.1668597) // Odds when black = 0

.20027803

di %6.4f .32975364 / .20027803 // Odds ratio

1.6465

# Compute OR for effect of race: predictions with other variables at means/modes
prs_black <- predictions(diabmod, newdata = datagrid(black = c("White", "Black")))
prs_black


 black Estimate Pr(>|z|)     S 2.5 % 97.5 %
 White    0.157   <0.001   Inf 0.147  0.168
 Black    0.235   <0.001 308.6 0.215  0.256

Type: invlink(link)

comparisons(diabmod, variables = "black", newdata = "mean") # MEM


 Estimate Std. Error    z Pr(>|z|)    S 2.5 % 97.5 %
   0.0778    0.00909 8.56   <0.001 56.3  0.06 0.0956

Term: black
Type: response
Comparison: Black - White

odds_black <- prs_black$estimate / (1 - prs_black$estimate)
odds_black                   # Odds when black = 0 and black = 1

[1] 0.1866140 0.3072561

odds_black[2] / odds_black[1] # Odds ratio

[1] 1.64648

Reminder—odds (Ω) are just a transformation of a probability (π):

\[\Omega = \frac{\pi}{1 - \pi}\]

Compute MEM and OR for effect of BMI

* Compute OR for effect of BMI
est restore diabmod

(results diabmod are active now)

sum bmi

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
         bmi |     17,123    27.83901    5.837704   10.57755    82.6728

mtable, at(bmi=(27.83901 28.83901)) decimal(7) atmeans post

Expression: Pr(diabetes), predict()

           |      bmi     Pr(y)
 ----------+-------------------
         1 |   27.839  0.1772895
         2 |   28.839  0.1924603

Specified values of covariates

           |                  1.        1.          2.        3.        4.
           |      age     black     woman      ed4cat    ed4cat    ed4cat
 ----------+-------------------------------------------------------------
   Current |   68.798   .147543   .583059    .3410534    .21185   .217023

mlincom 2 - 1 // MEM

             |   lincom    pvalue        ll        ul 
-------------+---------------------------------------
           1 |    0.015     0.000     0.014     0.016

di 0.1924603 / (1 - 0.1924603) // Odds when bmi = mean + 1

.23832921

di 0.1772895 / (1 - 0.1772895) // Odds when bmi = mean

.21549439

di %6.4f .23832921 / .21549439 // Odds ratio

1.1060

# Compute OR for effect of BMI: predictions with other variables at means/modes
summary(hrs$bmi)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max.     NAs 
  10.58   23.89   27.06   27.84   30.82   82.67     414

mean_bmi <- mean(hrs$bmi, na.rm = TRUE)

prs_bmi <- predictions(diabmod, newdata = datagrid(bmi = c(mean_bmi, mean_bmi + 1)))
prs_bmi


  bmi Estimate Pr(>|z|)   S 2.5 % 97.5 %
 27.8    0.157   <0.001 Inf 0.147  0.168
 28.8    0.171   <0.001 Inf 0.160  0.182

Type: invlink(link)

comparisons(diabmod, variables = list(bmi = 1), newdata = "mean") # MEM


 Estimate Std. Error    z Pr(>|z|)     S  2.5 % 97.5 %
   0.0138   0.000562 24.6   <0.001 440.6 0.0127 0.0149

Term: bmi
Type: response
Comparison: +1

odds_bmi <- prs_bmi$estimate / (1 - prs_bmi$estimate)
odds_bmi                 # Odds when bmi = mean and bmi = mean + 1

[1] 0.1865501 0.2063179

odds_bmi[2] / odds_bmi[1] # Odds ratio

[1] 1.105965

Given that the OR is simply a different way to present the same information a MEM summarizes—which would you prefer to interpret? Which is easier to understand?