Summaries By Group

This will produce summary statistics by Sex

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

Showing Statistical Tests

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)

Including Effect Sizes

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, Omega for ANOVA, Epsilon for Kruskal Wallis test, and Cramer’s V for categorical variables.

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

Parametric vs Non-parametric Comparisons

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)

Row vs Column Summaries

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')

Using Mean Summary Statistic Instead of Median

The default summary statistic for numerical variables is median. To specify which numerical variables should have mean displayed instead, change the use_mean argument. Unspecified xvars will use the default median

rm_compactsum(data=pembrolizumab, xvars=c('age','pdl1','change_ctdna_group'), grp = 'sex', use_mean = c("pdl1"), effSize = T, show.tests = T)

Custom Digits Argument

The digits and digits.cat arguments can be changed to a custom numerical value. The default digits is 1, the default digits.cat is 0

rm_compactsum(data=pembrolizumab, xvars=c('age','pdl1','change_ctdna_group'), grp = 'sex', use_mean = c("pdl1"), digits = 2, digits.cat = 1, effSize = T, show.tests = T)

To specify custom digit arguments for different numerical variables, change the digits argument. Unspecified variables will use the default value.

rm_compactsum(data=pembrolizumab, xvars=c('age','pdl1','l_size'), grp = 'sex', digits = c("age" = 3, "l_size" = 2), effSize = T, show.tests = T)

Viewing the Self-Generated Description

To view the self-generated description for the summary table:

summary_tab <- rm_compactsum(data=pembrolizumab, xvars=c('change_ctdna_group','orr', 'age'), grp = 'cohort', effSize = T, show.tests = T)
cat(attr(summary_tab, "description"))

Using Tidyselect

Both rm_covsum and rm_compactsum support tidyselect syntax, allowing you to use bare column names instead of character strings. This also enables piping with |>:

pembrolizumab |> rm_compactsum(xvars = c(age, pdl1, change_ctdna_group), grp = sex)

The xvars/grp aliases can also be used with rm_covsum for a consistent interface:

pembrolizumab |> rm_covsum(xvars = c(age, pdl1, change_ctdna_group), grp = sex)

Adding Header Rows

Header rows can be added to any table using the header_above argument in outTable(). This works for HTML, PDF and Word output. The argument takes a named vector where the names are the header labels and the values are the number of columns each label spans.

tab <- rm_covsum(data = pembrolizumab, maincov = cohort,
                  covs = c(age, pdl1, sex), full = FALSE,
                  pvalue = FALSE, tableOnly = TRUE)

outTable(tab, header_above = c(" " = 1, "Cohorts A-C" = 3, "Cohorts D,E" = 2))

Simple Linear Regression

If the response is continuous linear regression is the default. Using type = 'linear' will ensure linear regression.

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

Logistic Regression

If the response is binomial, logistic regression will be run (or specified with type = 'logistic').

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

Poisson Regression

If the response is integer, poisson regression will be run (or specified with type = 'poisson').

pembrolizumab$Counts <- rpois(nrow(pembrolizumab),lambda = 3)
rm_uvsum(data=pembrolizumab, response='Counts',covs=c('age','cohort'))

offset terms can be specified as well, but must correspond to a variable in the data set

pembrolizumab$length_followup <- rnorm(nrow(pembrolizumab),mean = 72,sd=3)
pembrolizumab$log_length_followup <- log(pembrolizumab$length_followup)
rm_uvsum(data=pembrolizumab, response='Counts',covs=c('age','cohort'),
         offset = "log_length_followup")

Negative Binomial Regression

To run negative binomial regression instead specify type = 'negbin'

rm_uvsum(data=pembrolizumab, response='Counts', type='negbin',
         covs=c('age','cohort'),
         offset = "log_length_followup")

Survival Analysis

If two response variables are specified and then survival analysis is run (specified with type='coxph').

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

Clustered analyses can also be run

data("diabetic", package = "survival")
rm_uvsum(data=diabetic, response=c('time','status'),
covs=c('trt','laser','age','eye'),id='id')

Competing Risk

Competing risk models need to be explicitly specified using type='crr'.

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

GEE Models

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)

MICE output

Outputs from multiple imputations can be summarised as follows (not run):

miss_data <- pembrolizumab
miss_data$age[sample(1:nrow(pembrolizumab),50,replace=F)] <- NA
miss_data$cohort[sample(1:nrow(pembrolizumab),50,replace=F)] <- NA
require(mice)
ini <- mice::mice(miss_data, maxit = 0)
pred_matrix <- ini$predictorMatrix
pred_matrix[1, ] <- 0
pred_matrix[, 1] <- 0
imp <- mice::mice(miss_data, m = 5, predictorMatrix = pred_matrix, printFlag = FALSE)
imputed_datasets <- purrr::map(1:5, \(m) mice::complete(imp, m))

