Difference between revisions of "Benefit-risk assessment of Baltic herring and salmon intake"

# This is code Op_en7748/hia on page [[Benefit-risk assessment of Baltic herring and salmon intake]]
#Code does not work yet.
library(OpasnetUtils)
library(ggplot2)

# amount does not know where to find jsp.
objects.latest("Op_en7749", code_name = "surveyjsp") # Uses jsp directly from survey data.

openv.setN(max(as.numeric(as.character(jsp$Iter)))) # Adjust N to data size

conc <- Ovariable(
  "conc",
  dependencies = data.frame(
    Name = c(
      "concentration", # [[EU-kalat]]
      "vit" # [[Concentrations of beneficial nutrients in fish]]
    ), 
    Ident = c("Op_en3104/initiate", "Op_en1838/initiate")
  ),
  formula = function(...) {
    colnames(concentration@output)[colnames(concentration@output)=="Compound"] <- "Exposure_agent"
    concentration$Fish <- as.factor(concentration$Fish) # This is character vector for some reason.
    out <- combine(concentration, vit)
    levels(out$Fish)[levels(out$Fish)=="Baltic herring"] <- "Herring"
    result(out) <- pmax(0,result(out))
    return(out)
  }
)

#objects.latest("Op_en4004", code_name = "initiate") # [[Mercury concentrations in fish in Finland]]
# ^ This code does not exist yet. Use EU-kalat as example.

# sumitem is needed by both expoRaw and bgexposure
sumitem <- function(
  ova, # ovariable that has locations to sum
  cond, # index column that contains the locations to sum
  condvalue, # vector of locations to sum
  sumvalue # location to be given to the rows with the sums
) {
  d <- ova
  d@output <- d@output[d@output[[cond]] %in% condvalue , ]
  d <- oapply(d, cols = cond, FUN = sum)
  d@output[[cond]] <- sumvalue
  ova@output <- orbind(ova, d)
  return(ova)
}

# Exposure with fish-related exposure and additional exposure from other sources.
# Does NOT contain infant's exposure from mother.

expoRaw <- Ovariable(
  "expoRaw",
  dependencies = data.frame(
    Name = c("conc", "amount", "addexposure"),
    Ident = c(NA, "Op_en7749/initiate", NA)
  ),
  formula = function(...) {
    amount <- unkeep(
      amount, 
      sources=TRUE, 
      prevresults=TRUE, 
      cols=c(
        "Eat.fish", 
        "How.often.fish",
        "Eat.salmon",
        "assumpUnit"
      )
    )
    
    # Make a summary for TEQ (Omega3 exists already)
    conc <- sumitem(conc, "Exposure_agent", c("PCDDF", "PCB"), "TEQ")
    levels(conc$Fish)[levels(conc$Fish) == "Baltic herring"] <- "Herring"
    expoRaw <- amount * conc
    
    return(expoRaw)
  }
)

## Background-exposure to vitamin D and omega-3
addexposure <- Ovariable(
  "addexposure",
  ddata = "Op_en7748",
  subset = "Background exposure"
)

# Should the background be specific for gender and country? At the moment it is.
bgexposure <- Ovariable(
  "bgexposure",
  dependencies = data.frame(Name="addexposure"),
  formula = function(...) {
    out <- addexposure
    # Empty values ("") in indices must be replaced by NA so that Ops works correctly.
    levels(out$Gender)[levels(out$Gender) == ""] <- NA
    levels(out$Country)[levels(out$Country) == ""] <- NA
    levels(out$Exposure_agent)[levels(out$Exposure_agent) == ""] <- NA
    out@output <- fillna(out@output, c("Country", "Gender", "Exposure_agent"))
    
    out <- sumitem(out, "Exposure_agent", c("PCDDF","PCB"), "TEQ")
    out <- sumitem(out, "Exposure_agent", c("EPA", "DHA"), "Omega3")
    out <- unkeep(out, prevresults = TRUE, sources = TRUE)
    
    # Make the background exposure uncertain rather than an index.
    out <- CollapseMarginal(out, cols = "Background", fun = "sample")
    
    return(out)
  }
)

