1.1 Rationale for the package

There has been a long relationship between the disciplines of Epidemiology / Public Health and Biostatistics. Students frequently find introductory textbooks explaining statistical methods and the maths behind them, but how to implement those techniques in a computer is most of the time not explained.

One of the most popular statistical software’s in public health settings is Stata. Stata has the advantage of offering a solid interface with functions that can be accessed via the use of commands or by selecting the proper input in the graphical unit interface (GUI). Furthermore, Stata offers “Grad Plans” to postgraduate students, making it relatively affordable from an economic point of view.

R is a free alternative to Stata. The use of R keeps growing; furthermore, with the relatively high number of packages and textbooks available, it’s popularity is also increasing.

In the case of epidemiology, there are already some good packages available for R, including: Epi, epibasix, epiDisplay, epiR and epitools. The pubh package does not intend to replace any of them, but to only provide a common syntax for the most frequent statistical analysis in epidemiology.

1.2 Syntax: the use of formulas

Most students and professionals from the disciplines of Epidemiology and Public Health analyse variables in terms of outcome, exposure and confounders. The following table shows the most common names used in the literature to characterise variables in a cause-effect relationships:

In R, formulas are used to declare relationships between variables. Formulas are common in classic statistical tests and in regression methods.

Formulas have the following standard syntax:

y ~ x, data = my_data

Where y is the outcome or response variable, x is the exposure or predictor of interest and my_data specifies the name of the data frame where x and y can be found.

The symbol ~ is used in R for formulas. It can be interpreted as depends on. In the most typical scenario, y ~ x means y depends on x or y is a function of x:

y = f(x)

Using epidemiology friendly terms:

Outcome = f(Exposure)

Is worth to mention that Stata requires for variables to be given in the same order as in formulas: outcome first, predictors next.

The pubh package integrates well with other packages of the tidyverse which use formulas and the pipe operator %>% to add layers over functions. In particular, pubh uses ggformula as a graphical interface for plotting and takes advantage of variable labels from sjlabelled. This versatility allows it to interact also with tables from moonBook and effect plots from sjPlot.

y ~ x|z, data = my_data

1.3 Pipes and the tidyverse

The tidyverse has become now the standard in data manipulation in R. The use of the pipe function %>% from magrittr allows for cleaner code. The principle of pipes is to change the paradigm in coding. In standard codding, when many functions are used, one goes from inner parenthesis to outer ones.

For example if we have three functions, a common syntax would look like:

f3(f2(f1(..., data), ...), ...)

With pipes, however, the code reads top to bottom and left to right:

data %>%
  f1(...) %>% 
  f2(...) %>% 
  f3(...)

Most of the functions from pubh are compatible with pipes and the tidyverse.

1.4 Contributions of the pubh package

The main contributions of the pubh package to students and professionals from the disciplines of Epidemiology and Public Health are:

1. Use of a common syntax for the most used analysis.
2. A function, glm_coef that displays coefficients from most common regression analysis in a way that can be easy interpreted and used for publications.
3. Integration with the ggformula package, introducing plotting functions with a relatively simple syntax.
4. Integration with the most common epidemiological analysis, using the standard formula syntax explained in the previous section.
5. Integration with the gtsummary and huxtable packages. The pubh use huxtables as standard as they can be easily exported to html, pdf and doc files when knitting. Tables of descriptive statistics and of regression coefficients are generated through functions from gtsummary and then converted to huxtable objects with convenient cosmetics.

2.1 Inferential statistics

From the descriptive statistics of Ratio we know that the relative dispersion in the Hodgkin group is almost as double as the relative dispersion in the Non-Hodgkin group.

For new users of R it helps to have a common syntax in most of the commands, as R could be challenging and even intimidating at times. We can test if the difference in variance is statistically significant with the var.test command, which uses can use a formula syntax:

var.test(Ratio ~ Group, data = Hodgkin)

    F test to compare two variances

data:  Ratio by Group
F = 0.6294, num df = 19, denom df = 19, p-value = 0.3214
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
 0.2491238 1.5901459
sample estimates:
ratio of variances 
         0.6293991 

What about normality? We can look at the QQ-plots against the standard Normal distribution.

Hodgkin %>%
  qq_plot(~ Ratio|Group) 

Let’s say we choose a non-parametric test to compare the groups:

wilcox.test(Ratio ~ Group, data = Hodgkin)

    Wilcoxon rank sum exact test

data:  Ratio by Group
W = 298, p-value = 0.007331
alternative hypothesis: true location shift is not equal to 0

2.2 Graphical output

For relatively small samples (for example, less than 30 observations per group) is a standard practice to show the actual data in dot plots with error bars. The pubh package offers two options to show graphically differences in continuous outcomes among groups:

• For small samples: strip_error
• For medium to large samples: bar_error

For our current example:

Hodgkin %>%
  strip_error(Ratio ~ Group)

The error bars represent 95% confidence intervals around mean values.

