Loading Required Packages

• This is an R Markdown presentation. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents.
• For more details on using R Markdown see http://rmarkdown.rstudio.com.
• Here we will use slidy to create slides using a markdwon file.

### Load required packages
library(arsenal) # for tableby
require(knitr)
require(survival)
require(kableExtra)
require(equatiomatic)
require(lme4)
require(gtsummary)
require(DiagrammeR)

Introduction

• This tutorial provides step-by-step directions for using the functions associated with tableby() to produce one of the most common tables in medical literature includes summary statistics for a set of variables, often stratified by some group (e.g. treatment arm).
• All functions presented here are available within the arsenal package.
• All the examples in this tutorial use a dataset called mockstudy made available by Paul Novotny which includes a variety of types of variables (character, numeric, factor, ordered factor, survival) to use as examples.

Utility Package KableExtra

• It is worthwhile mentioning another package i.e., kableExtra before getting started. It helps us build common complex tables and manipulate table styles.

• It imports the pipe %>% symbol from magrittr and verbalize all the functions, so basically you can add “layers” to a kable output in a way that is similar with ggplot2 and plotly. For more table styles refer to https://cran.r-project.org/web/packages/kableExtra/vignettes/awesome_table_in_html.html

• For users who are not very familiar with the pipe operator %>% in R, it passes the result along the chain for a more literal coding experience. Basically when we say A %>% B, technically it means sending the results of A to B as B’s first argument.

Bookstyle Table using KableExtra

Scroll Box table using KableExtra

-When you use scroll_box, you can specify either or both height or width to get a two-way scrollable box.

kbl(cbind(dt1, dt1)) %>%
  kable_paper() %>%
  scroll_box(width = "500px", height = "200px")

Pulling Underlying Equation Corresponding to a Statistical Model using Package Equatiomatic

lr <- glm(sex ~ species * bill_length_mm, data = penguins,family = binomial(link = "logit"))

extract_eq(lr, wrap = TRUE)

\[ \begin{aligned} \log\left[ \frac { P( \operatorname{sex} = \operatorname{male} ) }{ 1 - P( \operatorname{sex} = \operatorname{male} ) } \right] &= \alpha + \beta_{1}(\operatorname{species}_{\operatorname{Chinstrap}}) + \beta_{2}(\operatorname{species}_{\operatorname{Gentoo}}) + \beta_{3}(\operatorname{bill\_length\_mm})\ + \\ &\quad \beta_{4}(\operatorname{species}_{\operatorname{Chinstrap}} \times \operatorname{bill\_length\_mm}) + \beta_{5}(\operatorname{species}_{\operatorname{Gentoo}} \times \operatorname{bill\_length\_mm}) \end{aligned} \]

extract_eq(lr, wrap = TRUE, show_distribution = TRUE) # show the how the data are assumed to be distributed.

\[ \begin{aligned} \operatorname{sex} &\sim Bernoulli\left(\operatorname{prob}_{\operatorname{sex} = \operatorname{male}}= \hat{P}\right) \\ \log\left[ \frac { \hat{P} }{ 1 - \hat{P} } \right] &= \alpha + \beta_{1}(\operatorname{species}_{\operatorname{Chinstrap}}) + \beta_{2}(\operatorname{species}_{\operatorname{Gentoo}}) + \beta_{3}(\operatorname{bill\_length\_mm})\ + \\ &\quad \beta_{4}(\operatorname{species}_{\operatorname{Chinstrap}} \times \operatorname{bill\_length\_mm}) + \beta_{5}(\operatorname{species}_{\operatorname{Gentoo}} \times \operatorname{bill\_length\_mm}) \end{aligned} \]

## GLMM with individual-level variability (accounting for overdispersion)
## For this data set the model is the same as one allowing for a period:herd
## interaction, which the plot indicates could be needed.
cbpp$obs <- 1:nrow(cbpp)
gm2 <- glmer(cbind(incidence, size - incidence) ~ period +(1 | herd) +  (1|obs),family = binomial, data = cbpp)
extract_eq(gm2,wrap=TRUE)

