DV column in the output from predict function

Diagnostics
1

Hi,
I’m trying to get a diagnostic plot from modeling result.
When I use predict function, I can’t find dv column in the output. Because we usually compare dv vs pred or dv vs ipred so I think dv should be included in the output of predict function. Is there any way to handle this?

2

Hi.

I also have the same question. Do we need to merge observed DV with output from predict() to be able to make diagnostic plots?

3

Thanks for reporting @wpark. This has been reported and we are working on having this dv in the predict output in the next release, very soon.

4

Until the next release it out, a more complete example always makes it easier to help :slight_smile: A so called minimal working example (https://en.wikipedia.org/wiki/Minimal_working_example) makes it easier to show you the steps needed to do what you want with the current version. I will try to find time later today to a) find time to show how to do it currently, and b) make a PR to make sure that this is included in the next versions .

1 Like
5

I agree with you that I should have shown a working example. It’s a little late, but I made an example for the future.

#Source from https://tutorials.pumas.ai/html/exercises/workshop_solutions.html

using Pumas, Plots, CSV, Random
Random.seed!(0)

single_dose_regimen = DosageRegimen(100, time=0)
first(single_dose_regimen.data)
s1 = Subject(id=1, evs=single_dose_regimen,cvs=(Wt=70,))
choose_covariates() = (Wt = rand(55:80),)

pop = Population(map(i -> Subject(id = i,evs = single_dose_regimen, cvs =  choose_covariates()),1:24))

mymodel = @model begin
  @param   begin
    tvcl ∈ RealDomain(lower=0, init = 1.0)
    tvv ∈ RealDomain(lower=0, init = 20)
    tvka ∈ RealDomain(lower = 0, init= 1)
    Ω ∈ PDiagDomain(init=[0.09,0.09, 0.09])
    σ_prop ∈ RealDomain(lower=0,init=0.04)
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @pre begin
    CL = tvcl * (Wt/70)^0.75 * exp(η[1])
    V  = tvv * (Wt/70) * exp(η[2])
    Ka = tvka * exp(η[3])
  end
  @covariates Wt

  @dynamics Depots1Central1
    #@dynamics begin
    #    Depot' =  -Ka*Depot
    #    Central' =  Ka*Depot - (CL/V)*Central
    #end

  @derived begin
      cp = @. 1000*(Central / V)
      dv ~ @. Normal(cp, sqrt(cp^2*σ_prop))
    end
end

param = init_param(mymodel)
obs = simobs(mymodel, pop, param, obstimes=0:1:72)
#plot(obs)

simdf = DataFrame(obs)
simdf[!, :route] .= "ev"
#simdf = copy(obs)

simdf.cmt = ifelse.(ismissing.(simdf.cmt), 2, simdf.cmt)
est_df = simdf[.!((simdf.dv .== 0.0) .& (simdf.cmt .==2)),:]
#first(est_df,6)

data = read_pumas(est_df ,cvs = [:Wt], dvs=[:dv])
res = fit(mymodel,data,param,Pumas.FOCEI())
preds = DataFrame(predict(res))

After run this code, I got this result.

julia> first(preds, 6)
6×6 DataFrame
│ Row │ id     │ time    │ Wt    │ dv_pred │ dv_ipred │ pred_approx │
│     │ String │ Float64 │ Int64 │ Float64 │ Float64  │ Pumas.FOCEI │
├─────┼────────┼─────────┼───────┼─────────┼──────────┼─────────────┤
│ 1   │ 1      │ 1.0     │ 55    │ 3727.28 │ 5589.99  │ FOCEI()     │
│ 2   │ 1      │ 2.0     │ 55    │ 5212.33 │ 6498.21  │ FOCEI()     │
│ 3   │ 1      │ 3.0     │ 55    │ 5546.78 │ 6347.51  │ FOCEI()     │
│ 4   │ 1      │ 4.0     │ 55    │ 5472.67 │ 5976.31  │ FOCEI()     │
│ 5   │ 1      │ 5.0     │ 55    │ 5275.63 │ 5577.32  │ FOCEI()     │
│ 6   │ 1      │ 6.0     │ 55    │ 5048.29 │ 5193.63  │ FOCEI()     │

This is my package status.

(@v1.3) pkg> st
     Status `C:\Users\user\.juliapro\JuliaPro_v1.3.1-1\environments\v1.3\Project.toml`
   c52e3926 Atom v0.11.3 ⚲
   336ed68f CSV v0.5.22
   a93c6f00 DataFrames v0.20.0
   1313f7d8 DataFramesMeta v0.5.0
   82cc6244 DataInterpolations v1.3.1
   0c46a032 DifferentialEquations v6.10.1
   7073ff75 IJulia v1.20.2
   e5e0dc1b Juno v0.7.2 ⚲
   4722fa14 PkgAuthentication v0.1.2
   91a5bcdd Plots v0.27.1
   4f2c3c20 Pumas v0.9.0
   1a8c2f83 Query v0.12.2
   37e2e46d LinearAlgebra
1 Like
7

The dependent variables are included in the DataFrame in the next release of Pumas.