Hi Jay Shree,
I am sharing an example. I hope this will be helpful. In this example the PKPD estimation is divided into two parts - In the first part, the PK parameters are estimated. Individual parameters from the pk fit are used for PD estimation. Individual clearance and volume of distribution are used to calculate the concentration for the PD analysis. In the second part, the PD parameters are estimated using the PD model. Please let me know if this approach works for your model.
using Pumas
using CSV
using StatsBase
using Random
using DataFramesMeta
using StatsPlots
using PumasPlots
using Statistics
pkdata = CSV.File("iv_inf_pd_direct_add.csv"; missingstring="") |> DataFrame
pkdata.conc_obs = ifelse.(pkdata.time .== 0.0, missing, pkdata.conc_obs)
data4pkest = read_pumas(pkdata,
id = :id,
time = :time,
amt = :amt,
rate = :rate,
cmt = :cmt,
evid = :evid,
observations = [:conc_obs])
pk1cmt = @model begin
@param begin
tvcl ∈ RealDomain(lower=0.1)
tvvc ∈ RealDomain(lower=0.1)
Ω ∈ PDiagDomain(2)
σ_add_pk ∈ RealDomain(lower=0)
end
@random begin
η ~ MvNormal(Ω)
end
@pre begin
CL = tvcl * exp(η[1])
Vc = tvvc * exp(η[2])
end
@dynamics begin
Central' = - (CL/Vc)*Central
end
@derived begin
cp_ipred = @. (Central/Vc)
conc_obs ~ @. Normal(cp_ipred, sqrt(σ_add_pk))
end
end
pkparams = (tvcl = 1.5,
tvvc = 25.0,
Ω = Diagonal([0.09, 0.09]),
σ_add_pk = 2.0)
pk_1cmt_results = fit(pk1cmt, data4pkest, pkparams, Pumas.FOCE())
pk_1cmt_inspect = DataFrame(inspect(pk_1cmt_results))
pk_1cmt_icoef = DataFrame(icoef(pk_1cmt_results))
select!(pk_1cmt_icoef, ([:id,:CL,:Vc]))
rename!(pk_1cmt_icoef,:CL => :iCL, :Vc => :iVc)
pk_1cmt_icoef.id = map(x->parse(Int,x), pk_1cmt_icoef.id)
pkpddata_pkindvparams = innerjoin(pkpddata, pk_1cmt_icoef, on=:id)
# PD Modeling using Individual PK estimates
data4pdest = read_pumas(pkpddata_pkindvparams,
id = :id,
time = :time,
amt = :amt,
rate = :rate,
cmt = :cmt,
evid = :evid,
observations = [:resp_obs],
covariates = [:iCL,:iVc,:dose]
)
pd_add = @model begin
@param begin
tvbsl ∈ RealDomain(lower=0)
tvec50 ∈ RealDomain(lower=0)
tvemax ∈ RealDomain(lower=0)
Ω ∈ PDiagDomain(1)
σ_add_pd ∈ RealDomain(lower=0)
end
@random begin
η ~ MvNormal(Ω)
end
@covariates iCL iVc dose
@pre begin
CL = iCL
Vc = iVc
bsl = tvbsl * exp(η[1])
ec50 = tvec50
emax = tvemax
end
@dynamics begin
Central' = - (CL/Vc)*Central
end
@derived begin
cp_ipred = @. (Central/Vc)
resp_ipred = @. bsl + emax*cp_ipred/(ec50 + cp_ipred)
resp_obs ~ @. Normal(resp_ipred, sqrt(σ_add_pd))
end
end
pd_params = (tvbsl = 10,
tvec50 = 25,
tvemax = 10,
Ω = Diagonal([0.09]),
σ_add_pd = 2)
pd_fit = fit(pd_add, data4pdest, pd_params, Pumas.FOCE())
pd_fit_infer = infer(pd_fit)