Unable to specify a correct prior distribution

Estimation bayesian
1

Hi,
I am having difficulty specifying a prior distribution in Pumas with lower bound and initials, my codes is:

    @param   begin

      θ ~ MvNormal([16.6,13.1,4.05,4.89,0.623,0,0,0],
      Matrix{Float64}([0.0895	0.0167	0.00318	0.0033	0.000907	0	0	0;
                       0.0167	0.0655	-0.00765	0.00333	-0.000917	0	0	0;
                       0.00318	-0.00765	0.0238	0.0147	-0.000225	0	0	0;
                       0.0033	0.00333	0.0147	0.0161	-0.000243	0	0	0;
                       0.000907	-0.000917	-0.000225	-0.000243	0.0016	0	0	0;
                       0	0	0	0	0	1000000	0	0;
                       0	0	0	0	0	0	1000000	0;
                       0	0	0	0	0	0	0	1000000]),lower=[0,0,0,0,-Inf,-Inf,-Inf,-Inf],init=[18,14.2,2.13,4.29,1,1,1,1])

      Ω  ~ Wishart(3, Matrix{Float64}([0.282	0	0; 0	0.149	0; 0	0	0.0377]))

      σ_prop ∈ RealDomain(lower = 0)

    end

I mainly had 2 error message:
1.

LoadError: MethodError: no method matching MvNormal(::Array{Float64,1}, ::Array{Float64,2}; lower=[0.0, 0.0, 0.0, 0.0, -Inf, -Inf, -Inf, -Inf])

LoadError: ArgumentError: no Domain(::Wishart{Float64,PDMats.PDMat{Float64,Array{Float64,2}}}) constructor defined

I have google both of them and didn’t find anything matched…

Basically I am trying to specify a multivariate normal distribution for the Thetas with lower bound and initials, a Wishart distribution for the Omegas with initials.

May I have any recommendations on which part I missed in my codes?

Thanks!
Anqi

2

Hello,
Just follow up with the complete coding for this topic:

@time using Distributions, Pumas, CSV, TableView, StatsPlots, DataFrames, LinearAlgebra

ped_mod = @model begin

    @param   begin


      θ ~ Constrained(MvNormal([16.6,13.1,4.05,4.89,0.623,0,0,0],
      Matrix{Float64}([0.0895	0.0167	0.00318	0.0033	0.000907	0	0	0;
                       0.0167	0.0655	-0.00765	0.00333	-0.000917	0	0	0;
                       0.00318	-0.00765	0.0238	0.0147	-0.000225	0	0	0;
                       0.0033	0.00333	0.0147	0.0161	-0.000243	0	0	0;
                       0.000907	-0.000917	-0.000225	-0.000243	0.0016	0	0	0;
                       0	0	0	0	0	1000000	0	0;
                       0	0	0	0	0	0	1000000	0;
                       0	0	0	0	0	0	0	1000000])),lower=[0,0,0,0,-Inf,-Inf,-Inf,-Inf])


      Ω  ~ InverseWishart(3, diagm([0.282, 0.149, 0.0377]))


      σ_prop ∈ RealDomain(lower = 0)

    end

    @random begin
      η ~ MvNormal(Ω)
    end

    @covariates WTKG  CRCL

    @pre begin
    #--INSERT COVARIATE EFFECTS
    COV1=(CRCL/100)^θ[5]
    COV2=(WTKG/70)^θ[6]
    COV3=(WTKG/70)^θ[7]
    COV4=(WTKG/70)^θ[8]


    TVCLI = θ[1]*COV1*COV4
    TVCL  = TVCLI

    TVV1I = θ[2]*COV2
    TVV1  = TVV1I

    TVQI  = θ[3]
    TVQ   = TVQI

    TVV2I = θ[4]*COV3
    TVV2  = TVV2I

    CL  = TVCL*exp(η[1])
    V1  = TVV1*exp(η[2])
    Q   = TVQ
    V2  = TVV2*exp(η[3])

    Vc = V1
    Vp = V2

    S1 = V1

    #CALCULATION OF SECONDARY PARAMETERS
    KE  = CL/V1
    K12 = Q/V1
    K21 = Q/V2
    AA = KE+K12+K21
    ALPH = (AA+sqrt(AA*AA-4*KE*K21))/2
    BETA = (AA-sqrt(AA*AA-4*KE*K21))/2


    end

    @dynamics Central1Periph1 #a two compartment model


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

inputDataset = CSV.read("dat.csv")
df = read_pumas(inputDataset, id = :ID, dvs =[:DV],  cvs=[:WTKG, :CRCL], evid=:EVID, amt=:AMT,cmt=:CMT, rate=:RATE, time=:TIME)


