Reproducible epidemiology in R

What is `{gtsummary}`?

Create tables that are publication-ready
Highly customizable
Descriptive tables, regression tables, etc.

`gtsummary::tbl_summary()`

library(gtsummary)

tbl_summary(
  nlsy,
  by = sex_cat,
  include = c(sex_cat, race_eth_cat, region_cat,
              eyesight_cat, glasses, age_bir))

Characteristic	Male, N = 6,403¹	Female, N = 6,283¹
race_eth_cat
Hispanic	1,000 (16%)	1,002 (16%)
Black	1,613 (25%)	1,561 (25%)
Non-Black, Non-Hispanic	3,790 (59%)	3,720 (59%)
region_cat
Northeast	1,296 (21%)	1,254 (20%)
North Central	1,488 (24%)	1,446 (23%)
South	2,251 (36%)	2,317 (38%)
West	1,253 (20%)	1,142 (19%)
Unknown	115	124
eyesight_cat
Excellent	1,582 (38%)	1,334 (31%)
Very good	1,470 (35%)	1,500 (35%)
Good	792 (19%)	1,002 (23%)
Fair	267 (6.4%)	365 (8.5%)
Poor	47 (1.1%)	85 (2.0%)
Unknown	2,245	1,997
glasses	1,566 (38%)	2,328 (54%)
Unknown	2,241	1,995
age_bir	25 (21, 29)	22 (19, 27)
Unknown	3,652	3,091
¹ n (%); Median (IQR)

tbl_summary(
  nlsy,
  by = sex_cat,
  include = c(sex_cat, race_eth_cat, region_cat,
              eyesight_cat, glasses, age_bir),
  label = list(
    race_eth_cat ~ "Race/ethnicity",
    region_cat ~ "Region",
    eyesight_cat ~ "Eyesight",
    glasses ~ "Wears glasses",
    age_bir ~ "Age at first birth"
  ),
  missing_text = "Missing")

Characteristic	Male, N = 6,403¹	Female, N = 6,283¹
Race/ethnicity
Hispanic	1,000 (16%)	1,002 (16%)
Black	1,613 (25%)	1,561 (25%)
Non-Black, Non-Hispanic	3,790 (59%)	3,720 (59%)
Region
Northeast	1,296 (21%)	1,254 (20%)
North Central	1,488 (24%)	1,446 (23%)
South	2,251 (36%)	2,317 (38%)
West	1,253 (20%)	1,142 (19%)
Missing	115	124
Eyesight
Excellent	1,582 (38%)	1,334 (31%)
Very good	1,470 (35%)	1,500 (35%)
Good	792 (19%)	1,002 (23%)
Fair	267 (6.4%)	365 (8.5%)
Poor	47 (1.1%)	85 (2.0%)
Missing	2,245	1,997
Wears glasses	1,566 (38%)	2,328 (54%)
Missing	2,241	1,995
Age at first birth	25 (21, 29)	22 (19, 27)
Missing	3,652	3,091
¹ n (%); Median (IQR)

tbl_summary(
  nlsy,
  by = sex_cat,
  include = c(sex_cat, race_eth_cat,
              eyesight_cat, glasses, age_bir),
  label = list(
    race_eth_cat ~ "Race/ethnicity",
    eyesight_cat ~ "Eyesight",
    glasses ~ "Wears glasses",
    age_bir ~ "Age at first birth"
  ),
  missing_text = "Missing") |> 
  add_p(test = list(all_continuous() ~ "t.test", 
                    all_categorical() ~ "chisq.test")) |> 
  add_overall(col_label = "**Total**") |> 
  bold_labels() |> 
  modify_footnote(update = everything() ~ NA) |> 
  modify_header(label = "**Variable**", p.value = "**P**")

Variable	Overall, N = 12,686	Male, N = 6,403	Female, N = 6,283	P
Race/ethnicity				0.8
Hispanic	2,002 (16%)	1,000 (16%)	1,002 (16%)
Black	3,174 (25%)	1,613 (25%)	1,561 (25%)
Non-Black, Non-Hispanic	7,510 (59%)	3,790 (59%)	3,720 (59%)
Eyesight				<0.001
Excellent	2,916 (35%)	1,582 (38%)	1,334 (31%)
Very good	2,970 (35%)	1,470 (35%)	1,500 (35%)
Good	1,794 (21%)	792 (19%)	1,002 (23%)
Fair	632 (7.5%)	267 (6.4%)	365 (8.5%)
Poor	132 (1.6%)	47 (1.1%)	85 (2.0%)
Missing	4,242	2,245	1,997
Wears glasses	3,894 (46%)	1,566 (38%)	2,328 (54%)	<0.001
Missing	4,236	2,241	1,995
Age at first birth	23 (20, 28)	25 (21, 29)	22 (19, 27)	<0.001
Missing	6,743	3,652	3,091

`tbl_summary()`

Incredibly customizeable
Really helpful with Table 1
I often just view in the web browser and copy and paste into a Word document
Can also be used within quarto/R Markdown¹
If output is Word, I use as_flex_table()
Make even more customizeable with as_gt()
- Then can output to Word with gt::as_word()

Univariate regressions

Fit a series of univariate regressions of income on other variables.

tbl_uvregression(
  nlsy, 
  y = income,
  include = c(sex_cat, race_eth_cat,
              eyesight_cat, income, age_bir),
  method = lm)

