Different plots on same code

Simulation
1

On copying the same codes and running the simulation in another window, I tend to get different outputs.

using Pumas
using Plots
myfirstmodelinpumas = @model begin
    @param   begin
        tvka ∈ RealDomain(lower=0)
        tvcl ∈ RealDomain(lower=0)
        tvvc ∈ RealDomain(lower=0)
        tvvp ∈ RealDomain(lower=0)
        tvQ  ∈ RealDomain(lower=0)
        Ω ∈ PDiagDomain(5)
        σ_prop ∈ RealDomain(lower=0)
    end

    @random begin
      η ~ MvNormal(Ω)
    end

    @pre begin
        CL = tvcl * exp(η[1])
        Vc  = tvvc * exp(η[2])
        Ka = tvka * exp(η[3])
        Q = tvQ * exp(η[4])
        Vp = tvvp * exp(η[5])
    end

     @dynamics begin
         depot' = -Ka*depot
         central' =  Ka*depot - (CL/Vc)*central + (Q/Vp)*peripheral - (Q/Vc)*central
         peripheral' = -(Q/Vp)*peripheral + (Q/Vc)*central
     end

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

para = (tvka = 1,
          tvcl = 12.0,
          tvvc = 77.1,
          tvvp = 52.3,
          tvQ = 9.22,
          Ω = Diagonal([0.09,0.09,0.09,0.09,0.09]),
          σ_prop = 0.04)
          
evoral = DosageRegimen(120, time=0, cmt=1, ii=12, addl=5)
eviva = DosageRegimen(60, rate=40, time=1, cmt=2)
evcom = DosageRegimen(evoral, eviva)
ssotalol = Subject(id=1, evs=evcom)
obssotalol = simobs(myfirstmodelinpumas, ssotalol, para, obstimes=0:1:120)
plot(obssotalol)

pumas 1
pumas 11129×364 48.6 KB

pumas 2
pumas 21175×364 50.9 KB

why does that happen?

1 Like
2

This happens because you haven’t specified the random effects. In that case, the random effects are drawn randomly and the simulated will therefore be different for each time you invoke the simobs functions. Hence, it’s not related to extra window.

3

Got it.
Are there ways by which I can reduce the variability in the output.
Or including a function that could fix a range to my random effects.?

4

You can pass the exact random effects you want to use in the simulation. It’s the (optional) fourth argument to simobs. Hence, something like

obssotalol = simobs(myfirstmodelinpumas, ssotalol, para, (η=zeros(5),), obstimes=0:1:120)

would fix them all to zero and the top plots (cp) should be identical across invocations.

Alternatively, you can fix the global seed with

using Random
Random.seed!(SOMESEEDINTEGER)

Notice, that you’d also need to disable multithreading in simobs with the ensemblealg = EnsembleSerial() keyword argument if you want the simulation to be reproducible.

1 Like
5

That answered my question.!
Thanks Andreas.! :slightly_smiling_face: