Overview

reportRmd is a package designed to facilitate the reporting of common statistical outputs easily in RMarkdown documents. The package supports pdf, html and word output without any changes to the body of the report. The main features are Table 1 style summaries, combining multiple univariate regression models into a single table, tidy multivariable model output and combining univariate and multivariable regressions into a single table. Single table summaries of median survival times and survival probabilities are also provided. A highly customisable survival curve function, based on ggplot2 can be used to create publication-quality plots. Visualisation plots are also available for bivariate relationships and logistic regression models.

A word of caution:

The reportRmd package is designed to facilitate statistical reporting and does not provide any checks regarding the suitability of the models fit.

Summary statistics

Basic summary statistics

data("pembrolizumab")
rm_covsum(data=pembrolizumab, 
covs=c('age','sex'))

Set IQR = T for interquartile range instead of Min/Max

rm_covsum(data=pembrolizumab, 
covs=c('age','sex'),IQR=TRUE)

Or all.stats=T for both

rm_covsum(data=pembrolizumab, 
covs=c('age','sex'),all.stats = TRUE)

This will produce summary statistics by Sex

rm_covsum(data=pembrolizumab, maincov = 'sex',
covs=c('age','pdl1','change_ctdna_group'))

To indicate which statistical test was used use show.tests=TRUE

rm_covsum(data=pembrolizumab, maincov = 'sex',
covs=c('age','pdl1','change_ctdna_group'),
show.tests=TRUE)

Effect sizes can be added with effSize = TRUE. Effect size measures include the Wilcoxon r for Wilcoxon rank-sum test, Cohen’s d for t-test, Eta squared for Kruskal Wallis test and ANOVA, and Cramer’s V for categorical variables.

rm_covsum(data=pembrolizumab, maincov = 'sex',
covs=c('age','pdl1','change_ctdna_group'),
show.tests=TRUE, effSize=TRUE)

Group comparisons are non-parametric by default, specify testcont='ANOVA' for t-tests/ANOVA

rm_covsum(data=pembrolizumab, maincov = 'sex',
covs=c('age','pdl1'),
testcont='ANOVA',
show.tests=TRUE, effSize=TRUE)

The default is to indicate percentages by columns (ie percentages within columns add to 100)

rm_covsum(data=pembrolizumab, maincov = 'sex',
covs=c('cohort'),
pvalue = FALSE)

But you can also specify to show by row instead

rm_covsum(data=pembrolizumab, maincov = 'sex',
covs=c('cohort'),
pvalue = FALSE,
percentage='row')

Univariate regression

Combining multiple univariate regression analyses into a single table

rm_uvsum(data=pembrolizumab, response='orr',
covs=c('age','pdl1','change_ctdna_group'),p.adjust = 'holm')

If the response is continuous linear regression is the default

rm_uvsum(data=pembrolizumab, response='l_size',
covs=c('age','cohort'))

…unless two variables are specified and then survival analysis is run

rm_uvsum(data=pembrolizumab, response=c('os_time','os_status'),
covs=c('age','pdl1','change_ctdna_group'))

Correlated observations can be handled using GEE

data("ctDNA")
 rm_uvsum(response = 'size_change',
 covs=c('time','ctdna_status'),
 gee=TRUE,
 id='id', corstr="exchangeable",
 family=gaussian("identity"),
 data=ctDNA,showN=TRUE)

If you want to check the underlying models, set returnModels = TRUE

rm_uvsum(response = 'orr',
 covs=c('age','sex'),
 data=pembrolizumab,returnModels = TRUE)

## $age
## 
## Call:  glm(formula = orr ~ age, family = binomial, data = data)
## 
## Coefficients:
## (Intercept)          age  
##     4.12269     -0.04231  
## 
## Degrees of Freedom: 93 Total (i.e. Null);  92 Residual
## Null Deviance:       85.77 
## Residual Deviance: 82.53     AIC: 86.53
## 
## $sex
## 
## Call:  glm(formula = orr ~ sex, family = binomial, data = data)
## 
## Coefficients:
## (Intercept)      sexMale  
##      1.9859      -0.8873  
## 
## Degrees of Freedom: 93 Total (i.e. Null);  92 Residual
## Null Deviance:       85.77 
## Residual Deviance: 83.21     AIC: 87.21

