HPV disease model estimation (beginner issue)

Specific Domains Probabilistic programming
, ,
1

`Hi, I 'm trying to estimate HPV model parameters with Turing but I’m not sure about how to determine prior parameters especially with σ , so when I run the model it seems because of the bad selection of σ last parameter takes high values (it has to be between 0-1 but it seems chain stars near 1000.0 values for the last parameter).Lastly I want to compose my observation vector with two prevalence values(one for females and one for males) and two incidence values , so I guess I had to determine two different σ values, but I m not sure about σ has to be univariate or multivariate distribution.
Sorry for beginner code but I do not have advisor in my university , this will be my first attempt for feedback. I’m really stuck with that issue and I don’t even know where I lack knowledge any help will make my day.
Thanks.
Mevlüt.

using DifferentialEquations, LinearAlgebra, ModelingToolkit, Plots, Parameters, Turing, Distributions
using MCMCChains, StatsPlots
using Random

Random.seed!(14000);


using Logging
Logging.disable_logging(Logging.Warn)


const J = 3
const K = 5
const TF = 20
const sriskgr = range(0,stop=1,length=J)
const agegr = range(0,stop=1,length=K)
const timefinale = range(0,stop=1,length=TF)
const NF = 1500000
const NM = 1500000
const Rs = [6/10,3/10,1/10]
const Ra = [3/10,1/4,1/5,3/20,1/10]
const prv = 0.043


function init_nfemale(js,ka)
  J = length(js)
  K = length(ka)
  Nfemale = zeros(J,K)
  rsa = Rs*transpose(Ra)
  Nfemale = NF*rsa
  Nfemale
end

const Nfemale = init_nfemale(sriskgr,agegr)


function init_nmale(js,ka)
  J = length(js)
  K = length(ka)
  Nmale = zeros(J,K)
  rsa = Rs*transpose(Ra)
  Nmale = NM*rsa
  Nmale
end

const Nmale = init_nmale(sriskgr,agegr)

function init_smat(js,ka)
  J = length(js)
  L = length(js)
  K = length(ka)
  M = length(ka)
  cstar = 2

  actc = zeros(J,L,K,M)
  actcmult = zeros(J,K)
  actc[1,1,:,:] = I(K)
  actc[2,2,:,:] = I(K)
  actc[3,3,:,:] = I(K)
  actcnew = mapslices(sum, actc, dims = [2,4])
  actcnew2 = dropdims(actcnew; dims=2)
  actcnew3 = dropdims(actcnew2; dims=3)
  actcmult =  cstar * actcnew3

  actcmult
end

const actcmult = init_smat(sriskgr,agegr)


function ODEHPV(du,u,p,t)
  DIMM, DIP, DCN, DCC, DIPm, PCI, PCCN, PCC, BTRM, BTRF  = p





  @inbounds for j = 1, i in 1:J
    du[i,j,1] = u[i,j,5]/DIMM - u[i,j,1]*(BTRF*actcmult[i,j]*u[i,j,7]/Nmale[i,j]) - 0.1*u[i,j,1] + 0.1*Rs[i]*(u[5,j,2]+u[5,j,3]+u[5,j,4]+u[5,j,5])
    du[i,j,2] = u[i,j,1]*(BTRF*actcmult[i,j]*u[i,j,7]/Nmale[i,j]) - u[i,j,2]/DIP - 0.1*u[i,j,2]
    du[i,j,3] = (1-PCI)*u[i,j,2]/DIP - u[i,j,3]/DCN - 0.1*u[i,j,3]
    du[i,j,4] = (1-PCCN)*u[i,j,3]/DCN - u[i,j,4]/DCC - 0.1*u[i,j,4]
    du[i,j,5] = PCI*u[i,j,2]/DIP + PCCN*u[i,j,3]/DCN + PCC*u[i,j,4]/DCC - u[i,j,5]/DIMM - 0.1*u[i,j,5]
    du[i,j,6] = u[i,j,8]/DIMM - u[i,j,6]*(BTRM*actcmult[i,j]*u[i,j,2]/Nfemale[i,j]) - 0.1*u[i,j,6] + 0.1*Rs[i]*(u[5,j,7]+u[5,j,8])
    du[i,j,7] = u[i,j,6]*(BTRM*actcmult[i,j]*u[i,j,2]/Nfemale[i,j]) - u[i,j,7]/DIPm - 0.1*u[i,j,7]
    du[i,j,8] = u[i,j,7]/DIPm - u[i,j,8]/DIMM - 0.1*u[i,j,8]

  end

  @inbounds for j in 2:K, i in 1:J
    du[i,j,1] = u[i,j,5]/DIMM - u[i,j,1]*(BTRF*actcmult[i,j]*u[i,j,7]/Nmale[i,j]) - 0.1*u[i,j,1] + 0.1*u[i,j-1,1]
    du[i,j,2] = u[i,j,1]*(BTRF*actcmult[i,j]*u[i,j,7]/Nmale[i,j]) - u[i,j,2]/DIP - 0.1*u[i,j,2] + 0.1*u[i,j-1,2]
    du[i,j,3] = (1-PCI)*u[i,j,2]/DIP - u[i,j,3]/DCN - 0.1*u[i,j,3] + 0.1*u[i,j-1,3]
    du[i,j,4] = (1-PCCN)*u[i,j,3]/DCN - u[i,j,4]/DCC - 0.1*u[i,j,4] + 0.1*u[i,j-1,4]
    du[i,j,5] = PCI*u[i,j,2]/DIP + PCCN*u[i,j,3]/DCN + PCC*u[i,j,4]/DCC - u[i,j,5]/DIMM - 0.1*u[i,j,5] + 0.1*u[i,j-1,5]
    du[i,j,6] = u[i,j,8]/DIMM - u[i,j,6]*(BTRM*actcmult[i,j]*u[i,j,2]/Nfemale[i,j]) - 0.1*u[i,j,6] + 0.1*u[i,j-1,6]
    du[i,j,7] = u[i,j,6]*(BTRM*actcmult[i,j]*u[i,j,2]/Nfemale[i,j]) - u[i,j,7]/DIPm - 0.1*u[i,j,7] + 0.1*u[i,j-1,7]
    du[i,j,8] = u[i,j,7]/DIPm - u[i,j,8]/DIMM - 0.1*u[i,j,8] + 0.1*u[i,j-1,8]
  end

end

function init_mevotoar(js,ka)
  J = length(js)
  K = length(ka)
  u = zeros(J, K, 8)
  Rs = [6/10,3/10,1/10]
  Ra = [3/10,1/4,1/5,3/20,1/10]
  rsa = Rs*transpose(Ra)
  u[:,:,1] = (1-prv)* (NF * rsa)
  u[:,:,2] = prv* (NF * rsa)
  u[:,:,3] = zeros(3,5)
  u[:,:,4] = zeros(3,5)
  u[:,:,5] = zeros(3,5)
  u[:,:,6] = (1-prv)* (NM * rsa)
  u[:,:,7] = prv* (NM * rsa)
  u[:,:,8] = zeros(3,5)
  u
end


p = [8.0,6.0,8.0,5.0,5.0,0.88,0.65,0.6,0.1,0.1]

u0 = init_mevotoar(sriskgr,agegr)

tspan = (1.0,20.0)

prob = ODEProblem(ODEHPV,u0,tspan,p)

odep = modelingtoolkitize(prob)

jac = eval(ModelingToolkit.generate_jacobian(odep)[2])

f = ODEFunction(ODEHPV, jac=jac)

prob_jac = ODEProblem(f,u0,tspan,p)


serodataf = fill(0.043,(1,TF))
serodatam = fill(0.043,(1,TF))
datacnstage = fill(1430.0,(1,TF))
dataccstage = fill(85.0,(1,TF))

odedatadef = vcat(serodataf,serodatam,datacnstage,dataccstage)

odedatarndm = Array(odedatadef) + 0.01 * randn(size(Array(odedatadef)))


Turing.setadbackend(:forwarddiff)

@model function fitlv2(data,prob_jac) # data should be a Vector
    σ ~ InverseGamma(2, 3)
    DIMM ~ Uniform(4.0,12.0)
    DIP ~ Uniform(4.0,8.0)
    DCN ~ Uniform(6.0,10.0)
    DCC ~ Uniform(4.0,6.0)
    DIPm ~ Uniform(2.0,8.0)
    PCI ~ Beta(2,2)
    PCCN ~ Beta(2,2)
    PCC ~ Beta(2,2)
    BTRM ~ Beta(2,5)
    BTRF ~ Beta(2,5)


    p = [DIMM,DIP,DCN,DCC,DIPm,PCI,PCCN,PCC,BTRM,BTRF]
    prob = remake(prob_jac, p=p)
    predictednonesum = solve(prob,AutoTsit5(Rosenbrock23()), jac = true, saveat = 1.0)
    predicted0 = Array(predictednonesum)

    sumpredicted = mapslices(sum,predicted0, dims = 1)
    sumpredicted1 = dropdims(sumpredicted; dims=1)

    predInfectedFemale = Array(sumpredicted1[:,2,:])
    predInfectedMale = Array(sumpredicted1[:,7,:])

    predCNStagef = Array(sumpredicted1[:,3,:])
    predCCStagef = Array(sumpredicted1[:,4,:])

    pred_infectedfemale = mapslices(sum,predInfectedFemale, dims = 1)


    pred_infectedmale = mapslices(sum,predInfectedMale, dims = 1)


    predSeroPrevelancef = similar(pred_infectedfemale)
    predSeroPrevelancem = similar(pred_infectedmale)

    predSeroPrevelancef = pred_infectedfemale/NF
    predSeroPrevelancem = pred_infectedmale/NM

    predsumCNStagef = mapslices(sum, predCNStagef, dims = 1)

    predsumCCStagef = mapslices(sum, predCCStagef, dims = 1)

    predictedfnl = vcat(predSeroPrevelancef,predSeroPrevelancem,predsumCNStagef,predsumCCStagef)



    for i = 1:length(predictedfnl)
        data[i] ~ Normal(predictedfnl[i], σ) # predicted[i][2] is the data for y - a scalar, so we use Normal instead of MvNormal
    end
end

model2 = fitlv2(odedatarndm,prob_jac)


Threads.nthreads()


chain2 = sample(model2, NUTS(.65), MCMCThreads(), 5000, 3, progress=false)
2

Transform your values so that [-infty,infty] maps to [0,1]. You can do something like sigma = 1/(1+exp(-_sigma), and then you do parameter estimation on _sigma and transform inside of the model. In the end, it’ll estimate the correct _sigma and you then just have to transform it to interpret the results.

1 Like
3

Thanks a lot :relieved: