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

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

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

objects.latest("Op_en3104", code_name = "preprocess") # [[EU-kalat]]

# Hierarchical Bayes model.

# PCDD/F concentrations in fish.
# It uses the sum of PCDD/F (Pcdsum) as the total concentration of dioxin in fish.
# Cong_j is the fraction of a congener from pcdsum.
# pcdsum is log-normally distributed. cong_j follows Dirichlet distribution.
# pcdsum depends on age of fish, fish species and catchment area, but we only have species now so other variables are omitted.
# cong_j depends on fish species.

# Fish intake in humans
# Data from data.frame survey from page [[Goherr: Fish consumption study]]
# Start with salmon questions 46:49 (amounts eaten)
# Some preprocessing should be moved away from this code.
# Model assumes that questionnaire data follows binomial distribution with parameter p and max value 9.
# p depends on country (4 groups), sex (2), age (2) that are known individually.
# We should test whether sex makes a difference or whether it can be omitted.
# Predicted p for each question, determined by the group, are produced from the bayes model.
# p is used to sample answers, which are then combined with estimates of amounts related to each answer.
# Total amount eaten of each fish per modelled individual is finally calculated.

objects.latest("Op_en7749", "preprocess")
for (i in 46:49) {
  survey[[i]] <- as.character(survey[[i]])
}
survey[,46:49][is.na(survey[,46:49])] <- ""
# This row takes data and converts it to numbers. But what should be done is to take
# each variable, remove NA, order answers meaningfully and then convert to numbers.
answer <- as.data.frame(lapply(survey[46:49], FUN = function(x) as.numeric(as.factor(x))))

survey$Row <- 1:nrow(survey)
survey$Weighting <- as.double(levels(survey$Weighting))[survey$Weighting]
survey$Ages <- ifelse(as.numeric(as.character(survey$Age)) < 46, "18-45",">45")
agel <- unique(survey$Ages)
countryl <- unique(survey$Country)
genderl <- unique(survey$Gender)

conl <- unique(eu@output$Congener)
fisl <- unique(eu@output$Fish)
fishsamples <- reshape(
  eu@output, 
  v.names = "euResult", 
  idvar = "THLcode", 
  timevar = "Congener", 
  drop = c("Matrix", "euSource"), 
  direction = "wide"
)

# Find the level of quantification for dinterval function
LOQ <- unlist(lapply(fishsamples[3:ncol(fishsamples)], FUN = function(x) min(x[x!=0])))
names(LOQ) <- conl
cong <- data.matrix(fishsamples[3:ncol(fishsamples)])
cong[cong == 0] <- 0.01 # NA # Needed for dinterval

mod <- textConnection("
  model{
    for(i in 1:S) { # s = fish sample
      for(j in 1:C) { # C = congener
#        below.LOQ[i,j] ~ dinterval(-cong[i,j], -LOQ[j])
        cong[i,j] ~ dlnorm(mu[fis[i],j], tau[fis[i],j])
      }
    }
    for(i in 1:F) { # F = fish species
      for(j in 1:C) { # C = congener
        mu[i,j] ~ dunif(-3,3)#0.001,1000)
        tau[i,j] <- pow(sigma[i,j], -2)
        sigma[i,j] ~ dunif(0, 10)
        pcd.pred[i,j] ~ dlnorm(mu[i,j], tau[i,j]) # Model prediction
      }
    }
    
    for(i in 1:U) { # U = observation in the questionnaire
      for(j in 1:Q) { # Q = question in the questionnaire
        ans[i,j] ~ dbin(p[age[i],country[i],gender[i],j], 9)
      }
    }
    for(i in 1:A) { # Age groups
      for(j in 1:B) { # Countries
        for(k in 1:H) { #Genders
          for(l in 1:Q) { # Questions
            p[i,j,k,l] ~ dunif(0,1) 
            ans.pred[i,j,k,l] ~ dbin(p[i,j,k,l], 9)
          }
        }
      }
    }
  }
")

jags <- jags.model(
  mod,
  data = list(
    S = nrow(fishsamples),
    C = length(conl),
    F = length(fisl),
    cong = cong,
#    LOQ = LOQ,
#    below.LOQ = is.na(cong)*1,
    fis = match(fishsamples$Fish, fisl),
    U = nrow(answer),
    Q = ncol(answer),
    ans = answer,
    age = match(survey$Ages, agel),
    country = match(survey$Country, countryl),
    gender = match(survey$Gender, genderl),
    A = length(agel),
    B = length(countryl),
    H = length(genderl)
  ),
  n.chains = 4,
  n.adapt = 100
)

update(jags, 100)
samps <- jags.samples(jags, c('mu', 'pcd.pred', 'ans.pred'), 1000)
samps.coda <- coda.samples(jags, c('mu', 'pcd.pred', 'ans.pred'), 1000)
objects.store(samps)

#library(plyr)
#temp <- adply(samps$mu, c(1,2,3,4))

#if(TRUE){

  pl <- melt(
    array(samps$pcd.pred, 
     dim = c(length(conl), length(fisl), 1000, 4), 
     dimnames = list(
       Congener = conl, 
       Fish = fisl,
       Iter = 1:1000,
       Seed = c("S1","S2","S3","S4")
     )
  ))
  
  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,1E3))+
    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,1E3))+
  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") + coord_cartesian(ylim = c(1E-3,1E3))

#  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 = Fish))+geom_density(size = 1) + 
    facet_wrap(~ Congener, scale = "free_y") + coord_cartesian(xlim = c(-1.5,1.5))
  
  samps.coda <- coda.samples(jags, c('mu', 'pred'), 1000)
  
  out <- data.frame(
    Congener = conl,
    Fish = rep(fisl, each = length(conl)),
    Param = rep(c("mu", "pred"), each = length(conl)*length(fisl)),
    Value = (summary(samps.coda)[[1]][,1])
  )
  
  oprint(out)
  
  ggplot(out[out$Param == "mu", ], aes(x = Congener, y = Value, colour = Fish, group = Fish)) + 
    geom_line() + labs(title = "Congener amounts in fish compared with TCDD", y = "Log_10 ratio") 
  
  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
#} # IF FALSE