Characteristic	N	Beta	95% CI¹	p-value
age_bir	4,773	595	538, 652	<0.001
sex_cat	10,195
Male		—	—
Female		-358	-844, 128	0.15
race_eth_cat	10,195
Hispanic		—	—
Black		-1,747	-2,507, -988	<0.001
Non-Black, Non-Hispanic		3,863	3,195, 4,530	<0.001
eyesight_cat	6,789
Excellent		—	—
Very good		-578	-1,319, 162	0.13
Good		-1,863	-2,719, -1,006	<0.001
Fair		-4,674	-5,910, -3,439	<0.001
Poor		-6,647	-9,154, -4,140	<0.001
¹ CI = Confidence Interval

Can also do logistic regression

tbl_uvregression(
  nlsy, 
  y = glasses,
  include = c(sex_cat, race_eth_cat,
              eyesight_cat, glasses, age_bir),
  method = glm,
  method.args = list(family = binomial()),
  exponentiate = TRUE)

Characteristic	N	OR¹	95% CI¹	p-value
age_bir	5,813	1.02	1.01, 1.03	<0.001
sex_cat	8,450
Male		—	—
Female		1.97	1.81, 2.15	<0.001
race_eth_cat	8,450
Hispanic		—	—
Black		0.76	0.67, 0.86	<0.001
Non-Black, Non-Hispanic		1.34	1.19, 1.50	<0.001
eyesight_cat	8,444
Excellent		—	—
Very good		0.93	0.84, 1.03	0.2
Good		0.95	0.84, 1.07	0.4
Fair		0.81	0.68, 0.96	0.016
Poor		1.15	0.81, 1.63	0.4
¹ OR = Odds Ratio, CI = Confidence Interval

Customizable just like tbl_summary()

Some regressions

linear_model <- lm(income ~ sex_cat + age_bir + race_eth_cat, 
                   data = nlsy)

linear_model_int <- lm(income ~ sex_cat*age_bir + race_eth_cat, 
                   data = nlsy)

logistic_model <- glm(glasses ~ eyesight_cat + sex_cat + income, 
                      data = nlsy, family = binomial())

`gtsummary::tbl_regression()`

tbl_regression(
  linear_model, 
  intercept = TRUE,
  label = list(
    sex_cat ~ "Sex",
    race_eth_cat ~ "Race/ethnicity",
    age_bir ~ "Age at first birth"
  ))

Characteristic	Beta	95% CI¹	p-value
(Intercept)	2,147	493, 3,802	0.011
Sex
Male	—	—
Female	25	-654, 705	>0.9
Age at first birth	438	381, 495	<0.001
Race/ethnicity
Hispanic	—	—
Black	-772	-1,714, 171	0.11
Non-Black, Non-Hispanic	7,559	6,676, 8,442	<0.001
¹ CI = Confidence Interval

tbl_regression(
  logistic_model, 
  exponentiate = TRUE,
  label = list(
    sex_cat ~ "Sex",
    eyesight_cat ~ "Eyesight",
    income ~ "Income"
  ))

Characteristic	OR¹	95% CI¹	p-value
Eyesight
Excellent	—	—
Very good	0.92	0.82, 1.03	0.2
Good	0.92	0.80, 1.05	0.2
Fair	0.80	0.66, 0.98	0.028
Poor	1.03	0.69, 1.53	0.9
Sex
Male	—	—
Female	2.04	1.85, 2.25	<0.001
Income	1.00	1.00, 1.00	<0.001
¹ OR = Odds Ratio, CI = Confidence Interval

You could put several together

tbl_no_int <- tbl_regression(
  linear_model, 
  intercept = TRUE,
  label = list(
    sex_cat ~ "Sex",
    race_eth_cat ~ "Race/ethnicity",
    age_bir ~ "Age at first birth"
  ))

tbl_int <- tbl_regression(
  linear_model_int, 
  intercept = TRUE,
  label = list(
    sex_cat ~ "Sex",
    race_eth_cat ~ "Race/ethnicity",
    age_bir ~ "Age at first birth",
    `sex_cat:age_bir` ~ "Sex/age interaction"
  ))

You could put several together

tbl_merge(list(tbl_no_int, tbl_int), 
          tab_spanner = c("**Model 1**", "**Model 2**"))

Characteristic	Model 1			Model 2
Characteristic	Beta	95% CI¹	p-value	Beta	95% CI¹	p-value
(Intercept)	2,147	493, 3,802	0.011	4,064	1,884, 6,245	<0.001
Sex
Male	—	—		—	—
Female	25	-654, 705	>0.9	-3,635	-6,432, -838	0.011
Age at first birth	438	381, 495	<0.001	364	285, 443	<0.001
Race/ethnicity
Hispanic	—	—		—	—
Black	-772	-1,714, 171	0.11	-759	-1,701, 183	0.11
Non-Black, Non-Hispanic	7,559	6,676, 8,442	<0.001	7,550	6,668, 8,433	<0.001
Sex/age interaction
Female * Age at first birth				149	39, 260	0.008
¹ CI = Confidence Interval

Inline text

Instead of copying and pasting numbers from the table, pull them off automatically:

inline_text(tbl_no_int, sex_cat, level = "Female")

which produces: 25 (95% CI -654, 705; p>0.9)

and which can be further customized:

inline_text(tbl_no_int, sex_cat, level = "Female", pattern = "$\\beta = {estimate}$")

which produces:

bonus: `broom::tidy()`

bind_rows(
  broom::tidy(linear_model, conf.int = TRUE),
  broom::tidy(linear_model_int, conf.int = TRUE),
  .id = "model"
) |>
  mutate(model = factor(model, 
                        labels = c("main terms", "sex-age interaction"))) |> 
  ggplot(aes(x = term, y = estimate, 
             ymin = conf.low, ymax = conf.high,
             color = model)) +
  geom_point(position = position_dodge(width = .5)) +
  geom_errorbar(position = position_dodge(width = .5))