R script: Incorporating uncertainty in aggregate burden of disease measures.

Posted on June 10, 2011

1



The following script can be used for normal calculation of Disability Adjusted Life Years (DALYs) using data from life tables, but in addition enables the incorporation of the uncertainty of each parameter in the DALY calculation using Markov Chain Monte Carlo simulation. More detail can be found in the following paper (for which this script is also available as the Online Appendix):

Frank de Vocht, James Higgerson, Kathryn Oliver, Arpana Verma. Incorporating uncertainty in aggregate burden of
disease measures: an example of DALYs-averted by a smoking cessation campaign in the UK. J Epidemiol Community Health 2010; DOI   10.1136/jech.2010.119842.

The actual data from the script below are for the specific example described in the accompanying paper. It is straightforward to amend the script for specific research questions. After adjustment of the formulae this script can also be used for other aggregate burden of disease measures such as QALYs (quality adjusted life years), HALYs (health adjusted life years), DALEs (disability adjusted life expectancy), HALEs (health adjusted life expectancy)3 and PIMs (population impact measures)

***** Population burden of disease (DALYs) prior to intervention ***
for (i in 1:n){
#agegroup 1 (25-29)
rate.m.1[i]<-runif(1,0.0008,0.0009)
pop.m.1[i]<-runif(1,2068400,2068600)
no.deaths.m.1[i]<-rate.m.1[i]*pop.m.1a[i]
life.m.1[i]<-runif(1,49.48,53.32)
D.prev.m.1[i]<-rnorm(1,0.01,0.10)
dur.case.m.1[i]<-rnorm(1,10,3)
rate.f.1[i]<-runif(1,0.0003,0.0004)
pop.f.1[i]<-runif(1,2005700,2005900)
no.deaths.f.1[i]<-rate.f.1[i]*pop.f.1a[i]
life.f.1[i]<-runif(1,53.35,57.28)
D.prev.f.1[i]<-rnorm(1,0.01,0.12)
dur.case.f.1[i]<-rnorm(1,10,3)
disability.weight[i]<-rnorm(1,0.29,0.3606)
YLL.1[i]<-(no.deaths.m.1[i]*life.m.1[i])+(no.deaths.f.1[i]*life.f.1[i])
YLD.1[i]<-
((D.prev.m.1[i]/100)*pop.m.1[i]*dur.case.m.1[i]*disability.weight[i])+((D.prev.f.1[i]/100)*pop.f.1[i
]*dur.case.f.1[i]*disability.weight[i])
DALY.1[i]<-YLL.1[i]+YLD.1[i]
# similar calculations for agegroups 2(30-34) to 8 (60-64)
DALY.burden[i]<-
DALY.1[i]+DALY.2[i]+DALY.3[i]+DALY.4[i]+DALY.5[i]+DALY.6[i]+DALY.7[i]+DALY.8[i]
}
***** Population burden of disease (DALYs) after intervention ***
for (i in 1:n){
#agegroup 1 (25-29)
ces.reduction.m.1[i]<-(1-exp(rnorm(1,-0.3857,0.096)))
ces.rate.m.1[i]<-runif(1,0.0008,0.0009)
ces.pop.m.1[i]<-runif(1,2068400,2068600)
ces.no.deaths.m.1[i]<-rate.m.1[i]*pop.m.1[i]*(1-ces.reduction.m.1[i])
ces.life.m.1[i]<-runif(1,49.48,53.32)
ces.D.prev.m.1[i]<-rnorm(1,0.01,0.10)
ces.dur.case.m.1[i]<-rnorm(1,15,4)
ces.reduction.f.1[i]<-(1-exp(rnorm(1,-0.3857,0.096)))
ces.rate.f.1[i]<-runif(1,0.0003,0.0004)
ces.pop.f.1[i]<-runif(1,2005700,2005900)
ces.no.deaths.f.1[i]<-rate.f.1[i]*pop.f.1[i]*(1-ces.reduction.f.1[i])
ces.life.f.1[i]<-runif(1,53.35,57.28)
ces.D.prev.f.1[i]<-rnorm(1,0.01,0.12)
ces.dur.case.f.1[i]<-rnorm(1,15,4)
ces.disability.weight[i]<-rnorm(1,0.10,0.3606)
ces.YLL.1 [i]<-(ces.no.deaths.m.1[i]*ces.life.m.1[i])+(ces.no.deaths.f.1[i]*ces.life.f.1[i])
ces.YLD.1[i]<-
((ces.D.prev.m.1a[i]/100)*ces.pop.m.1a[i]*ces.dur.case.m.1a[i]*ces.disability.weight[i])+((ces.
D.prev.f.1a[i]/100)*ces.pop.f.1a[i]*ces.dur.case.f.1a[i]*ces.disability.weight[i])
De Vocht, Higgerson, Oliver, Verma. J Epi Com Health, Appendix
ces.DALY.1a[i]<-ces.YLL.1a[i]+ces.YLD.1a[i]
# similar calculations for agegroups 2(30-34) to 8 (60-64)
ces.DALY.burden[i]<-
ces.DALY.1[i]+ces.DALY.2[i]+ces.DALY.3[i]+ces.DALY.4[i]+ces.DALY.5[i]+ces.DALY.6[i]+ces
.DALY.7[i]+ces.DALY.8[i]
*** calculated of DALYs-averted ***
# 100% succesrate
DALY.averted[i]<-DALY.burden[i]-ces.DALY.burden[i]
# succesrate based on previous interventions (Table 1)
DALY.averted.real[i]<-((exp(rnorm(1,log(1.53),0.206)))-1)*(DALY.burden[i]-
ces.DALY.burden[i])
}

Advertisements