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:

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.