param = (
    θ=([18,14.2,2.13,4.29,1,1,1,1]),
    Ω  = Diagonal([0.25,0.25,0.25]),
    σ_prop = 0.04)

 pkres = fit(ped_mod, df, param, Pumas.FOCEI(),
              optimize_fn=Pumas.DefaultOptimizeFN(show_trace=true, extended_trace=false))

The error is saying no method matching inverse:

LoadError: MethodError: no method matching inverse!(::SubArray{Float64,1,Array{Float64,1},Tuple{UnitRange{Int64}},true}, ::Pumas.PSDTransform, ::Diagonal{Float64,Array{Float64,1}})

Closest candidates are:

inverse!(::AbstractArray{T,1} where T<:Real, !Matched::TransformVariables.ArrayTransform, ::AbstractArray) at C:\Users\.juliapro\JuliaPro_v1.4.2-2\packages\TransformVariables\a4AMY\src\aggregation.jl:79

inverse!(::AbstractArray, !Matched::Pumas.ElementArrayTransform, ::AbstractArray) at C:\Users\.juliapro\JuliaPro_v1.4.2-2\packages\Pumas\NCmSe\src\estimation\transforms.jl:25

inverse!(::AbstractArray, !Matched::Pumas.ConstantTransform, ::Any) at C:\Users\juliapro\JuliaPro_v1.4.2-2\packages\Pumas\NCmSe\src\estimation\transforms.jl:54

...

in expression starting at C:\Users\118

_inverse!_tuple(::Array{Float64,1}, ::Tuple{Pumas.ElementArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},Pumas.PSDTransform,TransformVariables.ShiftedExp{true,Int64}}, ::Tuple{Array{Float64,1},Diagonal{Float64,Array{Float64,1}},Float64}) at aggregation.jl:196

inverse!(::Array{Float64,1}, ::TransformVariables.TransformTuple{NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Pumas.ElementArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},Pumas.PSDTransform,TransformVariables.ShiftedExp{true,Int64}}}}, ::NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Array{Float64,1},Diagonal{Float64,Array{Float64,1}},Float64}}) at aggregation.jl:237

inverse(::TransformVariables.TransformTuple{NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Pumas.ElementArrayTransform{TransformVariables.ShiftedExp{true,Float64},1},Pumas.PSDTransform,TransformVariables.ShiftedExp{true,Int64}}}}, ::NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Array{Float64,1},Diagonal{Float64,Array{Float64,1}},Float64}}) at generic.jl:206

fit(::PumasModel{ParamSet{NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Constrained{MvNormal{Float64,PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}}},InverseWishart{Float64,PDMats.PDMat{Float64,Array{Float64,2}}},RealDomain{Int64,TransformVariables.Infinity{true},Float64}}}},var"#548#603",var"#549#604",var"#550#606",Central1Periph1,var"#551#607",var"#577#633"}, ::Array{Subject{NamedTuple{(:DV,),Tuple{Array{Union{Missing, Float64},1}}},Pumas.ConstantCovar{NamedTuple{(:WTKG, :CRCL),Tuple{Float64,Float64}}},Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1},Array{Float64,1},Float64},1}, ::NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Array{Float64,1},Diagonal{Float64,Array{Float64,1}},Float64}}, ::Pumas.FOCEI; optimize_fn::Pumas.DefaultOptimizeFN{Nothing,NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases),Tuple{Bool,Bool,Bool,Float64,Bool}}}, constantcoef::NamedTuple{,Tuple{}}, omegas::Tuple{}, ensemblealg::EnsembleSerial, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{,Tuple{}}}) at likelihoods.jl:1330

(::StatsBase.var"#fit##kw")(::NamedTuple{(:optimize_fn,),Tuple{Pumas.DefaultOptimizeFN{Nothing,NamedTuple{(:show_trace, :store_trace, :extended_trace, :g_tol, :allow_f_increases),Tuple{Bool,Bool,Bool,Float64,Bool}}}}}, ::typeof(fit), ::PumasModel{ParamSet{NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Constrained{MvNormal{Float64,PDMats.PDMat{Float64,Array{Float64,2}},Array{Float64,1}},VectorDomain{Array{Float64,1},Array{TransformVariables.Infinity{true},1},Array{Float64,1}}},InverseWishart{Float64,PDMats.PDMat{Float64,Array{Float64,2}}},RealDomain{Int64,TransformVariables.Infinity{true},Float64}}}},var"#548#603",var"#549#604",var"#550#606",Central1Periph1,var"#551#607",var"#577#633"}, ::Array{Subject{NamedTuple{(:DV,),Tuple{Array{Union{Missing, Float64},1}}},Pumas.ConstantCovar{NamedTuple{(:WTKG, :CRCL),Tuple{Float64,Float64}}},Array{Pumas.Event{Float64,Float64,Float64,Float64,Float64,Float64,Int64},1},Array{Float64,1},Float64},1}, ::NamedTuple{(:θ, :Ω, :σ_prop),Tuple{Array{Float64,1},Diagonal{Float64,Array{Float64,1}},Float64}}, ::Pumas.FOCEI) at likelihoods.jl:1328

