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)

# 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

info <- Ovariable(
  "info",
  dependencies = data.frame(Name = "jsp"),
  formula = function(...) {
    return(
      cbind(
        unique(
          jsp@output[c("Iter","Country","Gender","Ages")]
        ),
        Result = 1
      )
    )
  }
)

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(...) {
    expoRaw$Exposure <- "Other"
    dx.expo.child$Exposure <- "Breast feeding"
    out <- combine(expoRaw, dx.expo.child)
    out <- unkeep(out, sources=TRUE, prevresults=TRUE)
    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) * info + 1E-6 # Add minuscle amount to avoid log(0).
# iinfo brings back Gender, Ages, Country

objects.latest("Op_en6007", code_name="diagnostics")

objects.latest("Op_en2031", code_name = "initiate") # Default ERF but needs adjustment.
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.
objects.latest("Op_en2261", code_name = "dose")
dose <- unkeep(dose * 1, sources = TRUE, prevresults = TRUE)

######################################################### POPULATION

population <- Ovariable(
  "population", 
  data = data.frame(
    Country = rep(c("DK",  "EST", "FI",  "SWE"), each = 4),
    Gender = rep(c("Male","Female"), each=2, times=2),
    Ages = rep(c("18-45",">45"),4),
    Result = rep(c(1200000, 250000, 1350000, 2000000), each = 4)
  )
)

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.


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

RR <- EvalOutput(RR)
RR$Exposure <- NULL

# Merging of Responses cannot be done before the disease burdens or cases have been calculated
# and summed up. Otherwise the Responses will double match on Ops.
#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)
summary(BoDpaf,fun="length")
BoDcase <- EvalOutput(BoDcase)
summary(BoDcase,fun="length")
BoD <- EvalOutput(BoD)
oprint(summary(BoD))

#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