Difference between revisions of "Health impact assessment"

#	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))

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")

# 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")

# 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")

# This is code Op_en2261/exposure on page [[Health impact assessment]].
library(OpasnetUtils)

exposure <- Ovariable("exposure", data = data.frame(Result = 1))

objects.store(exposure)
cat("Ovariable exposure stored. page: Op_en2261, code_name: exposure.\n")

# 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")

# 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")

# 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")

# 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
      "bgexposure", # Background exposure (a level you use for comparison)
      "BW" # body weight
    ),
    Ident = c(
      "Op_en2261/exposure",
      "Op_en2261/bgexposure",
      "Op_en2261/BW"
    )
  ),
  formula = function(...) {
    
    ########### Create a single ovariable with exposure and background exposure.
    
    temp <- Ovariable( # Create alternative scenario with background exposure bgexposure.
      output = data.frame(Exposcen = c("BAU", "No exposure"), Result = c(1, 0)), 
      marginal = c(TRUE, FALSE)
    )
    
    out <- temp * exposure + (1 - temp) * bgexposure # Adds exposure and background to respective scenarios
    
    ######### 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.
    out2 <- out / BW
    out3 <- log10(out)
    
    out$Scaling <- "None"
    out2$Scaling <- "BW"
    out3$Scaling <- "Log10"
    
    out <- combine(out, out2, out3)
    
    return(out)
  }
)

objects.store(dose)
cat("Ovariable dose stored. page: Op_en2261, code_name: dose.\n")

# This is code Op_en2261/RR on page [[Health impact assessment]].
library(OpasnetUtils)

# Do we need testforrow any more, as RR does not use it?
testforrow <- function(x, y) {
  if (nrow(x@output) == 0 | nrow(y@output) == 0) return(FALSE)
  commons <- intersect(colnames(x@output), colnames(y@output))
  commons <- commons[!grepl("Result$", commons)]
  for (i in commons) {
    if (!any(unique(x@output[[i]]) %in% unique(y@output[[i]]))) return(FALSE)
  }
  return(TRUE)
}

# This code produces non-unique index combinations. They multiply when 
# formula is run and probably cause bias in results. Check!
RR <- Ovariable(
  "RR", # 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(
      "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.
      "frexposed", # fraction of population that is exposed
      "incidence", # This is only needed for OR and omitted otherwise.
      "mc2d"  # Function to run two-dimensional Monte Carlo
    ),
    Ident = c(
      "Op_en2261/dose",     # [[Health impact assessment]]
      "Op_en2031/initiate", # [[Exposure-response function]]
      "Op_en2031/initiate", # [[Exposure-response function]]
      "Op_en2261/frexposed",# [[Health impact assessment]]
      "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
    ERF <- ERF[ , !colnames(ERF@output) %in% c(
      "Source",
      "Exposure_unit"
    )]
    threshold <- threshold[ , !colnames(threshold@output) %in% c(
      "Source",
      "Exposure_unit"
    )]
    out <- NULL
    dose <- dose * ERF # Do the merge now to avoid redundant rows
    result(dose) <- dose$doseResult
    dose <- unkeep(dose, sources=TRUE)
    
    ####################################################################
    ####### This part is about risks relative to background.
    # Calcualte the risk ratio to each subgroup based on the exposure in that subgroup.
    # Combine pollutant-specific RRs by multiplying. For description, see [[Exposure-response function]].
    
    #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.
    #		}
    
    temp <- dose[dose$ER_function %in% c("RR", "OR") , ]
    if(nrow(temp@output)>0) {
      
      temp <- exp(log(ERF) * (temp - threshold)) # Actual function
      
      result(temp)[temp$doseResult < temp$thresholdResult] <- 1 # RR is 1 below threshold
      out <- temp
    }

    # Then take the relative Hill estimates
    
    temp <- dose[dose$ER_function %in% c("Relative Hill") , ]
    if(nrow(temp@output)>0) {

      temp <- 1 + (temp * ERF) / (temp + threshold)  # Actual function
      
      # ERF has parameter value for Imax. If Imax < 0, risk reduces.
      # threshold has parameter value for ED50.
      if(is.null(out)) out <- temp else out <- combine(out, temp)
    }
    temp <- NULL
    
    if(is.null(out)) {
      out <- oapply(dose, cols = c("Iter","Exposure_agent","Scaling","Exposure"), FUN=sum)
      result(out) <- 1
    } else {
      
      # 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"))
      }
    }
    if(mc2dparam$run2d) out <- mc2d(out) # Run two-dimensional Monte Carlo
    return(out)
  }
)

