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]]
Sys.time()
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

# For development and testing, use smaller ovariables
openv.setN(200) 
jsp <- jsp[as.numeric(as.character(jsp$Iter)) <= 200 , ]

conc <- Ovariable(
  "conc",
  dependencies = data.frame(
    Name = c(
      "conc_pcddf", # [[EU-kalat]]
      "conc_vit" # [[Concentrations of beneficial nutrients in fish]]
    ), 
    Ident = c("Op_en3104/initiate", "Op_en1838/initiate")
  ),
  formula = function(...) {
    colnames(conc_pcddf@output)[colnames(conc_pcddf@output)=="Compound"] <- "Exposure_agent"
    conc_pcddf$Fish <- as.factor(conc_pcddf$Fish) # This is character vector for some reason.
    levels(conc_vit$Exposure_agent)[levels(conc_vit$Exposure_agent) == "D_vitamin"] <- "Vitamin D"
    conc_vit <- conc_vit / 100 # From /100g to /g 
    result(conc_vit)[conc_vit$Exposure_agent == "Vitamin D"] <- result(conc_vit)[conc_vit$Exposure_agent == "Vitamin D"] * 1000 # from mg to ug
    out <- combine(conc_pcddf, conc_vit)
    levels(out$Fish)[levels(out$Fish)=="Baltic herring"] <- "Herring"
    result(out) <- pmax(0,result(out))
    return(out)
  }
)
conc@meta = c(conc@meta, list(units = "Vitamin D: ug/g; DHA, EPA, and omega3: mg/g; PCDD/F, PCB: pg/g f.w."))

#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", "bgexposure", "info"),
    Ident = c(NA, "Op_en7749/initiate", NA, 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 + bgexposure
    expoRaw <- oapply(expoRaw, cols = "Fish", FUN = sum) * info
    
    return(expoRaw)
  }
)

## Background-exposure to vitamin D and omega-3
addexposure <- Ovariable(
  "addexposure",
  ddata = "Op_en7748", # [[Benefit-risk assessment of Baltic herring and salmon intake]]
  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)  + 1E-6
    
    # Make the background exposure uncertain rather than an index.
    out <- CollapseMarginal(out, cols = "Background", fun = "sample")
    
    return(out)
  }
)

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

info <- EvalOutput(info)

exposure <- Ovariable(
  "exposure",
  dependencies = data.frame(
    Name = c(
      "expoRaw",
      "dx.expo.child",
      "info"
    ),
    Ident = c(
      NA,
      "Op_en7797/initiate", # [[Infant's dioxin exposure]]
      NA
    )
  ),
  formula = function(...) {
    expoRaw$Exposure <- NA
    dx.expo.child$Exposure <- "To child"
    out <- combine(expoRaw, dx.expo.child)
    out <- unkeep(out, sources=TRUE, prevresults=TRUE)
    return(out * info + 1E-6)
  }
)

# 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 = c("jsp","exposure")),
  formula = function(...) {
    out <- merge(
      unique(jsp@output[c("Gender", "Ages", "Country", "Iter")]),
      unique(exposure@output["Exposure"])
    )
    out$Result <- ifelse(
      out$Exposure == "To child",
      ifelse(
        out$Gender == "Female" & out$Ages == "18-45", 
        0.1, # Probability of birth during a year.
        0
      ),
      1
    )
    out$Result[is.na(out$Result)] <- 1
    return(Ovariable(output=out, marginal = !colnames(out) %in% "Result")) # Marginal setting does not work
  }
)
Sys.time()
frexposed <- EvalOutput(frexposed)

ERFchoice <- Ovariable(
  "ERFchoice",
  ddata = "Op_en7748", # [[Benefit-risk assessment of Baltic herring and salmon intake]]
  subset = "Exposure-response functions of interest"
)

objects.latest("Op_en2031", code_name = "initiate") # Default ERF but needs adjustment.

ERF <- EvalOutput(ERF)
Sys.time()
ERF <- ERF[ , !colnames(ERF@output) %in% c("Age","Sex","Source")]
levels(ERF$Exposure)[levels(ERF$Exposure) %in% c(
  "Ingestion etc. (as it was in Seveso) as log(TCDD serum concentration+1) in fat",
  "Maternal intake through placenta"
)] <- "To child"

threshold <- EvalOutput(threshold)
threshold <- threshold[ , !colnames(threshold@output) %in% c("Age","Sex","Source")]
levels(threshold$Exposure)[levels(threshold$Exposure) %in% c(
  "Ingestion etc. (as it was in Seveso) as log(TCDD serum concentration+1) in fat",
  "Maternal intake through placenta"
)] <- "To child"

cat("Exposure-responses of interest.\n")
oprint(ERFchoice@output)

# Default dose but adjustments are needed before use.
objects.latest("Op_en2261", code_name = "dose")
dose <- unkeep(EvalOutput(dose), sources = TRUE, prevresults = TRUE)
Sys.time()
######################################################### 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)
  )
)

mc2dparam<- list(
  N2 = 100, # Number of iterations in the new Iter
  run2d = TRUE, # Should the mc2d function be used or not?
  newmarginals = c("Gender", "Ages", "Country"), # Names of columns that are non-marginals but should be sampled enough to become marginals
  method = "bootstrap", # which method to use for 2D Monte Carlo? Currently bootsrap is the only option.
  fun = mean # Function for aggregating the first Iter dimension.
)

objects.latest("Op_en2261", code_name="RR") # Default RR but has to be adjusted before use.
RR <- EvalOutput(RR)
RR@marginal <- RR@marginal | colnames(RR@output) %in% mc2dparam$newmarginals

objects.latest("Op_en2261", code_name="casesabs") # Has to be adjusted before use.
casesabs <- EvalOutput(casesabs)
casesabs@marginal <- casesabs@marginal | colnames(casesabs@output) %in% mc2dparam$newmarginals

disabilityweight <- Ovariable(
  "disabilityweight",
  ddata = "Op_en7748",
  subset = "DALYs of responses"
)

duration <- 1

objects.latest("Op_en7422", code_name="BoD")

BoD <- EvalOutput(BoD, verbose = TRUE) # Fetches all objects

#BoDcase <- unkeep(EvalOutput(BoDcase), sources = TRUE) # Removes NA indices
# To avoid double counting, these responses must be removed.
#BoDcase <- BoDcase[!BoDcase$Resp %in% c("Stroke", "Heart (CHD)") , ]
#
#BoD <- EvalOutput(BoD) # Re-evaluate BoD without double counting

objects.latest("Op_en6007", code_name="diagnostics")
Sys.time()
oprint(showLoctable())
oprint(showind())

rm(wiki_username)

objects.store(list=ls())
cat("Ovariables", ls(), "stored.\n")
Sys.time()

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

objects.latest("Op_en7748", code_name="hia") #  [[Benefit-risk assessment of Baltic herring and salmon intake]]

cat("Units:", conc@meta$units, "\n")
oprint(summary(conc,fun=c("length","mean","sd")))
cat("Units: g/day\n")
oprint(summary(exposure, marginal=c("Exposure", "Exposure_agent", "Country","Background"), fun=c("length","mean","sd")))
oprint(summary(exposure[is.na(exposure$Exposure),], marginal=c("Exposure_agent", "Country","Background"), fun=c("length","mean","sd")))
oprint(summary(RR,fun=c("length","mean","sd")))
oprint(summary(casesrr,fun=c("length","mean","sd")))
oprint(summary(casesabs,fun=c("length","mean","sd")))
oprint(summary(BoDpaf,fun=c("length","mean","sd"))[-5])
oprint(summary(BoDcase,fun=c("length","mean","sd"))[-4])
oprint(summary(BoD,fun=c("length","mean","sd"))[-5])

ggplot(conc@output, aes(x=concResult, colour=Exposure_agent))+stat_ecdf()+
  facet_wrap(~Fish)+scale_x_log10()

ggplot(conc@output[conc$Fish %in% c("Herring", "Salmon"), ], aes(x=concResult, colour=Exposure_agent))+stat_ecdf()+
  facet_wrap(~Fish)+scale_x_log10()

ggplot(amount@output, aes(x=amountResult, colour=Country))+stat_ecdf()+
  facet_wrap(~Fish, scales="free")+scale_x_log10()

ggplot(exposure@output[is.na(exposure$Exposure) & result(exposure)>0.01,], aes(x = exposureResult, colour=Country))+
  stat_ecdf() + scale_x_log10() + theme_gray(base_size = 18)+
  facet_wrap(~ Exposure_agent, scales="free")+labs(title = "Exposure (DX: pg/d, O3: mg/d, VD: ug/d)")

ggplot(exposure@output[!is.na(exposure$Exposure) & result(exposure)>0.01,], aes(x = exposureResult, colour=Country))+
  stat_ecdf() + scale_x_log10() + theme_gray(base_size = 24)+
  facet_wrap(~ Exposure_agent, scales="free")+labs(title = "Exposure breast feeding (DX: log(pg/g) fat in child)")

ggplot(BoD@output, aes(x=Response, weight=BoDResult/openv$N, fill=Response))+geom_bar()+coord_flip()

ggplot(BoD@output, aes(x=BoDResult, colour=Response))+stat_ecdf()+
  scale_x_log10()

ggplot(ERF@output, aes(x=ERFResult, colour=Exposure_agent))+stat_ecdf()+
  facet_wrap(~ Response, scales = "free_x")

ggplot(casesabs@output, aes(x=casesabsResult, colour=Response))+stat_ecdf()+
  scale_x_log10()

ggplot(casesrr@output, aes(x=casesrrResult, colour=Response))+stat_ecdf()+
  scale_x_log10()

ggplot(BoDcase@output, aes(x=BoDcaseResult, colour=Response))+stat_ecdf()+
  scale_x_log10()

ggplot(BoDpaf@output, aes(x=BoDpafResult, colour=Response))+stat_ecdf()+
  scale_x_log10()

odag()

#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