Overlay `observations_vs_time` plot on `sim_plot`

1

Sometimes it is useful to have a scatter plot of the observations overlayed on the simulation results (not vpc) to get a quick idea of how a single simulation matches the data.

It can be done quickly using the following code

fig, ax, plt = sim_plot(model, simdata,  observations = [:Cvenn])
observations_vs_time!(ax, model, data ; loess = false )
display(fig)
1 Like
3

Hello @Vaibhavdixit02

This will certainly allow you to get a scatter plot of the observations overlayed on the simulation results. There is another way to do this using only sim_plot, so I will share an example of how you could do that.

We could use the nlme_sample dataset from PharmaDatasets for this example. You can load it and convert it to a population with the following code.

using Pumas, PumasUtilities
using PharmaDatasets # Get dataset
using CairoMakie # Customize plots and show legend

# Get dataset
pkdata = dataset("nlme_sample")

# Create population 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,
)

Now we can define a 1 compartment model to be used in our simulation:

# 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

Also, we should define the initial values for the model parameters. In this case, we can use the following.

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

Let’s choose one of the subjects from our population. I will be choosing one using the id (say, subject 8), but you could use different criteria depending on your needs.

subject = filter(i -> i.id == "8", population)

Now that we have our model, initial parameters, and subject, we can run the simulation using simobs.

simulation = simobs(model, subject, params)

Having done that, we can generate the plot using sim_plot.

sim_plot(
    model,
    simulation;
    observations=[:DV],
    axis=(;
        xlabel="Time (hr)",
        ylabel="Observed/Predicted \n Concentration (ng/mL)"
    ),
    color=:red
)

This will return a scatter plot of the observations and the simulated values as a red line on top. However, this plot will not have a legend, which I think would be very useful for visualizing the results.

To get a legend, we will take advantage of the fact that Makie.jl returns a tuple containing the figure, axis, and plot objects. We will store these values using the following syntax.

fig, ax, plt = sim_plot(
    model,
    simulation;
    observations=[:DV],
    axis=(;
        xlabel="Time (hr)",
        ylabel="Observed/Predicted \n Concentration (ng/mL)"
    ),
    color=:red
)

Now we can display the legend by using axislegend and then inspect fig to see the final result. We also use the keyword argument position to tell Makie to place the legend on the upper left portion of the figure.

axislegend(ax, position=:lt)
fig

This will produce the following result:

legend
legend800×600 30.8 KB

1 Like