top-level scope at tp4_5.jl:118

Thanks!
Anqi

3

Error 2. in your first post is a missing feature in Pumas but it’s tiny so we might be able to consider it a bugfix and include it Pumas 1.0. Right now, we only support the InverseWishart but extending the support to cover Wishart should be very easy.

The problem in your followup post is Ω = Diagonal([0.25,0.25,0.25]) because it restricts Ω do a diagonal matrix but the InverseWishart isn’t restricted to diagonal matrices so you should be able to make this work by setting Ω = diagm([0.25,0.25,0.25]) instead. We might be able to make your version to work as well. I’ve filed an issue.

4

Thank you! The error is fixed after I changed the initial to a diagm.

And may I just confirm that, earlier this morning you mentioned the prior is just a domain(at least for the current version), which sounds it won’t enter the posterior distribution like a real prior.
Is this the case for both FOCE and HMC-NUTS?
Or the distributional domain(prior) does enter the posterior as a real prior in HMC-NUTS?

Thanks!
Anqi

5

And may I just confirm that, earlier this morning you mentioned the prior is just a domain(at least for the current version), which sounds it won’t enter the posterior distribution like a real prior.
Is this the case for both FOCE and HMC-NUTS?

Yes. That is indeed the case for FOCE(I) (and LaplaceI) but for the NUTS sampler aka BayesMCMC the likelihood contribution of the prior is applied as you would expect.

6

I see, thank you for the confirmation!

7

I do have a further question regard the same model. Is there anyway that I could fix the off-diagonal elements of OMEGA matrix to 0?
I tried to use Diagonal() to force it as diagonal matrix both in @param and the initials, but it doesn’t run through.

Thank you!
Anqi

8

I think the only way to do that would be to specify each of the diagonal elements separately, i.e. something like

@random begin
  ηCL  ~ InverseGamma(...)
  ηV1  ~ InverseGamma(...)
  ηV2  ~ InverseGamma(...)
end

I chose InverseGamma here since it’s equivalent to a scalar InverseWishart but it might be better to use the Gamma since it has support at zero.

1 Like
9

This should be fixed on the new release @anqipan

10

@andreasnoack By the way for this issue, I had a confusion in Bayesian tutorial that:

Since the mode of InverseWishart is Ψ/(ν + p + 1) , the mode of our prior is diagm(fill(0.2, 3)) which corresponds to an inter-subject variability peaking at 20%.

In this, the IIV is 20% or sqrt(0.2)*100 % because the OMEGA is in the matrix.

Thank for your help.

11

Thanks. That is indeed correct. I should probably have been fill(0.2^2, 3). We’ll have it fixed.

1 Like
12

@andreasnoack In this case, could the code write that, the InverseGamma() distribution is for Omega or for the ETA?

And because I didn’t set the matrix for omega, then the IIV = 20℅ was represented as in SD = 0.2 rather than 0.2^2 in code block?

Thank you again for help!

param begin
    # Mode at [2.0, 0.2, 0.8, 2.0]
    θ ~ Constrained(
      MvNormal(
        [2.0, 0.2, 0.8, 2.0],
        Diagonal(ones(4))
      ),
      lower = zeros(4),
      upper = fill(10.0, 4),
      init  = [2.0, 0.2, 0.8, 2.0])

    # Mean at 0.2 and positive density at 0.0
    ωKa ~ InverseGamma(1.0, 0.2)
    ωCL ~ InverseGamma(1.0, 0.2)
    ωVc ~ InverseGamma(1.0, 0.2)

    # Mean at 0.5 and positive density at 0.0
    σ ~ Gamma(1.0, 0.5)
  end

  @random begin
    ηKa ~ Normal(0.0, ωKa)
    ηCL ~ Normal(0.0, ωCL)
    ηVc ~ Normal(0.0, ωVc)
  end
13

Dear Pumas team. if I want to impose a diagonal structure on Ω then should I use InverseGama rather than Gamma for each element of Ω because it’s equivalent to InverseWishart as mentioned above by @andreasnoack?

And here, should I use 0.2 (IIV = 20% for instance) rather than 0.2^2 in the InverseGamma because it’s not in the matrix structure?

