Benefit-risk assessment of Baltic herring and salmon intake

This page is a assessment. The page identifier is Op_en7748
Moderator:Arja (see all)
This page is a stub. You may improve it into a full page, and then a rating bar will appear here.
Upload data Show results

Scope

This assessment is part of the WP5 work in Goherr project. Purpose is to evaluate health benefits and risks caused of eating Baltic herring and salmon in four Baltic sea countries (Denmark, Estonia, Finland and Sweden). This assessment is currently on-going.

Question

What are the current population level health benefits and risks of eating Baltic herring and salmon in Finland, Estonia, Denmark and Sweden? How would the health effects change in the future, if consumption of Baltic herring and salmon changes due to actions caused by different management scenarios of Baltic sea fish stocks?

Intended use and users

Results of this assessment are used to inform policy makers about the health impacts of fish. Further, this assessment will be combined with the results of the other Goherr WPs to produce estimates of future health impacts of Baltic fish related to different policy options. Especially, results of this assessment will be used as input in the decision support model built in Goherr WP6.

Participants

National institute for health and welfare (THL)
Goherr project group

Boundaries

Four baltic sea countries (Denmark, Estonia, Finland, Sweden)
Current situation (fish use year 2016, pollutant levels in fish year 2010?)
Estimation for future (year 2020?)

Decisions and scenarios

Management scenarios developed in Goherr WP3 frames the following boundaries to the use and consumption of Baltic herring and salmon as human food. Effect of these scenarios to the dioxin levels and the human food use will be evalauted quantitatively and feed into the health benefit-risk model to assess the health effect changes.

Scenario 1: “Transformation to sustainability”
- Hazardous substances, including dioxins, are gradually flushed out and the dioxin levels in Baltic herring are below or close to the maximum allowable level.
- Fish stocks are allowed to recover to levels, which makes maximum sustainable yield possible and increases the total catches of wild caught fish. The catches of salmon by commercial fisheries has stabilized at low level, while the share of recreational catch increases slightly.
- The use of the Baltic herring catch for food increases. A regional proactive management plan for the use of catch has increased the capacity of the fishing fleets to fish herring for food and through product development and joint marketing, have increased consumer demand for Baltic herring.
Scenario 2: “Business-as-usual”
- The commercial catches of salmon continue to decrease. The demand for top predatory species, such as salmon and cod remains high, while the demand for herring decreased further as a result of demographic changes.
- Most of the herring catch are used for fish meal and oil production in the region.
- The use of Baltic herring from the southern parts of the Baltic Sea where the dioxin contents are not likely to exceed the maximum allowable level, are prioritised for human consumption. In the absence of the demand in many of the Baltic Sea countries, majority of the herring intended for direct human consumption are exported to Russia.
Scenario 3: “Inequality”
- The nutrient and dioxins levels continue to decrease slowly.
- The commercial catches of salmon have decreased further as the general attitudes favour recreational fishing, which has also resulted in decreased demand.
- The herring catches have increased slightly, but the availability of herring suitable for human consumption remains low due to both, dioxin levels that remain above the maximum allowable limit in the northern Baltic Sea and the poor capacity to fish for food.
- The use of the catch varies between countries. In Estonia, for example, where the whole catch has been traditionally used for human consumption, there is no significant change in this respect, but in Finland, Sweden and Denmark, herring fishing is predominantly feed directed.
Scenario 4: “Transformation to protectionism”
- The level of hazardous substances also increases as emission sources are not adequately addressed.
- Commercial salmon fisheries disappears almost completely from the Baltic Sea, although restocking keeps small scale fisheries going.
- Many of the Baltic herring stocks are also fished above the maximum sustainable yield and total catches are declining.
- Owing to the growing dioxin levels detected in herring, majority of the catch is used for aquaculture.

Timing

Model development during 2016 and 2017
First set of results in March 2017, draft publication in March 2018

Answer

This section will be updated as soon as preliminary results are available

Results

Conclusions

Rationale

Error creating thumbnail: Unable to save thumbnail to destination

Schematic picture of the health benefit-risk model for Baltic herring and salmon intake.

Stakeholders

Policy makers
- Food safety authorities
- Fisheries management
Researchers
- Food safety
- Health
NGO's
- WWF
- Active consumers
- Marine Stewardship Council
Baltic sea fishers and producers?

Dependencies

Content of all variable pages below are preliminary and will be updated when this assessment proceeds.

