Plot customization for Noncompartmental analysis - Single Ascending Dose

1 Libraries

We will load the NCA, NCAUtilities, and PumasUtilities packages in order to have access to Pumas’ core and utilities functions for Noncompartmental analysis:

using NCA
using NCAUtilities
using PumasUtilities

We will also use Makie.jl as the main plotting library and also AlgebraOfGraphics.jl.

CairoMakie.jl will act as our plotting backend and will also make possible the various customization options that we will go over. On the other hand, AlgebraOfGraphics.jl will provide a syntax inspired by the grammar of graphics and similar to R’s ggplot2 to create complex plots.

using CairoMakie
using AlgebraOfGraphics

Lastly, we will use the CSV.jl package to read external data and DataFramesMeta.jl for data wrangling.

using CSV
using DataFrames
using DataFramesMeta

2 Introduction

This tutorial will show you how to customize some of the plots available in the Pumas NCA workflow and create unique plots and layouts using Pumas and Makie.jl.

Note

We will be focusing on plotting and customization. Please refer to the tutorial Introduction to Noncompartmental Analysis (NCA) if anything feels unfamiliar or if you wish to explore the NCA capabilities that Pumas has to offer.

3 Motivation

Data visualization plays a central role in communicating NCA results. However, effective communication depends strongly on having control over how your data visualization tools will look like. By customizing the look and feel of your plots, you will be able to play emphasis on what you consider to be most the important results, create unique plots and layouts for your presentations, and follow the appropriate format for your tables, listings, and figures (TLF).

Luckily, Pumas and Makie.jl support a wide variety of plotting and customization options in order to fulfill your data visualization needs.

3.1 🎨 Color palette

In order to showcase the extensive support for customization in Pumas, we will be trying to apply a consistent color palette across the plots generated in this tutorial. As an example, consider the following two colors to be used as the main ones in the palette.

We will also use these three additional colors for highlighting and more complex plots.

These are the names and hex codes for the colors in the palette:

• Navy Blue: #0B1E3F
• Emerald Green: #50C878
• Coral: #FF7F50
• Mustard yellow: #FFDB58
• Deep purple: #6A5ACD

Note