exposure <- Ovariable(
  "exposure",
  dependencies = data.frame(
    Name = c(
      "expoRaw",
      "dx.expo.child"
    ),
    Ident = c(
      NA,
      "Op_en7797/initiate" # [[Infant's dioxin exposure]]
    )
  ),
  formula = function(...) {
    dx.expo.child <- unkeep(dx.expo.child, sources=TRUE, prevresults=TRUE)
    expoRaw$Exposure <- "Other"
    dx.expo.child$Exposure <- "Breast feeding"
    colnames(expoRaw@output)[colnames(expoRaw@output)=="expoRawResult"] <- "Result"
    colnames(dx.expo.child@output)[colnames(dx.expo.child@output)=="dx.expo.childResult"] <- "Result"
    
    out <- orbind(expoRaw, dx.expo.child) # Use orbind rather than combine because of non-marginals.
    out <- Ovariable(
      output=out, 
      marginal= colnames(out) %in% 
        c(
          colnames(expoRaw@output)[expoRaw@marginal],
          colnames(dx.expo.child@output)[dx.expo.child@marginal]
        ))
    return(out)
  }
)

# Limit the infant health responses to 10 % of females at age 18-45 a
# (assuming 10 % probability to give birth during a year)

frexposed <- Ovariable(
  "frexposed",
  dependencies = data.frame(Name = "exposure"),
  formula = function(...) {
    out <- unique(exposure@output[c("Gender", "Ages", "Exposure")])
    out$Result <- ifelse(
      out$Exposure == "Breast feeding",
      ifelse(
        out$Gender == "Female" & out$Ages == "18-45", 
        0.1, # Probability of birth during a year.
        0
      ),
      1
    )
    return(out)
  }
)

exposure <- EvalOutput(exposure)
objects.latest("Op_en6007", code_name="diagnostics")
oprint(showLoctable())
ovashapetest(exposure)
binoptest(concentration, vit)

objects.latest("Op_en2031", code_name = "initiate")
ERF <- EvalOutput(ERF)
threshold <- EvalOutput(threshold)
ERF <- ERF[ , !colnames(ERF@output) %in% c("Exposure", "Source")]
threshold <- threshold[ , !colnames(threshold@output) %in% c("Exposure", "Source")]

# Default dose but adjustments are needed before use.
# Log10 scaling changes some doses to -Inf and rows drop out in e.g. summary.
# Should we add some minuscule exposure or does that cause other problems?
objects.latest("Op_en2261", code_name = "dose")
dose <- unkeep(dose, sources = TRUE, prevresults = TRUE)

######################################################### VÄESTÖ

# Tarkastellaan totcases-laskennassa aluksi yksilöriskiä, ja vasta lopussa kerrotaan yksilötulokset yksilöiden edustamien ryhmien koolla.

population <- Ovariable(
  "population", 
  data = data.frame(
    Country = c("DK",  "EST", "FI",  "SWE"),
    Result = c(2400000, 500000, 2700000, 4000000) # Number for one gender
  )
)

dummy <- 0

BoDt <- Ovariable(
  "BoDt",
  dependencies = data.frame(Name="dummy"),
  formula = function(...) {
    BoDt <- opbase.data("Op_en5917", subset = "IHME_Goherr")# [[Disease risk]]
    BoDt <- BoDt[BoDt[[1]]=="DALY" , -1]
    # Measure has a special character in name, so use position.
    
    # levels(BoDt$Country)
    # [1] "Denmark" "Estonia" "Finland" "Sweden" 
    levels(BoDt$Country) <- c("DK", "EST", "FI", "SWE")
    
    # levels(BoDt$Age)
    #[1] "1 to 4"   "10 to 14" "15 to 19" "20 to 24" "25 to 29" "30 to 34" "35 to 39" "40 to 44"
    #[9] "45 to 49" "5 to 9"   "50 to 54" "55 to 59" "60 to 64" "65 to 69" "70 to 74" "75 to 79"
    #[17] "80 plus"  "All Ages"
    BoDt$Ages <- BoDt$Age
    levels(BoDt$Ages) <-  c(NA,NA,rep("18-45",6),">45","18-45",rep(">45",7),NA)
    BoDt <- BoDt[!is.na(BoDt$Ages) , ]
    
    #levels(BoDt$Response)
    #[1] "Hemorrhagic stroke"     "Ischemic heart disease" "Ischemic stroke"       
    #[4] "Permanent caries"      
    levels(BoDt$Response) <- c("Stroke","Ischaemic heart disease","Stroke","Tooth defect")
    
    BoDt$Result <- paste(BoDt$Lower, BoDt$Result, BoDt$Upper, sep=":")
    BoDt<- Ovariable(
      data = BoDt[!colnames(BoDt) %in% c("Upper", "Lower", "Year")]
    )
    BoDt <- EvalOutput(BoDt)
    BoDt$Source <- NULL
    BoDt@marginal <- colnames(BoDt@output) != "Result"
    BoDt <- oapply(BoDt, cols="Age", FUN=sum)
    return(BoDt)
  }
)
BoDt <- EvalOutput(BoDt)