\[ \begin{aligned} \operatorname{incidence}_{i} &\sim \operatorname{Binomial}(n = 1, \operatorname{prob}_{\operatorname{incidence} = 1} = \widehat{P}) \\ \log\left[\frac{\hat{P}}{1 - \hat{P}} \right] &=\alpha_{j[i],k[i]} \\ \alpha_{j} &\sim N \left(\gamma_{0}^{\alpha} + \gamma_{1}^{\alpha}(\operatorname{period}_{\operatorname{2}}) + \gamma_{2}^{\alpha}(\operatorname{period}_{\operatorname{3}}) + \gamma_{3}^{\alpha}(\operatorname{period}_{\operatorname{4}}), \sigma^2_{\alpha_{j}} \right) \text{, for obs j = 1,} \dots \text{,J} \\ \alpha_{k} &\sim N \left(\mu_{\alpha_{k}}, \sigma^2_{\alpha_{k}} \right) \text{, for herd k = 1,} \dots \text{,K} \end{aligned} \]

Using gtsummary to Create Publication Ready Output Tables from Fitted Models

-Start by creating a logistic regression model to predict tumor response using the variables age and grade from the trial data set.

tr<-head(trial)
tr %>%
  kbl(caption = "Simulated Trial Data") %>%
  kable_classic(full_width = F, html_font = "Cambria")

m1 <- glm(response ~ age + stage, trial, family = binomial)

tbl_regression(m1, exponentiate = TRUE)

Format Results into Data Frame with Global p-values

m1 %>%
  tbl_regression(
    exponentiate = TRUE, 
    pvalue_fun = ~style_pvalue(.x, digits = 2),
  ) %>% 
  add_global_p() %>%
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()

Adding False Dsicovery Rates

trial %>%
  select(response, age, grade) %>%
  tbl_uvregression(
    method = glm,
    y = response,
    method.args = list(family = binomial),
    exponentiate = TRUE,
    pvalue_fun = ~style_pvalue(.x, digits = 2)
  ) %>%
  add_global_p() %>%  # add global p-value 
  add_nevent() %>%    # add number of events of the outcome
  add_q() %>%         # adjusts global p-values for multiple testing
  bold_p() %>%        # bold p-values under a given threshold (default 0.05)
  bold_p(t = 0.10, q = TRUE) %>% # now bold q-values under the threshold of 0.10
  bold_labels()

## add_q: Adjusting p-values with
## `stats::p.adjust(x$table_body$p.value, method = "fdr")`

Create a Table Stratified by Treatment Arm

• We use a formula statement to specify the variables which we want summarized.
• The example below uses age (a continuous variable) and sex (a factor).
• Note results=“asis” when creating the R chunk below changes the layout slightly (compresses it) and bolds the variable names. If you want to use tabl as dataframe use as.data.frame.

tab1 <- tableby(arm ~ sex + age, data=mockstudy)
summary(tab1) # quick look at the table using  summary() on your tableby object, use text =true for pretty text version

Add Labels to Variables

• Either you could add labels to the variables in your dataset or you can use the built-in data.frame method for labels <-: as below

labels(mockstudy)  <- c(age = 'Age, yrs', sex = "Gender")

tab1 <- tableby(arm ~ sex + age, data=mockstudy)
summary(tab1)

Change Summary Statistics Globally

• By default it summarize continuous variables with: Number of missing values, Mean (SD), 25th - 75th quantiles, and Minimum-Maximum values with an ANOVA (t-test with equal variances) p-value.

• For categorical variables the default is to show: Number of missing values and count (column percent) with a chi-square p-value.

• This default can be modified using the tableby.control function. Note that test=FALSE and total=FALSE results in the total column and p-value column not being printed.

mycontrols  <- tableby.control(test=FALSE, total=FALSE,
                               numeric.test="kwt", cat.test="chisq",
                               numeric.stats=c("N", "median", "q1q3"),
                               cat.stats=c("countpct"),
                               stats.labels=list(N='Count', median='Median', q1q3='Q1,Q3'))
tab2 <- tableby(arm ~ sex + age, data=mockstudy, control=mycontrols)
summary(tab2)

Change Summary Statistics Within the Formula