Is relatively easy to add a line on top, to show that the two groups are significantly different. The function gf_star needs the reference point on how to draw an horizontal line to display statistical difference or to annotate a plot; in summary, gf_star:

• Draws an horizontal line from \((x1, y2)\) to \((x2, y2)\).
• Draws a vertical line from \((x1, y1)\) to \((x1, y1)\).
• Draws a vertical line from \((x2, y1)\) to \((x2, y1)\).
• Writes a character (the legend or “star”) at the mid point between \(x1\) and \(x2\) at high \(y3\).

Thus:

\[ y1 < y2 < y3 \]

In our current example:

Hodgkin %>%
  strip_error(Ratio ~ Group) %>%
  gf_star(x1 = 1, y1 = 4, x2 = 2, y2 = 4.05, y3 = 4.1, "**")

For larger samples we could use bar charts with error bars. For example:

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

birthwt %>%
  bar_error(bwt ~ smoke)

Quick normality check:

birthwt %>%
  qq_plot(~ bwt|smoke)

The (unadjusted) \(t\)-test:

birthwt %>% 
  t_test(bwt ~ smoke, detailed = TRUE) %>% 
  as.data.frame()

  estimate estimate1 estimate2 .y.     group1 group2  n1 n2 statistic     p
1 283.7767  3055.696  2771.919 bwt Non-smoker Smoker 115 74  2.729886 0.007
        df conf.low conf.high method alternative
1 170.1002 78.57486  488.9786 T-test   two.sided

The final plot with annotation to highlight statistical difference (unadjusted):

birthwt %>%
  bar_error(bwt ~ smoke) %>%
  gf_star(x1 = 1, x2 = 2, y1 = 3400, y2 = 3500, y3 = 3550, "**")

Both strip_error and bar_error can generate plots stratified by a third variable, for example:

birthwt %>%
  bar_error(bwt ~ smoke, fill = ~ Race) %>%
  gf_refine(ggsci::scale_fill_jama()) 

birthwt %>%
  bar_error(bwt ~ smoke|Race)

birthwt %>%
  strip_error(bwt ~ smoke, pch = ~ Race, col = ~ Race) %>%
  gf_refine(ggsci::scale_color_jama())

4.1 Contingency tables

Let’s say we want to look at the effect of ibuprofen on preventing death in patients with sepsis.

data(Bernard)
Bernard %>% head()

# A tibble: 6 × 9
     id treat     race             fate  apache o2del followup temp0 temp10
  <dbl> <fct>     <fct>            <fct>  <dbl> <dbl>    <dbl> <dbl>  <dbl>
1     1 Placebo   White            Dead      27  539.       50  35.2   36.6
2     2 Ibuprofen African American Alive     14   NA       720  38.7   37.6
3     3 Placebo   African American Dead      33  551.       33  38.3   NA  
4     4 Ibuprofen White            Alive      3 1376.      720  38.3   36.4
5     5 Placebo   White            Alive      5   NA       720  38.6   37.6
6     6 Ibuprofen White            Alive     13 1523.      720  38.2   38.2

Let’s look at the table:

Bernard %>%
  mutate(
    fate = relevel(fate, ref = "Dead"),
    treat = relevel(treat, ref = "Ibuprofen")
  ) %>%
  copy_labels(Bernard) %>%
  select(fate, treat) %>% 
  cross_tbl(by = "treat") %>%
  theme_pubh(2)

For epi.2by2 we need to provide the table of counts in the correct order, so we would type something like:

dat <- matrix(c(84, 140 , 92, 139), nrow = 2, byrow = TRUE)
epiR::epi.2by2(as.table(dat))

             Outcome +    Outcome -      Total        Inc risk *        Odds
Exposed +           84          140        224              37.5       0.600
Exposed -           92          139        231              39.8       0.662
Total              176          279        455              38.7       0.631

Point estimates and 95% CIs:
-------------------------------------------------------------------
Inc risk ratio                                 0.94 (0.75, 1.19)
Odds ratio                                     0.91 (0.62, 1.32)
Attrib risk in the exposed *                   -2.33 (-11.27, 6.62)
Attrib fraction in the exposed (%)            -6.20 (-33.90, 15.76)
Attrib risk in the population *                -1.15 (-8.88, 6.59)
Attrib fraction in the population (%)         -2.96 (-15.01, 7.82)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 0.260 Pr>chi2 = 0.610
Fisher exact test that OR = 1: Pr>chi2 = 0.631
 Wald confidence limits
 CI: confidence interval
 * Outcomes per 100 population units 

For contingency we only need to provide the information in a formula:

Bernard %>%
  contingency(fate ~ treat) 

           Outcome
Predictor   Dead Alive
  Ibuprofen   84   140
  Placebo     92   139

             Outcome +    Outcome -      Total        Inc risk *        Odds
Exposed +           84          140        224              37.5       0.600
Exposed -           92          139        231              39.8       0.662
Total              176          279        455              38.7       0.631