fits <- lapply(imputed_datasets, function(dat) {
  glm(orr ~ age+cohort, data = dat,family='binomial')
})
model <- mice::as.mira(fits)
rm_mvsum(model)           

Returning Model Objects

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

 rm_uvsum(response = 'orr',
 covs=c('age'),
 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

The data analysed can be examined by interrogating the data object appended to each model

mList <-  rm_uvsum(response = 'orr',
 covs=c('age'),
 data=pembrolizumab,returnModels = TRUE)
head(mList$age$data)

## # A tibble: 6 × 18
##   id       age sex   cohort l_size  pdl1   tmb baseline_ctdna change_ctdna_group
##   <fct>  <dbl> <fct> <fct>   <dbl> <dbl> <dbl>          <dbl> <fct>             
## 1 INS-A…  62.1 Male  A         140     2 0.574         114.   Increase from bas…
## 2 INS-A…  62.2 Male  A         108     0 0.921           4.26 Increase from bas…
## 3 INS-A…  70.9 Male  A          11   100 0.375         205.   Decrease from bas…
## 4 INS-A…  57.9 Male  A          39    45 0.824         169.   Increase from bas…
## 5 INS-A…  65.7 Fema… A          90     0 0.114         425.   Decrease from bas…
## 6 INS-A…  60.9 Male  A         124     2 0.717          15.8  Increase from bas…
## # ℹ 9 more variables: orr <fct>, cbr <fct>, os_status <dbl>, os_time <dbl>,
## #   pfs_status <dbl>, pfs_time <dbl>, Counts <int>, length_followup <dbl>,
## #   log_length_followup <dbl>

Adjusting p-values

Multiple comparisons can be controlled for with the p.adjust argument, which accepts any of the options from the p.adjust function.

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

Note: The raw p-value column is suppressed when there are categorical variables with >2 levels, to prevent three columns of p-values from appearing.

Including Unadjusted Estimates

As of version 0.1.3, you can include univariate (unadjusted) estimates alongside the multivariable estimates in a single table by setting include_unadjusted = TRUE. This provides a convenient way to compare unadjusted and adjusted estimates side-by-side.

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

The data are taken from the model itself, so the data argument is not needed. This includes Cox proportional hazards models, where the response is a Surv() object:

require(survival)
cox_fit <- coxph(Surv(os_time, os_status) ~ age + sex + baseline_ctdna,
                 data = pembrolizumab)
rm_mvsum(cox_fit, include_unadjusted = TRUE)

Transformed terms are summarised on the scale they were fit on, so the unadjusted estimate for log(l_size) below is the slope from glm(orr ~ log(l_size)), not from glm(orr ~ l_size):

rm_mvsum(glm(orr ~ log(l_size) + sex, family = 'binomial', data = pembrolizumab),
         include_unadjusted = TRUE)

Two points are worth knowing when reading these tables.

The univariate models are fit on the same complete-case sample as the multivariable model, so the two sets of estimates are directly comparable and the N column applies to both. This differs from calling rm_uvsum() yourself, which uses all rows available for each covariate.

A covariate that appears only inside an interaction term has no adjusted estimate of its own. Such covariates are appended to the end of the table with an unadjusted estimate and a blank adjusted column, rather than being dropped:

rm_mvsum(glm(orr ~ pdl1 + age:sex, family = 'binomial', data = pembrolizumab),
         include_unadjusted = TRUE)

Terms that expand to a matrix, such as poly(age, 2), can not be fit univariately; these are skipped with a warning and their unadjusted cells are left blank.

Scrollable Tables

If you are rendering to html format you can use the scrollTable function to create a scrollable table. This is useful for tables with many rows. If the output format is not html then a regular table will be displayed.

long_table <- rm_compactsum(data=pembrolizumab,xvars = c(age,sex,cohort,pdl1,tmb,baseline_ctdna,change_ctdna_group,orr,cbr))
scrolling_table(long_table,pixelHeight = 300)

Including Median Follow-up

The reverse Kaplan-Meier follow-up time can also be displayed.

rm_survsum(data=pembrolizumab,time='os_time',status='os_status',
 group="sex",survtimes=seq(12,24,12),survtimeunit='months',
 add.reverse.KM = T)

Survival Times With Covariate Adjustments

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')

Worked Example

Get some descriptive stats for the ctDNA data that comes with the package. The nicenames argument is TRUE by default so underscores are replaced by spaces

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

library(ggplot2)
p <- ggplot(data=ctDNA,aes(x=ctdna_status,y=size_change,colour=cohort))+
  geom_point()
replace_plot_labels(p)

ctDNA <- clear_labels(ctDNA)

Choosing Which Statistics to Report

The output argument selects statistics by name, and accepts several machine-learning measures that are not shown by default:

predsum(true = d$change_ctdna_group, predicted = d$predicted, cutoff = 0.5,
        output = c('accuracy', 'brier', 'f1', 'auc'))

Brier score, F1 and AUC-PR have no closed-form confidence interval, so predsum() obtains theirs by bootstrap. This is slower than the other statistics; nboot controls the number of replicates (default 2000).

Plotting bivariate relationships

These plots are designed for quick inspection of many variables, not for publication. This is the plotting version of rm_uvsum. As of 0.1.1 the variable names will be replaced by variable labels if they exist.

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

The plotuv function can also be used without a response variable to display summary of variables

plotuv(data = pembrolizumab, covs=c('age','cohort','pdl1','change_ctdna_group'), showN = T)

Plotting survival curves

Survival curves are now ggplot2-based, the older version, ggkmcif is deprecated from version 0.1.0

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

Forest Plots

Similar to rm_uvsum and rm_mvsum, forest plots can be created from univariate or multivariable models. forestplot2 is deprecated from version 0.1.0. Variable labels are not yet incorporated into the forest plots.

This will default to a log scale, but can be set to linear using logScale=FALSE

forestplotUV(response="orr", covs=c("change_ctdna_group", "sex", "age", "l_size"),
data=pembrolizumab, family='binomial')

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

glm_fit <- glm(orr~pdl1+age,
               family='binomial',
               data = pembrolizumab)
forestplotMV(glm_fit, include_unadjusted = TRUE)

excelCol

To identify the number of a column given the Excel column header

excelCol(G,AB,Az)

##  G AB AZ 
##  7 28 52

excelColLetters

To identify the Excel header given a columns number

excelColLetters(c(7,28,52))

##    7   28   52 
##  "G" "AB" "AZ"

extract_package_details

This function will search the current script or markdown file for R packages and record the package name, functions called, the package version and the citation(s) for the package. By default, function calls and packages that only appear in commented code will be ignored, but they can be included by setting ignore_comments = FALSE

extract_package_details() |> 
  outTable()

extract_package_details(ignore_comments = FALSE) |> 
  outTable()

Monte Carlo p-values

rm_covsum() and rm_compactsum() compare categorical variables with a chi-squared test, switching to Fisher’s exact test when any expected cell count falls below five. On a large or sparse table the exact test can exceed its workspace, in which case both functions fall back to a Monte Carlo approximation and say so:

Fisher's exact test could not be computed for site; using Monte Carlo
simulation with B = 1000 and seed = 1.

B sets the number of replicates and seed fixes the random state for the duration of the call; the session’s own random state is restored afterwards, so seeding here does not disturb anything else in your document. Setting seed = NULL restores the older behaviour in which the p-value varies between calls.

rm_covsum(data = pembrolizumab, maincov = 'cohort', covs = c('age', 'sex'),
          B = 5000, seed = 2024)

When a Monte Carlo p-value has been used, show.tests = TRUE records the replicate count alongside the test name, so the table itself documents how the p-value was obtained.

Because both functions now use the same expected-count rule, rm_covsum() and rm_compactsum() select the same test and report the same p-value for the same data.

Bootstrap confidence intervals

predsum() bootstraps the intervals for Brier score, F1 and AUC-PR. It takes the same seed argument, with the same default and the same restoration behaviour:

predsum(true = d$change_ctdna_group, predicted = d$predicted, cutoff = 0.5,
        output = c('brier', 'f1'), nboot = 5000, seed = 2024)

pembrolizumab

Survival status and ctDNA levels for patients receiving pembrolizumab

A data frame with 94 rows and 15 variables:

• id Patient ID
• age Age at study entry
• sex Patient Sex
• cohort Study Cohort: A = Squamous cell carcinoma of soft pallate, B = Triple negative breast cancer, C = Ovarian, high grade serous, D = Melanoma, E = Other Solid Tumor
• l_size Target lesion size at baseline
• pdl1 PD L1 percent
• tmb log of TMB
• baseline_ctdna Baseline ctDNA
• change_ctdna_group Did ctDNA increase or decrease from baseline to cycle 3
• orr Objective Response
• cbr Clinical Beneficial Response
• os_status Overall survival status, 0 = alive, 1 = deceased
• os_time Overall survival time in months
• pfs_status Progression free survival status, 0 = progression free, 1 = progressed
• pfs_time Progression free survival time in months

ctDNA

Longitudinal changes in tumour size since baseline for patients by changes in ctDNA status (clearance, decrease or increase) since baseline.

A data frame with 270 rows and 5 variables:

• id Patient ID
• cohort Study Cohort: A = Squamous cell carcinoma of soft pallate, B = Triple negative breast cancer, C = Ovarian, high grade serous, D = Melanoma, E = Other Solid Tumor
• ctdna_status Change in ctDNA since baseline
• time Number of weeks on treatment
• size_change Percentage change in tumour measurement