Add basic Pumas NLME tutorial scripts

Pumas_NLME/01-population.jl

@@ -0,0 +1,54 @@
+using DeepPumas
+using CairoMakie
+using PharmaDatasets
+# Set a nice plotting theme. There's also `deep_dark()`.
+set_theme!(deep_light()) 
+
+pkdata = dataset("nlme_sample")
+
+# check that we have the standard NM-TRAN column names such as
+# - id
+# - time
+# - dv
+# - amt
+# - evid
+# - cmt
+# - rate
+# We also have a covariates AGE and SEX
+pop = read_pumas(
+    pkdata;
+    id = :ID,
+    time = :TIME,
+    amt = :AMT,
+    covariates = [:AGE, :SEX],
+    observations = [:DV],
+    cmt = :CMT,
+    evid = :EVID,
+    rate = :RATE,
+)
+
+# The function that parses data into a Population is read_pumas
+# Check the docstring with ?read_pumas in your Julia REPL
+# Specifically, the keyword arguments
+
+#?read_pumas
+
+# Now let's read our DataFrame into a Population with read_pumas
+
+# A Population is simply a vector of Subjects
+pop[1]
+
+# You can also slice it same as with any vector
+pop[5:10]
+pop[begin:30]
+pop[80:end]
+
+# We can also convert back to a NM-TRAN DataFrame by using the DataFrame constructor
+reconstructed_pkdata = DataFrame(pop)
+
+# Or a single Subject of the Population
+reconstructed_subject = DataFrame(pop[1])
+
+# And, we can plot the individual timecourses. Slicing is useful to avoid getting too many
+# subplots at once.
+plotgrid(pop[1:12])

Pumas_NLME/02-model.jl

@@ -0,0 +1,46 @@
+include("01-population.jl")
+
+# Model definition
+model = @model begin
+    @metadata begin
+        desc = "1-compartment model"
+        timeu = u"hr"
+    end
+
+    @param begin
+        # here we define the parameters of the model
+        tvcl ∈ RealDomain(; lower = 0)
+        tvvc ∈ RealDomain(; lower = 0)
+        Ω ∈ PDiagDomain(2)
+        σ_add ∈ RealDomain(; lower = 0)
+        σ_prop ∈ RealDomain(; lower = 0)
+    end
+
+    @random begin
+        # here we define random effects
+        η ~ MvNormal(Ω)
+    end
+
+    @covariates AGE SEX
+
+    @pre begin
+        # pre computations and other statistical transformations
+        CL = tvcl * exp(η[1])
+        Vc = tvvc * exp(η[2])
+    end
+
+    @dynamics begin
+        # here we define compartments and dynamics
+        Central' = -CL / Vc * Central
+    end
+
+    @derived begin
+        # here is where we calculate concentration and add residual error
+        # tilde (~) means "distributed as"
+        CONC := @. Central / Vc
+        """
+        Drug Concentration (ng/mL)
+        """
+        DV ~ @. Normal(CONC, sqrt(CONC^2 * σ_prop^2 + σ_add^2))
+    end
+end

Pumas_NLME/03-fit.jl

@@ -0,0 +1,26 @@
+include("02-model.jl")
+
+# Initial parameter values
+params = (; tvvc = 5, tvcl = 0.2, Ω = Diagonal([0.01, 0.01]), σ_add = 0.1, σ_prop = 0.1)
+# Fit a the model with FOCE
+fit_foce = fit(model, pop, params, FOCE())
+
+# Fit a the model with NaivePooled
+fit_naivepooled = fit(model, pop, params, NaivePooled(); omegas = (:Ω,))
+
+# Fit a the model with LaplaceI
+fit_laplace = fit(model, pop, params, LaplaceI())
+
+# Fit a the model with FOCE and fixed parameters
+fit_foce_fixed = fit(model, pop, params, FOCE(); constantcoef = (; tvcl = 1.0))
+
+# Get a NamedTuple of the estimated parameter values
+coef(fit_foce)
+coef(fit_naivepooled)
+
+# Get a DataFrame of the estimated parameter values
+coeftable(fit_foce)
+coeftable(fit_naivepooled)
+
+# Get the icoefs from the model fit as a DataFrame
+DataFrame(icoef(fit_foce))

Pumas_NLME/04-predict.jl

@@ -0,0 +1,29 @@
+include("03-fit.jl")
+
+# Our original data has observations from 0 to 168 hours
+unique(pkdata.TIME)
+
+# predict generate pred/ipred and can take either:
+# - fit result
+# - fit result, population
+# - model, population, parameters
+preds = predict(fit_foce)
+preds = predict(fit_foce, pop)
+preds = predict(model, pop, coef(fit_foce))
+
+DataFrame(preds)
+
+# By default it will generate pred/ipred using the original data observation time profile
+# Suppose you want a more richer/denser pred/ipred time profile
+# You can do that with the keyword argument obstimes
+# it will "extend" the original observation profile to encompass the desired obstimes
+preds_custom = predict(fit_foce; obstimes = 168:172)
+DataFrame(preds_custom)
+
+# We can also plot the predictions
+plotgrid(preds[1:6])
+
+# Extending the observation times provides more detail on what happens between the 
+# timepoints for which we have data
+preds_dense = predict(fit_foce; obstimes = 0:172)
+plotgrid(preds_dense[1:6])