Thank you very much for helping!

   ωKa ~ InverseGamma(1.0, 0.2)
    ωCL ~ InverseGamma(1.0, 0.2)
    ωVc ~ InverseGamma(1.0, 0.2)

Then: 
@random begin
    ηKa ~ Normal(0.0, ωKa)
    ηCL ~ Normal(0.0, ωCL)
    ηVc ~ Normal(0.0, ωVc)
14

If you want to impose a diagonal structure on Ω then, indeed, you should define each scalar component the way you do it here.

Also, you are right that the InverseGamma is the scalar version of the InverseWishart.

However, that doesn’t mean that you should use that distribution. Actually, there are good reasons to avoid it so using Gamma with shape parameter equal to one might be a better choice since it has positive density at zero. In contrast, the InverseGamma forces the ω away from zero which might not always be reasonable.

Finally, if you decide to use the InverseGamma anyway then you’d have to choose the parameters differently. You can always check the mode and mean easily in Pumas (or Julia) with

julia> mean(InverseGamma(1.0, 0.2))
Inf

julia> mode(InverseGamma(1.0, 0.2))
0.1

I usually, look up the relationship between the parameters and the mean/variance on Wikipedia. Then you pick a mean value and a shape parameter value and then solve for the scale.

julia> mean(InverseGamma(2.0, 0.2))
0.2

julia> mode(InverseGamma(2.0, 0.2*3))
0.20000000000000004
1 Like
15

I am new to Pumas but I would like to point out that while using NUTS in a MCMC sampler a much faster and efficient prior would be a LKJ (2009) prior:

This is just a benchmark for random sampling using Distributions.jl:

using Distributions
using LinearAlgebra

function rand_invwishart()
           d = InverseWishart(6, diagm(fill(0.2, 3)))
           return rand(d)
end

function rand_lkj()
           d = LKJ(3, 1)
           return rand(d) * diagm(fill(0.2, 3))
end

function rand_lkj_corr()
           d = LKJCholesky(3, 1)
           A = rand(d)
           return A.L * A.U * diagm(fill(0.2, 3))
end

julia> @btime rand_invwishart()
@btime  1.718 μs (16 allocations: 1.38 KiB)
3×3 Matrix{Float64}:
 0.11615    0.0392431   0.0272846
 0.0392431  0.0425069   0.00645005
 0.0272846  0.00645005  0.0218819

julia> @btime rand_lkj()
  797.143 ns (15 allocations: 1.50 KiB)
3×3 Matrix{Float64}:
  0.2       0.046439   -0.150757
  0.046439  0.2         0.0195246
 -0.150757  0.0195246   0.2

julia> @btime rand_lkj_corr()
  529.849 ns (8 allocations: 944 bytes)
3×3 Matrix{Float64}:
  0.2        -0.0523386  0.136233
 -0.0523386   0.2        0.0646473
  0.136233    0.0646473  0.2

This is also a great paper that deals with the several options to specifying priors for covariance/correlation matrices in a Bayesian setting: New Prairie Press - Conference on Applied Statistics in Agriculture: BAYESIAN INFERENCE FOR A COVARIANCE MATRIX

Lewandowski, Daniel, Dorota Kurowicka, and Harry Joe. “Generating random correlation matrices based on vines and extended onion method.” Journal of multivariate analysis 100.9 (2009): 1989-2001.

One small caveat is that LKJCorr inside a Turing models will have to await for this issue to solve: Feature request: Add LKJCholesky · Issue #1629 · TuringLang/Turing.jl · GitHub

1 Like
16

@storopoli the LKJ prior is currently supported in the released Pumas. Here is an example of its usage:

theopmodel_bayes = @model begin
  @param begin
    θ ∈ VectorDomain(5,
      lower=[0.1,0.008,0.0004,0.1,0.0001],
      upper=[5,0.5,0.09,5,1.5],
      init=[1.9,0.0781,0.0463,1.5,0.4]
    )
    Ω  ~ LKJ(1, 2.0)
    σ² ~ Uniform(0.0, 2*0.388)
  end

  @random begin
    η ~ MvNormal(Ω)
  end

  @covariates SEX WT

  @pre begin
    Ka = (SEX == 1 ? θ[1] : θ[4]) + η[1]
    K  = θ[2]
    CL = θ[3]*(WT/70)^θ[5]
    Vc  = CL/K
    SC = Vc/(WT/70)
  end

  @vars begin
    conc = Central / SC
  end

  @dynamics Depots1Central1

  @derived begin
    dv ~ @. Normal(conc, sqrt(σ²))
  end
end

For those new to Pumas/Julia, to know what each parameter passed to LKJ above means, write ?LKJ in the Pumas/Julia REPL then press enter.

3 Likes