Toxic equivalency factor (TEF)
Consumption of fish
Pollutant and fatty acid levels
- Persistent organic pollutant levels in Baltic herring
- Persistent organic pollutant levels in Baltic salmon
- Nutrient levels in fish --# : Includes three tables, content modelled with Analytica years ago. How to update these? --Arja (talk) 07:20, 20 October 2016 (UTC)
Risk functions of health effects
Background health data (incidences and DALY)
- Burden of disease as DALY in Finland
- Disease risk: incidences in Finland and DALY in Europe
Disability weights of health effects
Length of disease

Analyses

Indices

Country (Denmark, Estonia, Finland, Sweden)
Year (current, future)
Sex
Age (categories?)
Fish species (Baltic herring, Baltic salmon)
Health end-point (ICD-10 code?, name? something else? --# : This need to be matched for disability weight, duration of disease, background incidence and disease burden data, dose-responses --Arja (talk) 09:39, 22 September 2016 (UTC)
Compound (pollutant, fatty acid) --# : This needs to be matched with dose-responses --Arja (talk) 09:39, 22 September 2016 (UTC)

Calculations

This section will have the actual health benefit-risk model (schematically described in the above figure) written with R. The code will utilise all variables listed in the above Dependencies section. Model results are presented as tables and figures when those are available.

Bayes model for intake

Model run 28.2.2017 [1]
Model run 28.2.2017 with corrected survey model [2]
Model run 28.2.2017 with Mu estimates [3]
Model run 1.3.2017 [4]

+ Show code - Hide code

# 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(fisl), length(conl), 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(fisl), length(conl), 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")

Exposure model

+ Show code - Hide code


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

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

objects.latest("Op_en7748", code_name = "bayes") #: pcd.pred, ans.pred, mu.pred

#### Dioxin exposure 

pl <- melt(pcd.pred)
levels(pl$Congener) <- gsub("HCDD", "HxCDD", levels(pl$Congener))
levels(pl$Congener) <- gsub("HCDF", "HxCDF", levels(pl$Congener))

colnames(pl)[colnames(pl) == "value"] <- "Result"
dioxconcsalmon <- subset(pl, pl$Fish == "Salmon" & pl$Seed == "S1")
dioxconcsalmon <- Ovariable("dioxconcsalmon", data = dioxconcsalmon)

tef <- opbase.data("Op_en4017.tef_values")
tef <-  Ovariable("tef", data = tef)

teqsalmon <- tef * dioxconcsalmon
teqsalmon <- aggregate(teqsalmon@output$Result, by = list(teqsalmon@output$Iter), sum)
colnames(teqsalmon) <- c("Iter","Result")
teqsalmon <- Ovariable("teqsalmon", data = teqsalmon)

dioxconcherring <- subset(pl, pl$Fish == "Baltic herring" & pl$Seed == "S1")
dioxconcherring <- Ovariable("dioxconcherring", data = dioxconcherring)

teqherring <- tef * dioxconcherring
teqherring <- aggregate(teqherring@output$Result, by = list(teqherring@output$Iter), sum)
colnames(teqherring) <- c("Iter","Result")
teqherring <- Ovariable("teqherring", data = teqherring)

objects.latest("Op_en7749", code_name = "calculate_amount") # amount of salmon and herring intake

expoherring <- teqherring * amountHerring
expoherring <- Ovariable("exposalmon", data = expoherring@output[c(1,5:7,29)])
exposalmon <- teqsalmon * amountSalmon
exposalmon <- Ovariable("exposalmon", data = exposalmon@output[c(1,5:7,29)])
expodiox <- exposalmon + expoherring

objects.latest("Op_en1838", code_name = "bayes") #nutrient concentrations

##Nutrient intake
nutconc <- melt(concddeo.pred)

Health impact model (Monte Carlo)

Model run 13.3.2017: a simple copy of op_fi:Silakan hyöty-riskiarvio [5]

+ Show code - Hide code

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

library(OpasnetUtils)
library(ggplot2)

exposure <- 1
BW <- 70

solet <- opbase.data("Op_fi3831.saantioletukset")

############# Tämän yläpuoliset eivät ole oikeaa koodia!!!

#!!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
## Tausta-altistus D-vitamiinille ja omega-3:lle
addexposure <- EvalOutput(Ovariable("addexposure", ddata = "Op_fi3831.tausta_altistus"))
#ii++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

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

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

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.

taustat <- Ovariable(
  output = data.frame(
    Iter = 1:get("N", envir = openv),
    Tausta = c("Kyllä", "Ei")[sample(2, get("N", envir = openv), replace = TRUE)],
    Result = 1
  ),
  marginal = c(TRUE, FALSE, FALSE)
)

addexposure <- addexposure * taustat

exposure <- exposure + addexposure

# Convert mother's dioxin intake (pg/d) into child's dioxin concentration (pg/g) in fat after 6 months of breast-feeding.
# For details, see [[ERF of dioxin#Dental defects]]

t0.5	<- Ovariable("t0.5", data = solet[solet$Muuttuja == "TEQ-puoliintumisaika" , ]["Result"])
f_ing	<- Ovariable("f_ing", data = solet[solet$Muuttuja == "TEQ-imeytymisosuus" , ]["Result"])
f_mtoc	<- Ovariable("f_mtoc", data = solet[solet$Muuttuja == "TEQ-osuus lapseen" , ]["Result"])
BF	<- Ovariable("BF", data = solet[solet$Muuttuja == "Imeväisen rasvamäärä" , ]["Result"])
BF <- EvalOutput(BF) # LISÄTTY 13.3.
temp <- exposure
temp@output <- temp@output[temp@output$Exposure_agent == "TEQ" , ]
temp <- log((temp * t0.5 * f_ing * f_mtoc / (log(2) * BF) ) + 1) # Actual conversion
temp@output$Exposure_agent <- "logTEQ"
temp <- unkeep(temp, prevresults = TRUE, sources = TRUE)
exposure@output <- orbind(exposure, temp)
exposure <- unkeep(exposure, prevresults = TRUE, sources = TRUE)

frexposed <- 1
bgexposure <- addexposure
date()
######################################################### 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(Result = 1))

# Väkimäärä. TÄTÄ KÄYTETÄÄN disincidence-ovariablen suhteuttamisessa CHD:n suhteen.

#!!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
pop <- opbase.data("Op_en2949") # [[Population of Finland]]
#ii+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

pop$Obs <- NULL
pop$Observation <- NULL
age <- pop$Result[pop$Age == "0-4"]
pop$Result[pop$Age == "0-4"] <- 4/5 * age
pop <- orbind(pop, data.frame(Age = "0", Result = 1/5 * age))
levels(pop$Age)[levels(pop$Age) == "0-4"] <- "1-4"
pop <- Ovariable("pop", data = pop)
date()
###############################################################################################
# ANNOSVASTEET JA TAUTIRISKIT 
# dioksiinille, 
# omega-3-rasvahapoille,
# vitamiineille (D-vitamiinille)
# ja metyylielohopealle ja niiden vaikutuksille
# Kuitenkaan metyylielohopeapitoisuuksia ei ole tässä malliversiossa, joten ne tippuvat pois.

# Annosvasteet

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

#!!++++++++++++++++++++++++++++++++++++++++++++++++++
#objects.latest("Op_en5823", code_name = "initiate") # [[ERF of dioxin]], ovariables ERF, threshold
#ii++++++++++++++++++++++++++++++++++++++++++++++++++

#temp1 <- ERF
#temp2 <- threshold

#!!++++++++++++++++++++++++++++++++++++++++++++++++++
#objects.latest("Op_en5830", code_name = "initiate") # [[ERF of omega-3 fatty acids]], ovariables ERF, threshold
#ii++++++++++++++++++++++++++++++++++++++++++++++++++

#temp1@data <- orbind(ERF@data, temp1@data)
#temp2@data <- orbind(threshold@data, temp2@data)

#!!++++++++++++++++++++++++++++++++++++++++++++++++++
#objects.latest("Op_en6866", code_name = "initiate") # [[ERFs of vitamins]], ovariables ERF, threshold
#ii++++++++++++++++++++++++++++++++++++++++++++++++++

#temp1@data <- orbind(ERF@data, temp1@data)
#temp2@data <- orbind(threshold@data, temp2@data)

#!!++++++++++++++++++++++++++++++++++++++++++++++++++
#objects.latest("Op_en5825", code_name = "initiate") # [[ERF of methylmercury]], ovariables ERF, threshold
#ii++++++++++++++++++++++++++++++++++++++++++++++++++

#ERF@data <- orbind(ERF@data, temp1@data)
#threshold@data <- orbind(threshold@data, temp2@data)

# Drop ERF rows that are not used in this assessment.

erfpois <- paste(ERF@data$Exposure_agent, ERF@data$Trait, ERF@data$ERF_parameter, ERF@data$Scaling) %in%
  c(
    "MeHg Child's IQ ERS BW",
    "DHA Child's intelligence ERS BW",
    "Omega3 Coronary heart disease RR None",
    "Omega3 CHD Relative Hill None",
    "logTEQ Developmental dental defects incl. agenesis ERS None",
    "TEQ Cancer, total CSF BW",
    "logTEQ Tooth defect ERS None" # Poistetaan Seveso koska suomalainen ERF uskottavampi
  )
ERF@data <- ERF@data[!erfpois , ]

# Drop nuisance indices because they use a lot of memory in oapply.
dropp <- c("Response_metric", "Exposure_route", "Exposure_metric", "Exposure_unit")
ERF <- unkeep(EvalOutput(ERF), dropp, sources = TRUE)
threshold <- unkeep(EvalOutput(threshold), dropp, sources = TRUE)


######################################## Tautiriski
#!!++++++++++++++++++++++++++++++++++++++++++++++++++
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++++++++++++++++++++++++++++++++++++++++++++++++++++
syy@data$Ikä <- gsub(" ", "", syy@data$Ikä)
colnames(syy@data)[colnames(syy@data) == "Ikä"] <- "Age"

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

pop <- EvalOutput(pop)

syyt <- syy / (pop / 2) * väli # Pop ei ole sukupuolen mukaan jaoteltu mutta syy on, joten pistetään kahtia.

disincidence <- disincidence * Ovariable(output = data.frame(Result = 1), marginal = FALSE)
marginals <- union(colnames(syyt@output)[syyt@marginal], colnames(disincidence@output)[disincidence@marginal])
syyt@output <- orbind(disincidence, syyt)
syyt@marginal <- colnames(syyt@output) %in% marginals
syyt <- unkeep(syyt, cols = c("Kuolemansyy", "Observation"), prevresults = TRUE, sources = TRUE)
syyt@output$Sukupuoli[syyt@output$Sukupuoli == "Miehet"] <- "Mies" # Ei ole faktori vaan character
syyt@output$Sukupuoli[syyt@output$Sukupuoli == "Naiset"] <- "Nainen"

disincidence <- syyt # Nyt sisältää sekä alkuperäisen disincidencen että ikäluokittaisen CHD:n
date()
###################################################################################
#### HIA-ovariablet

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

dose <- unkeep(EvalOutput(dose), prevresults = TRUE, sources = TRUE)

RR <- unkeep(EvalOutput(RR), prevresults = TRUE, sources = TRUE)

totcases <- EvalOutput(totcases)
date()

ggplot(totcases@output, aes(x = Response, weight = totcasesResult, fill = Sukupuoli))+
  geom_bar()

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

date()

exposure <- exposure * skaala * tiedot
amount <- amount * tiedot
concentration <- unkeep(concentration, prevresults = TRUE, sources = TRUE) * skaala

#Etsitään hedelmällisessä iässä olevat naiset ja lasten ÄO-vaste ja korjataan synnytyksen todennäköisyydellä.

indivrisk <- totcases

indivrisk <- indivrisk * tiedot

result(indivrisk) <- result(indivrisk) * ifelse(
  indivrisk@output$Trait %in% c("Child's IQ", "Tooth defect", "Dental defect"), 
  ifelse(
    indivrisk@output$Sukupuoli == "Nainen" & indivrisk@output$Hedelm == "TRUE", 
    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")
date()

####################################3

for(i in c("RR", "disincidence", "totcases", "ERF", "threshold")) {
  print(levels(get(i)@output$Response))
}
} # if(FALSE)

Plot concentrations and survey

Requires codes Op_en7748/bayes and indirectly Op_en7748/preprocess.
Model run 1.3.2017 [6]

+ Show code - Hide code

#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

References

Keywords

Related files

Benefit-risk assessment of Baltic herring and salmon intake

Contents

Scope

Question

Intended use and users

Participants

Boundaries

Decisions and scenarios

Timing

Answer

Results

Conclusions

Rationale

Stakeholders

Dependencies

Analyses

Indices

Calculations

Bayes model for intake

Exposure model

Health impact model (Monte Carlo)

Plot concentrations and survey

References

Keywords

See also

Related files

Navigation menu

Personal tools

Namespaces

Variants

Views

More

Search

Navigation

Page Tools

Tools

In other websites