Point estimates and 95% CIs:
-------------------------------------------------------------------
Inc risk ratio                                 0.94 (0.75, 1.19)
Odds ratio                                     0.91 (0.62, 1.32)
Attrib risk in the exposed *                   -2.33 (-11.27, 6.62)
Attrib fraction in the exposed (%)            -6.20 (-33.90, 15.76)
Attrib risk in the population *                -1.15 (-8.88, 6.59)
Attrib fraction in the population (%)         -2.96 (-15.01, 7.82)
-------------------------------------------------------------------
Uncorrected chi2 test that OR = 1: chi2(1) = 0.260 Pr>chi2 = 0.610
Fisher exact test that OR = 1: Pr>chi2 = 0.631
 Wald confidence limits
 CI: confidence interval
 * Outcomes per 100 population units 


    Pearson's Chi-squared test with Yates' continuity correction

data:  dat
X-squared = 0.17076, df = 1, p-value = 0.6794

Advantages of contingency:

1. Easier input without the need to create the table.
2. Displays the standard epidemiological table at the start of the output. This aids to check what are the reference levels on each category.
3. In the case that the \(\chi^2\)-test is not appropriate, contingency would show the results of the Fisher exact test at the end of the output.

outcome ~ stratum/exposure, data = my_data

data(oswego, package = "epitools")

oswego <- oswego %>%
  mutate(
    ill = factor(ill, labels = c("No", "Yes")),
    sex = factor(sex, labels = c("Female", "Male")),
    chocolate.ice.cream = factor(chocolate.ice.cream, labels = c("No", "Yes"))
  ) %>%
  var_labels(
    ill = "Developed illness",
    sex = "Sex",
    chocolate.ice.cream = "Consumed chocolate ice cream"
  )

oswego %>%
  mhor(ill ~ sex/chocolate.ice.cream)

                                   OR Lower.CI Upper.CI Pr(>|z|)
sexFemale:chocolate.ice.creamYes 0.47     0.11     2.06    0.318
sexMale:chocolate.ice.creamYes   0.24     0.05     1.13    0.072

                          Common OR Lower CI Upper CI Pr(>|z|)
Cochran-Mantel-Haenszel:       0.35     0.12     1.01    0.085

Test for effect modification (interaction): p =  0.5434 
 

4.2 Diagnostic tests

Example: We compare the use of lung’s X-rays on the screening of TB against the gold standard test.

Freq <- c(1739, 8, 51, 22)
BCG <- gl(2, 1, 4, labels=c("Negative", "Positive"))
Xray <- gl(2, 2, labels=c("Negative", "Positive"))
tb <- data.frame(Freq, BCG, Xray)
tb <- expand_df(tb)

diag_test(BCG ~ Xray, data = tb)

          Outcome +    Outcome -      Total
Test +           22           51         73
Test -            8         1739       1747
Total            30         1790       1820

Point estimates and 95% CIs:
--------------------------------------------------------------
Apparent prevalence *                  0.04 (0.03, 0.05)
True prevalence *                      0.02 (0.01, 0.02)
Sensitivity *                          0.73 (0.54, 0.88)
Specificity *                          0.97 (0.96, 0.98)
Positive predictive value *            0.30 (0.20, 0.42)
Negative predictive value *            1.00 (0.99, 1.00)
Positive likelihood ratio              25.74 (18.21, 36.38)
Negative likelihood ratio              0.27 (0.15, 0.50)
False T+ proportion for true D- *      0.03 (0.02, 0.04)
False T- proportion for true D+ *      0.27 (0.12, 0.46)
False T+ proportion for T+ *           0.70 (0.58, 0.80)
False T- proportion for T- *           0.00 (0.00, 0.01)
Correctly classified proportion *      0.97 (0.96, 0.98)
--------------------------------------------------------------
* Exact CIs

Using the immediate version (direct input):

diag_test2(22, 51, 8, 1739)

          Outcome +    Outcome -      Total
Test +           22           51         73
Test -            8         1739       1747
Total            30         1790       1820

Point estimates and 95% CIs:
--------------------------------------------------------------
Apparent prevalence *                  0.04 (0.03, 0.05)
True prevalence *                      0.02 (0.01, 0.02)
Sensitivity *                          0.73 (0.54, 0.88)
Specificity *                          0.97 (0.96, 0.98)
Positive predictive value *            0.30 (0.20, 0.42)
Negative predictive value *            1.00 (0.99, 1.00)
Positive likelihood ratio              25.74 (18.21, 36.38)
Negative likelihood ratio              0.27 (0.15, 0.50)
False T+ proportion for true D- *      0.03 (0.02, 0.04)
False T- proportion for true D+ *      0.27 (0.12, 0.46)
False T+ proportion for T+ *           0.70 (0.58, 0.80)
False T- proportion for T- *           0.00 (0.00, 0.01)
Correctly classified proportion *      0.97 (0.96, 0.98)
--------------------------------------------------------------
* Exact CIs