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]]

library(OpasnetUtils)
library(ggplot2)

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

#objects.latest("Op_en7749", code_name = "initiate") # [[Goherr: Fish consumption study]]
objects.latest("Op_en7749", code_name = "surveyjsp") # Uses jsp directly from survey data.

openv.setN(max(as.numeric(as.character(jsp@data$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"
    out <- combine(concentration, vit)
    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.

#objects.latest("Op_en4017", code_name = "initiate") # [[TEF]]

conc <- EvalOutput(conc)

exposure <- Ovariable(
  "exposure",
  dependencies = data.frame(
    Name = c("conc", "amount"),
    Ident = c(NA, "Op_en7749/initiate")
  ),
  formula = function(...) {
    # sumitem only works in this code
    conc <- sumitem(conc, "Exposure_agent", c("TEQdx", "TEQpcb"), "TEQ")
    levels(conc$Fish)[levels(conc$Fish) == "Baltic herring"] <- "Herring"
    exposure <- amount * conc
    return(exposure)
  }
)

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



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

# Empty values ("") in indices must be replaced by NA so that Ops works correctly.
levels(addexposure$Gender)[levels(addexposure$Gender) == ""] <- NA
levels(addexposure$Exposure_agent)[levels(addexposure$Exposure_agent) == ""] <- NA
addexposure@output <- fillna(addexposure@output, c("Gender", "Exposure_agent"))

addexposure <- sumitem(addexposure, "Exposure_agent", c("PCDDF","PCB"), "TEQ")
addexposure <- sumitem(addexposure, "Exposure_agent", c("EPA", "DHA"), "Omega3")
addexposure <- unkeep(addexposure, prevresults = TRUE, sources = TRUE)

# Make the background exposure uncertain rather than an index.
addexposure <- CollapseMarginal(addexposure, cols = "Background", fun = "sample")

exposure <- exposure + addexposure

#ii++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

objects.latest("Op_en7797", code_name="initiate") # [[Infant's dioxin exposure]]

#########################

bgexposure <- addexposure

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

#!!+++++++++++++++++++++++ THIS CODE REPLACES ALL POLLUTANT-SPECIFIC CODES:
objects.latest("Op_en2031", code_name = "initiate") # [[Exposure-response function]]

################################### MAKE A SINGLE BoD OVARIABLE FROM THE CODE BELOF UNTIL HIA.
##### THIS SHOULD BE BASED ON BOD RATHER THAN INCIDENCES
#!!++++++++++++++++++++++++++++++++++++++++++++++++++
objects.latest("Op_en5917", code_name = "initiate") # [[:op_en:Disease risk]] ovariable disincidence
#ii++++++++++++++++++++++++++++++++++++++++++++++++++

disincidence <- EvalOutput(disincidence) # TEMPORARY 13.3.

disincidence <- disincidence / 100000 # Change from 1/100000py to 1/py.
disincidence <- unkeep(disincidence, cols = c("Age", "Sex", "Population", "Unit", "Response_metric")) 
# Näitä indeksejä ei tarvita tässä arvioinnissa.
disincidence@output <- disincidence@output[disincidence@output$Response != "CHD" , ] # Tulee ikävakioidusta alta.

############################################################################
# Ikävakioidut kuolemansyyt

# Sairastuvuus (tarkemmin sanottuna [[Kuolemansyyt Suomessa]])
#!!++++++++++++++++++++++++++++++++++++++++++++++++++++
syy <- Ovariable(
  "syy", 
  data = opbase.data(
    "Op_fi4558", 
    subset = "2012", 
    include = list(Kuolemansyy = c(
      "04-21 Syövät (C00-C97)",
      "27 Iskeemiset sydäntaudit (I20-I25)",
      "29 Aivoverisuonien sairaudet (I60-I69)"
    )
    )
  )
)

#ii++++++++++++++++++++++++++++++++++++++++++++++++++++
# Syöpää ei haluta linkata ikäryhmittäin ERF:iin, vaan sitä käytetään muodostamaan diskontattu elinaikainen
# odotettu eliniän menetys jos kuolee syöpään tietyn ikäisenä. Ks. myöhemmin.

väli <- Ovariable(output = data.frame(
  Kuolemansyy = c(
    rep("27 Iskeemiset sydäntaudit (I20-I25)", 2), 
    "04-21 Syövät (C00-C97)", 
    "29 Aivoverisuonien sairaudet (I60-I69)"
  ),
  Response = c("CHD", "CHD2", "Syövät", "Stroke"), # Syövät nimetään hassusti, koska sen EI haluta linkkautuvan ERF:iin.
  # Lisäksi Tehdään kaksinkertainen sydäntaululistaus, koska on kaksi ERFiä: Mozaffarian ja Cohen.
  Result = 1
))

###################################################################################
#### HIA-ovariablet

#!!++++++++++++++++++++++++++++++++++++++++++++++++++
objects.latest("Op_en2261", code_name = "totcases") # [[:op_en:Health impact assessment]] 
# ovariable totcases and recursive ovariables dose, RR, and threshold
#ii++++++++++++++++++++++++++++++++++++++++++++++++++

totcases <- EvalOutput(totcases)

# WHY ARE THERE RESPONSES THAT DO NOT MATCH WITH EXPOSURES?

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


colnames(exposure@output)

temp <- oapply(totcases, cols = "Iter", FUN = mean)
oprint(temp@output)  

Loctable <- showLoctable()
oprint(showLoctable())
oprint(Loctable[Loctable$Index == "Age",])
##################################################################################
# Tapauskohtaiset postprosessoinnit

if(FALSE) {
  
  ### 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)
  )
  
  tiedot <- rivit * background * Ovariable(
    output = data.frame(
      silakka[c("Rivi", "Sukupuoli", "Age", "Hedelm", "Ikä", "Paino", "Silakkamäärä")],
      Result = 1
    ),
    marginal = c(FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE)
  )
  tiedot <- unkeep(tiedot, cols = "Rivi", prevresults = TRUE, sources = TRUE)
  
  exposure <- exposure * skaala * tiedot
  amount <- amount * tiedot
  concentration <- unkeep(concentration, prevresults = TRUE, sources = TRUE) * skaala
  
  indivrisk <- totcases
  
  indivrisk <- indivrisk * tiedot
  
  ################################3
  # Limit the infant health responses to 10 % of females at age 18-45 a
  # (assuming 10 % probability to give birth during a year)
  result(indivrisk) <- result(indivrisk) * ifelse(
    indivrisk$Response %in% c("Child's IQ", "Tooth defect", "Dental defect"), 
    ifelse(
      indivrisk$Gender == "Female" & indivrisk$Ages == "18-45", 
      0.1, # Probability of birth during a year.
      0
    ),
    1
  )
  ###############################
  
  #ages <- c(
  #	"0", "1-4", "5-9", "10-14", "15-19", "20-24",
  #	"25-29", "30-34", "35-39", "40-44", "45-49", "50-54", 
  #	"55-59", "60-64", "65-69", "70-74", "75-79"
  #)
  indivrisk@output$Age <- factor(indivrisk@output$Age, levels = ages, ordered = TRUE)
  
  indivrisk@output$Iter <- as.numeric(as.character(indivrisk@output$Iter))
  
  #!!+++++++++++++++++++++++++++++++++++++++++++++++++
  #  objects.store(indivrisk, pop, exposure, amount, ERF, threshold, silakka, disincidence, concentration, rivit)
  #ii+++++++++++++++++++++++++++++++++++++++++++++++++
  cat("indivrisk, pop, exposure, amount, ERF, threshold, silakka, disincidence, concentration, rivit stored. \n")
  
  ####################################3
  
  for(i in c("RR", "disincidence", "totcases", "ERF", "threshold")) {
    print(levels(get(i)@output$Response))
  }
} # if(FALSE)

######################################## TESTIALUE

totcases2 <- 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(FALSE) {
        print("RR")                           
        print(head(dose@output))                           
        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
        print("pci1")                           
        print(head(pci@output))                           
        # 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)
        print("temp")
        print(head(temp@output))                           
        pci <- pci / temp
        temp <- NULL
      } # if(FALSE)       
      # 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 * population * (RR-1)/RR #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.
      #print("out1")
      #print(head(out1@output))
      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
      print("UR")                           
      print(head(dose@output))                           
      
      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)
  }
)

totcases2 <- EvalOutput(totcases2)
unique(totcases2@output$Response)
nrow(RR@output)

oprint(head(totcases2@output))
ggplot(totcases2@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