LloydChapman / VLSpatiotemporalModelling

Code for 'Inferring transmission trees to guide targeting of interventions against visceral leishmaniasis and post–kala-azar dermal leishmaniasis'
GNU General Public License v3.0
4 stars 1 forks source link

VLStm.jl observations #2

Open t-pollington opened 3 years ago

t-pollington commented 3 years ago

Hi @LloydChapman,

Please pardon the brevity:

CalcHHDists()

KnlHH()

Tim.

t-pollington commented 3 years ago

Note also in newer Julia, the following (which appears twice) as CSV.read("data_final2.csv") should be replaced as CSV.read("data_final2.csv", DataFrame).

Also join() on a DataFrame no longer works, but assuming that innerjoin() is fine (when this occurs twice) in PlotSims.jl.

t-pollington commented 3 years ago

In functions.jl what does IM_IN do on line 102?

t-pollington commented 3 years ago

In functions.jl line 100 why is 8-status assigned to Event? ie why is Event treated like a counter in a relative sense when it could just be set to different values that describe what that event was?

t-pollington commented 3 years ago

In functions.jl line 103, why is only the first person sharing the internal migration of IM_IN used? Is this because multiple rows can account for the same person?

t-pollington commented 3 years ago

In tI[i] = t+rand(MersenneTwister(),NegativeBinomial(r1,p1))+1 why not allow infection to start at exposure? Is it better to use a truncation on the NB to restrict the support to >=1 instead, than doing this +1 correction?

t-pollington commented 3 years ago

"I think there's an issue here, as coded like this their PKDL onset could be after their death. Simulate birth and death times? Otherwise I'm assuming death rate is unaffected by KA." I agree. Also natural death from non-VL causes as a function of age and sex not included. Best for those who do get infected, for risk of death to also be modelled but only for the infectious part rather than the asyx part.

t-pollington commented 3 years ago

Why no +1 on this line tP[i] = tDor[i] + rand(MersenneTwister(),NegativeBinomial(r3,p3))

t-pollington commented 3 years ago

How did you obtain r0 and p0 since these are output as vectors from the following code:

    OTparams = readdlm(joinpath(@__DIR__,"OTparams.csv"),',',Float64)
    r0 = OTparams[:,1]
    p0 = OTparams[:,2]
    # r0 = 1.369639756715030*ones(9,1)
    # p0 = 0.380308479571719*ones(9,1)

Also 1.36 and 0.38 don't match with mean() or median() of the vector, nor the 1.34 value (0.38 was correct) in the SI just after Eqn S5.

t-pollington commented 3 years ago

🐛 In tI[i] = t + rand(MersenneTwister(),NegativeBinomial(r1,p1)) + 1 I think r1 and p1 are incorrect params for the exposure to onset time. Shouldn't this be r=3 and p=1.34ish instead?

t-pollington commented 3 years ago
tR[Asx] = t .+ rand(MersenneTwister(),Geometric(p2),A[1]) .+ 1
    tI[Exp] = t .+ randSizeBiasedNegativeBinomial(r1,p1,E[1])
    tR[Exp] = tI[Exp] + rand(MersenneTwister(),NegativeBinomial(r0[1],p0[1]),E[1]) .+ 1
    tDor[infDor] = t .+ randSizeBiasedNegativeBinomial(r0[1],p0[1],sum(infDor))
    tP[infDor] = tDor[infDor] + rand(MersenneTwister(),NegativeBinomial(r3,p3),sum(infDor))
    tR[infDor] = tP[infDor] + rand(MersenneTwister(),NegativeBinomial(r5,p5),sum(infDor)) .+ 1
    tR[infRec] = t .+ randSizeBiasedNegativeBinomial(r0[1],p0[1],sum(infRec))
    tP[Dor] = t .+ randSizeBiasedNegativeBinomial(r3,p3,D[1])
    tR[Dor] = tP[Dor] + rand(MersenneTwister(),NegativeBinomial(r5,p5),D[1]) .+ 1
    tR[PKDL] = t .+ randSizeBiasedNegativeBinomial(r5,p5,P[1])

Some missing +1. Is this OK?

t-pollington commented 3 years ago

🐛 Need to consider when at some times the tR[infRec] = t .+ randSizeBiasedNegativeBinomial(r0[1],p0[1],sum(infRec)) will evaluate to 0×1 Array{Int64,2} if randSizeBiasedNegativeBinomial(...,n=0). Better for these functions not to be performed and alternative control structures to guide the program away from using them, as 0×1 Array{Int64,2} may lead to instability when given as later function input.

t-pollington commented 3 years ago

🐛 As model can only work with integer months, randSizeBiasedNegativeBinomial() shouldn't accept non-integer r.

t-pollington commented 3 years ago

According to the parameters in the SI which are named as in Table S2 on p8, there is not a parameter for the success probability of the KA IP NB distribution. It is loaded in RunSims.jl as p1_str = "p1_DELTA_HGHR_PRESX_ASX_INFCTS_AllParas.csv" and is called p1 in the code, not to be confused with p1 in Table S2.

t-pollington commented 3 years ago

Why isn't r0[1] and p0[1] used when medians could be taken of these and used instead? It appears r0[1] & p0[1] are used initially for t=w=12 since r0[n] & p0[n] relate to the estimate in the study year n. Is this right? If so why is it that r0[9]=1 exactly and p0[9]=0.5 exactly, and why doesn't SI p5 section D mention the year-specific estimates for r0 & p0? Note these appear as r1 and p1 in the SI text.

t-pollington commented 3 years ago

Need to understand Event codes better and why they are used to advance individuals through Status using Status + Event.

t-pollington commented 3 years ago

Why is recovery from KA, and, AsxInf -> AsxDor grouped together for event 3?

elseif (inf[i] & infRec[i]) | (Asx[i] & AsxDor[i]) # add 3 for direct recovery from KA, or for progression from asx infection to dormant infection

since those who have recovered from KA are not infectious. Is it because AsxDor aren't infectious either but could become inf if they progress to PKDL?

t-pollington commented 3 years ago

I don't understand:

Event[i] = Status[IM_OUT[IM_IN.==i][1]] - Status[i]
t-pollington commented 3 years ago
elseif !(t<tIM[i]) & !(t>tD[i]) & !(t>tEM[i]) # individual didn't internally migrate in and hasn't already died/emigrated
            Event[i] = 1 # add 1 for birth, progression to KA, progression from KA to dormant infection, progression to PKDL, and recovery from PKDL
        end

Doesn't the above also add +1 if you are in Status=7 ie Rec and will move you to Status=8 ie death?