### Relation to traditional cost-effectiveness plane

A traditional PSA output is a two-dimensional black and white scatter plot presented on a CE-plane (Fig. 1a). The PSA-ReD plot (Fig. 1b) combines a multi-colored density plot (Additional file 1: Fig. S1a) and a contour plot (Additional file 1: Fig. S1b). The combination of these two plots allows the reader to identify and distinguish high density areas using a color scale, as well as a quantification of the point estimate density within the CE-plane, thus visualizing the information that remains hidden in the traditional scatter plot. This increases the information that can be understood from the scatter plot and improves understanding of the parameter uncertainty which a PSA is aimed to address.

To explain how the PSA-ReD plot works, we provide a simulated example in Fig. 1b. It assumes a model with only two standard normally distributed parameters (mean = 0 and standard deviation = 1) that define the incremental costs and incremental QALYs.

In Fig. 1a, the base case would be at zero incremental costs and zero incremental QALYs. Though we can see that the borders of the area are less densely populated, it is unclear how the density of iterations is spread over the populated area. If instead we look at the PSA-ReD plot in Fig. 1b, it becomes clear that the density is evenly spread around the base case, as would be expected for this normally distributed data. Additionally, the contours give insight into the spread of the iterations. In this case, the area containing 95% of the iterations will approximate that of a 95% confidence interval because we used normal distributions. A bivariate normal distribution is distributed according to the χ^{2}-distribution with two degrees of freedom [14]. Taking the square root of the critical value for the 95% confidence interval (5.99), results in the area borders of the confidence interval (2.45). This is clearly shown by the contours in the PSA-ReD plot. In general, the probability that the values for two variables within a joint distribution together fall in any area of their two dimensions is given by the volume (or cumulative probability) under the density function above that area. This is exactly what the PSA-ReD method calculates, as is explained in the next section.

### The density plot

The density plot (Additional file 1: Fig. S1a) could be interpreted as a two-dimensional histogram. Like a traditional one-dimensional histogram, the two axes are divided in sections. These sections on both axes divide the two-dimensional space in distinct rectangular regions. Then, as in a one-dimensional histogram, the number of data points per region is counted and transformed to present the relative frequency using a color scale. Low density is presented by a green to blue scale and high density is presented by a yellow to red scale.

When using a histogram, the choice of the anchor point of the plot area (i.e., the range and starting points of the axes) has influence on the graphical outcome [15]. This effect would be most pronounced in a one-dimensional histogram with relatively few datapoints: Shifting the x-axis would change the histogram as, by chance, the number of datapoints falling within each bin would differ with each x-axis shift. Crucially, the underlying data remains the same and this bias would also be present in two-dimensional histograms [15]. To overcome this, bivariate kernel density estimation (kde) is used as it provides a more accurate representation of the probability density [15]. Instead of counting the number of data points per rectangular section, each data point is surrounded by a kernel which are summed to yield the kernel density estimate. Each data point is thereby smoothed over a small surrounding area (data kernel) instead of being a single data point [15]. The size of this area is determined by the data as explained in the ‘technical aspects’ paragraph. Consequently, the way the resulting plot looks does not depend anymore on the size of any bins over the x-axis or y-axis.

### The contour plot

The PSA-ReD plot combines the density plot with a contour plot (Additional file 1: Fig. S1b). The contours indicate the boundaries of regions with a specific cumulative probability. This cumulative probability is calculated by summing the area density estimates (retrieved from the kernel estimation method). For this summing of density estimates to cumulative probabilty, densities are sorted from high density to low density with the summing starting from the highest density. These cumulative probabilities are then mapped to a range of 0 to 1. In this way, the total density in the plot area sums to one, reflecting the cumulative probability. A contour line is then drawn joining areas with specific pre-specified values of cumulative probability. For each individual plot area in the PSA-ReD, two values are now available: the density per area and the cumulative probability that is reached per area. A contour line is then drawn joining areas with specific pre-specified values of cumulative probability. These values can be chosen by the user (e.g. 0.1, 0.5, 0.95).

When there are multiple, disjointed high (or low) density-areas, it is possible that separate contours with equal cumulative probability values are drawn over these separate areas. For example, when there are separate high-density areas that together amount to 50% of the data points, seperate contour-lines with the value of 0.5 would be drawn around both high density values.

### Hardware and software

The script for the PSA-ReD plot was developed and tested using R version 3.5.1 and Rstudio version 1.1.453 (R Core team [11, 16]). For our analyses, we used a standard consumer grade personal computer (Dell OptiPlex 9020). In the technical appendix (Additional file 2), we provide detailed information on the hardware and software used.