Pumas_NLME/05-compare.jl

@@ -0,0 +1,136 @@
+using PumasUtilities
+
+include("04-predict.jl")
+
+# Create a second model for comparison, this time a two-compartment PK model
+model_2cmt = @model begin
+    @metadata begin
+        desc = "2-compartment model"
+        timeu = u"hr"
+    end
+
+    @param begin
+        """
+        Clearance (L/hr)
+        """
+        tvcl ∈ RealDomain(; lower = 0)
+        """
+        Volume Central Compartment (L)
+        """
+        tvvc ∈ RealDomain(; lower = 0)
+        """
+        Intercompartmental Clearance (L/hr)
+        """
+        tvq ∈ RealDomain(; lower = 0)
+        """
+        Volume Peripheral Compartment (L)
+        """
+        tvvp ∈ RealDomain(; lower = 0)
+        """
+          - ΩC
+          - ΩV
+        """
+        Ω ∈ PDiagDomain(2)
+        """
+        Additive RUV
+        """
+        σ_add ∈ RealDomain(; lower = 0)
+        """
+        Proportional RUV
+        """
+        σ_prop ∈ RealDomain(; lower = 0)
+    end
+
+    @random begin
+        η ~ MvNormal(Ω)
+    end
+
+    @pre begin
+        CL = tvcl * exp(η[1])
+        Vc = tvvc * exp(η[2])
+        Q = tvq
+        Vp = tvvp
+    end
+
+    @dynamics Central1Periph1
+
+    @derived begin
+        CONC := @. Central / Vc
+        DV ~ @. Normal(CONC, sqrt(CONC^2 * σ_prop^2 + σ_add^2))
+    end
+end
+
+params_2cmt = (;
+    tvvc = 5,
+    tvcl = 0.02,
+    tvq = 0.01,
+    tvvp = 10,
+    Ω = Diagonal([0.01, 0.01]),
+    σ_add = 0.1,
+    σ_prop = 0.1,
+)
+
+fit_2cmt = fit(model_2cmt, pop, params_2cmt, FOCE())
+
+
+fit_1cmt = fit_foce # Make the name more clear
+
+# You can get all model metrics from a fit result with metrics_table
+# It returns a DataFrame
+metrics_table(fit_1cmt)
+metrics_table(fit_2cmt)
+
+# Additionally, everything in metrics_table can be individually retrieved with specialized functions
+loglikelihood(fit_1cmt)
+aic(fit_1cmt)
+bic(fit_1cmt)
+ηshrinkage(fit_1cmt)
+ϵshrinkage(fit_1cmt)
+
+# One highlight is the log-likelihood function
+# It can compute a loglikelihood for any model given any population, parameter values, and estimation method
+# This is helpful for model conversions from other software/tools
+loglikelihood(
+    model,
+    pop,
+    (; # NamedTuple of parameter values
+        tvcl = 0.2,
+        tvvc = 5,
+        Ω = Diagonal([0.1, 0.1]),
+        σ_add = 0.01,
+        σ_prop = 0.05,
+    ),
+    FOCE(),
+)
+
+# You can use the inspect function to calculate in one go:
+# - population predictions (pred) and individual predictions (ipred)
+# - weighted residuals (wres)
+# - empirical bayes estimates (ebe)
+# - individual parameters (icoef)
+# - dose control parameters (dcp), if applicable
+# It takes as input a single fitted Pumas model
+inspect_1cmt = inspect(fit_1cmt)
+inspect_2cmt = inspect(fit_2cmt)
+
+# There are several plotting functions available in Pumas
+# We will not cover all of them but one deserves a highlight: goodness_of_fit
+# It takes a Pumas result from inspect and outputs a 4-panel plot with:
+# 1. observations versus pred
+# 2. observation versus ipred
+# 3. wres versus time
+# 4. wres versus ipred
+goodness_of_fit(inspect_1cmt)
+goodness_of_fit(inspect_2cmt)
+
+
+# We can overlay individual plots for additional comparisons
+plt = plotgrid(preds_dense[1:6])
+predict_2cmt_dense = predict(fit_2cmt; obstimes=0:0.5:180)
+# The ! in a function name indicates that something is being modified rather than created.
+plotgrid!(
+    plt,
+    predict_2cmt_dense[1:6]; 
+    pred = (; color=:green, label = "2cmt pred", linestyle=:dash),
+    ipred = (; color=:red, label = "2cmt ipred", linestyle=:dash),
+    )

Pumas_NLME/06-vpc.jl

@@ -0,0 +1,16 @@
+include("05-compare.jl")
+
+# You can do visual predictive checks (VPC) with any fitted Pumas model
+# First, use the vpc function
+vpc_1cmt = vpc(fit_1cmt)
+vpc_2cmt = vpc(fit_2cmt)
+
+# Then you can plot the vpc result with vpc_plot
+vpc_plot(vpc_1cmt)
+vpc_plot(vpc_2cmt)
+
+# You can also use a prediction-corrected VPC with the keyword argument prediction_correction
+vpc_1cmt_pc = vpc(fit_1cmt; prediction_correction = true)
+vpc_2cmt_pc = vpc(fit_2cmt; prediction_correction = true)
+vpc_plot(vpc_1cmt_pc)
+vpc_plot(vpc_2cmt_pc)