NPDE vs prediction plots

1

Hi all,
I wanted to generate NPDE_vs_predictions and NPDE_vs_time plots. Can anyone help me out in inserting the the percentiles for the predicted and observed data and the prediction bands.
The npde_vs_predictions (inspect) code doesn’t provide the keyword to insert these intervals.

Regards,
Sucharitha

2

Hello @sucharitha

Here is a code example of how you could add percentiles to an NPDE vs time plot for a dataset from PharmaDatasets.jl.

We first load the data, define the model, fit it, and run the inspection:

using Pumas, PumasUtilities
using PharmaDatasets # Get dataset
using CairoMakie # Customize plots and show legend
using PlottingUtilities # generate the pair plot
using DataFramesMeta # Data wrangling
using AlgebraOfGraphics # Abbreviated to AoG

# Get dataset
pkdata = dataset("nlme_sample")

# Create popoulation with read_pumas
population = read_pumas(
    pkdata;
    id = :ID,
    time = :TIME,
    amt = :AMT,
    covariates = [:WT, :AGE, :SEX, :CRCL, :GROUP],
    observations = [:DV],
    cmt = :CMT,
    evid = :EVID,
    rate = :RATE,
)

# Create the model
model = @model begin
    @metadata begin
        desc = "base model: 1comp"
        timeu = u"hr"
    end

    @param begin
        "Clearance (L/hr)"
        tvcl ∈ RealDomain(; lower = 0.0001)
        "Volume (L)"
        tvvc ∈ RealDomain(; lower = 0.0001)
        """
        - ΩCL
        - ΩVc
        """
        Ω ∈ PDiagDomain(2) #covariance matrix 
        "Additive RUV"
        σ²_add ∈ RealDomain(; lower = 0.0001) # variance 
        "Proportional RUV"
        σ²_prop ∈ RealDomain(; lower = 0.0001) # variance  
    end

    @random begin
        η ~ MvNormal(Ω)
    end

    @covariates WT AGE SEX CRCL GROUP

    @pre begin
        CL = tvcl * exp(η[1])
        Vc = tvvc * exp(η[2])
    end

    @dynamics Central1

    @derived begin
        CONC := @. Central / Vc
        """
        DrugY Concentration (ng/mL)
        """
        DV ~ @. Normal(CONC, sqrt(CONC^2 * σ²_prop + σ²_add)) # using variance
    end
end

params = (;
    tvvc = 5, 
    tvcl = 0.2, 
    Ω = Diagonal([0.01, 0.01]), 
    σ²_add = 0.01, 
    σ²_prop = 0.01
)

# Fit the model
fitted_model = fit(model, population, params, FOCE())

# Inspect the model
inspection = inspect(fitted_model; nsim=1000)

Now that we have our inspection results (inspection), we can create our NPDE vs time plot using npde_vs_time as you have probably already done:

fig, ax, plt = npde_vs_time(
    inspection;
    loess=false,
    zeroline=false,
    ols=false,
    markercolor=:grey
)

Notice that we have saved the results from npde_vs_time since we will be adding the percentiles plot on top of it. In order to do this, we will convert our inspection results to a DataFrame so that we can extract the NPDE values and create the plots:

inspection_data = DataFrame(inspection)
@rsubset! inspection_data !ismissing(:DV)

df = @by inspection_data :time begin
    :q5 = quantile(:DV_npde, 0.05)
    :q50 = quantile(:DV_npde, 0.5)
    :q95 = quantile(:DV_npde, 0.95)
end

We calculated the percentiles by grouping our data by time (or predictions) and using the quantile function from Statistics. Here we chose the 5th, 50th, and 95th percentiles, but you can easily adjust these values in the quantile calls.

Now we change our DataFrame to a long format for plotting and we define our plot using AlgebraOfGraphics.jl:

df = stack(df, Not(:time))

quantiles_plt = data(df) *
                    mapping(
                        :time,
                        :value;
                        color=:variable
                    ) *
                    visual(Lines; linewidth=3)

To calculate the bands for our percentiles we use a similar procedure, this time grouping by our percentile value:

df_bands = @by df :variable begin
    :ql = quantile(:value, 0.025)
    :qh = quantile(:value, 0.975)
end

bands_plt = data(df_bands) *
                mapping(
                    :ql,
                    :qh;
                    color=:variable
                ) *
                visual(HSpan; alpha=0.2)

We chose a 95% CI for this case, but that could be modified in the quantile calls. Now we can add our plots to our Axis (ax) with draw!:

plt = draw!(
    ax,
    quantiles_plt + plt_bands;
    palettes=(;
        color=[:blue, :red, :blue]
    )
)

To see the results, inspect the Figure variable (fig in this case). You should get the following output:

plot_50
plot_50800×600 55.9 KB

I hope you find this useful. Please let me know if this is what you were expecting or if you have any questions about the code.

3 Likes