objects.store(RR, testforrow)
cat("Ovariables RR, testforrow saved. page: Op_en2261, code_name: RR.\n")

# 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")

# This is code Op_en2261/casesabs on page [[Health impact assessment]].
library(OpasnetUtils)

casesabs <- Ovariable(
  "casesabs", # This calculates the burden of disease for background-independent endpoints.
  dependencies = data.frame(
    Name = c(
      "population", # Population divided into subgroups as necessary
      "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.
      "frexposed", # fraction of population that is exposed
      "sumExposcen", # function that calculates difference between exposure scenarios
      "mc2d"  # Function to run two-dimensional Monte Carlo
    ), 
    Ident = c(
      "Op_en2261/population", # [[Health impact assessment]]
      "Op_en2261/dose",       # [[Health impact assessment]]
      "Op_en2031/initiate",   # [[Exposure-response function]]
      "Op_en2031/initiate",   # [[Exposure-response function]]
      "Op_en2261/frexposed",  # [[Health impact assessment]]
      "Op_en2261/sumExposcen",# [[Health impact assessment]]
      "Op_en7805/mc2d"        # [[Two-dimensional Monte Carlo]]
    )
  ),
  formula = function(...) {
    out <- NULL
    dose2 <- dose * ERF # Do the merge first
    result(dose2) <- dose2$doseResult

    temp <- dose2[dose2$ER_function %in% c("UR", "CSF", "ERS"), ]
    
    # Dose could be simplified with combine.
    if(nrow(temp@output)>0) {
      out <- (threshold + temp * ERF * frexposed) * population # 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
      # But this interpretation is problematic. Threshold should have same units as exposure, not response.
    }

    # Step estimates: value is 1 below threshold and above ERF, and 0 in between.
    
    temp <- dose2[dose2$ER_function %in% c("Step", "ADI", "TDI", "RDI", "NOAEL") , ]
    if(nrow(temp@output)>0) {
      
      temp <- (1 - (temp >= threshold) * (temp <= ERF)) * frexposed * population # Actual equation
      
      if(is.null(out)) out <- temp else out <- combine(out, temp)
    }
    out <- oapply(out, NULL, sum, c("Exposure_agent", "Exposure", "ER_function", "Scaling"))
    if(mc2dparam$run2d) out <- mc2d(out) # Run two-dimensional Monte Carlo
    
    return(sumExposcen(out))
  }
)

objects.store(casesabs)
cat("Ovariable casesabs stored. page: Op_en2261, code_name: casesabs.\n")

# This is code Op_en2261/casesrr on page [[Health impact assessment]].
library(OpasnetUtils)

casesrr <- Ovariable(
  "casesrr", # This calculates the burden of disease for endpoints using RR.
  dependencies = data.frame(
    Name = c(
      "population", # Population divided into subgroups as necessary
      "RR", # Relative risks for the given exposure
      "incidence", # incidence of responses
      "sumExposcen" # function that calculates difference between exposure scenarios
    ), 
    Ident = c(
      "Op_en2261/population", # [[Health impact assessment]]
      "Op_en2261/RR",         # [[Health impact assessment]]
      "Op_en5917/initiate",   # [[Disease risk]]
      "Op_en2261/sumExposcen" # [[Health impact assessment]]
    )
  ),
  formula = function(...) {
    AF <- (RR > 1) * (1 - 1/RR) + (RR <= 1) * (RR - 1)
    out <- population * incidence * AF

    return(sumExposcen(out))
  }
)

objects.store(casesrr)
cat("Ovariable casesrr stored. page: Op_en2261, code_name: casesrr.\n")

# 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")