if(FALSE) {
  # WHY ARE THERE RESPONSES THAT DO NOT MATCH WITH EXPOSURES?
  
  ggplot(concentration@output, aes(x=concentrationResult, Colour=Compound))+geom_density()+
    facet_wrap(~Fish, scales="free_y")+scale_x_log10()
  
  ggplot(amount@output, aes(x=amountResult, Colour=Country))+stat_ecdf()+
    facet_wrap(~Fish, scales="free_y")+scale_x_log10()
  
  ggplot(totcases@output, aes(x = Response, weight = totcasesResult, fill = Gender))+
    geom_bar()+coord_flip() + theme_gray(base_size = 24)
  
  ggplot(exposure@output, aes(x = Result, colour = Gender))+
    geom_density() + scale_x_log10() + theme_gray(base_size = 24)+
    facet_wrap(~ Country)+labs(title = "Dioxin exposure (pg/d)")
  
  
  ### Muutetaan exposure yksikköön g /d kaikkien Exposure_agentien osalta.
  skaala <- Ovariable(
    output = data.frame(
      Exposure_unit = "g /d", 
      Exposure_agent = c("Vitamin_D", "EPA", "DHA", "Omega3", "PCDDF", "PCB", "TEQ", "MeHg"), 
      Result = c(1E-6, 1E-3, 1E-3, 1E-3, 1E-12, 1E-12, 1E-12, 1E-6)),
    marginal = c(FALSE, TRUE, FALSE)
  )
  
  exposure <- exposure * skaala
  
  concentration <- unkeep(concentration, prevresults = TRUE, sources = TRUE) * skaala
}

############# This part is about absolute risks (i.e., risk is not affected by background rates).
# Unit risk is essentially an individual function and cannot be used for BoD.??    
# Unit risk (UR), cancer slope factor (CSF), and Exposure-response slope (ERS) estimates.
# threshold is not used for RR. And the interpretation is strange for ERS (use as X intercept instead).

# These have NA and are fillna'ed wrong
# UR concSource, expoRawSource, BWSource.
# Step concSource, expoRawSource
# out  Country, populationSource, frexposedSource

# Why is Country a marginal in casesabs? It should not be.
# RR does not have Country as index, so it has disappeared. Why?
# After casesabs, there should be some sorting to choose the best ER_functions and Responses for the purpose.

# 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)
    )
    out <- oapply(out, NULL, sum, "Exposcen")
  }
  return(out)
}

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
    ), 
    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]]
    )
  ),
  formula = function(...) {
    out <- NULL
    UR <- ERF[ERF$ER_function %in% c("UR", "CSF", "ERS") , ]
    
    # Dose could be simplified with combine.
    if(testforrow(UR, dose)) {
      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
      # But this interpretation is problematic. Threshold should have same units as exposure, not response.
      
      UR <- oapply(UR, NULL, sum, "Exposure_agent")
      out <- population * 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 function is essentially an individual function and cannot be used for BoD calculations???
    
    Step <- ERF[ERF$ER_function %in% c("Step", "ADI", "TDI", "RDI", "NOAEL") , ]
    if(testforrow(Step, dose)) { # See RR for explanation.
      
      Step <- (1 - (dose >= threshold) * (dose <= Step)) * frexposed # Actual equation
      
      Step <- oapply(Step, NULL, sum, "Exposure_agent")
      Step <- population * Step
      
      if(is.null(out)) out <- Step else out <- combine(out, Step)
    }
    Step <- NULL
    
    return(sumExposcen(out))
  }
)

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(...) {
    
    out <- population * incidence * (RR - 1)/RR

    return(sumExposcen(out))
  }
)

