Health impact assessment

This page is a knowledge crystal of subtype method. The page identifier is Op_en2261
Moderator:Nobody (see all) Click here to sign up.
Upload data {{#opasnet_base_link:Op_en2261}}

<section begin=glossary />

Health impact assessment is an assessment method that is used to estimate the health impacts of a particular event or policy. In Europe, it is most widely used in UK, Finland, and the Netherlands.

<section end=glossary />

Question

How to calculate health impacts based on information about exposure, population, disease, and exposure-response function?

Answer

For simple calculations, you can use the concept of attributable fraction. This is presented here. For more complex and comprehensive methods, you may want to consider these:

• Multistage model
• Life table

An example model run by the model below [1].

Number of iterations: + Show code- Hide code # name:exposuresource|description:Which exposure data do you want to use?|type:selection|options:'Op_en5918';Exposures in Finland| library(OpasnetUtils) library(ggplot2) objects.latest("Op_en5917", "initiate") # [[Disease risk]] #objects.latest(exposuresource, "initiate") # [[Exposures in Finland]] objects.latest("Op_en2261", "initiate") # [[Health impact assessment]] dose, RR, totcases, AF objects.latest('Op_en5827', code_name = 'initiate') # [[ERFs of environmental pollutants]] ERF, threshold openv.setN(N) exposure <- 1 bgexposure <- 0 # background exposure frexposed <- 1 # fraction exposed BW <- 70 # body weight population <- 100000 # population size ls() disincidence <- disincidence / 100000 # # /100000py -> # /py cat("Exposure-response functions used:\n") oprint(summary(EvalOutput(ERF))) totcases <- EvalOutput(totcases) cat("Variable totcases for 100000 population and nominal 1 unit exposure.\n") oprint(summary(totcases), digits = 4) ggplot(totcases@output, aes(x = Trait, weight = result(totcases)/get("N", envir=openv), fill = Exposure_agent)) + geom_bar() + theme_grey(base_size = 24) + theme(axis.text.x = element_text(angle = 90, hjust = 1))

Inputs

If you are able to describe your data in the format similar to the tables below, you can use ready-made tools in Opasnet and things are quite straightforward. The example tables show data about radon in indoor air.

Exposure

• The table has an index Observation with four locations: Exposed fraction, Background, Exposure, and Description.

Pollutant	Exposure route	Exposure metric	Exposure parameter	Population	Exposure unit	Exposed fraction	Background	Exposure	Description
Radon	Inhalation	Annual average concentration	Population average	Finland	Bq/m3	1	5	100	Kurttio Päivi, 2006: STUK otantatutkimus 100 (95 – 105); background 5 (4 – 9)
Radon	Inhalation	Annual average concentration	Guidance value for new apartments	Finland	Bq/m3	1	0	200	STM decision 944/92 for new apartments [2]
Radon	Inhalation	Annual average concentration	Guidance value for old apartments	Finland	Bq/m3	1	0	400	STM decision 944/92 for old apartments [3]

Disease response

• The table has index Observation with two locations: Response and Description.

Disease	Response metric	Population	Unit	Response	Description
Lung cancer	Incidence	Finland	1/100000 py	38.058	2020/5307690*100000 [4]
Lung cancer	Burden of disease	Finland	DALY	14000	Olli Leino 2010: Includes trachea, bronchus, and lung cancers. [5]

Exposure-response function

Pollutant	Disease	Response metric	Exposure route	Exposure metric	Exposure unit	Threshold	ERF parameter	ERF	Description
Radon	Lung cancer	Incidence	Inhalation	Annual average concentration	Bq/m3	0	RR	1.0016	Darby 2004: 1.0016 (1.0005 – 1.0031)

Population

Population	Year	Sex	Age	Amount	Description
Finland	2010	Total	All	5307690	[6]

Rationale

Structure of the health impact assessment model. Each node is an ovariable. Arrows depict dependencies.

These are the equations you should use:

RR for exposure = EXP(LN(RR)*(Exposure Result - MAX(Exposure Background, Exposure-response function threshold)))

Attributable fraction in the whole population = Exposed fraction * (RR for exposure – 1) / (Exposed fraction *(RR for exposure – 1)+1)

Extra cases per year =Disease incidence * Population * attributable fraction

Burden of disease of exposure = Burden of disease of the disase * attributable fraction

Personal lifetime risk = Extra cases per year * life expectancy * population

Attributable fraction is (RR-1)/RR=1-1/RR if RR>1. If smaller, you must compare the other way round: control group is considered an exposure to lack of a protective agent and thus the exposure group is the reference. In this comparison, the attributable fraction of lack of protection (AF_lp) is calculated from a new rate ratio RR_lp = 1/RR and

Failed to parse (SVG (MathML can be enabled via browser plugin): Invalid response ("Math extension cannot connect to Restbase.") from server "https://wikimedia.org/api/rest_v1/":): {\displaystyle AF_{lp} = 1 - \frac{1}{RR_{lp}} = 1 - \frac{1}{1/RR} = 1 - RR}

When multiplied by the number of cases, we get the number of excess cases (that would not have occurred if the population had not been exposed to lack of protection). This comparison is symmetric and we can use either counterfactual situation as the reference just by calculating the difference the other way round, i.e. changing the sign of the value. Therefore, the number of cases avoided with exposure to a protective agent is -AF_lp = RR - 1. So, AF is calculated as 1-1/RR or RR-1 depending on whether RR>1 or not, respectively.

Calculations

Ovariables for calculating up to RR

+ Show code- Hide code # This is code Op_en2261/BW on page [[Health impact assessment]]. library(OpasnetUtils) BW <- Ovariable("BW", data = data.frame(Result = 70)) # 70 kg objects.store(BW) cat("Ovariable BW (body weight) stored. page: Op_en2261, code_name: BW.\n")

+ Show code- Hide code # This is code Op_en2261/frexposed on page [[Health impact assessment]]. library(OpasnetUtils) frexposed <- Ovariable("frexposed", data = data.frame(Result = 1)) objects.store(frexposed) cat("Ovariable frexposed stored. page: Op_en2261, code_name: frexposed.\n")

+ Show code- Hide code # This is code Op_en2261/expo_dir on page [[Health impact assessment]]. library(OpasnetUtils) expo_dir <- Ovariable( "expo_dir", dependencies=data.frame(Name=c("amount","conc","expo_background","ERFchoice"), Ident=c(NA,NA,NA,"Op_en2031/ERFchoice")), formula = function(...) { out <- amount * conc * Ovariable( output = data.frame(Exposcen = c("BAU","No exposure"), Result=c(1,0)), marginal=c(TRUE,FALSE)) colnames(out@output)[colnames(out)=="Compound"] <- "Exposure_agent" out <- OpasnetUtils::combine(out, expo_background) out <- out * oapply(ERFchoice,c("Exposure_agent","Scaling"),mean) return(out) }, unit = "?g/day" ) objects.store(expo_dir) cat("Ovariable expo_dir stored.\n")

+ Show code- Hide code # This is code Op_en2261/exposure on page [[Health impact assessment]]. library(OpasnetUtils) exposure <- Ovariable( "exposure", dependencies = data.frame( Name = c( "expo_dir", # direct exposure, i.e. the person eats or breaths the exposure agent themself "expo_indir", # indirect exposure, i.e. the person (typically fetus or infant) is exposed via someone else (mother) "mc2d" # 2D Monte Carlo function ), Ident = c( "Op_en2261/expo_dir", # [[Health impact assessment]] "Op_en7797/expo_indir2", # [[Infant's dioxin exposure]] # expo_indir "Op_en7805/mc2d") # [[Two-dimensional Monte Carlo]] ), formula = function(...) { expo_dir$Exposure <- "Direct" if(nrow(expo_indir@output)>0) out <- OpasnetUtils::combine(expo_dir, expo_indir) else out <- expo_dir out <- mc2d(out) return(out) } ) objects.store(exposure) cat("Ovariable exposure stored. page: Op_en2261, code_name: exposure.\n")

+ Show code- Hide code # This is code Op_en2261/bgexposure on page [[Health impact assessment]]. library(OpasnetUtils) bgexposure <- Ovariable("bgexposure", data = data.frame(Result = 0)) objects.store(bgexposure) cat("Ovariable bgexposure stored. page: Op_en2261, code_name: bgexposure.\n")

+ Show code- Hide code # This is code Op_en2261/population on page [[Health impact assessment]]. library(OpasnetUtils) population <- Ovariable("population", data = data.frame(Result = 1)) objects.store(population) cat("Ovariable population stored. page: Op_en2261, code_name: population.\n")

+ Show code- Hide code # This is code Op_en2261/incidence on page [[Health impact assessment]]. library(OpasnetUtils) incidence <- Ovariable("incidence", data = data.frame(Result = 0.1)) objects.store(incidence) cat("Ovariable incidence stored. page: Op_en2261, code_name: incidence.\n")

+ Show code- Hide code # This is code Op_en2261/dose on page [[Health impact assessment]]. library(OpasnetUtils) dose <- Ovariable( "dose", # This calculates the body-weight-scaled exposure or "dose" to be used with ERFs. dependencies = data.frame( Name = c( "exposure", # Exposure to the pollutants "BW" # body weight ), Ident = c( "Op_en2261/exposure", "Op_en2261/BW" ) ), formula = function(...) { ######### Body weight scaling: In some cases, exposure is given as per body weight and in some cases as absolute amounts. # Here we add one index to define for this difference. out <- Ovariable(output=data.frame(Exposcen="Remove",Result=0),marginal=FALSE) # must contain > 0 rows for(i in unique(exposure$Scaling)) { if(i=="BW") out <- OpasnetUtils::combine(out, exposure[exposure$Scaling==i,] / BW) if(i=="Log10") out <- OpasnetUtils::combine(out, log10(exposure[exposure$Scaling==i,])) if(i=="None") out <- OpasnetUtils::combine(out, exposure[exposure$Scaling==i,]) } out <- out[out$Exposcen!="Remove",] return(out) } ) objects.store(dose) cat("Ovariable dose stored. page: Op_en2261, code_name: dose.\n")

+ Show code- Hide code # This is code Op_en2261/RR2 on page [[Health impact assessment]]. library(OpasnetUtils) RR <- Ovariable( "RR", # This calculates the relative risk at given exposure level. dependencies = data.frame( Name = c( "dose", # Exposure to the pollutants "ERF" # Exposure-response function of the pollutants or agents (RR for unit exposure) ), Ident = c( "Op_en2261/dose", # [[Health impact assessment]] "Op_en2031/ERF2" # [[Exposure-response function]] ) ), formula = function(...) { t1 <- t2 <- Ovariable() # This function cleans the ERF ovariable when it is split into parameters # Otherwise there is a risk that two ovariables of the same name are opsed and a result is taken from a wrong column. # ova is the sliced ERF ovariable # obs is the number of the parameter, currently 1 for ERF and 2 for Threshold. Numbering makes it easier to change if needed. clean <- function(ova, obs) { parameter <- c("ERF", "Threshold") ova <- ova[ova$Observation==parameter[obs] , colnames(ova@output) != "Observation"] colnames(ova@output)[colnames(ova@output)==paste0(ova@name,"Result")] <- "Result" ova@name <- character() return(ova) } #################################################################### ####### This part is about risks relative to background. # Calcualte the risk ratio to each subgroup based on the exposure in that subgroup. # 1. First take the relative risk estimates. Convert ORs to RRs by using incidence. # We need OR but not yet crucial, so let's postpone this. See [[:op_en:Converting between exposure-response parameters]] # #Then take the odds ratio estimates # OR <- ERF[ERF$ERF_parameter %in% c("OR", "OR bw") , ] # if(testforrow(OR, dose)) { # See ERFrr for explanation # out <- OR / (1 - incidence + OR*incidence) # Actual function with background incidence. # } tmp <- ERF[ERF$ER_function %in% c("RR", "OR") , ] if(nrow(tmp@output)>0) { rr <- clean(tmp, 1) threshold <- clean(tmp, 2) if("OR" %in% tmp$ER_function) warning("OR is simply assumed to equal RR, which is fair approximation with small OR and low incidence.\n") dose2 <- dose - threshold result(dose2) <- pmax(0,result(dose2)) t1 <- tryCatch(exp(log(rr) * dose2), error = function(e) return(Ovariable())) # Actual function } # 2. Then take the relative Hill estimates tmp <- ERF[ERF$ER_function %in% c("Relative Hill") , ] if(nrow(tmp@output)>0) { Imax <- clean(tmp, 1) ed50 <- clean(tmp, 2) t2 <- tryCatch(1 + (dose * Imax) / (dose + ed50), error = function(e) return(Ovariable())) # Actual function # ERF has parameter value for Imax. If Imax < 0, risk reduces. # threshold has parameter value for ED50. } if(nrow(t1@output)>0) { if(nrow(t2@output)>0) { out <- OpasnetUtils::combine(t1, t2) } else out <- t1 } else out <- t2 if (nrow(out@output) == 0) { out <- Ovariable(output=data.frame(Result = 1)) } return(out) } ) objects.store(RR) cat("Ovariable RR saved.\n")

+ Show code- Hide code # This is code Op_en2261/RR on page [[Health impact assessment]]. library(OpasnetUtils) RR <- Ovariable( "RR", # This calculates the relative risk at given exposure level. dependencies = data.frame( Name = c( "dose", # Exposure to the pollutants "ERF", # Exposure-response function of the pollutants or agents (RR for unit exposure) "threshold" # exposure level below which the agent has no impact. # "incidence", # This is only needed for OR and omitted otherwise. # "mc2d" # Function to run two-dimensional Monte Carlo # This is now done in PAF ), Ident = c( "Op_en2261/dose", # [[Health impact assessment]] "Op_en2031/initiate", # [[Exposure-response function]] "Op_en2031/initiate" # [[Exposure-response function]] # "Op_en2261/incidence", # [[Health impact assessment]] # "Op_en7805/mc2d" # [[Two-dimensional Monte Carlo]] ) ), formula = function(...) { # Make sure that these are not marginals # ERF@marginal[colnames(ERF@output) %in% c("ERF_parameter", "Scaling")] <- FALSE # Remove redundant columns. No need to remove. If they cause problems, use CollapseMarginals. # ERF <- ERF[ , !colnames(ERF@output) %in% c( # "Source", # "Exposure_unit" # )] # threshold <- threshold[ , !colnames(threshold@output) %in% c( # "Source", # "Exposure_unit" # )] #################################################################### ####### This part is about risks relative to background. # Calcualte the risk ratio to each subgroup based on the exposure in that subgroup. # 1. First take the relative risk estimates. Convert ORs to RRs by using incidence. # We need OR but not yet crucial, so let's postpone this. See [[:op_en:Converting between exposure-response parameters]] # #Then take the odds ratio estimates # OR <- ERF[ERF$ERF_parameter %in% c("OR", "OR bw") , ] # if(testforrow(OR, dose)) { # See ERFrr for explanation # out <- OR / (1 - incidence + OR*incidence) # Actual function with background incidence. # } t1 <- ERF[ERF$ER_function %in% c("RR", "OR") , ] if("OR" %in% t1$ER_function) warning("OR is simply assumed to equal RR, which is fair approximation with small OR and low incidence.\n") tmp <- dose - threshold result(tmp) <- pmax(0,result(tmp)) t1 <- tryCatch(exp(log(t1) * tmp), error = function(e) return(Ovariable())) # Actual function # 2. Then take the relative Hill estimates t2 <- ERF[ERF$ER_function %in% c("Relative Hill") , ] t2 <- tryCatch(1 + (dose * t2) / (dose + threshold), error = function(e) return(Ovariable())) # Actual function # ERF has parameter value for Imax. If Imax < 0, risk reduces. # threshold has parameter value for ED50. out <- orbind(t1@output, t2@output) if(nrow(out)==0) { out <- data.frame(Result=1) } # This is done in PAF nowadays. # # Dilute the risk in the population if not all are exposed i.e. frexposed < 1. # out <- frexposed * (out - 1) + 1 # out <- unkeep(out, prevresults = TRUE, sources = TRUE, cols = c("Scaling", "ER_function")) # if(length(unique(out$Exposure_agent)) > 1) { # Could we just oapply everything? # out <- oapply(out, cols = c("Exposure_agent", "Exposure"), FUN = prod) # } else { # out <- unkeep(out, cols = c("Exposure_agent", "Exposure")) # } # } return(out) } ) objects.store(RR) cat("Ovariable RR saved. page: Op_en2261, code_name: RR.\n")

Ovariables for calculating PAF and BoD

These codes calculate population attributable fraction and burden of disease.

NOTE! Ovariables casesabs and casesrr used to be here but they were replaced by PAF and archived.

+ Show code- Hide code # This is code Op_en2261/sumExposcen on page [[Health impact assessment]]. library(OpasnetUtils) # sumExposcen calculates the difference between scenarios BAU and No exposure. sumExposcen <- function (out) { if ("Exposcen" %in% colnames(out@output)) { out <- out * Ovariable( output = data.frame(Exposcen = c("BAU", "No exposure"), Result = c(1, -1)), marginal = c(TRUE, FALSE) ) # Remove ERF-related indices as they are no longer needed. # out <- oapply(out, NULL, sum, c("Exposcen","Exposure","ER_function","Exposure_unit","Scaling")) } return(out) } objects.store(sumExposcen) cat("Function sumExposcen dose stored. page: Op_en2261, code_name: sumExposcen.\n")

The version below uses ERF ovariable that combines the previous ERF and threshold ovariables. This is newer and should be implemented with all exposure-response functions.

+ Show code- Hide code # This is code Op_en2261/PAF2 on page [[Health impact assessment]]. library(OpasnetUtils) PAF <- Ovariable( "PAF", # This calculates the burden of disease for both background-dependent and independent endpoints. dependencies = data.frame( Name = c( "dose", # Exposure to the pollutants "ERF", # Other ERFs than those that are relative to background. "RR", # relative risk at given exposure level. "incidence", # background incidence of the disease. It cancels out if response not dependent on this "frexposed", # fraction of population that is exposed "P_illness", # for microbial ERFs (codes 4-6), probability of illness given infection "sumExposcen", # function that calculates difference between exposure scenarios "mc2d", # Function to run two-dimensional Monte Carlo "mc2dparam" # Parameters for mc2d function. Default: don't run. ), Ident = c( "Op_en2261/dose", # [[Health impact assessment]] "Op_en2031/ERF2", # [[Exposure-response function]] "Op_en2261/RR2", # [[Health impact assessment]] "Op_en2261/incidence", # [[Health impact assessment]] "Op_en2261/frexposed", # [[Health impact assessment]] "Op_en7948/P_illness", # [[ERF of waterborne microbes]] "Op_en2261/sumExposcen",# [[Health impact assessment]] "Op_en7805/mc2d", # [[Two-dimensional Monte Carlo]] "Op_en7805/mc2dparam" # [[Two-dimensional Monte Carlo]] ) ), formula = function(...) { out1 <- out2 <- out3 <- out4 <- out5 <- out6 <- Ovariable() # This function cleans the ERF ovariable when it is split into parameters # Otherwise there is a risk that two ovariables of the same name are opsed and a result is taken from a wrong column. # ova is the sliced ERF ovariable # obs is the number of the parameter, currently 1 for ERF and 2 for Threshold. Numbering makes it easier to change if needed. clean <- function(ova, obs) { parameter <- c("ERF", "Threshold") ova <- ova[ova$Observation==parameter[obs] , colnames(ova@output) != "Observation"] colnames(ova@output)[colnames(ova@output)==paste0(ova@name,"Result")] <- "Result" ova@name <- character() return(ova) } # 1. Linear dose-responses y = k * (dose - threshold) tmp <- ERF[ERF$ER_function %in% c("UR", "CSF", "ERS") , ] if(nrow(tmp@output)>0) { k <- clean(tmp, 1) threshold <- clean(tmp, 2) dose2 <- dose - threshold result(dose2) <- pmax(0,result(dose2)) out1 <- tryCatch(k * dose2 * frexposed / incidence, error = function(e) return(NULL)) # Actual equation if(is.null(out1)) { out1 <- Ovariable() } else { colnames(out1@output)[colnames(out1@output)==paste0(out1@name,"Result")] <- "ERF_linResult" out1@name <- "ERF_lin" } } # 2. Step estimates: value is 1 below lower (threshold) and above upper (ERF), and 0 in between. tmp <- ERF[ERF$ER_function %in% c("Step", "ADI", "TDI", "TWI", "RDI", "NOAEL") , ] if(nrow(tmp@output)>0) { upper <- clean(tmp, 1) lower <- clean(tmp, 2) out2 <- tryCatch((1 - (dose >= lower) * (dose <= upper)) * frexposed / incidence, error = function(e) return(NULL)) # Actual equation if(is.null(out2)) { out2 <- Ovariable() } else { colnames(out2@output)[colnames(out2@output)==paste0(out2@name,"Result")] <- "ERF_stepResult" out2@name <- "ERF_step" } } # 3. RR-based PAF: (rr-1)/rr if rr>1; rr-1 if rr<1 if(nrow(RR@output)==1 & ncol(RR@output)==2 & result(RR)[1]==1) { out3 <- Ovariable() # No data so pass empty ovariable } else { out3 <- frexposed * (RR - 1) # RR diluted to get PAF instead of AF r <- result(out3) result(out3) <- (r > 0) * (r / (r + 1)) + (r <= 0) * r # Negative if RR < 1 colnames(out3@output)[colnames(out3@output)==paste0(out3@name,"Result")] <- "ERF_relResult" out3@name <- "ERF_rel" } # 4. beta poisson approximation: 1-(1+Dose/Param2)^-Param1 where Param1=ERF, Param2=threshold tmp <- ERF[ERF$ER_function %in% c("beta poisson approximation") , ] if(nrow(tmp@output)>0) { param1 <- clean(tmp, 1) param2 <- clean(tmp, 2) out4 <- tryCatch((1-(1+dose/param2)^(-1 * param1)) * frexposed / incidence, error = function(e) return(NULL)) # Actual equation if(is.null(out4)) { out4 <- Ovariable() } else { out4 <- out4 * P_illness # The original column Response is replaced by Illness from P_illness out4$Response <- NULL colnames(out4@output)[colnames(out4@output)=="Illness"] <- "Response" colnames(out4@output)[colnames(out4@output)==paste0(out4@name,"Result")] <- "ERF_bpaResult" out4@name <- "ERF_bpa" } } # 5. exact beta poisson: 1-exp(-(Param1/(Param1+Param2))*Dose) tmp <- ERF[ERF$ER_function %in% c("exact beta poisson") , ] if(nrow(tmp@output)>0) { param1 <- clean(tmp, 1) param2 <- clean(tmp, 2) out5 <- tryCatch((1 - exp(-1 * param1 / (param1 + param2) * dose)) * frexposed / incidence, error = function(e) return(NULL)) # Actual equation if(is.null(out5)) { out5 <- Ovariable() } else { out5 <- out5 * P_illness # The original column Response is replaced by Illness from P_illness out5$Response <- NULL colnames(out5@output)[colnames(out5@output)=="Illness"] <- "Response" colnames(out5@output)[colnames(out5@output)==paste0(out5@name,"Result")] <- "ERF_ebpResult" out5@name <- "ERF_ebp" } } # 6. exponential: 1-exp(-k*Dose) where k is rate tmp <- ERF[ERF$ER_function %in% c("exponential") , ] if(nrow(tmp@output)>0) { k <- clean(tmp, 1) out6 <- tryCatch((1 - exp(-1 * k * dose)) * frexposed / incidence, error = function(e) return(NULL)) # Actual equation if(is.null(out6)) { out6 <- Ovariable() } else { out6 <- out6 * P_illness # The original column Response is replaced by Illness from P_illness out6$Response <- NULL colnames(out6@output)[colnames(out6@output)=="Illness"] <- "Response" colnames(out6@output)[colnames(out6@output)==paste0(out6@name,"Result")] <- "ERF_expoResult" out6@name <- "ERF_expo" } } # Combine results from different ERFs. out <- Ovariable() if(nrow(out1@output)>0) out <- out1 if(nrow(out2@output)>0) {if(nrow(out@output)>0) out <- OpasnetUtils::combine(out,out2) else out <- out2} if(nrow(out3@output)>0) {if(nrow(out@output)>0) out <- OpasnetUtils::combine(out,out3) else out <- out3} if(nrow(out4@output)>0) {if(nrow(out@output)>0) out <- OpasnetUtils::combine(out,out4) else out <- out4} if(nrow(out5@output)>0) {if(nrow(out@output)>0) out <- OpasnetUtils::combine(out,out5) else out <- out5} if(nrow(out6@output)>0) {if(nrow(out@output)>0) out <- OpasnetUtils::combine(out,out6) else out <- out6} out <- mc2d(out) # Run two-dimensional Monte Carlo if needed return(sumExposcen(out)) } ) objects.store(PAF) cat("Ovariable PAF stored.\n")

The code below is old and does use threshold variable.

+ Show code- Hide code # This is code Op_en2261/PAF on page [[Health impact assessment]]. library(OpasnetUtils) PAF <- Ovariable( "PAF", # This calculates the burden of disease for both background-dependent and independent endpoints. dependencies = data.frame( Name = c( "dose", # Exposure to the pollutants "ERF", # Other ERFs than those that are relative to background. "threshold", # exposure level below which the agent has no impact. "RR", # relative risk at given exposure level. "incidence", # background incidence of the disease. It cancels out if response not dependent on this "frexposed", # fraction of population that is exposed "sumExposcen", # function that calculates difference between exposure scenarios "mc2d", # Function to run two-dimensional Monte Carlo "mc2dparam" # Parameters for mc2d function. Default: don't run. ), Ident = c( "Op_en2261/dose", # [[Health impact assessment]] "Op_en2031/initiate", # [[Exposure-response function]] "Op_en2031/initiate", # [[Exposure-response function]] "Op_en2261/RR", # [[Health impact assessment]] "Op_en2261/incidence", # [[Health impact assessment]] "Op_en2261/frexposed", # [[Health impact assessment]] "Op_en2261/sumExposcen",# [[Health impact assessment]] "Op_en7805/mc2d", # [[Two-dimensional Monte Carlo]] "Op_en7805/mc2dparam" # [[Two-dimensional Monte Carlo]] ) ), formula = function(...) { # 1. Linear dose-responses out1 <- ERF[ERF$ER_function %in% c("UR", "CSF", "ERS") , ] if(nrow(out1@output)>0) { tmp <- dose - threshold result(tmp) <- pmax(0,result(tmp)) out1 <- tryCatch(out1 * tmp * frexposed / incidence, error = function(e) return(NULL)) # Actual equation if(is.null(out1)) { out1 <- Ovariable() } else { colnames(out1@output)[colnames(out1@output)==paste0(out1@name,"Result")] <- "ERF_linResult" out1@name <- "ERF_lin" } } # 2. Step estimates: value is 1 below threshold and above ERF, and 0 in between. out2 <- ERF[ERF$ER_function %in% c("Step", "ADI", "TDI", "TWI", "RDI", "NOAEL") , ] if(nrow(out2@output)>0) { out2 <- tryCatch(1 - (dose >= threshold) * (dose <= out2) * frexposed / incidence, error = function(e) return(NULL)) # Actual equation if(is.null(out2)) { out2 <- Ovariable() } else { colnames(out2@output)[colnames(out2@output)==paste0(out2@name,"Result")] <- "ERF_stepResult" out2@name <- "ERF_step" } } # 3. RR-based PAF if(nrow(RR@output)==1 & ncol(RR@output)==2 & result(RR)[1]==1) { out3 <- Ovariable() # No data so pass empty ovariable } else { out3 <- frexposed * (RR - 1) # RR diluted to get PAF instead of AF r <- result(out3) result(out3) <- (r > 0) * (r / (r + 1)) + (r <= 0) * r # Negative if RR < 1 colnames(out3@output)[colnames(out3@output)==paste0(out3@name,"Result")] <- "ERF_relResult" out3@name <- "ERF_rel" } # 4. beta poisson approximation: 1-(1+Dose/Param2)^-Param1 where Param1=ERF, Param2=threshold out4 <- ERF[ERF$ER_function %in% c("beta poisson approximation") , ] if(nrow(out4@output)>0) { out4 <- tryCatch((1-(1+dose/threshold)^(-1 * out4)) * frexposed / incidence, error = function(e) return(NULL)) # Actual equation if(is.null(out1)) { out4 <- Ovariable() } else { colnames(out4@output)[colnames(out4@output)==paste0(out4@name,"Result")] <- "ERF_bpaResult" out4@name <- "ERF_bpa" } } # 5. exact beta poisson: 1-exp(-(Param1/(Param1+Param2))*Dose) out5 <- ERF[ERF$ER_function %in% c("exact beta poisson") , ] if(nrow(out5@output)>0) { out5 <- tryCatch((1 - exp(-1 * out5 / (out5 + threshold) * dose)) * frexposed / incidence, error = function(e) return(NULL)) # Actual equation if(is.null(out5)) { out5 <- Ovariable() } else { colnames(out5@output)[colnames(out5@output)==paste0(out5@name,"Result")] <- "ERF_ebpResult" out5@name <- "ERF_ebp" } } # 6. exponential: 1-exp(-Param1*Dose) out6 <- ERF[ERF$ER_function %in% c("exponential") , ] if(nrow(out6@output)>0) { out6 <- tryCatch((1 - exp(-1 * out6 * dose)) * frexposed / incidence, error = function(e) return(NULL)) # Actual equation if(is.null(out6)) { out6 <- Ovariable() } else { colnames(out6@output)[colnames(out6@output)==paste0(out6@name,"Result")] <- "ERF_expoResult" out6@name <- "ERF_expo" } } # Combine results from different ERFs. out <- Ovariable() if(nrow(out1@output)>0) out <- out1 if(nrow(out2@output)>0) {if(nrow(out@output)>0) out <- OpasnetUtils::combine(out,out2) else out <- out2} if(nrow(out3@output)>0) {if(nrow(out@output)>0) out <- OpasnetUtils::combine(out,out3) else out <- out3} if(nrow(out4@output)>0) {if(nrow(out@output)>0) out <- OpasnetUtils::combine(out,out4) else out <- out4} if(nrow(out5@output)>0) {if(nrow(out@output)>0) out <- OpasnetUtils::combine(out,out5) else out <- out5} if(nrow(out6@output)>0) {if(nrow(out@output)>0) out <- OpasnetUtils::combine(out,out6) else out <- out6} out <- mc2d(out) # Run two-dimensional Monte Carlo if needed return(sumExposcen(out)) } ) objects.store(PAF) cat("Ovariable PAF stored. page: Op_en2261, code_name: PAF.\n")

BoD is the total burden of selected diseases, BoDattr is the fraction attributed to selected exposure agents.

+ Show code- Hide code # This is code Op_en2261/BoD on page [[Health impact assessment]]. library(OpasnetUtils) BoD <- Ovariable( "BoD", dependencies=data.frame( Name=c( "incidence", "case_burden", "population" ), Ident=c( NA,NA,NA # There are no good default values available atm. )), formula=function(...) { out <- incidence * population * case_burden return(out) } ) objects.store(BoD) cat("Ovariable BoD stored. page: Op_en2261, code_name: BoD.\n")

+ Show code- Hide code # This is code Op_en2261/BoDattr2 on page [[Health impact assessment]]. library(OpasnetUtils) BoDattr <- Ovariable( "BoDattr", dependencies = data.frame( Name=c( "BoD", # disease burden of selected causes due to all risk factors "PAF" # population attributable fraction of given risk factors ), Ident=c( "Op_en2261/BoD", # [[Health impact assessment]] "Op_en2261/PAF2" ) ), formula = function(...) { out <- BoD * PAF return(out) } ) objects.store(BoDattr) cat("Ovariable BoDattr stored.\n")

+ Show code- Hide code # This is code Op_en2261/BoDattr on page [[Health impact assessment]]. library(OpasnetUtils) BoDattr <- Ovariable( "BoDattr", dependencies = data.frame( Name=c( "BoD", # disease burden of selected causes due to all risk factors "PAF" # population attributable fraction of given risk factors ), Ident=c( "Op_en2261/BoD", # [[Health impact assessment]] "Op_en2261/PAF" ) ), formula = function(...) { out <- BoD * PAF return(out) } ) objects.store(BoDattr) cat("Ovariable BoDattr stored. page: Op_en2261, code_name: BoDattr.\n")

Depreciated code

• The code was restructured and old code archived 13 June, 2015.
• Code AF (attributable fraction) was moved to page Population attributable fraction.
• See also Seturi: Excel file, [7]

+ Show code- Hide code library(OpasnetUtils) dummy <- 0 HIA <- Ovariable("HIA", dependencies = data.frame(Name = "dummy"), formula = function(...) { cat("This code is outdated. Instead, use Op_en2261/totcases on page Health impact assessment.\n") } ) totcases <- Ovariable("totcases", dependencies = data.frame(Name = "dummy"), formula = function(...) { cat("This code is outdated. Instead, use Op_en2261/totcases on page Health impact assessment.\n") } ) AF <- Ovariable("AF", dependencies = data.frame(Name = "dummy"), formula = function(...) { cat("This code is outdated. Instead, use Op_en6211/AF on page Population attributable fraction.\n") } ) objects.store(HIA, totcases, AF, dummy) cat("Warnings created about old method.\n")

Totcases (old version)

+ Show code- Hide code # This is code Op_en2261/totcases on page [[Health impact assessment]]. library(OpasnetUtils) totcases <- Ovariable("totcases", # This calculates the total number of cases in each population subgroup. # The cases are calculated for specific (combinations of) causes. However, these causes are NOT visible in the result. dependencies = data.frame( Name = c( "population", # Population divided into subgroups as necessary "dose", # Exposure to the pollutants "disincidence", # Incidence of the disease of interest "RR", # Relative risks for the given exposure "ERF", # Other ERFs than those that are relative to background. "threshold", # exposure level below which the agent has no impact. "frexposed" # fraction of population that is exposed ), Ident = c( "Op_en2261/population", # [[Health impact assessment]] "Op_en2261/dose", # [[Health impact assessment]] "Op_en5917/initiate", # [[Disease risk]] "Op_en2261/RR", # [[Health impact assessment]] "Op_en2031/initiate", # [[Exposure-response function]] "Op_en2031/initiate", # [[Exposure-response function]] "Op_en2261/frexposed" # [[Health impact assessment]] ) ), formula = function(...) { test <- list() ERF@marginal[colnames(ERF@output) %in% c("ERF_parameter", "Scaling")] <- FALSE # Make sure that these are not marginals ############### First look at the relative risks based on RR if(testforrow(RR, dose)) { # If an ovariable whose nrow(ova@output) == 0 # is used in Ops, it is re-EvalOutput'ed, and therefore ERFrr*dose may have rows even if ERFrr doesn't. # takeout is a vector of column names of indices that ARE in population but NOT in the disease incidence. # However, populationSource is kept because oapply does not run if there are no indices. if(class(population) == "ovariable") { takeout <- setdiff(colnames(population@output)[population@marginal], colnames(disincidence@output)[disincidence@marginal] ) if(length(takeout) > 0) {# Aggregate to larger subgroups. pop <- oapply(population, NULL, sum, takeout) } else { pop <- population } } else { takeout <- character() pop <- population } # pci is the proportion of cases across different population subroups based on differential risks and # population sizes. pci sums up to 1 for each larger subgroup found in disincidence. # See [[Population attributable fraction]]. pci <- population * RR # Divide pci by the values of the actually exposed group (discard nonexposed) # The strange Ovariable thing is needed to change the name of temp to avoid problems later. temp <- pci * Ovariable(data = data.frame(Result = 1)) if ("Exposcen" %in% colnames(temp@output)) { temp@output <- temp@output[temp@output$Exposcen == "BAU" , ] temp <- unkeep(temp, cols = "Exposcen", prevresults = TRUE, sources = TRUE) } temp <- unkeep(temp, prevresults = TRUE, sources = TRUE) if(length(takeout) > 0) temp <- oapply(temp, NULL, sum, takeout) #if(length(takeout) > 0) temp <- ooapply(temp, cols = takeout, FUN = "sum", use_plyr = TRUE) # if(length(takeout) > 0) temp <- osum(temp, cols = takeout) pci <- pci / temp temp <- NULL # The cases are divided into smaller subgroups based on weights in pci. # This is why the larger groups of population are used (pop instead of population). out1 <- disincidence * pop * unkeep(pci, prevresults = TRUE, sources = TRUE) out1 <- unkeep(out1, cols = "populationResult") # populationResult comes from pop and not from pci that actually contains # the population weighting for takeout indices. Therefore it would be confusing to leave it there. test <- c(test, out1) } out1 <- NULL ########################################################################## ############# This part is about absolute risks (i.e., risk is not affected by background rates). # Unit risk (UR), cancer slope factor (CSF), and Exposure-response slope (ERS) estimates. UR <- ERF UR@output <- UR@output[UR@output$ERF_parameter %in% c("UR", "CSF", "ERS") , ] if(testforrow(UR, dose)) { # See RR for explanation. UR <- threshold + UR * dose * frexposed # Actual equation # threshold is here interpreted as the baseline response (intercept of the line). It should be 0 for # UR and CSF but it may have meaningful values with ERS UR <- oapply(UR, NULL, sum, "Exposure_agent") UR <- population * UR test <- c(test, UR) } UR <- NULL # Step estimates: value is 1 below threshold and above ERF, and 0 in between. # frexposed cannot be used with Step because this may be used at individual and maybe at population level. Step <- ERF Step@output <- Step@output[Step@output$ERF_parameter %in% c("Step", "ADI", "TDI", "RDI", "NOAEL") , ] if(testforrow(Step, dose)) { # See RR for explanation. Step <- 1 - (dose >= threshold) * (dose <= Step) # Actual equation # Population size should be taken into account here. Otherwise different population indices may go unnoticed.(?) Step <- oapply(Step, NULL, sum, "Exposure_agent") test <- c(test, Step) } Step <- NULL ##################################################################### # Combining effects if(length(test) == 0) return(data.frame()) if(length(test) == 1) out <- test[[1]]@output if(length(test) == 2) out <- orbind(test[[1]], test[[2]]) if(length(test) == 3) out <- orbind(orbind(test[[1]], test[[2]]), test[[3]]) # Find out the right marginals for the output marginals <- character() nonmarginals <- character() for(i in 1:length(test)) { marginals <- c(marginals, colnames(test[[i]]@output)[test[[i]]@marginal]) nonmarginals <- c(nonmarginals, colnames(test[[i]]@output)[!test[[i]]@marginal]) } test <- NULL out <- out[!colnames(out) %in% c("populationSource", "populationResult")] # These are no longer needed. out <- Ovariable(output = out, marginal = colnames(out) %in% setdiff(marginals, nonmarginals)) if("Exposcen" %in% colnames(out@output)) { out <- out * Ovariable( output = data.frame(Exposcen = c("BAU", "No exposure"), Result = c(1, -1)), marginal = c(TRUE, FALSE) ) out <- oapply(out, NULL, sum, "Exposcen") } return(out) } ) objects.store(totcases) cat("Ovariable totcases saved. page: Op_en2261, code_name: totcases.\n")

NOTE! These ovariables used to utilise ooapply function, but it was archived after improved oapply.

The codes above are based on these input variables:

• Exposure-response function
• Disease risk (case-spacific data)
• Exposures in Finland (case-specific data)
• Burden of disease in Finland
• Disability weights (not used yet)
• Duration of morbidity (not used yet)
• heande:File:Impact Calculation Tool.ana (not used, but functionalities should be merged to this page)

See also related page: ISTE EBD.

See also

• Population attributable fraction
• WHO Health impact assessment tools
• WHO: HIA tools
• DALY calculator in R
• Converting between exposure-response parameters
• The effectiveness of health impact assessment. WHO 2007
• Ihmisiin kohdistuvien vaikutusten arviointi (käsikirja) (in Finnish)
• Ihmisiin kohdistuvien vaikutusten arviointi (sivusto) (in Finnish)
• An article about HIA in Finland
• Health impact assessment suggested to be released in Intarese
• Intarese Health Effects Methodology (D13 Final, July 2007)
• Health impact assessment in Wikipedia
• Impact assessment in Wikipedia
• IngentaConnect: Impact assessment and Project Appraisal ISSN 1471-5465 21: 4 (2003).
• Scott-Samuel: Health impact assessment BMJ.1996; 313: 183-184
• Assessing health impact assessment: multidisciplinary and international perspectives. J Epidemiol Community Health 2003;57:659-662
• Health Impact Assessment in Urban Settings New South Wales Public Health Bulletin 18: 9 & 10, 2007.
• Health Impact Assessment. Bulletin of the World Health Organization (BLT): 81: 6: 387-472, 2003.
• Health Impact Assessment. Community knowledge wiki
• IMPACT - International Health Impact Assessment Consortium
• Health Impact Assessment page by WHO
• ISO 31000:2009 Risk management -- Principles and guidelines
• Risk Observatory, based in the European Agency for Safety and Health at Work
• Health-EU Portal
• STM: Hyvinvointi (Well-being, in Finnish)
• John Kemm, Jayne Parry, Stephen Palmer: Health impact assessment: concepts, theory, techniques, and applications‎
• Interdepartmental Liaison Group on Risk Assessment in UK
• Life cycle assessment
• Four-step impact assessment by HSPH.
• OpasnetUtils/Drafts

---

Further reading

The text below is a description of HIA by A. Knol and B. Staatsen from RIVM. It was originally written for use in Intarese project.

Health Impact Assessment

One way to compare different policy options is by carrying out a health impact assessment (HIA). HIA is a combination of procedures, methods and instruments used for assessing the potential health impacts of certain matters. These can vary from a single environmental factor to a more complicated set of factors, for instance in an infrastructural or industrial project. For quantifying health impacts, the following steps can be distinguished (Hertz-Picciotto, 1998):

• Selection of health endpoints with sufficient proof (based on expert judgements) of a causal relationship with the risk factor
• Assessment of population exposure (combination of measurements, models and demographic data)
• Identification of exposure-response relations (relative risks, threshold values) based on (meta) analyses and epidemiological and toxicological research.
• Estimation of the (extra) number of cases with the specific health state, attributable to exposure to the risk factor. This is a function of the population distribution, exposure-response relation and base prevalence of the health state in the population.
• Computation of the total health burden, or costs to society of all risk factors (if wanted/necessary)

A common problem is that the health effects of environmental factors can vary considerably with regard to their severity, duration and magnitude. These differences hamper the comparison of policies (comparative risk assessment) or the costs of policy measures (cost effectiveness analysis). An integrated health measure, using the same denominator for all health effects, can help with interpretation and comparison of health problems and policies.

Integrated health measures

Common health measures include mortality, morbidity, healthy life expectancy, attributable burden of disease measures, and monetary valuation. Some of these measures will be further described below. All methods have several associated difficulties, such as imprecision of the population exposure assessment; uncertain shapes of the exposure-response curves for the low environmental exposure levels; insufficient (quality of) epidemiological data; extrapolation from animal to man or from occupational to the general population; generalisation of exposure-response relations from locally collected data for use on regional, national or global scale; combined effects in complex mixtures, etc.

Mortality figures The annual mortality risk or the number of deaths related to a certain (environment-related) disease can be compared with this risk or number in another region or country, or with data from another period in time. Subsequently, different policies can be compared and policies that do or do not work can be identified. Within a country, time trends can be analyzed. This method is easy to comprehend. No ethical questions are attached; everyone is treated equal. Since this method only includes mortality, it is not suitable for assessing factors with less severe consequences (morbidity). Also, it is difficult to attribute mortality to specific environmental causes.

Morbidity figures Similar to mortality figures, morbidity numbers (prevalences or incidences based on hospital admissions or doctor visits) can be used to evaluate a (population) health state. Advantages and drawbacks are comparable to those applying to using mortality figures. The use of morbidity numbers is therefore similarly limited, especially when (environmental) causes of the diseases vary.

Healthy life expectancy Using mortality tables, one can calculate the total average life expectancy for different age groups in a population, subdivided into years with good and years with less-than-good health. This measure is especially useful to review the generic health state in a country for the long term, but it doesn’t give insight into specific health effects, effects of specific policy interventions, or trends in certain subgroups.

Attributable burden of disease Health impact assessments can also be executed by calculating the attributable burden of disease. There are several ways to assess the burden of disease attributable to an (environmental) factor, such as the QALY and the DALY. Quality Adjusted Life Years, QALYs, capture both the quality and quantity elements of health in one indicator. Essentially, time spent in ill health (measured in years) is multiplied by a weight measuring the relative (un)desirability of the illness state. Thereby a number is obtained which represents the equivalent number of years with full health. QALYs are commonly used for cost-utility analysis and to appraise different forms of health care. To do that, QALYs combine life years gained as a result of these health interventions/health care programs with a judgment about the quality of these life years. Disability adjusted life years, DALYs, are comparable to QALYs in that they both combine information on quality and quantity of life. However, contrary to QALYs, DALYs give an indication of the (potential) number of healthy life years lost due to premature mortality or morbidity and are estimated for particular diseases, instead of a health state. Morbidity is weighted for the severity of the disorder.

With QALY, the focus is on assessing individual preference for different non-fatal health outcomes that might result from a specific intervention, whereas the DALY was developed primarily to compare relative burdens among different diseases and among different populations (Morrow and Bryant, 1995). DALYs are suitable for analyzing particular disorders or specific factors that influence health. Problems associated with the DALY approach include the difficulty of estimating the duration of the effects (which have hardly been studied) and the severity of a disease; and allowing for combined effects in the same individual (first you have symptoms, then you go to a hospital and then you may die). The DALY concept, which has been used in our study, will be further described in the next chapter. More information on the drawbacks of the method can be found in Chapter 6.4.

Monetary valuation Another approach to health impact assessment is monetary valuation. In this measure, money is used as a unit to express health loss or gain, thereby facilitating the comparison of policy costs and benefits. It can help policy makers in allocating limited (health care) resources and setting priorities. There are different approaches to monetary valuation such as ‘cost of illness’ and ‘willingness to pay/accept’.

The cost of illness (COI) approach estimates the material costs related to mortality and morbidity. It includes the costs for the whole society and considers loss of income, productivity and medical costs. This approach does not include immaterial costs, such as impact of disability (pain, fear) or decrease in quality of life. This could lead to an underestimation of the health costs. Furthermore, individual preferences are not considered.

The willingness to pay (WTP) approach measures how much money one would be willing to pay for improvement of a certain health state or for a reduction in health risk. The willingness to accept (WTA) approach measures how much money one wants to receive to accept an increased risk. WTP and WTA can be estimated by observing the individual’s behaviour and expenditures on related goods (revealed preference). For example, the extra amount of money people are willing to pay for safer or healthier products (e.g. cars with air bags), or the extra salary they accept for compensation of a risky occupation (De Hollander, 2004). Another similar method is contingent valuation (CV), in which people are asked directly how much money they would be willing to pay (under hypothetical circumstances) for obtaining a certain benefit (e.g. clean air or good health).

Source: Knol, A.B. en Staatsen, B.A.M. (2005). Trends in the environmental burden of disease in the Netherlands, 1980-2020. Rapport 500029001, RIVM, Bilthoven. Downloadable at http://www.rivm.nl/bibliotheek/rapporten/500029001.html