1 Introduction

The aim of this vignette is to illustrate the use of pubh functions for common regression analysis in Public Health. In particular, the vignette shows the use of the following functions from pubh:

1. cosm_reg for reporting tables of coefficients.
2. cross_tbl for reporting tables of descriptive statistics by exposure of interest.
3. multiple for adjusting confidence intervals and p-values for post-hoc analysis.
4. get_r2 for accessing \(R^2\) or pseudo-\(R^2\) from regression models.
5. glm_coef for some special cases of regression models.

The advantages and limitations of glm_coef are:

1. Recognises the main models used in epidemiology/public health.
2. Automatically back-transforms estimates and confidence intervals, when the model requires it.
3. Can use robust standard errors for the calculation of confidence intervals.
   • Standard errors are used by default.
   • The use of standard errors is restricted by the following classes of objects (models): gee, glm and survreg.
4. Can display nice labels for the names of the parameters.
5. Returns a data frame that can be modified and/or exported as tables for publications (with further editing).

We start by loading relevant packages and setting the theme for the plots (as suggested in the Template of this package):

rm(list = ls())
library(tidyverse)
library(rstatix)
library(parameters)
library(performance)
library(jtools)
library(pubh)
library(sjlabelled)
library(sjPlot)

theme_set(sjPlot::theme_sjplot2(base_size = 10))
theme_update(legend.position = "top")
options('huxtable.knit_print_df' = FALSE)
options('huxtable.autoformat_number_format' = list(numeric = "%5.2f"))
knitr::opts_chunk$set(comment = NA)

2 Multiple Linear Regression

For continuous outcomes there is no need of exponentiating the results unless the outcome was fitted in the log-scale. In our first example we want to estimate the effect of smoking and race on the birth weight of babies.

We can generate factors and assign labels in the same pipe stream:

data(birthwt, package = "MASS")
birthwt <- birthwt %>%
  mutate(
    smoke = factor(smoke, labels = c("Non-smoker", "Smoker")),
    race = factor(race, labels = c("White", "African American", "Other"))
    ) %>%
  var_labels(
    bwt = 'Birth weight (g)',
    smoke = 'Smoking status',
    race = 'Race'
    )

Is good to start with some basic descriptive statistics, so we can compare the birth weight between groups.

Graphical analysis:

birthwt %>%
  box_plot(bwt ~ smoke, fill = ~ race)

Another way to compare the means between the groups is with gen_bst_df which estimates means with corresponding bootstrapped CIs.

birthwt %>%
  gen_bst_df(bwt ~ race|smoke) %>%
  as_hux() %>% theme_pubh()

We fit a linear model.

model_norm <- lm(bwt ~ smoke + race, data = birthwt)

Note: Model diagnostics are not be discussed in this vignette.

Traditional output from the model:

model_norm %>% Anova()

Anova Table (Type II tests)

Response: bwt
            Sum Sq  Df F value    Pr(>F)    
smoke      7322575   1 15.4588 0.0001191 ***
race       8712354   2  9.1964 0.0001557 ***
Residuals 87631356 185                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_norm %>% parameters()

Parameter               | Coefficient |     SE |             95% CI | t(185) |      p
-------------------------------------------------------------------------------------
(Intercept)             |     3334.95 |  91.78 | [3153.89, 3516.01] |  36.34 | < .001
smoke [Smoker]          |     -428.73 | 109.04 | [-643.86, -213.60] |  -3.93 | < .001
race [African American] |     -450.36 | 153.12 | [-752.45, -148.27] |  -2.94 | 0.004 
race [Other]            |     -452.88 | 116.48 | [-682.67, -223.08] |  -3.89 | < .001

Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
  using a Wald t-distribution approximation.

model_norm %>% performance()

# Indices of model performance

AIC      |     AICc |      BIC |    R2 | R2 (adj.) |    RMSE |   Sigma
----------------------------------------------------------------------
3012.223 | 3012.551 | 3028.432 | 0.123 |     0.109 | 680.924 | 688.246

Table of coefficients for publication:

model_norm %>%
  tbl_regression() %>% 
  cosm_reg() %>% theme_pubh() %>% 
  add_footnote(get_r2(model_norm), font_size = 9)

model_norm %>% 
  glm_coef(labels = model_labels(model_norm)) %>%
  as_hux() %>% set_align(everywhere, 2:3, "right") %>%
  theme_pubh() %>%
  add_footnote(get_r2(model_norm), font_size = 9)

Function glm_coef allows the use of robust standard errors.

model_norm %>%
  glm_coef(se_rob = TRUE, labels = model_labels(model_norm)) %>%
  as_hux() %>% set_align(everywhere, 2:3, "right") %>%
  theme_pubh() %>%
  add_footnote(paste(
    get_r2(model_norm), "\n",
    "CIs and p-values estimated with robust standard errors."),
    font_size = 9)

To construct the effect plot, we can use plot_model from the sjPlot package. The advantage of plot_model is that recognises labelled data and uses that information for annotating the plots.

model_norm %>%
  plot_model("pred", terms = ~race|smoke, dot.size = 1.5, title = "")

When the explanatory variables are categorical, another option is emmip from the emmeans package. We can include CIs in emmip but as estimates are connected, the resulting plots look more messy, so I recommend emmip to look at the trace.

emmip(model_norm, smoke ~ race) %>%
  gf_labs(y = get_label(birthwt$bwt), x = "", col = "Smoking status")

3 Logistic Regression

For logistic regression we are interested in the odds ratios. We will look at the effect of amount of fibre intake on the development of coronary heart disease.

data(diet, package = "Epi")
diet <- diet %>%
  mutate(
    chd = factor(chd, labels = c("No CHD", "CHD"))
  ) %>%
  var_labels(
    chd = "Coronary Heart Disease",
    fibre = "Fibre intake (10 g/day)"
    )

We start with descriptive statistics:

diet %>% estat(~ fibre|chd) %>%
  as_hux() %>% theme_pubh()

diet %>% na.omit() %>%
  copy_labels(diet) %>%
  box_plot(fibre ~ chd)

We fit a linear logistic model:

model_binom <- glm(chd ~ fibre, data = diet, family = binomial)

Summary:

model_binom %>% parameters(exponentiate = TRUE)

Parameter   | Odds Ratio |   SE |       95% CI |     z |     p
--------------------------------------------------------------
(Intercept) |       0.95 | 0.58 | [0.30, 3.16] | -0.08 | 0.936
fibre       |       0.33 | 0.12 | [0.15, 0.66] | -2.94 | 0.003

Uncertainty intervals (profile-likelihood) and p-values (two-tailed)
  computed using a Wald z-distribution approximation.

model_binom %>% performance()

# Indices of model performance

AIC     |    AICc |     BIC | Tjur's R2 |  RMSE | Sigma | Log_loss | Score_log | Score_spherical |   PCP
--------------------------------------------------------------------------------------------------------
257.535 | 257.571 | 265.151 |     0.028 | 0.337 | 1.000 |    0.381 |    -6.359 |           0.017 | 0.773

Table of coefficients for publication:

model_binom %>% 
  tbl_regression(exponentiate = TRUE) %>% 
  cosm_reg() %>% theme_pubh() %>% 
  add_footnote(get_r2(model_binom), font_size = 9)

Effect plot:

model_binom %>%
  plot_model("pred", terms = "fibre [all]", title = "")

data(bdendo, package = "Epi") 
bdendo <- bdendo %>%
  mutate(
    cancer = factor(d, labels = c('Control', 'Case')),
    gall = factor(gall, labels = c("No GBD", "GBD")),
    est = factor(est, labels = c("No oestrogen", "Oestrogen"))
  ) %>%
  var_labels(
    cancer = 'Endometrial cancer',
    gall = 'Gall bladder disease',
    est = 'Oestrogen'
  )