BoDpaf <- Ovariable(
  "BoDpaf", # This calculates the burden of disease for endpoints using PAF.
  dependencies = data.frame(
    Name = c(
      "population", # Population relative to the scale (often 1 or 1M) of BoDt
      "BoDt", # Total burden of disease of the studied responses per a defined time interval
      "RR", # Relative risks for the given exposure
      "sumExposcen" # function that calculates difference between exposure scenarios
    ), 
    Ident = c(
      "Op_en2261/population", # [[Health impact assessment]]
      "Op_en5917/BoDt",       # [[Disease risk]]
      "Op_en2261/RR",        # [[Health impact assessment]]
      "Op_en2261/sumExposcen" # [[Health impact assessment]]
    )
  ),
  formula = function(...) {
    
    out <- BoDt * population * (RR-1)/RR

    return(sumExposcen(out))
    
  }
)

BoDcase <- Ovariable(
  "BoDcase", # This calculates the burden of disease for endpoints using numbers of cases.
  dependencies = data.frame(
    Name = c(
      "casesrr", # Number of cases from relative endpoints with RR.
      "casesabs", # Number of cases from absolute ERFs.
      "disabilityweight", # Disability weights for each response.
      "duration" # Duration of a response case.
    ), 
    Ident = c(
      "Op_en2261/casesrr", # [[Health impact assessment]]
      "Op_en2261/casesabs", # [[Health impact assessment]]
      "Op_en2307/disabilityweight", # [[Disability weights]]
      "Op_en2307/duration" # [[Disability weights]]
    )
  ),
  formula = function(...) {
    
    out <- combine(casesrr, casesabs)
    out <- out * disabilityweight * duration
    return(out)
    
  }
)

# Note! There may double counting if there are responses that have both BoDt and incidence.
# Some warning against double counting should be created.
BoD <- Ovariable(
  "BoD", # This calculates the total burden of disease.
  dependencies = data.frame(
    Name = c(
      "BoDpaf", # Burden of disease as calculated from PAF
      "BoDcase" # Burden of disease as calculated from number of cases
    ), 
    Ident = c(
      "Op_en2261/BoDpaf", # [[Health impact assessment]]
      "Op_en2261/BoDcase" # [[Health impact assessment]]
    )
  ),
  formula = function(...) {
    
    out <- combine(BoDpaf, BoDcase)
    doublecount <- intersect(unique(BoDpaf$Response), unique(BoDcase$Response))
    if(length(doublecount>0)) warning(paste(
      "The following response(s) are probably double counted from both disease burdens and incidences:",
      paste(doublecount, collapse = ", ")
    ))
    return(out)
  }
)

objects.latest("Op_en2261", code_name="RR") # Default RR but has to be adjusted before use.

# Why do RR have interations for each combination of Gender and Ages although they are not marginals?

RR <- EvalOutput(RR)
#levels(RR$Response)
#[1] "CHD2 mortality"                   "Coronary heart disease mortality"
#[3] "Stroke mortality"                 "CHD arrythmia mortality"         
levels(RR$Response) <- c(NA,"Ischaemic heart disease", "Stroke",NA)
RR <- RR[!is.na(RR$Response) , ]

disabilityweight <- 0.1
duration <- 1

BoDpaf <- EvalOutput(BoDpaf)
oprint(summary(BoDpaf,fun="length"))
BoDcase <- EvalOutput(BoDcase)
oprint(summary(BoDcase,fun="length"))
Bod <- EvalOutput(BoD)
oprint(summary(BoD,fun="length"))

oprint(showLoctable())

odag()