• We can also modify the test and summary statistics for specific variables within the formula statement. For example, below we use the kwt (Kruskal-Wallis rank-based) for age and anova for BMI, and no test for the variable “ast”.

• The tests function can do a quick check on what tests were performed on each variable in tableby.

tab.test <- tableby(arm ~ kwt(age) + anova(bmi) + notest(ast), data=mockstudy)
tt<-tests(tab.test)
tt %>%
  kbl(caption = "Tests function for a quick check on what tests were performed") %>%
  kable_classic(full_width = F, html_font = "Cambria")

Modifying Summary Statistics

• Summary statistics for any individual variable can also be modified, but it must be done as secondary arguments to the test function.
• The function names must be strings that are functions already written for tableby.

summary(tab.test, pfootnote=TRUE) # default summary table with footnote display option TRUE

1. Kruskal-Wallis rank sum test
2. Linear Model ANOVA

tab.test1 <- tableby(arm ~ kwt(ast, "Nmiss2","median") + anova(age, "N","mean") +
                    notest(bmi, "Nmiss","median"), data=mockstudy)
summary(tab.test1)

Categorical Tests (Chisq and Fisher’s)

• The formal tests for categorical variables against the levels of the by variable, chisq and FE, have options to simulate p-values. We show how to turn on the simulations for these with 500 replicates for the Fisher’s test (FE).

• The chi-square test on 2x2 tables applies Yates’ continuity correction by default,with and without the correction that is applied to treatment arm by sex, if we use subset to ignore one of the three treatment arms.

set.seed(100)
tab.catsim <- tableby(arm ~ sex + race, cat.test="fe", simulate.p.value=TRUE, B=500, data=mockstudy)
t1<-tests(tab.catsim)
t1 %>%
  kbl(caption = "Fisher's Exact Test for Count Data with simulated p-value") %>%
  kable_classic(full_width = F, html_font = "Cambria")

cat.correct <- tableby(arm ~ sex + race, cat.test="chisq", subset = !grepl("^F", arm), data=mockstudy)
t2<-tests(cat.correct)

t2 %>%
  kbl(caption = "Chi-Square Test with correction for Count Data with simulated p-value") %>%
  kable_classic(full_width = F, html_font = "Cambria")

cat.nocorrect <- tableby(arm ~ sex + race, cat.test="chisq", subset = !grepl("^F", arm),chisq.correct=FALSE, data=mockstudy)
t3<-tests(cat.nocorrect)
t3 %>%
  kbl(caption = "Chi-Square Test without correction for Count Data with simulated p-value") %>%
  kable_classic(full_width = F, html_font = "Cambria")

Summarize Without a Group Variable

tab.noby <- tableby(~ bmi + sex + age, data=mockstudy)
summary(tab.noby)

Summarize an Ordered Factor

When comparing groups of ordered data there are two options.

• The default uses a general independence test available from the coin package. For two-group comparisons, this is essentially the Armitage trend test.

• The other option is to specify the Kruskal Wallis test.

summary(tableby(sex ~ age.ord, data = mockstudy), pfootnote = TRUE)

1. Trend test for ordinal variables

Summarize a Survival Variable

• Information presented by the survfit() function,can be seen with tableby.
• The default is to show the median survival (time at which the probability of survival = 50%).
• It is also possible to obtain summaries of the % survival at certain time points (say the probability of surviving 1-year).

survfit(Surv(fu.time, fu.stat)~sex, data=mockstudy) # first look at info from survfit func 

## Call: survfit(formula = Surv(fu.time, fu.stat) ~ sex, data = mockstudy)
## 
##              n events median 0.95LCL 0.95UCL
## sex=Male   916    829    550     515     590
## sex=Female 583    527    543     511     575

survdiff(Surv(fu.time, fu.stat)~sex, data=mockstudy) 

## Call:
## survdiff(formula = Surv(fu.time, fu.stat) ~ sex, data = mockstudy)
## 
##              N Observed Expected (O-E)^2/E (O-E)^2/V
## sex=Male   916      829      830  0.000370  0.000956
## sex=Female 583      527      526  0.000583  0.000956
## 
##  Chisq= 0  on 1 degrees of freedom, p= 1