bdendo %>% 
  mutate(
    cancer = relevel(cancer, ref = "Case"),
    gall = relevel(gall, ref = "GBD"),
    est = relevel(est, ref = "Oestrogen")
  ) %>%
  copy_labels(bdendo) %>%
  select(cancer, gall, est) %>% 
  tbl_strata(
    strata = est,
    .tbl_fun = ~ .x %>%
      tbl_summary(by = gall)
  ) %>% 
  cosm_sum(bold = TRUE, head_label = " ") %>% 
  theme_pubh(2) %>% 
  set_align(1, everywhere, "center")

bdendo %>%
  gf_percents(~ cancer|gall, fill = ~ est, position = "dodge", alpha = 0.6) %>%
  gf_labs(
    y = "Percent",
    x = "",
    fill = ""
  )

require(survival, quietly = TRUE)
model_clogit <- clogit(cancer == 'Case' ~ est * gall + strata(set), data = bdendo)

model_clogit %>%
  glm_coef(labels = c("Oestrogen/No oestrogen", "GBD/No GBD",
                      "Oestrogen:GBD Interaction")) %>%
  as_hux() %>% set_align(everywhere, 2:3, "right") %>%
  theme_pubh() %>%
  add_footnote(get_r2(model_clogit), font_size = 9)

#require(MASS, quietly = TRUE)
data(housing, package = "MASS")
housing <- housing %>%
  var_labels(
    Sat = "Satisfaction",
    Infl = "Perceived influence",
    Type = "Type of rental",
    Cont = "Contact"
    )

model_ord <- MASS::polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing, Hess = TRUE)

model_ord %>% 
  tbl_regression(exponentiate = TRUE) %>% 
  cosm_reg() %>% theme_pubh() %>% 
  add_footnote(get_r2(model_ord), font_size = 9)

model_ord %>%
  plot_model(type = "pred", terms = c("Infl", "Cont"), 
           dot.size = 1, title = "") %>%
  gf_theme(axis.text.x = element_text(angle = 45, hjust = 1))

model_ord %>%
  plot_model(type = "pred", terms = c("Infl", "Type"), 
           dot.size = 1, title = "") %>%
  gf_theme(axis.text.x = element_text(angle = 45, hjust = 1))

emmip(model_ord, Infl ~ Type |Cont) %>%
  gf_labs(x = "Type of rental", col = "Perceived influence")

4 Poisson Regression

For Poisson regression we are interested in incidence rate ratios. We will look at the effect of sex, ethnicity and age group on number of absent days from school in a year.

data(quine, package = "MASS")
levels(quine$Eth) <- c("Aboriginal", "White")
levels(quine$Sex) <- c("Female", "Male")

quine <- quine %>%
  var_labels(
    Days = "Number of absent days",
    Eth = "Ethnicity",
    Age = "Age group"
    )

Descriptive statistics:

quine %>% 
  cross_tbl(by = "Eth") %>% 
  theme_pubh(2) %>% 
  add_footnote("n (%); Median (IQR)", font_size = 9)

quine %>%
  box_plot(Days ~ Age|Sex, fill = ~ Eth)

We start by fitting a standard Poisson linear regression model:

model_pois <- glm(Days ~ Eth + Sex + Age, family = poisson, data = quine)

model_pois %>% 
  tbl_regression(exponentiate = TRUE) %>% 
  cosm_reg() %>% theme_pubh() %>% 
  add_footnote(get_r2(model_pois), font_size = 9)

model_pois %>% performance()

# Indices of model performance

AIC      |     AICc |      BIC | Nagelkerke's R2 |   RMSE | Sigma | Score_log | Score_spherical
-----------------------------------------------------------------------------------------------
2342.982 | 2343.586 | 2360.883 |           0.896 | 14.940 | 1.000 |    -7.983 |           0.049

model_pois %>% check_overdispersion()

# Overdispersion test

       dispersion ratio =   13.890
  Pearson's Chi-Squared = 1944.553
                p-value =  < 0.001

Overdispersion detected.

model_negbin <- MASS::glm.nb(Days ~ Eth + Sex + Age, data = quine)

model_negbin %>% 
  tbl_regression(exponentiate = TRUE) %>% 
  cosm_reg() %>% theme_pubh() %>% 
  add_footnote(get_r2(model_negbin), font_size = 9)