Multivariable analysis

To create a nice display for multivariable models the model needs to first be fit. By default, the variance inflation factor will be shown to check for multicollinearity. To suppress this column set vif=FALSE. Note: variance inflation factors are not computed (yet) for multilevel or GEE models.

glm_fit <- glm(orr~change_ctdna_group+pdl1+age,
               family='binomial',
               data = pembrolizumab)
rm_mvsum(glm_fit,p.adjust = 'holm')

p-values can be adjusted for multiple comparisons using any of the options available in the p.adjust function. This argument is also available for univariate models run with rm_uvsum.

rm_mvsum(glm_fit, showN = TRUE, vif=TRUE,p.adjust = 'holm')

Combining univariate and multivariable models

To display both models in a single table run the rm_uvsum and rm_mvsum functions with tableOnly=TRUE and combine.

uvsumTable <- rm_uvsum(data=pembrolizumab, response='orr',
covs=c('age','sex','pdl1','change_ctdna_group'),tableOnly = TRUE)

glm_fit <- glm(orr~change_ctdna_group+pdl1,
               family='binomial',
               data = pembrolizumab)
mvsumTable <- rm_mvsum(glm_fit, showN = TRUE,tableOnly = TRUE)

rm_uv_mv(uvsumTable,mvsumTable)

This can also be done with adjusted p-values, but when combined the raw p-values are dropped

uvsumTable <- rm_uvsum(data=pembrolizumab, response='orr',
covs=c('age','sex','pdl1','change_ctdna_group'),tableOnly = TRUE,p.adjust='holm')

glm_fit <- glm(orr~change_ctdna_group+pdl1,
               family='binomial',
               data = pembrolizumab)
mvsumTable <- rm_mvsum(glm_fit,tableOnly = TRUE,p.adjust='holm')

rm_uv_mv(uvsumTable,mvsumTable)

Changing the output

If you need to make changes to the tables, setting tableOnly=TRUE will return a data frame for any of the rm_ functions. Changes can be made, and the table output using outTable()

mvsumTable <- rm_mvsum(glm_fit, showN = TRUE,tableOnly = TRUE)
names(mvsumTable)[1] <-'Predictor'
outTable(mvsumTable)

Combining tables

Tables can be nested with the nestTable() function

cohortA <- rm_uvsum(data=subset(pembrolizumab,cohort=='A'), 
                     response = 'pdl1',
                     covs=c('age','sex'),
                     tableOnly = T)
cohortA$Cohort <- 'Cohort A'
cohortE <- rm_uvsum(data=subset(pembrolizumab,cohort=='E'), 
                     response = 'pdl1',
                     covs=c('age','sex'),
                     tableOnly = T)
cohortE$Cohort <- 'Cohort E'
nestTable(rbind(cohortA,cohortE),head_col = 'Cohort',to_col = 'Covariate')

Simple Survival Summaries

Displaying survival probabilities at different times by sex using Kaplan Meier estimates

rm_survsum(data=pembrolizumab,time='os_time',status='os_status',
 group="sex",survtimes=seq(12,36,12),survtimeunit='months')

Survival Times in Long Format

Displaying survival probabilities at different times by sex using Cox PH estimates

rm_survtime(data=pembrolizumab,time='os_time',status='os_status',
 strata="sex",survtimes=c(12,24),survtimeunit='mo',type='PH')

Displaying survival probabilities at different times by sex, adjusting for age using Cox PH estimates

rm_survtime(data=pembrolizumab,time='os_time',status='os_status', covs='age',
 strata="sex",survtimes=c(12,24),survtimeunit='mo',type='PH')

Stratified Survival Summary

rm_survdiff(data=pembrolizumab,time='os_time',status='os_status', 
            covs='sex',strata='cohort',digits=1)

Working with Labels

Variable labels will be shown in the nicenames argument is set to TRUE (the default). Variable labels are set using the label attribute of individual variables (assigned using reportRmd or another package like haven).

