This tutorial teaches how to create side by side univariate and multivariate table, such as the one reported in the linked paper below, using simulated data.

https://www.aaojournal.org/action/showFullTableHTML?isHtml=true&tableId=tbl2&pii=S0161-6420%2825%2900410-5

Install Packages

# Install packages if needed
# install.packages(c("gtsummary", "geepack", "dplyr", "tidyr", "gt"))

library(gtsummary)
library(geepack)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
library(gt)

Create simulated eye-level data

set.seed(1234)

n_patients <- 500

# Create 1 or 2 eyes per patient
df <- data.frame(
  Patient_ID = 1:n_patients,
  n_eyes = sample(
    c(1, 2),
    n_patients,
    replace = TRUE,
    prob = c(0.30, 0.70)
  )
) |>
  uncount(n_eyes, .id = "Eye_ID")



# Create patient/eye characteristics


df <- df |>
  mutate(

    # Continuous variable
    Age = round(
      pmax(
        pmin(rnorm(n(), mean = 68, sd = 10), 95),
        40
      )
    ),

    # Categorical variables
    BMI_Category = factor(
      sample(
        c("Normal", "Overweight", "Obese"),
        n(),
        replace = TRUE,
        prob = c(0.30, 0.40, 0.30)
      ),
      levels = c("Normal", "Overweight", "Obese")
    ),

    SEX = factor(
      sample(
        c("Male", "Female"),
        n(),
        replace = TRUE,
        prob = c(0.45, 0.55)
      ),
      levels = c("Male", "Female")
    ),

    Race = factor(
      sample(
        c("White", "Black", "Other"),
        n(),
        replace = TRUE,
        prob = c(0.65, 0.20, 0.15)
      ),
      levels = c("White", "Black", "Other")
    ),

    Surgery_Type = factor(
      sample(
        c(
          "Cataract",
          "Cataract with glaucoma",
          "Glaucoma",
          "Other"
        ),
        n(),
        replace = TRUE,
        prob = c(0.55, 0.20, 0.15, 0.10)
      ),
      levels = c(
        "Cataract",
        "Cataract with glaucoma",
        "Glaucoma",
        "Other"
      )
    ),

    ASA = factor(
      sample(
        c("1 and 2", "3 and 4", "Not available"),
        n(),
        replace = TRUE,
        prob = c(0.50, 0.45, 0.05)
      ),
      levels = c(
        "1 and 2",
        "3 and 4",
        "Not available"
      )
    ),

    Insurance = factor(
      sample(
        c(
          "Medicare",
          "Commercial",
          "Medicaid",
          "Uninsured"
        ),
        n(),
        replace = TRUE,
        prob = c(0.55, 0.30, 0.12, 0.03)
      ),
      levels = c(
        "Medicare",
        "Commercial",
        "Medicaid",
        "Uninsured"
      )
    )
  )



#  Generate binary outcome

# Y = "Pre"  = did not pass
# Y = "Pass" = passed
#
# The probabilities are generated from a logistic model.

linear_predictor <-
  -0.5 +
  (-0.03 * (df$Age - 68)) +
  (-0.10 * (df$BMI_Category == "Overweight")) +
  (-0.25 * (df$BMI_Category == "Obese")) +
  (-0.15 * (df$SEX == "Female")) +
  (-0.50 * (df$Race == "Black")) +
  (-0.10 * (df$Race == "Other")) +
  (-0.35 * (df$Surgery_Type == "Cataract with glaucoma")) +
  (-0.80 * (df$Surgery_Type == "Glaucoma")) +
  (-0.10 * (df$Surgery_Type == "Other")) +
  (-0.40 * (df$ASA == "3 and 4")) +
  (-0.10 * (df$ASA == "Not available")) +
  ( 0.10 * (df$Insurance == "Commercial")) +
  (-0.20 * (df$Insurance == "Medicaid")) +
  ( 0.20 * (df$Insurance == "Uninsured"))

prob_pass <- plogis(linear_predictor)

# Check probabilities
summary(prob_pass)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 0.04189 0.17365 0.24049 0.24616 0.30576 0.55971
# Generate binary outcome
df$Y <- rbinom(
  n = nrow(df),
  size = 1,
  prob = prob_pass
)

# Convert to factor
df$Y <- factor(
  df$Y,
  levels = c(0, 1),
  labels = c("Pre", "Pass")
)

# Set reference level
df$Y <- relevel(df$Y, ref = "Pre")


#Check the outcome

table(df$Y)
## 
##  Pre Pass 
##  644  208
prop.table(table(df$Y))
## 
##       Pre      Pass 
## 0.7558685 0.2441315
# Check for missing values
sum(is.na(df$Y))
## [1] 0
# Check that there are no invalid values
stopifnot(
  !anyNA(df$Y)
)

Convert Outcome to numeric for GEEglm fitting and test for one variable

library(geepack)

df <- df |>
  mutate(
    Y_numeric = if_else(Y == "Pass", 1, 0)
  )