model_negbin %>% performance()

# Indices of model performance

AIC      |     AICc |      BIC | Nagelkerke's R2 |   RMSE | Sigma | Score_log | Score_spherical
-----------------------------------------------------------------------------------------------
1109.653 | 1110.464 | 1130.538 |           0.210 | 15.087 | 1.000 |    -3.908 |           0.067

model_negbin %>%
  plot_model(type = "pred", terms = c("Age", "Eth"), 
           dot.size = 1.5, title = "") 

emmip(model_negbin, Eth ~ Age|Sex) %>%
  gf_labs(y = "Number of absent days", x = "Age group", col = "Ethnicity")

multiple(model_negbin, ~ Age|Eth)$df

   contrast        Eth ratio   SE null z.ratio p.value lower.CL upper.CL
1   F1 / F0 Aboriginal  0.69 0.16    1   -1.57    0.40     0.38     1.26
2   F2 / F0 Aboriginal  1.20 0.28    1    0.77    0.86     0.66     2.17
3   F2 / F1 Aboriginal  1.73 0.35    1    2.65    0.04     1.02     2.92
4   F3 / F0 Aboriginal  1.29 0.31    1    1.04    0.72     0.69     2.40
5   F3 / F1 Aboriginal  1.86 0.40    1    2.89    0.02     1.07     3.22
6   F3 / F2 Aboriginal  1.08 0.23    1    0.34    0.99     0.62     1.88
7   F1 / F0      White  0.69 0.16    1   -1.57    0.40     0.38     1.26
8   F2 / F0      White  1.20 0.28    1    0.77    0.86     0.66     2.17
9   F2 / F1      White  1.73 0.35    1    2.65    0.04     1.02     2.92
10  F3 / F0      White  1.29 0.31    1    1.04    0.72     0.69     2.40
11  F3 / F1      White  1.86 0.40    1    2.89    0.02     1.07     3.21
12  F3 / F2      White  1.08 0.23    1    0.34    0.99     0.62     1.88

multiple(model_negbin, ~ Age|Eth)$fig_ci %>%
  gf_labs(x = "IRR")

5 Survival Analysis

We will use an example on the effect of thiotepa versus placebo on the development of bladder cancer.

data(bladder)

Warning in data(bladder): data set 'bladder' not found

bladder <- bladder %>%
  mutate(times = stop,
         rx = factor(rx, labels=c("Placebo", "Thiotepa"))
         ) %>%
  var_labels(times = "Survival time",
             rx = "Treatment")

model_surv <- survreg(Surv(times, event) ~ rx, data = bladder)

model_surv %>%
  glm_coef(labels = c("Treatment: Thiotepa/Placebo", "Scale"), se_rob = TRUE) %>%
  as_hux() %>% set_align(everywhere, 2:3, "right") %>%
  theme_pubh(1) %>%
  add_footnote(get_r2(model_surv), font_size = 9)

model_exp <- survreg(Surv(times, event) ~ rx, data = bladder, dist = "exponential")

model_exp %>%
  glm_coef(labels = c("Treatment: Thiotepa/Placebo"), se_rob = TRUE) %>%
  as_hux() %>% set_align(everywhere, 2:3, "right") %>%
  theme_pubh() %>%
  add_footnote(get_r2(model_exp), font_size = 9)

model_exp %>%
  glm_coef(labels = c("Treatment: Thiotepa/Placebo")) %>%
  as_hux() %>% set_align(everywhere, 2:3, "right") %>%
  theme_pubh() %>%
  add_footnote(get_r2(model_exp), font_size = 9)

model_exp %>%
  plot_model(type = "pred", terms = ~ rx, dot.size = 1.5, title = "") %>%
  gf_labs(y = "Survival time", x = "Treatment", title = "")

model_cox <-  coxph(Surv(times, event) ~ rx, data = bladder)

model_cox %>% 
  tbl_regression(exponentiate = TRUE) %>% 
  cosm_reg() %>% theme_pubh() %>% 
  add_footnote(get_r2(model_cox), font_size = 9)