The R script that we used is available in a GitHub repository [17]. We adhered to Google’s R Style guide and provide step-by-step guidance using comments embedded in the script [18]. The R script is licensed under the GNU General Public License v3.0 [19]. In short, this means that users are free to run, study, share and modify the software. The license dictates, among other things, that the software (or derivative work) must be open source and that derivative work must be published using the same license [19]. This guarantees that our project can be used and optimized by anyone whilst ensuring that it remains open to all.

### Technical aspects of plot generation in R

In R, we use the kde2d function from the Modern Applied Statistics with S (MASS) package to perform the aforementioned kernel density estimation [20]. In essence, the outcome of kde is a density value per area of a prespecified size, comparable to the number of data points within each bin in histograms. Detailed information is provided in the work by Silverman and in the documentation of the MASS package [20, 21]. As these density values are very small and hard to interpret, we normalize these values by taking the reciprocal of the maximum density value to yield values ranging from 0.0 to 1.0. With these, we generate an easy to interpret plot with a scale from 1.0 (highest density) to 0.0 (lowest density).

The kde2d function has, besides the x and y values, two arguments that influence the kde. These are n (the number of bins in each dimension) and h (the bandwidth that determines the level of smoothing). The number of bins defines the number of sections on each axis. The total number of areas within the resulting plots is therefore horizontal bins * vertical bins (e.g. 100*100 = 10,000). An easy analogy of these areas would be to regard them as pixels, the bins then determine the resolution in both directions. This pixel-analogy only reflects to the number of underlying bins. As we outline in the Additional file 2 (Technical Appendix) regarding the saving of plots, the actual resolution of the figures can be specified and is irrespective of the number of bins used. As the number of bins can be interpreted as the resolution of the figure, a larger number of bins produces a more precise figure. However, increasing the number of bins also increases computation time which means a balance must be struck.

In Additional file 1: Fig. S2, we present the influence of different bin sizes. Using 50-500 bins (Additional file 1: Fig. S2a, b), yields a density gradient that is not smooth and may appear like the image is pixelated. With 1000 bins (Additional file 1: Fig. S2c), the image is smooth, no pixilation can be identified and all the computation is performed within 1 minute on the aforementioned consumer grade computer. With more bins (2000, Additional file 1: Fig. S2d), the image does not get smoother but it does lead to increased RAM usage and computation time. Another consequence of a lower number of bins (as in Additional file 1: Fig. S2a) is that due to the lack of smoothing, the contours are placed different from when more bins are used. Essentially, with fewer bins, the contours are very rough. We therefore recommend using 1000 bins and have used this number of bins in all figures throughout the manuscript, unless otherwise stated.

The h argument of the kde2d function determines the bandwidth of the kernel areas. It can be interpreted as the size of the kernels that is applied when converting each data point to a data kernel. We have chosen to leave this at the default setting where the bandwidth is automatically selected based on the data by the well-established MASS package (specifically, the bandwidth.nrd function) [20]. This guarantees generalizability of results.

### Number of PSA iterations

As in any PSA, it is preferred to run as many iterations as necessary to reach model convergence [6]. We explored the influence of the number of iterations used by varying this between 1000 and 100,000 iterations, as presented in Additional file 1: Fig. S3. As RAM usage and computation time increases when more iterations are used, we recommend using a maximum of 10,000 iterations. Running the script with 10,000 iterations takes a maximum of 1 min. In all figures throughout this manuscript, we have used 10,000 iterations unless otherwise stated.

### User modifications

Other parameters that can be altered by the user are contour levels, axis-, legend- and plot titles, font sizes and font types. In the supplied script, it is explained how and where this can be done. Apart from these cosmetic changes, we provide means to zoom on a particular area of the plot and generate a new plot from that specific area. Additional file 1: Fig. S4 displays this zooming capability. We also provide a feature that allows users to plot willingness-to-pay (WTP) thresholds in the PSA-ReD plot, as well as plotting the average incremental costs and incremental effects for the PSA and the results of the deterministic base case scenario. The technical appendix (Additional file 2) provides in-depth explanations on the use of the various features described above.

### Case study demonstration

To demonstrate the novel graphical presentation, the concept was applied to two exemplary case studies. These real-world case studies were selected as a convenience sample as we needed access to the raw PSA results. To increase transparency, we opted for published case studies. The two selected cases each show a different pattern within the PSA results. Both patterns are commonly seen in economic evaluations. The first real-world case study assesses the influence of three characteristics (cost, specificity and sensitivity) on cost-effectiveness of a hypothetical pharmacogenomic test for prevention of angiotensin-converting enzyme inhibitor induced angioedema (denoted as ‘eHTA study’) [22]. The second real-world case study used a three-state partitioned survival model to investigate cost-effectiveness of periodic therapeutic drug monitoring of endoxifen levels in breast cancer patients (denoted as ‘TDM study’) [23].