Makie uses Colors.jl to define colors. This means that you can pass a string representing any [CSS color] (https://developer.mozilla.org/en-US/docs/Web/CSS/color) to Pumas and Makie’s plotting functions.

For instance, you could define the color red with "rgb(255,0,0)", "#FF0000", or "red". You can also use one of the named colors from Colors.jl by passing them as a Symbol (e.g :navyblue)

We will store the colors in a Vector so that we can access them later more easily. Let’s begin by defining the main and complementary colors. Notice that we are creating strings with the hex codes from above:

main_colors = ["#0B1E3F", "#50C878"];
complementary = ["#FF7F50", "#FFDB58", "#6A5ACD"];

We can join both lists of colors using vcat, which will concatenate main_colors and complementary_colors to create a single color palette vector.

palette = vcat(complementary, main_colors);

We would also like to use a shade of gray for the background color instead of the default white color

background = "#FAF9F6";

4 Load Data

We will be using simulated data for an orally administered drug through a single dose ranging from 100 to 1,600 mg.

Note

This tutorial is based on Alison Margolskee and Andy Stein’s PK - Single Ascending Dose workflow. You can get the dataset used in this tutorial here.

Read the data using CSV.jl and convert it to a DataFrame

df = CSV.read("Single_Ascending_Dose_Dataset2.csv", DataFrame; missingstring = ["NA"])
first(df, 5)

5×15 DataFrame

Row	ID	TIME	NOMTIME	TIMEUNIT	AMT	LIDV	CMT	NAME	EVENTU	CENS	EVID	WEIGHTB	SEX	TRTACT	DOSE
	Int64	Float64	Float64	String7	Int64	Float64?	Int64	String31	String7	Int64	Int64	Float64	String7	String7	Int64
1	1	-0.089	-0.1	Hours	0	missing	2	PK Concentration	ng/mL	0	0	77.9	Male	100 mg	100
2	1	0.0	0.0	Hours	100	missing	1	Dosing	mg	0	1	77.9	Male	100 mg	100
3	1	0.17	0.1	Hours	0	0.665	2	PK Concentration	ng/mL	0	0	77.9	Male	100 mg	100
4	1	0.565	0.5	Hours	0	0.997	2	PK Concentration	ng/mL	0	0	77.9	Male	100 mg	100
5	1	1.145	1.0	Hours	0	1.35	2	PK Concentration	ng/mL	0	0	77.9	Male	100 mg	100

4.1 Data wrangling

Tip

We have a series of tutorials on Data Wrangling in Julia, where you can learn more about reading and manipulating data.

Before defining a population with read_nca function, we need to do some transformations on the dataset. First, we will switch the time units to days and create a column for the route of administration as extravascular ("ev").

@rtransform! df begin
    :NOMTIME = :NOMTIME / 24
    :ROUTE = "ev"
end;

We also have to get rid of the rows where the time is negative:

@rsubset! df begin
    :NOMTIME >= 0
end;

Finally, we have to define the units for the time (days), concentration (ng/mL), and dose amount (mg)

timeu = u"d";
concu = u"ng/mL";
amtu = u"mg";

Note

Notice that we have defined units using the u string literal (as in u"..."). This syntax corresponds to Unitful.jl, which is the Julia package from which Pumas obtains its unit support.

Now we are ready to run read_nca. Notice that we are setting the keyword argument group as [:DOSE] because we want to be able to account for different dosing levels in our analysis.

nca_data = read_nca(
    df;
    id = :ID,
    time = :NOMTIME,
    observations = :LIDV,
    amt = :AMT,
    route = :ROUTE,
    group = [:DOSE],
    timeu = timeu,
    concu = concu,
    amtu = amtu,
);

5 Exploring concentration time profiles

In order to get the concentration time profile for each subject in the population, we use the function observations_vs_time.

Tip

We recommend setting the keyword argument paginate to true for datasets containing many subjects (more than 9). By doing this, observations_vs_time will return a vector of plots instead of trying to cram all of them into the limited space of a single figure.

Default

This is the default look for the plots.

obs_default = observations_vs_time(nca_data; paginate = true);
obs_default[1]

Marker/Lines Customizations

Let’s try to customize it to follow our color palette. We will start by changing the color of the markers and the lines, which can be done using the keyword arguments markercolor and color.

obs_color = observations_vs_time(
    nca_data;
    paginate = true,
    markercolor = main_colors[1],
    color = main_colors[2],
);
obs_color[1]

Labels Customizations

Now let’s change the background to our color of choice and customize the axis labels. In order to do this, we have to use the axis and figure keyword arguments. We will also use the facet keyword to combine labels where possible and avoid cluttering the plots:

obs_final = observations_vs_time(
    nca_data;
    paginate = true,
    markercolor = main_colors[1],
    color = main_colors[2],
    axis = (;
        backgroundcolor = background,
        xlabel = "Time (days)",
        ylabel = "Concentration (ng/mL)",
    ),
    figure = (; backgroundcolor = background, resolution = (1000, 800), fontsize = 18),
    facet = (; combinelabels = true),
);
obs_final[1]

Note

We have also changed the figure’s resolution and fontsize. Makie.jl offers a very wide range of customization options, so please make sure to check our tutorial on Customization of AlgebraOfGraphics.jl Plots and Makie’s tutorials for more information.

Lastly, Pumas’ documentation on plotting is also a great place to learn more about the customization options available for functions like observations_vs_time.

With just a few keyword arguments, we have changed look and feel of our plot, and we have made it our own by applying our custom color palette.

6 Exploring the effect of dosing

To explore the effect of dosing on drug concentration we will generate the average concentration time profile for each dose level. This can be accomplished by using the function summary_observations_vs_time.

Color and Labels Customizations

Once again, we will use the axis, figure, and color keyword arguments to make the result look the way we want.

f1, ax1, p1 = summary_observations_vs_time(
    nca_data;
    separate = false,
    group = [:DOSE],
    figure = (; backgroundcolor = background, resolution = (900, 500), fontsize = 16),
    axis = (;
        backgroundcolor = background,
        xlabel = "Time (days)",
        ylabel = "Concentration (ng/mL)",
    ),
    color = :Set1_5,
)

axislegend(ax1; bgcolor = background)
f1

Note

Notice that we had to run the function axislegend which accepts an Axis (in our case ax1) from Makie.jl in order to get a legend for the plot. We also had to store the figure, axis, and plot objects when calling the function summary_observations_vs_time in order to use the axislegend function and display the result.

Tip

You probably noticed that we assigned Set1_5 to the color keyword argument. Currently, summary_observations_vs_time only accepts a single color or a predefined colormap from ColorSchemes.jl. We are currently working to add support for a custom colormap and expect this feature to be available soon.

We recommend using a colormap if you are setting separate to false since this will allow you to get a different color for each line.

Y-axis Scale

The last plot is very useful to observe the behavior of the concentration and determine the maximum concentration for each dosage regimen, but other interesting information such as the number of compartments that would be needed for a PK model, half-life, and nonlinearity of clearance is better observed on a logarithmic scale. To do this, simply add yscale=log10 to axis.

f2, ax2, p2 = summary_observations_vs_time(
    nca_data;
    separate = false,
    group = [:DOSE],
    figure = (; backgroundcolor = background, resolution = (900, 500), fontsize = 16),
    axis = (;
        backgroundcolor = background,
        xlabel = "Time (days)",
        ylabel = "Concentration (ng/mL)",
        yscale = log10,
    ),
    color = :Set1_5,
)

axislegend(ax2; bgcolor = background)
f2

Tip

You can set yscale (or xscale) to be any function you like. For instance, you could set the scale to be the natural logarithm (log) or the logarithm to the base 2 (log2). Finally, you can also pass any custom user-defined function and not be restricted by Julia’s or Pumas’ existing functions.

6.1 Dose normalized concentration

In addition to the concentration time profile, it could also be interesting to analyze the dose-normalized concentration profile. Normalizing the concentration will allow us to assess the dose linearity of exposure and we will use the normalized concentration for the subsequent analysis.

We begin by creating a new column, NLIDV, which will be the concentration (LIDV) divided by the corresponding dose (DOSE):

@rtransform! df begin
    :NLIDV = :LIDV / :DOSE
end;

Now we create another NCAPopulation, but this time the observations will be the normalized concentrations.

nca_data_norm = read_nca(
    df;
    id = :ID,
    time = :NOMTIME,
    observations = :NLIDV,
    amt = :AMT,
    route = :ROUTE,
    group = [:DOSE],
);

In order to get the normalized concentration profile for each dose, we use the same code as before, just changing the NCAPopulation for the one that we just computed (this time we only show you the plot on a logarithmic scale):

f3, ax3, p3 = summary_observations_vs_time(
    nca_data_norm;
    separate = false,
    group = [:DOSE],
    figure = (; backgroundcolor = background, resolution = (900, 500), fontsize = 16),
    axis = (;
        backgroundcolor = background,
        xlabel = "Time (days)",
        ylabel = "Normalized Concentration (ng/mL)/mg",
        yscale = log10,
    ),
    color = :Set1_5,
)
axislegend(ax3; bgcolor = background)
f3

Now that we have normalized the concentration, we will calculate the AUC and Cₘₐₓ metrics. To do this, we must first run the analysis on the population by using the function run_nca.

nca_results_norm = run_nca(nca_data_norm);

In order to visualize these results, we will generate violin plots for both the AUC and the Cₘₐₓ at each dose level. Also, we would like to observe them side by side.

To do this, use the function parameters_vs_group to generate each plot, and we will put them together using the Figure object from CairoMakie.jl.

Figure

Let’s start by creating the figure:

fAUC_Cmax = Figure(; backgroundcolor = background, resolution = (1000, 500));

Plot 1

Now we will add the plots to the figure by indexing it and passing it as the first argument of the parameters_vs_group function. We will start with the AUC.

parameters_vs_group(
    fAUC_Cmax[1, 1],
    nca_results_norm,
    group = "DOSE",
    parameter = :aucinf_obs;
    axis = (;
        backgroundcolor = background,
        xlabel = "Dose (mg)",
        ylabel = "AUC (ng/mL)*day",
        title = "AUC",
    ),
    color = main_colors[1],
);

Tip

Notice that this will not return the plot. If you ever want to check how the plot is looking, inspect the variable associated with the figure (in this case AUC_Cmax).

Plot 2

Now let’s add the plot corresponding to Cₘₐₓ. In this case, we want to place it on the right of the AUC plot, so we index the figure at [1, 2]:

parameters_vs_group(
    fAUC_Cmax[1, 2],
    nca_results_norm;
    group = "DOSE",
    parameter = :cmax,
    axis = (;
        backgroundcolor = background,
        xlabel = "Dose (mg)",
        ylabel = "Cₘₐₓ (ng/mL)/mg",
        title = "Cₘₐₓ",
    ),
    color = main_colors[2],
);

And now we have completed our custom layout. Again, we can see our results by inspecting the figure.

fAUC_Cmax

Tip

This procedure can be generalized to be used with any other plotting functions, just pass the axis in which you want to place the plot (e.g figure[2, 1]) as the first argument of the function. Make sure to check our documentation section on customizing plots for a more detailed explanation.

We would also like to explore the effect of dosing on these metrics. One way in which we might go about this is by performing a linear regression with respect to dose levels. This can be achieved by using the function power_model:

First, we have to create a DoseLinearityPowerModel object. We will create one for both Cₘₐₓ and AUC, which we will name dpcmax and dpauc.

dpcmax = NCA.DoseLinearityPowerModel(nca_results_norm, :cmax);

dpauc = NCA.DoseLinearityPowerModel(nca_results_norm, :aucinf_obs);

Now we will apply the same procedure as with the violin plots in order to put the two visualizations next to each other:

Figure

fAUC_Cmax_reg =
    Figure(resolution = (1500, 800), backgroundcolor = background, fontsize = 20);

Plot 1

Let’s determine the unique dose levels in order to use them as our x-axis ticks:

doses = unique(df.DOSE);

Now we can add our plot to the Figure:

power_model(
    fAUC_Cmax_reg[1, 1],
    dpauc;
    axis = (;
        backgroundcolor = background,
        xlabel = "Dose (mg)",
        ylabel = "Dose-normalized AUC (ng/mL)*day",
        title = "AUC",
        xticks = doses,
        xscale = log10,
    ),
    color = main_colors[1],
    markercolor = main_colors[2],
);

Plot 2

Now we can add our plot to the Figure:

power_model(
    fAUC_Cmax_reg[1, 2],
    dpcmax;
    axis = (;
        backgroundcolor = background,
        xlabel = "Dose (mg)",
        ylabel = "Dose-normalized Cₘₐₓ (ng/mL)",
        title = "Cₘₐₓ",
        xticks = doses,
        xscale = log10,
    ),
    color = main_colors[1],
    markercolor = main_colors[2],
);

Legend

We could stop here, but let’s add a figure legend with figurelegend:

figurelegend(
    fAUC_Cmax_reg;
    bgcolor = background,
    position = :b,
    nbanks = 2,
    tellheight = true,
);

Again, we can see our results by inspecting our Figure (fAUC_Cmax_reg):

fAUC_Cmax_reg

7 Exploring the effects of covariates

As you may have noticed, the dataset that we are using includes information about the subject’s sex and weight. We will explore the effect of these two covariates on the concentration profiles, AUC and Cₘₐₓ

7.1 Sex

We will generate the mean concentration profile for each sex. In order to do this, we will have to create a new NCAPopulation with read_nca. However, this time we’ll specify SEX as the group variable

nca_data_sex = read_nca(
    df;
    id = :ID,
    time = :NOMTIME,
    observations = :NLIDV,
    amt = :AMT,
    route = :ROUTE,
    group = [:SEX],
);

Note

We could add any number of columns to the group keyword argument. For instance, we could have done group = [:SEX, :DOSE] to let Pumas know that we want to group subjects both by sex and by dose level. However, this will result in 10 groups (cartesian product for a set of 2 categories of sex and 5 dose levels), and we want to limit this analysis to only the sex covariate.

Y-axis Scale

Now let’s generate the concentration profile on a logarithmic scale following the same procedure as described above.

f4, ax4, p4 = summary_observations_vs_time(
    nca_data_sex;
    separate = false,
    group = "SEX",
    figure = (; backgroundcolor = background, resolution = (900, 500), fontsize = 16),
    axis = (;
        backgroundcolor = background,
        xlabel = "Time (days)",
        ylabel = "Normalized Concentration (ng/mL)/mg",
        yscale = log10,
    ),
    color = :Set1_5,
)

axislegend(ax4; bgcolor = background)
f4

Here we notice that, on average, female subjects exhibit a higher concentration profile. Let’s also compute the values for AUC and Cₘₐₓ in order to examine more closely this difference. Just like before, we must use the function run_nca on our new NCA population.

nca_results_sex = run_nca(nca_data_sex);

Violin Plots

Now we can generate violin plots for each of these metrics by using the same method as when we were exploring the effect of dosing:

fAUC_Cmax_sex = Figure(; backgroundcolor = background, resolution = (1000, 500))

parameters_vs_group(
    fAUC_Cmax_sex[1, 1],
    nca_results_sex;
    group = "SEX",
    parameter = :aucinf_obs,
    axis = (;
        xlabel = "Sex",
        ylabel = "AUC (ng/mL)*day",
        title = "AUC",
        backgroundcolor = background,
    ),
    color = main_colors[1],
)

parameters_vs_group(
    fAUC_Cmax_sex[1, 2],
    nca_results_sex;
    group = "SEX",
    parameter = :cmax,
    axis = (;
        backgroundcolor = background,
        xlabel = "Sex",
        ylabel = "Cₘₐₓ (ng/mL)/mg",
        title = "Cₘₐₓ",
    ),
    color = main_colors[2],
)

fAUC_Cmax_sex

7.2 Weight

To explore the effects of weight, we will divide subjects into those under 70 kg and those over. This means that we have to create the following column in our dataset.

@rtransform! df begin
    :WEIGHT = :WEIGHTB > 70 ? "Over 70 kg" : "Under 70 kg"
end;

Having done that, we can create an NCAPopulation setting the group as :WEIGHT.

nca_data_weight = read_nca(
    df;
    id = :ID,
    time = :NOMTIME,
    observations = :NLIDV,
    amt = :AMT,
    route = :ROUTE,
    group = [:WEIGHT],
);

f5, ax5, p5 = summary_observations_vs_time(
    nca_data_weight;
    separate = false,
    group = [:WEIGHT],
    figure = (; backgroundcolor = background, resolution = (900, 500), fontsize = 16),
    axis = (;
        backgroundcolor = background,
        xlabel = "Time (days)",
        ylabel = "Normalized Concentration (ng/mL)/mg",
        yscale = log10,
    ),
    color = :Set1_5,
)

axislegend(ax5; bgcolor = background)
f5

This indicates that drug concentration is higher for individuals weighing less than 70 kg. We could also create the violin plots for AUC and Cₘₐₓ by following the same procedure as before, but we’ll leave it here since we already went over the procedure

8 Exploring both covariates simultaneously

For our last visualization, let’s combine both covariates in a single analysis. To do this, we will plot AUC and Cₘₐₓ against body weight, and we will also examine if there are differences related to sex.

Let’s start by creating an NCA population specifying the two groups we are interested in: SEX and WEIGHT.

nca_data_sex_weight = read_nca(
    df;
    id = :ID,
    time = :NOMTIME,
    observations = :NLIDV,
    amt = :AMT,
    route = :ROUTE,
    group = [:SEX, :WEIGHTB],
);

Now that we have created the population, we can run the analysis:

nca_results_sex_weight = run_nca(nca_data_sex_weight).reportdf;

Tip

You can get the results of run_nca as a DataFrame by inspecting the field reportdf.

Scatter Plot AUC

Let’s begin by generating a scatter plot of AUC vs weight and perform a linear regression on top of it.

Note

We are going to be using AlgebraOfGraphics.jl for plotting. If you are unfamiliar with its syntax and use, feel free to refer to our series of tutorials on data visualization in Julia, starting with Introduction to AlgebraOfGraphics.jl.

auc_plot =
    data(nca_results_sex_weight) *
    mapping(:WEIGHTB, :aucinf_obs; color = :SEX => "Sex") *
    (visual(Scatter) + AlgebraOfGraphics.linear());

We can check our mapping and visuals by using the draw function

draw(auc_plot)

Scatter Plot Cₘₐₓ

It seems like we got the result we wanted. Let’s do the same for Cₘₐₓ:

cmax_plot =
    data(nca_results_sex_weight) *
    mapping(:WEIGHTB, :cmax; color = :SEX => "Sex") *
    (visual(Scatter) + AlgebraOfGraphics.linear());

We can check our mapping and visuals by using the draw function

draw(cmax_plot)

Combined Scatter Plots

This time we can trust it will look the way we want. Let’s now create a figure where both plots are shown next to each other and we apply our custom color palette:

fig_two_cov = Figure(; backgroundcolor = background, resolution = (1000, 600));
axs_auc = fig_two_cov[1, 1];
axs_cmax = fig_two_cov[1, 2];

auc = draw!(
    axs_auc,
    auc_plot,
    axis = (;
        backgroundcolor = background,
        xlabel = "Weight (kg)",
        ylabel = "AUC (ng/mL)*day",
        title = "AUC",
    ),
    palettes = (; color = main_colors),
)
cmax = draw!(
    axs_cmax,
    cmax_plot,
    axis = (;
        backgroundcolor = background,
        xlabel = "Weight (kg)",
        ylabel = "Cₘₐₓ (ng/mL)/mg",
        title = "Cₘₐₓ",
    ),
    palettes = (; color = main_colors),
)
legend!(
    fig_two_cov[end+1, 1:2],
    cmax;
    orientation = :horizontal,
    tellheight = true,
    bgcolor = background,
)
fig_two_cov

This plot lets us see that, on average, AUC and Cₘₐₓ decrease as the person’s weight increases, and that these metrics generally have a higher value for female subjects than for male subjects of the same weight.

9 Conclusions

We have given you a very broad overlook of Pumas’ plot customization capabilities for Noncompartmental analysis. We hope that you feel more comfortable with applying your custom color palette to your data visualization outputs and creating your custom layouts by mixing different plotting functions and creating your own with AlgebraOfGraphics.jl and Makie.jl.

We would like to end with an invitation to check our documentation on plotting. You’ll find that this is the most complete source of information for Pumas’ plotting capabilities. Don’t forget to check the documentation for Makie.jl and our tutorials on data visualization in Julia to get more information on the available customization options for things like figures and axes.