Visualising results from rsimsum
Alessandro Gasparini
2024-03-03
This vignette requires the following R packages:
library(rsimsum)
library(ggplot2)
Data
We use data from a simulation study on misspecification of the
baseline hazard in survival models. This dataset is included in
rsimsum and can be loaded with:
data("relhaz", package = "rsimsum")
Inspecting the structure of the dataset and the first 15 rows of
data:
str(relhaz)
#> 'data.frame': 1200 obs. of 6 variables:
#> $ dataset : int 1 2 3 4 5 6 7 8 9 10 ...
#> $ n : num 50 50 50 50 50 50 50 50 50 50 ...
#> $ baseline: chr "Exponential" "Exponential" "Exponential" "Exponential" ...
#> $ theta : num -0.88 -0.815 -0.143 -0.333 -0.483 ...
#> $ se : num 0.333 0.325 0.305 0.314 0.306 ...
#> $ model : chr "Cox" "Cox" "Cox" "Cox" ...
head(relhaz, n = 15)
#> dataset n baseline theta se model
#> 1 1 50 Exponential -0.88006151 0.3330172 Cox
#> 2 2 50 Exponential -0.81460242 0.3253010 Cox
#> 3 3 50 Exponential -0.14262887 0.3050516 Cox
#> 4 4 50 Exponential -0.33251820 0.3144033 Cox
#> 5 5 50 Exponential -0.48269940 0.3064726 Cox
#> 6 6 50 Exponential -0.03160756 0.3097203 Cox
#> 7 7 50 Exponential -0.23578090 0.3121350 Cox
#> 8 8 50 Exponential -0.05046332 0.3136058 Cox
#> 9 9 50 Exponential -0.22378715 0.3066037 Cox
#> 10 10 50 Exponential -0.45326446 0.3330173 Cox
#> 11 11 50 Exponential -0.71402510 0.3251902 Cox
#> 12 12 50 Exponential -0.32956944 0.3073481 Cox
#> 13 13 50 Exponential -0.15351788 0.3056453 Cox
#> 14 14 50 Exponential -0.82742207 0.3283561 Cox
#> 15 15 50 Exponential -0.14594648 0.3255636 Cox
Summarise results
We use the simsum function to summarise results:
s1 <- simsum(
data = relhaz, estvarname = "theta", se = "se", true = -0.50,
methodvar = "model", by = c("n", "baseline"), x = TRUE
)
#> 'ref' method was not specified, Cox set as the reference
s1
#> Summary of a simulation study with a single estimand.
#> True value of the estimand: -0.5
#>
#> Method variable: model
#> Unique methods: Cox, Exp, RP(2)
#> Reference method: Cox
#>
#> By factors: n, baseline
#>
#> Monte Carlo standard errors were computed.
We call simsum with x = TRUE as that is
required for some types of plots (e.g. zip plots, scatter plots,
etc.).
summary(s1)
#> Values are:
#> Point Estimate (Monte Carlo Standard Error)
#>
#> Non-missing point estimates/standard errors:
#> n baseline Cox Exp RP(2)
#> 50 Exponential 100 100 100
#> 50 Weibull 100 100 100
#> 250 Exponential 100 100 100
#> 250 Weibull 100 100 100
#>
#> Average point estimate:
#> n baseline Cox Exp RP(2)
#> 50 Exponential -0.4785 -0.4761 -0.4817
#> 50 Weibull -0.5282 -0.3491 -0.5348
#> 250 Exponential -0.5215 -0.5214 -0.5227
#> 250 Weibull -0.5120 -0.3518 -0.5139
#>
#> Median point estimate:
#> n baseline Cox Exp RP(2)
#> 50 Exponential -0.4507 -0.4571 -0.4574
#> 50 Weibull -0.5518 -0.3615 -0.5425
#> 250 Exponential -0.5184 -0.5165 -0.5209
#> 250 Weibull -0.5145 -0.3633 -0.5078
#>
#> Average variance:
#> n baseline Cox Exp RP(2)
#> 50 Exponential 0.1014 0.0978 0.1002
#> 50 Weibull 0.0931 0.0834 0.0898
#> 250 Exponential 0.0195 0.0191 0.0194
#> 250 Weibull 0.0174 0.0164 0.0172
#>
#> Median variance:
#> n baseline Cox Exp RP(2)
#> 50 Exponential 0.1000 0.0972 0.0989
#> 50 Weibull 0.0914 0.0825 0.0875
#> 250 Exponential 0.0195 0.0190 0.0194
#> 250 Weibull 0.0174 0.0164 0.0171
#>
#> Bias in point estimate:
#> n baseline Cox Exp RP(2)
#> 50 Exponential 0.0215 (0.0328) 0.0239 (0.0326) 0.0183 (0.0331)
#> 50 Weibull -0.0282 (0.0311) 0.1509 (0.0204) -0.0348 (0.0311)
#> 250 Exponential -0.0215 (0.0149) -0.0214 (0.0151) -0.0227 (0.0149)
#> 250 Weibull -0.0120 (0.0133) 0.1482 (0.0093) -0.0139 (0.0137)
#>
#> Relative bias in point estimate:
#> n baseline Cox Exp RP(2)
#> 50 Exponential -0.0430 (0.0657) -0.0478 (0.0652) -0.0366 (0.0662)
#> 50 Weibull 0.0564 (0.0623) -0.3018 (0.0408) 0.0695 (0.0622)
#> 250 Exponential 0.0430 (0.0298) 0.0427 (0.0301) 0.0455 (0.0298)
#> 250 Weibull 0.0241 (0.0267) -0.2963 (0.0186) 0.0279 (0.0274)
#>
#> Empirical standard error:
#> n baseline Cox Exp RP(2)
#> 50 Exponential 0.3285 (0.0233) 0.3258 (0.0232) 0.3312 (0.0235)
#> 50 Weibull 0.3115 (0.0221) 0.2041 (0.0145) 0.3111 (0.0221)
#> 250 Exponential 0.1488 (0.0106) 0.1506 (0.0107) 0.1489 (0.0106)
#> 250 Weibull 0.1333 (0.0095) 0.0929 (0.0066) 0.1368 (0.0097)
#>
#> % gain in precision relative to method Cox:
#> n baseline Cox Exp RP(2)
#> 50 Exponential -0.0000 (0.0000) 1.6773 (3.2902) -1.6228 (1.7887)
#> 50 Weibull -0.0000 (0.0000) 132.7958 (16.4433) 0.2412 (3.7361)
#> 250 Exponential 0.0000 (0.0000) -2.3839 (3.0501) -0.1491 (0.9916)
#> 250 Weibull -0.0000 (0.0000) 105.8426 (12.4932) -4.9519 (2.0647)
#>
#> Mean squared error:
#> n baseline Cox Exp RP(2)
#> 50 Exponential 0.1073 (0.0149) 0.1056 (0.0146) 0.1089 (0.0154)
#> 50 Weibull 0.0968 (0.0117) 0.0640 (0.0083) 0.0970 (0.0117)
#> 250 Exponential 0.0224 (0.0028) 0.0229 (0.0028) 0.0225 (0.0028)
#> 250 Weibull 0.0177 (0.0027) 0.0305 (0.0033) 0.0187 (0.0028)
#>
#> Model-based standard error:
#> n baseline Cox Exp RP(2)
#> 50 Exponential 0.3185 (0.0013) 0.3127 (0.0010) 0.3165 (0.0012)
#> 50 Weibull 0.3052 (0.0014) 0.2888 (0.0005) 0.2996 (0.0012)
#> 250 Exponential 0.1396 (0.0002) 0.1381 (0.0002) 0.1394 (0.0002)
#> 250 Weibull 0.1320 (0.0002) 0.1281 (0.0001) 0.1313 (0.0002)
#>
#> Relative % error in standard error:
#> n baseline Cox Exp RP(2)
#> 50 Exponential -3.0493 (6.9011) -4.0156 (6.8286) -4.4305 (6.8013)
#> 50 Weibull -2.0115 (6.9776) 41.4993 (10.0594) -3.6873 (6.8549)
#> 250 Exponential -6.2002 (6.6679) -8.3339 (6.5160) -6.4133 (6.6528)
#> 250 Weibull -0.9728 (7.0397) 37.7762 (9.7917) -4.0191 (6.8228)
#>
#> Coverage of nominal 95% confidence interval:
#> n baseline Cox Exp RP(2)
#> 50 Exponential 0.9500 (0.0218) 0.9400 (0.0237) 0.9500 (0.0218)
#> 50 Weibull 0.9700 (0.0171) 0.9900 (0.0099) 0.9500 (0.0218)
#> 250 Exponential 0.9300 (0.0255) 0.9200 (0.0271) 0.9300 (0.0255)
#> 250 Weibull 0.9400 (0.0237) 0.8500 (0.0357) 0.9400 (0.0237)
#>
#> Bias-eliminated coverage of nominal 95% confidence interval:
#> n baseline Cox Exp RP(2)
#> 50 Exponential 0.9500 (0.0218) 0.9500 (0.0218) 0.9500 (0.0218)
#> 50 Weibull 0.9500 (0.0218) 1.0000 (0.0000) 0.9500 (0.0218)
#> 250 Exponential 0.9400 (0.0237) 0.9400 (0.0237) 0.9400 (0.0237)
#> 250 Weibull 0.9500 (0.0218) 0.9900 (0.0099) 0.9400 (0.0237)
#>
#> Power of 5% level test:
#> n baseline Cox Exp RP(2)
#> 50 Exponential 0.3600 (0.0480) 0.3800 (0.0485) 0.3700 (0.0483)
#> 50 Weibull 0.4300 (0.0495) 0.0900 (0.0286) 0.4700 (0.0499)
#> 250 Exponential 0.9800 (0.0140) 0.9900 (0.0099) 0.9900 (0.0099)
#> 250 Weibull 0.9700 (0.0171) 0.8600 (0.0347) 0.9700 (0.0171)
rsimsum implements the autoplot method for
objects of classes: simsum, summary.simsum,
multisimsum, summary.multisimsum.
See ?ggplot2::autoplot() for details on the S3 generic
function.
Scatter plots
Scatter plots allow to assess serial trends in estimates and standard
errors.
For instance, if we want to compare the point estimates from
different methods (across data-generating mechanisms):
autoplot(s1, type = "est")
#> `geom_smooth()` using formula = 'y ~ x'
Analogously, if we want to compare standard errors:
autoplot(s1, type = "se")
#> `geom_smooth()` using formula = 'y ~ x'
These two plot types allow comparing estimates (and standard errors)
obtained from different methods with ease; in ideal settings, the points
of the scatterplots should lay on the diagonal (dashed line). An
estimated regression line (of X vs Y, blue
line) is superimposed by default to ease the comparison even more.
In addition to plots comparing estimates and standard errors,
Bland-Altman-type plots are supported as well:
autoplot(s1, type = "est_ba")
#> `geom_smooth()` using formula = 'y ~ x'
Bland-Altman plots compare the difference between estimates from two
competing methods (on the y-axis) with the mean of the estimates from
the two methods (x-axis). In the ideal scenario, there should be no
trend and all the points from the scatter plot should lay around the
horizontal dashed line. To ease comparison, a regression line is
included here as well. Bland-Altman plots for standard errors could be
obtained by setting the argument type = "se_ba".
Ridgeline plots
Another way of visually comparing estimates from different methods is
given by ridgeline plots. According to the documentation of the
ggridges package (source),
Ridgeline plots are partially overlapping line plots that create the
impression of a mountain range. They can be quite useful for visualizing
changes in distributions over time or space.
In the settings of simulation studies, we aim to visualise changes in
distribution over data-generating mechanisms.
For instance, say we want to compare estimates across data-generating
mechanisms and methods:
autoplot(s1, type = "est_ridge")
#> Picking joint bandwidth of 0.077
This allows us to see how the estimates from the Exp
method are different from the other methods in two of the four
data-generating mechanisms: 50, Weibull and
250, Weibull.
To obtain a similar plot for standard errors, call the
autoplot method with type = "se_ridge
instead.
Lolly plots
Lolly plots are used to present estimates for a given summary
statistic with confidence intervals based on Monte Carlo standard errors
(if calling the autoplot method on summary
objects). They allow to easily compare methods.
Say we are interested in bias:
autoplot(summary(s1), type = "lolly", stats = "bias")
It is straightforward to identify the exponential model as yielding
biased results when the true baseline hazard is Weibull, irrespectively
of the sample size.
On a relative scale, we can plot relative bias as well:
autoplot(summary(s1), type = "lolly", stats = "rbias")
If confidence intervals based on Monte Carlo errors are not required,
it is sufficient to call the autoplot method on the
simsum object:
autoplot(s1, type = "lolly", stats = "bias")
Analogously, for coverage:
autoplot(summary(s1), type = "lolly", stats = "cover")
Forest plots
Forest plots could be an alternative to lolly plots, with similar
interpretation:
autoplot(s1, type = "forest", stats = "bias")
autoplot(summary(s1), type = "forest", stats = "bias")
Zipper plots
Zipper plots (or zip plots), introduced in Morris et al.
(2019), help understand coverage by visualising the confidence intervals
directly. For each data-generating mechanism and method, the confidence
intervals are centile-ranked according to their significance against the
null hypothesis \(H_0: \theta\) = \(\theta_{\text{true}}\), assessed via a
Wald-type test. This ranking is used for the vertical axis and is
plotted against the intervals themselves.
When a method has 95% coverage, the colour of the intervals switches
at 95 on the vertical axis. Finally, the horizontal lines represent
confidence intervals for the estimated coverage based on Monte Carlo
standard errors.
autoplot(s1, type = "zip")
The zipper plot for the exponential model with n = 50 and a true
Weibull baseline hazard shows that coverage is approximately 95%;
however, there are more intervals to the right of \(\theta\) = -0.50 than to the left: this
indicates that the model standard errors must be overestimating the
empirical standard error, because coverage is appropriate despite
bias.
It is also possible to zoom on the top x% of the zip plot to
increase readability, e.g. on the top 30%:
autoplot(s1, type = "zip", zoom = 0.3)
Heat plots
Heat plots are a new visualisation that we suggest and include here
for the first time. With heat plots, we produce a heat-map-like plot
where the filling of each tile represents a given summary statistic, say
bias:
autoplot(s1, type = "heat", stats = "bias")
This visualisation type automatically puts all data-generating
mechanisms on the y-axis; by default, the methods are included on the
x axis. Therefore, this plot is most useful when a
simulation study includes different methods to be compared and many
data-generating mechanisms. Using a heat plot, it is immediate to
identify visually which method performs better and under which
data-generating mechanisms.
By default, heat plots use the default ggplot scale for
the filling aesthetic. It is recommended to use a different colour
palette with better characteristics, e.g. the viridis colour palette
from matplotlib; see the next section for details on how to
do this, and here
for details on the viridis colour palette.
Contour plots and hexbin plots
Individual point estimates and standard errors could also be plotted
using contour plots or hexbin plots.
Contour plots represent a 3-dimensional surface by plotting constant
z slices (called contours) on a 2-dimensional format. That is, given a
value for z, lines are drawn for connecting the (x, y) coordinates where
the value of z is (relatively) homogenous.
Hexbin plots are useful to represent the relationship of 2 numerical
variables when you have a lot of data points: instead of overlapping,
the plotting window is split in several hexbins, and the number of
points per hexbin is counted. The colour filling denotes then the number
of points.
Both plots provide an alternative to scatter plots when there is a
large number of data points that overlap. Contour plots and hexbin plots
can be easily obtained using the autoplot method once
again, and using the argument type = "est_density",
type = "se_density", type = "est_hex", or
type = "se_hex". For instance, focussing on point
estimates:
autoplot(s1, type = "est_density")
#> `geom_smooth()` using formula = 'y ~ x'
autoplot(s1, type = "est_hex")
#> `geom_smooth()` using formula = 'y ~ x'
Of course, analogous plots could be obtained for standard errors.
Custom plotting
All plots produced by rsimsum are meant to be quick
explorations of results from Monte Carlo simulation studies: they are
not meant to be final manuscript-like-quality plots (although they can
be useful as a starting point).
Generally, the output of all types of autoplot calls are
ggplot objects; hence, it should be generally
straightforward to customise all plots. For instance, say we want to add
a custom theme:
autoplot(summary(s1), type = "lolly", stats = "bias") +
ggplot2::theme_bw()
Compare to the default plot:
autoplot(summary(s1), type = "lolly", stats = "bias")
We also mentioned before that the colour palette of heat plots should
be customised. Say we want to use the default viridis colour
palette:
autoplot(s1, type = "heat", stats = "bias") +
ggplot2::scale_fill_viridis_c()
Analogously, say we want to customise the colour palette of ridgeline
plots, again using the viridis colour palette:
autoplot(s1, type = "est_ridge") +
scale_fill_viridis_d() +
scale_colour_viridis_d()
#> Picking joint bandwidth of 0.077