# now apply tableby
summary(tableby(sex ~ Surv(fu.time, fu.stat), data=mockstudy))

summary(survfit(Surv(fu.time/365.25, fu.stat)~sex, data=mockstudy), times=1:5) # using survfit

## Call: survfit(formula = Surv(fu.time/365.25, fu.stat) ~ sex, data = mockstudy)
## 
##                 sex=Male 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    626     286   0.6870  0.0153       0.6576       0.7177
##     2    309     311   0.3437  0.0158       0.3142       0.3761
##     3    152     151   0.1748  0.0127       0.1516       0.2015
##     4     57      61   0.0941  0.0104       0.0759       0.1168
##     5     24      16   0.0628  0.0095       0.0467       0.0844
## 
##                 sex=Female 
##  time n.risk n.event survival std.err lower 95% CI upper 95% CI
##     1    380     202   0.6531  0.0197       0.6155        0.693
##     2    190     189   0.3277  0.0195       0.2917        0.368
##     3     95      90   0.1701  0.0157       0.1420        0.204
##     4     51      32   0.1093  0.0133       0.0861        0.139
##     5     18      12   0.0745  0.0126       0.0534        0.104

summary(tableby(sex ~ Surv(fu.time/365.25, fu.stat), data=mockstudy, times=1:5, surv.stats=c("NeventsSurv","NriskSurv"))) # using tableby

Summarising Subset data using the Subset Argument within tableby

summary(tableby(sex ~ ps + hgb + bmi, subset=age>50 & arm=="F: FOLFOX", data=mockstudy))

Create Combinations of Variables

## create a variable combining the levels of mdquality.s and sex
with(mockstudy, table(interaction(mdquality.s,sex)))

0.Male 1.Male 0.Female 1.Female 77 686 47 437

summary(tableby(arm ~ interaction(mdquality.s,sex), data=mockstudy))

## create a new grouping variable with combined levels of arm and sex
summary(tableby(interaction(mdquality.s, sex) ~  age + bmi, data=mockstudy, subset=arm=="F: FOLFOX"))

Variable Transformations

• If the transformation includes any of the symbols / - + ^ * then surround the new variable by I(), so that R knows to treat it as a variable transformation and not some special model feature.

trans <- tableby(arm ~ I(age/10) + log(bmi) + factor(mdquality.s, levels=0:1, labels=c('N','Y')),
                 data=mockstudy)
summary(trans)

• Change lables to desired ones

labels(trans)[2:4] <- c('Age per 10 yrs', 'log(BMI)', 'MD Quality')
labels(trans)

##                                                     arm 
##                                         "Treatment Arm" 
##                                               I(age/10) 
##                                        "Age per 10 yrs" 
##                                                log(bmi) 
##                                              "log(BMI)" 
## factor(mdquality.s, levels = 0:1, labels = c("N", "Y")) 
##                                            "MD Quality"

summary(trans)

Subsetting Tables

mytab <- tableby(arm ~ sex + alk.phos + age, data=mockstudy)
mytab2 <- mytab[c('age','sex','alk.phos')] #Change the ordering of the variables
summary(mytab2)

summary(mytab[c(3,1)], digits = 3) # delete a variable

summary(mytab[mytab < 0.5]) # filter p-value<0.5

summary(sort(mytab, decreasing = TRUE)) # sort p-value

Add a Title to the Table

• When creating a pdf the tables are automatically numbered and the title appears below the table. In Word and HTML, the titles appear un-numbered and above the table.

t1 <- tableby(arm ~ sex + age, data=mockstudy)
summary(t1, title='Demographics')

#With multiple left-hand sides, you can pass a vector or list to determine labels for each table:
summary(tableby(list(arm, sex) ~ age, data = mockstudy), title = c("arm table", "sex table"))

The End!!

• R Markdown contains a combination of YAML (metadata), text (narratives), and code chunks.To understand and/or change structure of a Markdown document an excellent source to refer to is https://bookdown.org/yihui/rmarkdown-cookbook/rmarkdown-process.html

• DiagrammeR provides methods to build graphs for a number of different graphing languages. We will present a Graphviz example in this section,6 but you can also use pure R code to create graphs and diagrams with DiagrammeR