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

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

library(OpasnetUtils)
library(reshape2)
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") # [[Goherr: Fish consumption study]]: survey, surcol

# The following three rows should be in the preprocessing code.
survey$Row <- 1:nrow(survey)
survey$Weighting <- as.double(as.character(survey$Weighting))
survey$Ages <- ifelse(as.numeric(as.character(survey$Age)) < 46, "18-45",">45")

# Interesting fish eating questions
surv <- survey[c(1,3,4,16,29,30,31,46:49,95:98,158)]

agel <- unique(survey$Ages)
countryl <- unique(survey$Country)
genderl <- unique(survey$Gender)
ans <- data.matrix(surv[-c(1:3, ncol(surv))]) - 1 # Zero is the first option
qlen <-  unlist(lapply(
  surv, 
  FUN = function(x) {length(levels(x))}))[-c(1:3, ncol(surv))]
questionl <- colnames(surv)
agel
countryl
genderl

conl <- as.character(unique(eu@output$Congener))
fisl <- sort(as.character(unique(eu@output$Fish)))
conl
fisl
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) # Why does this not work with dnorm(0, 0.001)?
        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], qlen[j])
      }
    }
    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], qlen[l])
          }
        }
      }
    }
  }
")

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(ans),
    Q = ncol(ans),
    ans = ans,
    qlen = qlen-1, # Subtract 1 because starts from 0.
    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))

pcd.pred <- array(
  samps$pcd.pred, 
   dim = c(length(conl), length(fisl), 1000, 4), 
   dimnames = list(
     Fish = fisl,
     Congener = conl, 
     Iter = 1:1000,
     Seed = c("S1","S2","S3","S4")
   )
)

mu.pred <- array(
  samps$mu, 
  dim = c(length(conl), length(fisl), 1000, 4), 
  dimnames = list(
    Fish = fisl,
    Congener = conl, 
    Iter = 1:1000,
    Seed = c("S1","S2","S3","S4")
  )
)

ans.pred <- array(
  samps$ans.pred, 
  dim = c(
    length(agel),
    length(countryl),
    length(genderl), 
    ncol(ans),
    1000, # iterations
    4 # Seeds
  ), 
  dimnames = list(
    Age = agel,
    Country = countryl, 
    Gender = genderl,
    Question = colnames(ans),
    Iter = 1:1000,
    Seed = c("S1","S2","S3","S4")
  )
)

objects.store(pcd.pred, mu.pred, ans.pred)
cat("Arrays pcd.pred, mu.pred, ans.pred stored.\n")

#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 colnames(ans)) {
  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