reportRmd supports the addition, removal and export of labels using the following functions:

• set_labels will set labels for a data frame from a lookup table
• set_var_labels allows you to set individual variable labels to a data frame
• clear_labels removes all labels from a data frame
• export_labels extracts variable labels from a data frame and returns a data frame of variable names and variable labels

data(ctDNA)
rm_covsum(data=ctDNA,
          covs=c('cohort','ctdna_status','size_change'))

ctDNA_names <- data.frame(var=names(ctDNA),
                          label=c('Patient ID',
                                  'Study Cohort',
                                  'Change in ctDNA since baseline',
                                  'Number of weeks on treatment',
                                  'Percentage change in tumour measurement'))
ctDNA <- set_labels(ctDNA,ctDNA_names)

rm_covsum(data=ctDNA,
          covs=c('cohort','ctdna_status','size_change'))

ctDNA <- set_var_labels(ctDNA,
                        cohort="A new cohort label")
rm_covsum(data=ctDNA,
          covs=c('cohort','ctdna_status','size_change'))

var_labels <- extract_labels(ctDNA)
var_labels

##       variable                                   label
## 1           id                              Patient ID
## 2       cohort                      A new cohort label
## 3 ctdna_status          Change in ctDNA since baseline
## 4         time            Number of weeks on treatment
## 5  size_change Percentage change in tumour measurement

ctDNA <- clear_labels(ctDNA)

Plotting Functions

plotuv(data=pembrolizumab, response='orr',
covs=c('age','cohort','pdl1','change_ctdna_group'))

forestplotUV(data=pembrolizumab, response='orr',
covs=c('age','sex','pdl1','change_ctdna_group'))

require(ggplot2)
glm_fit <- glm(orr~change_ctdna_group+pdl1,
               family='binomial',
               data = pembrolizumab)
forestplotMV(glm_fit)

Plotting combined univariable and multivariable model odd ratios

To display odds ratios from univariate and multivariable models in a single forest plot, run the forestplotUV and forestplotMV functions and combine.

uvFP <- forestplotUV(data=pembrolizumab, response='orr',
covs=c('age','sex','pdl1','change_ctdna_group'))

glm_fit <- glm(orr~change_ctdna_group+pdl1,
               family='binomial',
               data = pembrolizumab)
mvFP <- forestplotMV(glm_fit)

forestplotUVMV(uvFP,mvFP,showN=T,showEvent=T)

This can also be done with log scale odds ratios (default is linear scale). Number of subjects and/or number of events can also be turned off, as well as different colours used.

uvFP <- forestplotUV(data=pembrolizumab, response='orr',
covs=c('age','sex','pdl1','change_ctdna_group'))

glm_fit <- glm(orr~change_ctdna_group+pdl1,
               family='binomial',
               data = pembrolizumab)
mvFP <- forestplotMV(glm_fit)

forestplotUVMV(uvFP,mvFP,showN=F,showEvent=F,colours=c("orange","black","blue"),logScale=T)

ggkmcif(response = c('os_time','os_status'),
cov='cohort',
data=pembrolizumab)

Options

The following options can be set:

• reportRmd.digits, detault is 2. Sets the default number of digits in output functions. Does not affect p-values.
• reportRmd.forceWald, default is FALSE. If set to TRUE then Wald CI, as opposed to likelihood profile CI will be used in the rm_uvsum function.
• reportRmd.logScale, default is TRUE. The scale of the forest plots will default to log unless otherwise specified.

Example:

 rm_uvsum(response = 'baseline_ctdna',
 covs=c('age','sex','l_size','pdl1','tmb'),
 data=pembrolizumab)

 options('reportRmd.digits'=1) 
 
rm_uvsum(response = 'baseline_ctdna',
 covs=c('age','sex','l_size','pdl1','tmb'),
 data=pembrolizumab)

PDF Output

For pdf to be correctly generate when using survival curves it is recommended that the cairo format be used. This can be specified with the following command in the setup code chunk:

knitr::opts_chunk$set(message = FALSE, warning = FALSE,dev="cairo_pdf")

Data Sets