#ggplot(BoD@output, aes(x = Response, weight = totcasesResult, fill = Gender))+geom_bar(position="dodge")

#This is code Op_en7748/ on page [[Benefit-risk assessment of Baltic herring and salmon intake]]

library(OpasnetUtils)
library(ggplot2)
library(reshape2)

objects.latest("Op_en7748", code_name = "bayes") #: pcd.pred, ans.pred, mu.pred
objects.latest("Op_en3104", code_name = "preprocess") # [[EU-kalat]]: eu

pl <- melt(pcd.pred)
mul <- melt(mu.pred)
ql <- melt(ans.pred)

############## PCDDF concentrations plotted

ggplot(eu@output, aes(x = euResult, colour = Fish))+geom_density(size = 1) +
  scale_x_log10() +
  facet_wrap(~ Congener, scales = "free_y") + coord_cartesian(xlim = c(1E-3,1E2))+
  labs(title = "Probability densities of PCDD/F concentrations directly from data")

ggplot(eu@output, aes(x = euResult, colour = Congener))+geom_density(size = 1) +
  scale_x_log10() +
  facet_wrap(~ Fish, scales = "free_y") + coord_cartesian(xlim = c(1E-3,1E2))+
  labs(title = "Probability densities of PCDD/F concentrations directly from data")

ggplot(pl, aes(x = value, colour = Fish))+geom_density(size = 1) + 
  scale_x_log10() +
  facet_wrap(~ Congener, scales = "free_y") + coord_cartesian(xlim = c(1E-3,1E2))+
labs(title = "Probability densities of PCDD/F concentrations from bayes model")

ggplot(pl, aes(x = value, colour = Congener))+geom_density(size = 1) + 
  scale_x_log10() +
  facet_wrap(~ Fish, scales = "free_y") + coord_cartesian(xlim = c(1E-3,1E2))+
  labs(title = "Probability densities of PCDD/F concentrations from bayes model")

ggplot(pl[pl$Fish == "Baltic herring",], aes(x = Iter, y = value,  group = Seed, colour = Seed))+
  geom_line() + scale_y_log10() +
  facet_wrap(~ Congener, scales = "free_y") +
  labs(title = "Mixing of model. Runs 1..1000, no thinning")

#  aggregate(pl["value"], by = pl[c("Congener", "Fish")], FUN = mean)

ggplot(pl[pl$Fish %in% c("Baltic herring", "Rainbow trout", "Pike-perch", "Pike", "Salmon", "Sprat") , ],
       aes(x = value, colour = Congener))+geom_density(size = 1) + scale_x_log10() +
  facet_wrap(~ Fish, scales = "free_y") + coord_cartesian(xlim = c(1E-3,1E2))+
  labs(title = "Selected fish species")

########## Questionnaire plotted

for(i in unique(ql$Question)) {
  print(ggplot(ql[ql$Question == i,], aes(x = value, fill = Gender))+
    geom_bar(position = "dodge") +
    facet_grid(Country ~ Age)+
    labs(title = paste(i, "from bayes survey model")))
}

##################### Mu plotted

ggplot(mul, aes(x = value, colour = Congener))+stat_ecdf(size = 1) + 
  scale_x_log10() +
  facet_wrap(~ Fish, scales = "free_y") + #coord_cartesian(xlim = c(1E-3,1E2))+
  labs(title = "Mu of PCDD/F concentrations from bayes model")

ggplot(mul[mul$Fish == "Baltic herring",], aes(x = Iter, y = value,  group = Seed, colour = Seed))+
  geom_line() + scale_y_log10() +
  facet_wrap(~ Congener, scales = "free_y") +
  labs(title = "Mixing of model. Runs 1..1000, no thinning")

#tef <- Ovariable("tef", ddata = "Op_en4017", subset = "TEF values")
#tef <- EvalOutput(tef)

#levels(eu$Congener) <- gsub("HCDD", "HxCDD", levels(eu$Congener))
#levels(eu$Congener) <- gsub("HCDF", "HxCDF", levels(eu$Congener))
#levels(eu$Congener) <- gsub("CoPCB", "PCB", levels(eu$Congener))
#euteq <- eu * tef