reportRmd Package
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'))
|
|
n=94
|
|
Age at study entry
|
|
|
Mean (sd)
|
57.9 (12.8)
|
|
Median (Min,Max)
|
59.1 (21.1, 81.8)
|
|
Patient Sex
|
|
|
Female
|
58 (62)
|
|
Male
|
36 (38)
|
Set IQR = T for interquartile range instead of
Min/Max
rm_covsum(data=pembrolizumab,
covs=c('age','sex'),IQR=TRUE)
|
|
n=94
|
|
Age at study entry
|
|
|
Mean (sd)
|
57.9 (12.8)
|
|
Median (Q1,Q3)
|
59.1 (49.5, 68.7)
|
|
Patient Sex
|
|
|
Female
|
58 (62)
|
|
Male
|
36 (38)
|
Or all.stats=T for both
rm_covsum(data=pembrolizumab,
covs=c('age','sex'),all.stats = TRUE)
|
|
n=94
|
|
Age at study entry
|
|
|
Mean (sd)
|
57.9 (12.8)
|
|
Median (Q1,Q3)
|
59.1 (49.5, 68.7)
|
|
Range (min, max)
|
(21.1, 81.8)
|
|
Patient Sex
|
|
|
Female
|
58 (62)
|
|
Male
|
36 (38)
|
Summaries By Group
This will produce summary statistics by Sex
rm_covsum(data=pembrolizumab, maincov = 'sex',
covs=c('age','pdl1','change_ctdna_group'))
|
|
Full Sample (n=94)
|
Female (n=58)
|
Male (n=36)
|
p-value
|
|
Age at study entry
|
|
|
|
0.30
|
|
Mean (sd)
|
57.9 (12.8)
|
56.9 (12.6)
|
59.3 (13.1)
|
|
|
Median (Min,Max)
|
59.1 (21.1, 81.8)
|
56.6 (34.1, 78.2)
|
61.2 (21.1, 81.8)
|
|
|
PD L1 percent
|
|
|
|
0.76
|
|
Mean (sd)
|
13.9 (29.2)
|
15.0 (30.5)
|
12.1 (27.3)
|
|
|
Median (Min,Max)
|
0 (0, 100)
|
0.5 (0.0, 100.0)
|
0 (0, 100)
|
|
|
Missing
|
1
|
0
|
1
|
|
|
Did ctDNA increase or decrease from
baseline to cycle 3
|
|
|
|
0.84
|
|
Decrease from baseline
|
33 (45)
|
19 (48)
|
14 (42)
|
|
|
Increase from baseline
|
40 (55)
|
21 (52)
|
19 (58)
|
|
|
Missing
|
21
|
18
|
3
|
|
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)
|
|
Full Sample (n=94)
|
Female (n=58)
|
Male (n=36)
|
p-value
|
StatTest
|
|
Age at study entry
|
|
|
|
0.30
|
Wilcoxon Rank Sum
|
|
Mean (sd)
|
57.9 (12.8)
|
56.9 (12.6)
|
59.3 (13.1)
|
|
|
|
Median (Min,Max)
|
59.1 (21.1, 81.8)
|
56.6 (34.1, 78.2)
|
61.2 (21.1, 81.8)
|
|
|
|
PD L1 percent
|
|
|
|
0.76
|
Wilcoxon Rank Sum
|
|
Mean (sd)
|
13.9 (29.2)
|
15.0 (30.5)
|
12.1 (27.3)
|
|
|
|
Median (Min,Max)
|
0 (0, 100)
|
0.5 (0.0, 100.0)
|
0 (0, 100)
|
|
|
|
Missing
|
1
|
0
|
1
|
|
|
|
Did ctDNA increase or decrease from
baseline to cycle 3
|
|
|
|
0.84
|
Chi Sq
|
|
Decrease from baseline
|
33 (45)
|
19 (48)
|
14 (42)
|
|
|
|
Increase from baseline
|
40 (55)
|
21 (52)
|
19 (58)
|
|
|
|
Missing
|
21
|
18
|
3
|
|
|
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)
|
|
Full Sample (n=94)
|
Female (n=58)
|
Male (n=36)
|
p-value
|
Effect Size
|
|
Age at study entry
|
|
|
|
0.30
|
0.11
|
|
Mean (sd)
|
57.9 (12.8)
|
56.9 (12.6)
|
59.3 (13.1)
|
|
|
|
Median (Min,Max)
|
59.1 (21.1, 81.8)
|
56.6 (34.1, 78.2)
|
61.2 (21.1, 81.8)
|
|
|
|
Did ctDNA increase or decrease from
baseline to cycle 3
|
|
|
|
0.84
|
0.020
|
|
Decrease from baseline
|
33 (45)
|
19 (48)
|
14 (42)
|
|
|
|
Increase from baseline
|
40 (55)
|
21 (52)
|
19 (58)
|
|
|
|
Missing
|
21
|
18
|
3
|
|
|
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)
|
|
Full Sample (n=94)
|
Female (n=58)
|
Male (n=36)
|
p-value
|
Effect Size
|
StatTest
|
|
Age at study entry
|
|
|
|
0.39
|
0.18
|
t-test, Cohen’s d
|
|
Mean (sd)
|
57.9 (12.8)
|
56.9 (12.6)
|
59.3 (13.1)
|
|
|
|
|
Median (Min,Max)
|
59.1 (21.1, 81.8)
|
56.6 (34.1, 78.2)
|
61.2 (21.1, 81.8)
|
|
|
|
|
PD L1 percent
|
|
|
|
0.63
|
0.100
|
t-test, Cohen’s d
|
|
Mean (sd)
|
13.9 (29.2)
|
15.0 (30.5)
|
12.1 (27.3)
|
|
|
|
|
Median (Min,Max)
|
0 (0, 100)
|
0.5 (0.0, 100.0)
|
0 (0, 100)
|
|
|
|
|
Missing
|
1
|
0
|
1
|
|
|
|
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)
|
|
Full Sample (n=94)
|
Female (n=58)
|
Male (n=36)
|
|
Study Cohort
|
|
|
|
|
A
|
16 (17)
|
3 (5)
|
13 (36)
|
|
B
|
18 (19)
|
18 (31)
|
0 (0)
|
|
C
|
18 (19)
|
18 (31)
|
0 (0)
|
|
D
|
12 (13)
|
7 (12)
|
5 (14)
|
|
E
|
30 (32)
|
12 (21)
|
18 (50)
|
But you can also specify to show by row instead
rm_covsum(data=pembrolizumab, maincov = 'sex',
covs=c('cohort'),
pvalue = FALSE,
percentage='row')
|
|
Full Sample (n=94)
|
Female (n=58)
|
Male (n=36)
|
|
Study Cohort
|
|
|
|
|
A
|
16
|
3 (19)
|
13 (81)
|
|
B
|
18
|
18 (100)
|
0 (0)
|
|
C
|
18
|
18 (100)
|
0 (0)
|
|
D
|
12
|
7 (58)
|
5 (42)
|
|
E
|
30
|
12 (40)
|
18 (60)
|
Compact summary statistics
Basic summary statistics
rm_compactsum(data = pembrolizumab, xvars = c("change_ctdna_group", "l_size"))
|
|
Full Sample (n=94)
|
Missing
|
|
Did ctDNA increase or decrease from
baseline to cycle 3 - Increase from baseline
|
40 (55%)
|
21
|
|
Target lesion size at baseline
|
73.5 (49.2-108.8)
|
0
|
Set iqr = T for interquartile range instead of
Min/Max
rm_compactsum(data=pembrolizumab, xvars = c("change_ctdna_group", "l_size"), iqr=TRUE)
|
|
Full Sample (n=94)
|
Missing
|
|
Did ctDNA increase or decrease from
baseline to cycle 3 - Increase from baseline
|
40 (55%)
|
21
|
|
Target lesion size at baseline
|
73.5 (49.2-108.8)
|
0
|
Or all.stats=T for both
rm_compactsum(data=pembrolizumab, xvars = c("change_ctdna_group", "l_size"), all.stats = T)
|
|
Full Sample (n=94)
|
Missing
|
|
Did ctDNA increase or decrease from
baseline to cycle 3 - Increase from baseline
|
40 (55%)
|
21
|
|
Target lesion size at baseline
|
|
0
|
|
Mean (sd)
|
87.9 (59.6)
|
|
|
Median (Q1-Q3)
|
73.5 (49.2-108.8)
|
|
|
Range (min-max)
|
(11.0-387.0)
|
|
Summaries By Group
This will produce summary statistics by Sex:
rm_compactsum(data=pembrolizumab, xvars=c('age','pdl1','change_ctdna_group'), grp = 'sex')
|
|
Full Sample (n=94)
|
Female (n=58)
|
Male (n=36)
|
p-value
|
Missing
|
|
Age at study entry
|
59.1 (49.5-68.7)
|
56.6 (45.8-67.8)
|
61.2 (52.0-69.4)
|
0.30
|
0
|
|
PD L1 percent
|
0.0 (0.0-10.0)
|
0.5 (0.0-13.8)
|
0.0 (0.0-4.5)
|
0.76
|
1
|
|
Did ctDNA increase or decrease from
baseline to cycle 3 - Increase from baseline
|
40 (55%)
|
21 (52%)
|
19 (58%)
|
0.84
|
21
|
Showing Statistical Tests
To indicate which statistical test was used use
show.tests=TRUE
rm_compactsum(data=pembrolizumab, xvars=c('age','pdl1','change_ctdna_group'), grp = 'sex', show.tests=TRUE)
|
|
Full Sample (n=94)
|
Female (n=58)
|
Male (n=36)
|
p-value
|
Missing
|
pTest
|
|
Age at study entry
|
59.1 (49.5-68.7)
|
56.6 (45.8-67.8)
|
61.2 (52.0-69.4)
|
0.30
|
0
|
Wilcoxon Rank Sum
|
|
PD L1 percent
|
0.0 (0.0-10.0)
|
0.5 (0.0-13.8)
|
0.0 (0.0-4.5)
|
0.76
|
1
|
Wilcoxon Rank Sum
|
|
Did ctDNA increase or decrease from
baseline to cycle 3 - Increase from baseline
|
40 (55%)
|
21 (52%)
|
19 (58%)
|
0.84
|
21
|
ChiSq
|
Including Effect Sizes
Effect sizes can be added with effSize = TRUE. If
show.tests = TRUE as well, the effStat will also be
shown:
rm_compactsum(data=pembrolizumab, xvars=c('age','pdl1','change_ctdna_group'), grp = 'sex', effSize = T, show.tests = T)
|
|
Full Sample (n=94)
|
Female (n=58)
|
Male (n=36)
|
p-value
|
Effect Size (95% CI)
|
Missing
|
pTest
|
effStat
|
|
Age at study entry
|
59.1 (49.5-68.7)
|
56.6 (45.8-67.8)
|
61.2 (52.0-69.4)
|
0.30
|
0.11 (-0.09, 0.21)
|
0
|
Wilcoxon Rank Sum
|
Wilcoxon r
|
|
PD L1 percent
|
0.0 (0.0-10.0)
|
0.5 (0.0-13.8)
|
0.0 (0.0-4.5)
|
0.76
|
0.03 (-0.18, 0.06)
|
1
|
Wilcoxon Rank Sum
|
Wilcoxon r
|
|
Did ctDNA increase or decrease from
baseline to cycle 3 - Increase from baseline
|
40 (55%)
|
21 (52%)
|
19 (58%)
|
0.84
|
0.02 (-0.19, 0.05)
|
21
|
ChiSq
|
Cramer’s V
|
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)
|
|
Full Sample (n=94)
|
Female (n=58)
|
Male (n=36)
|
p-value
|
Effect Size (95% CI)
|
Missing
|
pTest
|
effStat
|
|
Age at study entry
|
59.12 (49.54-68.71)
|
56.64 (45.80-67.83)
|
61.16 (51.99-69.44)
|
0.30
|
0.11 (-0.09, 0.21)
|
0
|
Wilcoxon Rank Sum
|
Wilcoxon r
|
|
PD L1 percent
|
13.90 (29.24)
|
15.02 (30.55)
|
12.06 (27.26)
|
0.63
|
0.10 (-0.38, 0.20)
|
1
|
t-test
|
Cohen’s d
|
|
Did ctDNA increase or decrease from
baseline to cycle 3 - Increase from baseline
|
40 (54.8%)
|
21 (52.5%)
|
19 (57.6%)
|
0.84
|
0.02 (-0.17, 0.05)
|
21
|
ChiSq
|
Cramer’s V
|
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)
|
|
Full Sample (n=94)
|
Female (n=58)
|
Male (n=36)
|
p-value
|
Effect Size (95% CI)
|
Missing
|
pTest
|
effStat
|
|
Age at study entry
|
59.121 (49.542-68.708)
|
56.643 (45.800-67.834)
|
61.164 (51.992-69.438)
|
0.30
|
0.11 (-0.08, 0.21)
|
0
|
Wilcoxon Rank Sum
|
Wilcoxon r
|
|
PD L1 percent
|
0.0 (0.0-10.0)
|
0.5 (0.0-13.8)
|
0.0 (0.0-4.5)
|
0.76
|
0.03 (-0.16, 0.06)
|
1
|
Wilcoxon Rank Sum
|
Wilcoxon r
|
|
Target lesion size at baseline
|
73.50 (49.25-108.75)
|
68.00 (44.25-97.75)
|
93.00 (65.50-121.00)
|
0.066
|
0.19 (0.00, 0.36)
|
0
|
Wilcoxon Rank Sum
|
Wilcoxon r
|
Row vs Column Summaries
The default is to indicate percentages by columns (ie. percentages
within columns add to 100). But you can also specify to show by row
instead
rm_compactsum(data=pembrolizumab, xvars=c('change_ctdna_group','orr'), grp = 'cohort', effSize = T, show.tests = T, percentage = "row")
|
|
Full Sample (n=94)
|
A (n=16)
|
B (n=18)
|
C (n=18)
|
D (n=12)
|
E (n=30)
|
p-value
|
Effect Size (95% CI)
|
pTest
|
effStat
|
Missing
|
|
Did ctDNA increase or decrease from
baseline to cycle 3 - Increase from baseline
|
40
|
8 (20%)
|
7 (18%)
|
5 (12%)
|
2 (5%)
|
18 (45%)
|
0.18
|
0.30 (0.20, 0.51)
|
Fisher Exact
|
Cramer’s V
|
21
|
|
Objective Response - SD/PD
|
78
|
13 (17%)
|
18 (23%)
|
18 (23%)
|
4 (5%)
|
25 (32%)
|
<0.001
|
0.55 (0.76, 1.03)
|
Fisher Exact
|
Cramer’s V
|
0
|
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"))
Univariate regression
Combining multiple univariate regression analyses into a single
table. The function will try to determine the most appropriate model
from the data.
rm_uvsum(data=pembrolizumab, response='orr',
covs=c('age','pdl1','change_ctdna_group'))
|
|
OR(95%CI)
|
p-value
|
N
|
Event
|
|
Age at study entry
|
0.96 (0.91, 1.00)
|
0.089
|
94
|
78
|
|
PD L1 percent
|
0.97 (0.95, 0.98)
|
<0.001
|
93
|
77
|
|
Did ctDNA increase or decrease from
baseline to cycle 3
|
|
|
73
|
58
|
|
Decrease from baseline
|
Reference
|
|
33
|
19
|
|
Increase from baseline
|
28.74 (5.20, 540.18)
|
0.002
|
40
|
39
|
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'))
|
|
Estimate(95%CI)
|
p-value
|
N
|
|
Age at study entry
|
-0.58 (-1.54, 0.38)
|
0.23
|
94
|
|
Study Cohort
|
|
|
94
|
|
A
|
Reference
|
|
16
|
|
B
|
-38.04 (-74.95, -1.13)
|
0.044
|
18
|
|
C
|
20.35 (-16.56, 57.26)
|
0.28
|
18
|
|
D
|
-24.79 (-65.82, 16.23)
|
0.23
|
12
|
|
E
|
31.69 (-1.56, 64.95)
|
0.062
|
30
|
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'))
|
|
OR(95%CI)
|
p-value
|
N
|
Event
|
|
Age at study entry
|
0.96 (0.91, 1.00)
|
0.089
|
94
|
78
|
|
Study Cohort
|
|
|
94
|
78
|
|
A
|
Reference
|
|
16
|
13
|
|
B
|
7.3e+07 (8.5e-76, NA)
|
0.99
|
18
|
18
|
|
C
|
7.3e+07 (7.7e-74, NA)
|
0.99
|
18
|
18
|
|
D
|
0.12 (0.02, 0.60)
|
0.015
|
12
|
4
|
|
E
|
1.15 (0.21, 5.49)
|
0.86
|
30
|
25
|
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'))
|
|
RR(95%CI)
|
p-value
|
N
|
|
Age at study entry
|
1.00 (0.99, 1.01)
|
0.51
|
94
|
|
Study Cohort
|
|
|
94
|
|
A
|
Reference
|
|
16
|
|
B
|
0.89 (0.62, 1.28)
|
0.53
|
18
|
|
C
|
0.97 (0.68, 1.38)
|
0.85
|
18
|
|
D
|
0.94 (0.63, 1.40)
|
0.77
|
12
|
|
E
|
0.83 (0.60, 1.16)
|
0.26
|
30
|
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")
|
|
RR(95%CI)
|
p-value
|
N
|
|
Age at study entry
|
1.00 (0.99, 1.01)
|
0.53
|
94
|
|
Study Cohort
|
|
|
94
|
|
A
|
Reference
|
|
16
|
|
B
|
0.89 (0.62, 1.28)
|
0.52
|
18
|
|
C
|
0.99 (0.69, 1.42)
|
0.97
|
18
|
|
D
|
0.95 (0.64, 1.42)
|
0.81
|
12
|
|
E
|
0.83 (0.60, 1.16)
|
0.27
|
30
|
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")
|
|
RR(95%CI)
|
p-value
|
N
|
|
Age at study entry
|
1.00 (0.99, 1.01)
|
0.54
|
94
|
|
Study Cohort
|
|
|
94
|
|
A
|
Reference
|
|
16
|
|
B
|
0.89 (0.61, 1.29)
|
0.53
|
18
|
|
C
|
0.99 (0.69, 1.43)
|
0.97
|
18
|
|
D
|
0.95 (0.63, 1.42)
|
0.82
|
12
|
|
E
|
0.83 (0.60, 1.17)
|
0.28
|
30
|
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")
|
|
HR(95%CI)
|
p-value
|
N
|
Event
|
|
Age at study entry
|
0.99 (0.97, 1.01)
|
0.16
|
94
|
64
|
|
PD L1 percent
|
0.99 (0.98, 1.00)
|
0.026
|
93
|
63
|
|
Did ctDNA increase or decrease from
baseline to cycle 3
|
|
|
73
|
46
|
|
Decrease from baseline
|
Reference
|
|
33
|
14
|
|
Increase from baseline
|
3.06 (1.62, 5.77)
|
<0.001
|
40
|
32
|
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')
|
|
HR(95%CI)
|
p-value
|
N
|
Event
|
|
Age at study entry
|
0.99 (0.97, 1.00)
|
0.15
|
94
|
64
|
|
PD L1 percent
|
0.99 (0.98, 1.00)
|
0.017
|
93
|
63
|
|
Did ctDNA increase or decrease from
baseline to cycle 3
|
|
|
73
|
46
|
|
Decrease from baseline
|
Reference
|
|
33
|
14
|
|
Increase from baseline
|
3.06 (1.64, 5.69)
|
<0.001
|
40
|
32
|
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)
|
|
Estimate(95%CI)
|
p-value
|
N
|
|
Number of weeks on treatment
|
-0.12 (-0.44, 0.19)
|
0.44
|
262
|
|
Change in ctDNA since baseline
|
|
|
262
|
|
Clearance
|
Reference
|
|
134
|
|
No clearance, decrease from baseline
|
61.29 (37.49, 85.09)
|
<0.001
|
42
|
|
No clearance, increase from baseline
|
82.52 (64.75, 100.28)
|
<0.001
|
86
|
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')
|
|
OR(95%CI)
|
p-value
|
N
|
Event
|
|
Age at study entry
|
0.96 (0.91, 1.00)
|
0.11
|
94
|
78
|
|
Patient Sex
|
|
|
94
|
78
|
|
Female
|
Reference
|
|
58
|
51
|
|
Male
|
0.41 (0.13, 1.22)
|
0.11
|
36
|
27
|
|
PD L1 percent
|
0.97 (0.95, 0.98)
|
<0.001
|
93
|
77
|
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.
Multivariable analysis
To create a nice display for multivariable models the multivariable
model first needs to 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, showN = TRUE, vif=TRUE)
|
|
OR(95%CI)
|
p-value
|
N
|
Event
|
VIF
|
|
Did ctDNA increase or decrease from
baseline to cycle 3
|
|
|
73
|
58
|
1.03
|
|
Decrease from baseline
|
Reference
|
|
33
|
19
|
|
|
Increase from baseline
|
23.92 (3.69, 508.17)
|
0.006
|
40
|
39
|
|
|
PD L1 percent
|
0.97 (0.95, 0.99)
|
0.011
|
73
|
58
|
1.24
|
|
Age at study entry
|
0.94 (0.87, 1.00)
|
0.078
|
73
|
58
|
1.23
|
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')
|
|
OR(95%CI)
|
p-value
|
N
|
Event
|
VIF
|
|
Did ctDNA increase or decrease from
baseline to cycle 3
|
|
|
73
|
58
|
1.03
|
|
Decrease from baseline
|
Reference
|
|
33
|
19
|
|
|
Increase from baseline
|
23.92 (3.69, 508.17)
|
0.018
|
40
|
39
|
|
|
PD L1 percent
|
0.97 (0.95, 0.99)
|
0.022
|
73
|
58
|
1.24
|
|
Age at study entry
|
0.94 (0.87, 1.00)
|
0.078
|
73
|
58
|
1.23
|
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)
|
|
Unadjusted OR(95%CI)
|
p
|
Adjusted OR(95%CI)
|
p (adj)
|
|
Age at study entry
|
0.96 (0.91, 1.00)
|
0.089
|
|
|
|
Patient Sex
|
|
|
|
|
|
Female
|
Reference
|
|
|
|
|
Male
|
0.41 (0.13, 1.22)
|
0.11
|
|
|
|
PD L1 percent
|
0.97 (0.95, 0.98)
|
<0.001
|
0.98 (0.95, 1.00)
|
0.024
|
|
Did ctDNA increase or decrease from
baseline to cycle 3
|
|
|
|
|
|
Decrease from baseline
|
Reference
|
|
Reference
|
|
|
Increase from baseline
|
28.74 (5.20, 540.18)
|
0.002
|
24.71 (4.19, 479.13)
|
0.004
|
Note: 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)
|
|
Unadjusted OR(95%CI)
|
p
|
Adjusted OR(95%CI)
|
p (adj)
|
|
Age at study entry
|
0.96 (0.91, 1.00)
|
0.18
|
|
|
|
Patient Sex
|
|
|
|
|
|
Female
|
Reference
|
|
|
|
|
Male
|
0.41 (0.13, 1.22)
|
0.18
|
|
|
|
PD L1 percent
|
0.97 (0.95, 0.98)
|
<0.001
|
0.98 (0.95, 1.00)
|
0.024
|
|
Did ctDNA increase or decrease from
baseline to cycle 3
|
|
|
|
|
|
Decrease from baseline
|
Reference
|
|
Reference
|
|
|
Increase from baseline
|
28.74 (5.20, 540.18)
|
0.005
|
24.71 (4.19, 479.13)
|
0.007
|
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)
|
Predictor
|
OR(95%CI)
|
p-value
|
N
|
Event
|
VIF
|
|
Did ctDNA increase or decrease from
baseline to cycle 3
|
|
|
73
|
58
|
1.01
|
|
Decrease from baseline
|
Reference
|
|
33
|
19
|
|
|
Increase from baseline
|
24.71 (4.19, 479.13)
|
0.004
|
40
|
39
|
|
|
PD L1 percent
|
0.98 (0.95, 1.00)
|
0.024
|
73
|
58
|
1.01
|
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')
|
|
Estimate(95%CI)
|
p-value
|
N
|
|
Cohort A
|
|
|
|
|
Age at study entry
|
2.94 (-0.70, 6.58)
|
0.10
|
15
|
|
Patient Sex
|
|
|
15
|
|
Female
|
Reference
|
|
3
|
|
Male
|
-40.25 (-96.25, 15.75)
|
0.14
|
12
|
|
Cohort E
|
|
|
|
|
Age at study entry
|
-0.44 (-1.02, 0.15)
|
0.14
|
30
|
|
Patient Sex
|
|
|
30
|
|
Female
|
Reference
|
|
12
|
|
Male
|
-14.86 (-32.57, 2.85)
|
0.097
|
18
|
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')
|
Group
|
Events/Total
|
Median (95%CI)
|
12months (95% CI)
|
24months (95% CI)
|
36months (95% CI)
|
|
Female
|
39/58
|
14.29 (9.69, 23.82)
|
0.55 (0.44, 0.69)
|
0.34 (0.24, 0.50)
|
0.29 (0.18, 0.45)
|
|
Male
|
25/36
|
11.24 (6.14, 25.33)
|
0.50 (0.36, 0.69)
|
0.31 (0.18, 0.52)
|
0.27 (0.15, 0.48)
|
|
|
|
|
Log Rank Test
|
ChiSq
|
0.5 on 1 df
|
|
|
|
|
|
p-value
|
0.46
|
Stratified Survival Summary
To combine estimates across strata
rm_survdiff(data=pembrolizumab,time='os_time',status='os_status',
covs='sex',strata='cohort',digits=1)
|
group
|
N
|
Observed
|
Expected
|
Median (95%CI)
|
|
Overall
|
94
|
64
|
|
14.0 (9.0, 20.1)
|
|
Female
|
58
|
39
|
43.0
|
14.3 (9.7, 23.8)
|
|
Male
|
36
|
25
|
21.0
|
11.2 (6.1, 25.3)
|
|
Log Rank Test
|
|
|
|
ChiSq = 1.9 on 1 df
|
|
stratified by cohort
|
|
|
|
p-value = 0.17
|
Working with Labels
New in version 0.1.0 is the ability to incorporate variable labels
into summary tables (but not yet all plots). If variables contain a
label attribute this will be displayed automatically, to
disable this set nicenames=F
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 tableset_var_labels allows you to set individual variable
labels to a data frameclear_labels removes all labels from a data frameexport_labels extracts variable labels from a data
frame and returns a data frame of variable names and variable
labels
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'))
|
|
n=270
|
|
Study Cohort
|
|
|
A
|
50 (19)
|
|
B
|
14 (5)
|
|
C
|
18 (7)
|
|
D
|
88 (33)
|
|
E
|
100 (37)
|
|
Change in ctDNA since baseline
|
|
|
Clearance
|
137 (51)
|
|
No clearance, decrease from baseline
|
44 (16)
|
|
No clearance, increase from baseline
|
89 (33)
|
|
Percentage change in tumour
measurement
|
|
|
Mean (sd)
|
-29.7 (52.8)
|
|
Median (Min,Max)
|
-32.5 (-100.0, 197.1)
|
|
Missing
|
8
|
set_labels
If we have a lookup table of variable names and labels that we
imported from a data dictionary we can set the variable labels for the
data frame and these will be used in the rm_ functions
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'))
|
|
n=270
|
|
Study Cohort
|
|
|
A
|
50 (19)
|
|
B
|
14 (5)
|
|
C
|
18 (7)
|
|
D
|
88 (33)
|
|
E
|
100 (37)
|
|
Change in ctDNA since baseline
|
|
|
Clearance
|
137 (51)
|
|
No clearance, decrease from baseline
|
44 (16)
|
|
No clearance, increase from baseline
|
89 (33)
|
|
Percentage change in tumour
measurement
|
|
|
Mean (sd)
|
-29.7 (52.8)
|
|
Median (Min,Max)
|
-32.5 (-100.0, 197.1)
|
|
Missing
|
8
|
set_var_labels
Individual labels can be changed with with the
set_var_labels command
ctDNA <- set_var_labels(ctDNA,
cohort="A new cohort label")
rm_covsum(data=ctDNA,
covs=c('cohort','ctdna_status','size_change'))
|
|
n=270
|
|
A new cohort label
|
|
|
A
|
50 (19)
|
|
B
|
14 (5)
|
|
C
|
18 (7)
|
|
D
|
88 (33)
|
|
E
|
100 (37)
|
|
Change in ctDNA since baseline
|
|
|
Clearance
|
137 (51)
|
|
No clearance, decrease from baseline
|
44 (16)
|
|
No clearance, increase from baseline
|
89 (33)
|
|
Percentage change in tumour
measurement
|
|
|
Mean (sd)
|
-29.7 (52.8)
|
|
Median (Min,Max)
|
-32.5 (-100.0, 197.1)
|
|
Missing
|
8
|
replace_plot_labels
This function will accept a ggplot plot and replace the variable
names (x-axis, y-axis and legend) with the variable labels. This is
useful for more professional looking plots.
library(ggplot2)
p <- ggplot(data=ctDNA,aes(x=ctdna_status,y=size_change,colour=cohort))+
geom_point()
replace_plot_labels(p)
clear_labels
Or clear them all
ctDNA <- clear_labels(ctDNA)
Plotting Functions
{width=100%}
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')
###
Multivariable Model Forest Plot
glm_fit <- glm(orr~change_ctdna_group+pdl1+age,
family='binomial',
data = pembrolizumab)
forestplotMV(glm_fit)
Combining univariable and multivariable models into a single
plot
UVp = forestplotUV(response="orr", covs=c("change_ctdna_group", "sex", "age",
"l_size"), data=pembrolizumab, family='binomial')
MVp = forestplotMV(glm(orr~change_ctdna_group+sex+age+l_size,
data=pembrolizumab,family = 'binomial'))
forestplotUVMV(UVp, MVp)
This can also be done with linear scale odds ratios. 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=F)
Miscellaneous Functions
excelCol
To identify the number of a column given the Excel column header
## 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"
Options
The following options can be set:
- reportRmd.digits, default 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)
|
|
Estimate(95%CI)
|
p-value
|
N
|
|
Age at study entry
|
0.82 (-10.13, 11.76)
|
0.88
|
94
|
|
Patient Sex
|
|
|
94
|
|
Female
|
Reference
|
|
58
|
|
Male
|
56.61 (-228.71, 341.93)
|
0.69
|
36
|
|
Target lesion size at baseline
|
1.21 (-1.12, 3.54)
|
0.31
|
94
|
|
PD L1 percent
|
-3.50 (-8.27, 1.27)
|
0.15
|
93
|
|
log of TMB
|
18.78 (-125.18, 162.74)
|
0.80
|
94
|
options('reportRmd.digits'=1)
rm_uvsum(response = 'baseline_ctdna',
covs=c('age','sex','l_size','pdl1','tmb'),
data=pembrolizumab)
|
|
Estimate(95%CI)
|
p-value
|
N
|
|
Age at study entry
|
0.8 (-10.1, 11.8)
|
0.88
|
94
|
|
Patient Sex
|
|
|
94
|
|
Female
|
Reference
|
|
58
|
|
Male
|
56.6 (-228.7, 341.9)
|
0.69
|
36
|
|
Target lesion size at baseline
|
1.2 (-1.1, 3.5)
|
0.31
|
94
|
|
PD L1 percent
|
-3.5 (-8.3, 1.3)
|
0.15
|
93
|
|
log of TMB
|
18.8 (-125.2, 162.7)
|
0.80
|
94
|
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
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