table(df$Y, df$Y_numeric, useNA = "ifany")
##       
##          0   1
##   Pre  644   0
##   Pass   0 208
gee_age <- geeglm(
  Y_numeric ~ Age,
  data = df,
  family = binomial(link = "logit"),
  id = Patient_ID,
  corstr = "exchangeable"
)

summary(gee_age)
## 
## Call:
## geeglm(formula = Y_numeric ~ Age, family = binomial(link = "logit"), 
##     data = df, id = Patient_ID, corstr = "exchangeable")
## 
##  Coefficients:
##              Estimate   Std.err   Wald Pr(>|W|)    
## (Intercept)  0.821639  0.559390  2.157 0.141884    
## Age         -0.029172  0.008306 12.336 0.000444 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation structure = exchangeable 
## Estimated Scale Parameters:
## 
##             Estimate Std.err
## (Intercept)        1 0.06197
##   Link = identity 
## 
## Estimated Correlation Parameters:
##       Estimate Std.err
## alpha -0.02717  0.0412
## Number of clusters:   500  Maximum cluster size: 2
library(gtsummary)

tbl_age <-
  gee_age |>
  tbl_regression(
    exponentiate = TRUE,
    pvalue_fun = label_style_pvalue(digits = 2)
  )

tbl_age
Characteristic OR 95% CI p-value
Age 0.97 0.96, 0.99 <0.001
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

Now Fit Univariate Regression for all Variables of interest

table_vars <- c(
  "Age",
  "BMI_Category",
  "SEX"
)

tbl_unadjusted <-
  df |>
  tbl_uvregression(
    method = geeglm,
    y = Y_numeric,
    include = all_of(table_vars),
    method.args = list(
      family = binomial(link = "logit"),
      id = Patient_ID,
      corstr = "exchangeable"
    ),
    exponentiate = TRUE,
    pvalue_fun = label_style_pvalue(digits = 2)
  )

tbl_unadjusted
Characteristic N OR 95% CI p-value
Age 852 0.97 0.96, 0.99 <0.001
BMI_Category 852


    Normal
— —
    Overweight
0.72 0.49, 1.05 0.090
    Obese
0.60 0.40, 0.90 0.013
SEX 852


    Male
— —
    Female
0.70 0.52, 0.96 0.025
Abbreviations: CI = Confidence Interval, OR = Odds Ratio 
## Add Formatting
tbl_unadjusted <-
  tbl_unadjusted |>
  bold_p(t = 0.05) |>
  bold_labels()

tbl_unadjusted
Characteristic N OR 95% CI p-value
Age 852 0.97 0.96, 0.99 <0.001
BMI_Category 852


    Normal
— —
    Overweight
0.72 0.49, 1.05 0.090
    Obese
0.60 0.40, 0.90 0.013
SEX 852


    Male
— —
    Female
0.70 0.52, 0.96 0.025
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

Create Adjusted GEE

gee_adjusted <- geeglm(
  Y_numeric ~
    Age +
    BMI_Category +
    SEX,
  data = df,
  family = binomial(link = "logit"),
  id = Patient_ID,
  corstr = "exchangeable"
)

tbl_adjusted <-
  gee_adjusted |>
  tbl_regression(
    exponentiate = TRUE,
    pvalue_fun = label_style_pvalue(digits = 2)
  ) |>
  bold_p(t = 0.05) |>
  bold_labels()

tbl_adjusted
Characteristic OR 95% CI p-value
Age 0.97 0.95, 0.99 <0.001
BMI_Category


    Normal — —
    Overweight 0.70 0.48, 1.03 0.071
    Obese 0.58 0.38, 0.88 0.010
SEX


    Male — —
    Female 0.68 0.49, 0.93 0.015
Abbreviations: CI = Confidence Interval, OR = Odds Ratio

Now Merge them

table2 <-
  tbl_merge(
    tbls = list(
      tbl_unadjusted,
      tbl_adjusted
    ),
    tab_spanner = c(
      "**Unadjusted**",
      "**Adjusted**"
    )
  ) |>
  modify_header(
    estimate_1 = "**OR (95% CI)**",
    p.value_1 = "***P* Value**",
    estimate_2 = "**OR (95% CI)**",
    p.value_2 = "***P* Value**"
  )

table2
Characteristic
Unadjusted
Adjusted
N OR (95% CI) 95% CI P Value OR (95% CI) 95% CI P Value
Age 852 0.97 0.96, 0.99 <0.001 0.97 0.95, 0.99 <0.001
BMI_Category 852





    Normal
— —
— —
    Overweight
0.72 0.49, 1.05 0.090 0.70 0.48, 1.03 0.071
    Obese
0.60 0.40, 0.90 0.013 0.58 0.38, 0.88 0.010
SEX 852





    Male
— —
— —
    Female
0.70 0.52, 0.96 0.025 0.68 0.49, 0.93 0.015
Abbreviations: CI = Confidence Interval, OR = Odds Ratio