Benefit-risk assessment of Baltic herring and salmon intake

Assessment presentation · show ready-made model results (without salmon)

# This is code Op_en7748/ on page [[Benefit-risk assessment of Baltic herring and salmon intake]]
# This code simply runs ggplot objects from the code further down on this page at "Second part"

library(OpasnetUtils)
objects.latest("Op_en7748", code_name="plots")

plot1
plot2
plot3
plot4
plot5
plot6
plot7
plot8
plot9
plot10
plot11

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 (not year specific)

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

This page is a knowledge crystal of subtype assessment. The page identifier is Op_en7748
Moderator:Arja (see all)
This page is a stub. You may improve it into a full page.
Upload data {{#opasnet_base_link:Op_en7748}}

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

Calculation of cases of disease

totcases (Op_en2261/totcases on page Health impact assessment) with dose and RR, automatic intermediate variables. Indices: Age, Gender, Country, Response.
- population, case-specific, from main model. Indices: Age, Gender, Country.
- disincidence,Incidence of the disease of interest. Op_en5917/initiate Disease risk. Indices: Age, Gender, Country, Response.
- ERF and threshold, exposure-response functions. Op_en2031/initiate Exposure-response function. Indices: Exposure_agent, Response, Exposure, Exposure_unit, ER_function, Scaling.
- frexposed, fraction of population that is exposed from Goherr: Fish consumption study. Indices: Age, Gender, Country.
- exposure, from the main model because has case-specific adjustments. Indices: Age, Gender, Country, Compound.
  - amount, Consumption of fish. Indices: Gender, Age, Country, Fish.
  - EU-kalat: concentration, Pollutant and fatty acid concentrations in fish. Indices: Fish_species → Fish, POP → Compound, Catch_square (Catch_site, Catch_location) → Area, Length_mean_mm → Length.
    - Concentrations of beneficial nutrients in fish. Indices: Fish, Compound. ----#: . Good data for Baltic herring, Salmon not taken into account yet --Arja (talk) 13:27, 13 March 2017 (UTC) (type: truth; paradigms: science: comment) ←--#: . Take data from Fineli. --Jouni (talk) 16:03, 13 April 2017 (UTC) (type: truth; paradigms: science: defence) ----#: . Omega-3 content in salmon: update and change the answer to point to the main page. --Jouni (talk) 16:03, 13 April 2017 (UTC) (type: truth; paradigms: science: comment)
    - Mercury concentrations in fish in Finland. Suggested indices: Fish, Location, Size, Year, Compound. ←--#: . update code by using Bayesian model on Kerty database. --Jouni (talk) 16:03, 13 April 2017 (UTC) (type: truth; paradigms: science: defence)
  - Toxic equivalency factor (TEF). Indices: TEFversion, Compound. ←--#: . Done. --Jouni (talk) 14:22, 19 April 2017 (UTC) (type: truth; paradigms: science: defence)

Calculation of DALYs:

Decisions(-)
Obs	Decision maker	Decision	Option	Variable	Cell	Change	Result	Description
1	Food safety authority	Background	Yes	bgexposure		Identity		Consider background from other sources
2	Food safety authority	Background	No	bgexposure		Replace	0	Do not consider background from other sources

Background exposure(-)
Obs	Country	Gender	Exposure_agent	Result	Unit	Description
1	FI	Male	Vitamin D	11.7	µg /d	Finriski 12 - 0.3 silakasta
2	SWE	Male	Vitamin D	11.7	µg /d	Finriski 12 - 0.3 silakasta
3	EST	Male	Vitamin D	11.7	µg /d	Finriski 12 - 0.3 silakasta
4	DK	Male	Vitamin D	11.7	µg /d	Finriski 12 - 0.3 silakasta
5		Female	Vitamin D	8.5	µg /d	Finriski 8.7 - 0.2 silakasta
6		Male	EPA	120	mg /d	Finriski 125 - 4.6 silakasta
7		Female	EPA	96	mg /d	Finriski 100 - 3.9 silakasta
8		Male	DHA	118	mg /d	Finriski 125 - 6.7 silakasta
9		Female	DHA	94	mg /d	Finriski 100 - 5.4 silakasta
10			PCDDF	0	pg /d (TEQ)
11			PCB	0	pg /d (TEQ)
12			MeHg	0	µg /d
13			logTEQ	0	log(pg /g)

Exposure-response functions of interest(-)
Obs	Exposure_agent	Resp	Response	ER_function	Scaling	Dummy
1	TEQ	Tooth defect	Yes or no dental defect	ERS	None	1
2	TEQ	Cancer	Cancer morbidity	CSF	BW	1
3	TEQ	Dioxin TDI	Dioxin recommendation tolerable daily intake	TDI	BW	1
4	DHA	Child's IQ	Loss in child's IQ points	ERS	None	1
5	Omega3	Heart (CHD)	CHD2 mortality	Relative Hill	None	1
6	Omega3	Stroke	Stroke mortality	Relative Hill	None	1
7	Vitamin D	Vitamin D intake	Vitamin D recommendation	Step	None	1
8	MeHg	Child's IQ	Loss in child's IQ points	ERS	BW	1

DALYs of responses(DALY /case)
Obs	Resp	DALY	Description
1	Heart (CHD)	5 - 15	Assumes DW 1 and D 10 U 50%
2	Stroke	5 - 15	Assumes DW 1 and D 10 U 50 %
3	Tooth defect	0 - 0.12	DW 0.001 D 60 U 100 %. Or should we use this: Developmental defect: caries or missing tooth, 0.008 (0.003 - 0.017), Periodontitis weight from IHME. D: 1. U from IHME?
4	Cancer	0 - 0.28	DW 0.1 D 20, in addition loss of life expectancy 5 a. This comes from a lifetime exposure, so it is (linearly( assumed that 1/50 of this is caused by one-year exposure. U 100 %
5	Vitamin D intake	0.0001 - 0.0101	DW 0.001 D 1 U 101x. 1 if not met
6	Dioxin TDI	0.0001 - 0.0101	DW 0.001 D 1 U 101x. 1 if not met
7	Child's IQ	0.0517 (0.03 - 0.0817)	Intellectual disability, mild (IQ<70): 0.031/IQ point (0.018-0.049) From IHME. D 50 U from IHME.

DW = disability weight

D = duration (a)

U = uncertainty

Diox conc fishing areas(TEQ)
Obs	Area	Fish	Result	Concentrations used for countries
1	24	Herring	3.08	DK
2	25-29	Herring	1.95	SE
3	28	Herring	2.84
4	30-31	Herring	2.29	FI (reference: reflects EU-kalat study)
5	32	Herring	2.45	EE
6	24	Salmon	16.0	DK
7	25-29	Salmon	10.9
8	28	Salmon	10.9
9	30-31	Salmon	9.3	FI, SE (reference: reflects EU-kalat study)
10	32	Salmon	10.7	EE

PCDD/F and PCB concentration distributions for each fish species and country are estimated in the following way:

EU-kalat study is used as the reference, because it has a large number of measurements. They come mostly from Bothnian Sea (subdivision 30).
The distributions are like those from the EU-kalat study, except that the means are scaled based on the area of interest relative to the reference (see table above).
The relative standard deviations are assumed to be the same for each fish population.
Baltic herring (>= 17 cm) and Small Baltic herring (<17 cm) are treated in the model as if they were separate species. Estonians are assumed to eat only small, others only large Baltic herring. ----#: . It should be checked how realistic assumption this is. --Jouni (talk) 13:40, 10 October 2017 (UTC) (type: truth; paradigms: science: comment)

Population data from Eurostat database

Population(n)
Obs	Country	Gender	Age	Result
1	DK	Male	18-45	1032266
2	DK	Male	>45	1207260
3	DK	Female	18-45	1003587
4	DK	Female	>45	1296678
5	EST	Male	18-45	252485
6	EST	Male	>45	237413
7	EST	Female	18-45	239174
8	EST	Female	>45	339881
9	FI	Male	18-45	970016
10	FI	Male	>45	1182808
11	FI	Female	18-45	921887
12	FI	Female	>45	1339537
13	SWE	Male	18-45	1823042
14	SWE	Male	>45	2065334
15	SWE	Female	18-45	1738408
16	SWE	Female	>45	2199156

**Landing of fish from different areas in Baltic Sea by country.**
Fish species and area	Country
Fish species and area	DK	FI	EE	SE
Tot herring (1000t)	4.8	132.4	35.3	66.5
Total salmon (t)	122	367	9	315
Herring 24	73 %			1 %
Herring 25-29	6 %	20 %	21 %	76 %
Herring 28.1			47 %
Herring 30		73 %		21 %
Herring 31		3 %		0 %
Herring 32	0 %	3 %	33 %	0 %
Salmon 25-29		7 %	11 %	6 %
Salmon 30-31		78 %		94 %
Salmon 32		15 %	67 %	0 %
Salmon 28			22 %
Salmon 24	92 %

**ICES subdivisions**
Subdivision	Name
21	Kattegat
22	Great Belt
23	Öresund
24	Baltic West of Bornholm
25	Southern Central Baltic West
26	Southern Central Baltic East
27	North West of Gotland
28	Eastern Main Basin
28.1	Gulf of Riga
28.2	Eastern Main Basin
29	Åland Sea
30	Bothnian Sea
31	Bothnian Bay
32	Gulf of Finland

**Number of samples in EU-kalat study for PCDD/F and PCB measurements.**
Catch location	Subdivision	Fish_species	Number of samples
Hanko	29	Small Baltic herring (< 17 cm)	2
Kotka	32	Small herring	2
Oulu	31	Small herring	2
Pori	30	Small herring	32
Turku	29	Small herring	2
Vaasa	30	Small herring	2
Hanko	29	Large Baltic herring (>= 17 cm)	3
Kotka	32	Large herring	3
Oulu	31	Large herring	2
Pori	30	Large herring	42
Turku	29	Large herring	3
Vaasa	30	Large herring	3
Hanko	29	Salmon	1
Kotka	32	Salmon	2
Oulu	31	Salmon	2
Pori	30	Salmon	2
Turku	29	Salmon	2

Analyses

Indices

Country (Denmark, Estonia, Finland, Sweden)
Year (current, future)
Gender
Age: 18-45 years or >45 years
Fish species (Baltic herring, Baltic salmon)
Health end-point, specified by name
Compound: TEQ (PCDD/F and PCB), Vitamin D, Omega3 (includes EPA and DHA), MeHg

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.

18.5.2017: Archived exposure model Op7748/exposure by Arja (used separate ovariables for salmon and herring) [1]

Health impact model (Monte Carlo)

Model run 13.3.2017: a simple copy of op_fi:Silakan hyöty-riskiarvio [2]
Model run 13.3.2017 with showLocations function [3]
Model run 13.3.2017 produces totcases results but are not meaningful yet [4]
Model run 14.3.2017 with exposure graph [5]
Model run 14.3.2017 bugs not fixed [6]
Model run 30.5.2017 [7]
Model run 12.6.2017 with 2D Monte Carlo [8]
Model run 8.9.2017 with known bugs fixed [9]
Model run 6.10.2017: added updated mc2d function, better exposure$Exposure (To eater or To child), Background as marginal after mc2d, MeHg concentrations; bugs fixed with conc_mehg and IQ disabilityweight [10]
Model run 6.10.2017: 1000 iterations [11]
Model run 11.10.2017: country-specific concentrations [12]
Model run 12.10.2017: Background bug fixed [13] objects.get("150780268107")
Model run 12.10.2017: with marginal background and new Limit decision [14] objects.get("150783160307")
Model run 19.10.2017: n = 2117 and mc2d$run2d = FALSE [15]
Model run 27.10.2017: salmon removed from amount. [16]
Model run 15.11.2017: salmon included [17]
Model run 21.11.2017: using EU-kalat concentration model from May 2017 [18]
Model run 22.11.2017: correcting for conc and Large herring, also debugged mc2d. [19]

⇤--#: . We must check that Baltic herring, Large herring and Small herring are used in all variables coherently. There may be problems in conc. --Jouni (talk) 09:39, 21 November 2017 (UTC) (type: truth; paradigms: science: attack)

----#: . This must be corrected back when the bayesian conc_pcddf works again. --Jouni (talk) 14:11, 22 November 2017 (UTC) (type: truth; paradigms: science: comment)

⇤--#: . Country, Ages, and Gender are marginals in both RR an casesabs although mc2dparam$run2d==FALSE. Yet, in casesabs all combinations have 2117 rows, while in RR all combinations have 80-190 rows and only the combinations Exposcen*Limit*Background*Resp have 2117 rows each. So, RR works as expected with run2d==FALSE but casesabs does not. What is wrong? --Jouni (talk) 14:11, 22 November 2017 (UTC) (type: truth; paradigms: science: attack)

+ Show code - Hide code

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

# Decisions
DecisionTableParser(opbase.data("Op_en7748", subset = "Decisions"))

# amount does not know where to find jsp.
objects.latest("Op_en7749", code_name = "surveyjsp") # Uses jsp directly from survey data.

# Assumes that Estonians eat only small herring while all others eat only large herring.
jsp$Fish <- as.factor(
  ifelse(
    jsp$Fish == "Herring",
    ifelse(
      jsp$Country == "EST",
      "Small herring",
      "Large herring"
    ),
    as.character(jsp$Fish)
  )
)

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_mehg <- Ovariable(
  "conc_mehg",
  dependencies = data.frame(
    Name = c("Hg"),
    Ident = "Op_en4004/Hg" # Contains Kerty measurements
  ),
  formula = function(...) {
    
    test <- paste(Hg$Parameter, Hg$Näytelaji, Hg$Helcom.alue)
    meanSalm  <- Hg$Result[test == "Mean lohi "]
    meanHerrN <- Hg$Result[test == "Mean silakka C-Selkämeri"]
    meanHerrS <- Hg$Result[test == "Mean silakka H-Pohjoinen Itämeri"]
    sdSalm    <- Hg$Result[test == "SD lohi "]
    sdHerrN   <- Hg$Result[test == "SD silakka C-Selkämeri"]
    sdHerrS   <- Hg$Result[test == "SD silakka H-Pohjoinen Itämeri"]
    
    N <- openv$N
    
    out <- data.frame(
      Exposure_agent = "MeHg", # Assuming all mercury in fish is methylmercury
      Fish = rep(c("Herring", "Salmon"), each = N),
      Iter = rep(1:N, 2),
      Country = rep(c("FI", "SWE", "DK", "EST"), each = 2*N),
      Result = c(
        rnorm(N, meanHerrN, sdHerrN),
        rnorm(N, meanSalm, sdSalm),
        rnorm(N, meanHerrS, sdHerrS),
        rnorm(N, meanSalm, sdSalm),
        rnorm(N, meanHerrS, sdHerrS),
        rnorm(N, meanSalm, sdSalm),
        rnorm(N, meanHerrS, sdHerrS),
        rnorm(N, meanSalm, sdSalm)
      )
    )
    return(out)
  }
)
conc_mehg <- EvalOutput(conc_mehg)
conc_mehg@marginal <- c(rep(TRUE,4), FALSE, TRUE)

#Adjust PCDD/F measurements from the Bothnian Sea [[EU-kalat]] with country-specific coefficients
objects.latest("Op_en3104", code_name = "initiate")

# Fetch a model run by hand from May 2017, because the JAGS model does not currently work in Opasnet.
objects.get("149562691490") # [[EU-kalat]] code bayes. Model run http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=kNNzEMTSD4N2f0Yy

conc <- Ovariable(
  "conc",
  dependencies = data.frame(
    Name = c(
      "conc_pcddf", # [[EU-kalat]]
      "conc_vit", # [[Concentrations of beneficial nutrients in fish]]
      "conc_mehg" # data from [[Mercury concentrations in fish in Finland]]
    ), 
    Ident = c("Op_en3104/initiate", "Op_en1838/initiate", NA)
  ),
  formula = function(...) {
    
    ############ Adjust PCDD/F concentrations with country
    ## coeff: country-specific coefficients
    d <- opbase.data("Op_en7748", subset = "Diox conc fishing areas")

    coeff <- Ovariable(output = data.frame(
      Fish = rep(c("Small herring", "Large herring", "Salmon"), each = 4),
      Country = c("DK", "EST", "FI", "SWE"),
      Result = c(
        d$Result[d$Fish=="Herring"&d$Area=="24"],
        d$Result[d$Fish=="Herring"&d$Area=="32"],
        d$Result[d$Fish=="Herring"&d$Area=="30-31"],
        d$Result[d$Fish=="Herring"&d$Area=="25-29"],
        d$Result[d$Fish=="Herring"&d$Area=="24"],
        d$Result[d$Fish=="Herring"&d$Area=="32"],
        d$Result[d$Fish=="Herring"&d$Area=="30-31"],
        d$Result[d$Fish=="Herring"&d$Area=="25-29"],
        d$Result[d$Fish=="Salmon"&d$Area=="24"],
        d$Result[d$Fish=="Salmon"&d$Area=="32"],
        d$Result[d$Fish=="Salmon"&d$Area=="30-31"],
        d$Result[d$Fish=="Salmon"&d$Area=="30-31"]
      )
    )) / Ovariable(output = data.frame(
      Fish = c("Small herring", "Large herring", "Salmon"),
      Result = c(
        d$Result[d$Fish=="Herring"&d$Area=="30-31"],
        d$Result[d$Fish=="Herring"&d$Area=="30-31"],
        d$Result[d$Fish=="Salmon"&d$Area=="30-31"]
      )
    ))
    
    ## This fish name change is temporarily needed because conc_pcddf data from May 2017 is used
    conc_pcddf$Fish <- as.factor(conc_pcddf$Fish)
    levels(conc_pcddf$Fish)[levels(conc_pcddf$Fish)=="Baltic herring"] <- "Large herring"
    conc_pcddf <- conc_pcddf * coeff
    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.
    conc_pcddf <- unkeep(conc_pcddf, sources = TRUE)
    
    ################### Adjust conc_vit and conc_mehg
    
    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_vit, conc_mehg)
    out@output <- fillna(out@output, marginals = "Country")
    temp <- out[out$Fish == "Herring",] # Have to adjust conc_vit and conc_mehg
    levels(temp$Fish)[levels(temp$Fish)=="Herring"] <- "Large herring"
    levels(out$Fish)[levels(out$Fish)=="Herring"] <- "Small herring"

    out <- combine(out, temp, conc_pcddf)
    result(out) <- pmax(0,result(out))
    return(out)
  }
)
conc@meta <- c(
  conc@meta,
  list(units = "PCDDF, PCB, TEQ: pg /g; Vitamin D, MeHg: µg /g; DHA, EPA, Omega3: mg /g fresh weight")
)

# Exposure with fish-related exposure and additional exposure from other sources.
# Does NOT contain infant's exposure from mother.

######## Standard amount but needs a Limit scenario

objects.latest("Op_en7749", code_name = "initiate") # [[Goherr: Fish consumption study]]
amount <- amount * Ovariable(
  output = data.frame(Limit = c("No limit", "Max 3 g/d"), Result = 1),
  marginal = c(TRUE, FALSE)
)
result(amount) <- ifelse(result(amount) > 3 & amount$Limit != "No limit", 3, result(amount))

########## For herring-specific scenario, remove salmon
# amount <- amount[amount$Fish != "Salmon" , ]

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"
      )
    )
    
    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"))
    
    temp1 <- out[out$Exposure_agent %in% c("PCDDF","PCB") , ]
    temp1 <- oapply(temp1, cols = "Exposure_agent", FUN = sum)
    temp1$Exposure_agent <- "TEQ"

    temp2 <- out[out$Exposure_agent %in% c("EPA", "DHA") , ]
    temp2 <- oapply(temp2, cols = "Exposure_agent", FUN = sum)
    temp2$Exposure_agent <- "Omega3"

    out <- combine(out, temp1, temp2)
    out <- unkeep(out, prevresults = TRUE, sources = TRUE)
    
    return(out)
  }
)

bgexposure <- CheckDecisions(EvalOutput(bgexposure))

info <- Ovariable(
  "info",
  dependencies = data.frame(Name = c("jsp")),
  formula = function(...) {
    out <- unique(jsp@output[c("Iter","Country","Gender","Ages")])
    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 <- as.factor(
      ifelse(
        expoRaw$Exposure_agent %in% c("DHA", "MeHg"),
        "To child",
        "To eater"
      )
    )
    dx.expo.child$Exposure <- "To child"
    out <- combine(expoRaw, dx.expo.child)
    out <- unkeep(out, sources=TRUE, prevresults=TRUE)
    return(out * info + 0.01) # Roughtly 98 % of all exposures for Exposure_agents is > 0.01.
  }
)
exposure@meta <- c(
	exposure@meta, 
	list(units = "PCDDF, PCB, TEQ: pg /day; Vitamin D, MeHg: µg /day; DHA, EPA, Omega3: mg /day")
)

# 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","info")),
  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(info * Ovariable(output=out))
  }
)
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")]
ERF$Exposure <- as.factor(ifelse(
  ERF$Exposure %in% c(
    "Maternal intake through placenta",
    "Maternal ingested intake",
    "Intake through placenta and mother's milk"
  ),
  "To child",
  "To eater"
))

threshold <- EvalOutput(threshold)
threshold <- threshold[ , !colnames(threshold@output) %in% c("Age","Sex","Source")]
threshold$Exposure <- as.factor(ifelse(
  threshold$Exposure %in% c(
    "Maternal intake through placenta",
    "Maternal ingested intake",
    "Intake through placenta and mother's milk"
  ),
  "To child",
  "To eater"
))

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 = 1000, # Number of iterations in the new Iter
  run2d = FALSE, # Should the mc2d function be used or not?
  info = info, # Ovariable that contains additional indices, e.g. newmarginals.
  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

BoD$Group <- paste(BoD$Gender, BoD$Ages)
BoD$Group <- factor(BoD$Group, levels = c(
  "Female 18-45",
  "Male 18-45",
  "Female >45",
  "Male >45"
))
BoD@meta <- c(BoD@meta, list(units = "DALY /year with continuous exposure"))

amount@meta <- c(amount@meta, list(units = "g /day"))

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

Second part

Model run 4.6.2017 [20]
Model run 11.6.2017 with 2D Monte Carlo [21]
Model run 6.10.2017 [22] barchart shows blocks wrong in Opasnet but correctly on own computer
Model run 12.10.2017 [23]
Model run 16.10.2017: ggplot objects stored [24]
Model run 21.11.2017: plots stored [25]

←--#: . These are the equal sizes for different graphics settings. A typically good base_size is 24:

Opasnet graphics

png(width = 1024, height=768) # (in pixels)

pdf(width = 14, height=10.5) # (in inches) --Jouni (talk) 15:52, 21 November 2017 (UTC) (type: truth; paradigms: science: defence)

+ Show code - Hide code

# This is code Op_en7748/plots 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(RR,fun=c("length","mean","sd")))[-5]
#oprint(summary(ERF,fun=c("length","mean"))[-2])
#oprint(summary(casesrr,fun=c("length","mean","sd")))
#oprint(summary(casesabs,fun=c("length","mean","sd")))[-1]
#oprint(summary(BoDt,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])

BS <- 24

plot1 <- ggplot(conc@output, aes(x=concResult, colour=Fish))+stat_ecdf()+
  theme_grey(base_size=BS)+
  facet_wrap(~Exposure_agent, scales="free_x")+scale_x_log10()

plot2 <- ggplot(amount@output, aes(x=amountResult, colour=Country))+stat_ecdf()+
  theme_grey(base_size=BS)+
  facet_wrap(~Fish, scales="free_x")+scale_x_log10()

plot3 <- ggplot(exposure@output[result(exposure)>0.01 & exposure$Background == "No",], aes(x = exposureResult, colour=Country))+
  stat_ecdf() + scale_x_log10() + theme_gray(base_size = BS)+
  facet_wrap(~ Exposure_agent, scales="free_x")+labs(title = "Exposure (DX: pg/d, O3: mg/d, VD: ug/d)")

plot4 <- ggplot(exposure@output[exposure$Background == "No",], aes(x = exposureResult, colour=Country))+
  stat_ecdf() + scale_x_log10() + theme_gray(base_size = BS)+
  facet_wrap(~ Exposure_agent, scales="free_x")+labs(title = "Exposure (DX: pg/d, O3: mg/d, VD: ug/d)")

exposure$Group <- paste(exposure$Gender, exposure$Ages)
exposure$Group <- factor(exposure$Group, levels = c(
  "Female 18-45",
  "Male 18-45",
  "Female >45",
  "Male >45"
))

plot5 <- ggplot(exposure@output[exposure$Background == "No" & exposure$Exposure_agent == "TEQ",], aes(x = exposureResult, colour=Country))+
  stat_ecdf(size=2) + scale_x_log10() + theme_gray(base_size = BS)+
  facet_grid(Ages ~ Gender, scales="free_x")+labs(title = "Exposure (DX: pg/d, O3: mg/d, VD: ug/d)")

plot6 <- ggplot(ERF@output, aes(x=ERFResult, colour=Exposure_agent))+stat_ecdf()+
  theme_grey(base_size=BS)+
  facet_wrap(~ Resp, scales = "free_x")

plot7 <- ggplot(casesabs@output, aes(x=casesabsResult, colour=Background))+stat_ecdf()+
  theme_grey(base_size=BS)+
  facet_wrap(~ Resp, scales = "free_x")

plot8 <- ggplot(casesrr@output, aes(x=casesrrResult, colour=Background))+stat_ecdf()+
  theme_grey(base_size=BS)+
  facet_wrap(~ Resp, scales = "free_x")

plot9 <- ggplot(BoD@output[BoD$Background=="No" & BoD$Limit == "No limit" , ],
       aes(x=BoDResult, colour=Group))+stat_ecdf()+
  theme_grey(base_size=BS)+
  facet_wrap(~ Resp, scales = "free_x")

temp <- oapply(BoD, cols = "Iter", FUN = mean)@output

if(FALSE) BS <- 42

varit7 <- c(
  "slateblue",
  "mediumseagreen",
  "violetred",
  "greenyellow",
  "darkgoldenrod1",
  "firebrick1",
  "dodgerblue2"
)

plot10 <- ggplot(temp[temp$Limit == "No limit" , ],  aes(x=Group, weight=BoDResult, fill=Resp))+
  geom_bar()+facet_grid(Background ~ Country)+
  theme_grey(base_size=BS)+
  theme(
#    legend.position = "bottom",
    axis.text.x = element_text(angle = -90)
  )+
  labs(
    title="Burden of disease by country and background nutrient intake",
    y = "Disability-adjusted life years"
  )+
  scale_fill_manual(values=varit7)
#  coord_cartesian(ylim=c(-1000,1000))

plot11 <- ggplot(temp[temp$Country == "DK",],  aes(x=Group, weight=BoDResult, fill=Resp))+
  geom_bar()+facet_grid(Background ~ Limit)+
  theme_grey(base_size=BS)+
  theme(
#    legend.position = "bottom",
    axis.text.x = element_text(angle = -90)
  )+
  labs(
    title="Burden of disease by intake limit and background nutrient intake in Denmark",
    y = "Disability-adjusted life years"
  )+
  scale_fill_manual(values=varit7)

if(FALSE) {

png(file="BoD by country and background.png", height=1200, width=2400)
plot10
dev.off()

png(file="BoD in DK by limit and background.png", height=1200, width=2400)
plot11
dev.off()

} else {

print(plot1)
print(plot2)
print(plot3)
print(plot4)
print(plot5)
print(plot6)
print(plot7)
print(plot8)
print(plot9)
print(plot10)
print(plot11)
}

#objects.store(plot1, plot2, plot3, plot4, plot5, plot6, plot7, plot8, plot9, plot10, plot11)
#cat("ggplots plot1 - plot11 stored.\n")

odag()

Plot concentrations and survey

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

+ 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

Interface for BBN model

Health risk-benefit assessment model (BRA) implemented on this page produces also data for the overall Goherr Bayesian belief network model (BBN). In brief, the models are built in a way that they share the important nodes. This section describes how data from BRA is managed to fit BBN.

Data files and indices needed:

Consumer country (out.population from population)
- Country (DK, EST, FI, SWE): total population size for each country
Consumer gender (out.population from population)
- Gender (Female, Male): gender-specific number of people for each country
Consumer age group (out.population from population)
- Age (18-45, 45>): gender and age-specific number of people for each country
Effect of improved information
- Information improvements (policy options to be determined)
- Country, age, and gender-specific distributions ----#: . Is it realistic be so specific? --Jouni (talk) 12:47, 22 November 2017 (UTC) (type: truth; paradigms: science: comment)
Human consumption of salmon (out.salmonintake from amount[Fish=="Salmon",])
- Country, age, and gender-specific consumption (g/day)
- The same with or without salmon recommendation policy given improved information (to be determined)
Human consumption of herring: same as for salmon (out.herringintake from amount[Fish!="Salmon",])
Omega3 intake from salmon and herring (mg/day) (out.omegaintake from exposure[Exposure_agent=="Omega3"&Background=="No",])
Human intake of dioxin from herring (out.herringdioxinexp from amount * conc)
Human intake of dioxin from salmon (out.salmondioxinexp from amount * conc)
Dioxin intake total (out.dioxinintake from exposure[Exposure_agent=="TEQ",])
Other intake of omega3 (out.omegatotal$Background from exposure[Exposure_agent=="Omega3",])
Omega3 intake total (out.omegatotal from exposure[Exposure_agent=="Omega3",])
Net burden of disease (out.BoD from BoD summed over Response)

+ Show code - Hide code

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

library(OpasnetUtils)

objects.latest("Op_en7748", code_name="hia")

#d <- readLines("//cesium/jtue$/_Downloads/GOHERR_model_March17_Oregrund_testing.oobn")
#nam <- d[grepl("label =", d)]

out.population <- unkeep(population, sources=TRUE, prevresults=TRUE)@output
temp <- unkeep(amount, sources=TRUE, prevresults=TRUE)@output
out.salmonintake  <- temp[temp$Fish=="Salmon",]
out.herringintake <- temp[temp$Fish!="Salmon",]

temp <- unkeep(exposure, sources=TRUE, prevresults=TRUE)@output
temp <- temp[temp$Exposure=="To eater", colnames(temp)!="Exposure"]
out.omegaintake <- temp[temp$Exposure_agent=="Omega3"&temp$Background=="No",colnames(temp)!="Background"]
out.dioxinintake <- temp[temp$Exposure_agent=="TEQ"&temp$Background=="No",colnames(temp)!="Background"]
out.omegatotal <- temp[temp$Exposure_agent=="Omega3",]

temp <- unkeep(conc, sources=TRUE, prevresults=TRUE)@output
out.herringconc <- temp[temp$Exposure_agent=="TEQ" & temp$Fish!="Salmon",]
out.salmonconc <- temp[temp$Exposure_agent=="TEQ" & temp$Fish=="Salmon",]

temp <- unkeep(amount * conc[conc$Exposure_agent=="TEQ",], sources=TRUE, prevresults=TRUE)
out.herringdioxinexp <- temp@output[temp$Fish!="Salmon",]
out.salmondioxinexp <- temp@output[temp$Fish=="Salmon",]

out.BoD <- oapply(BoD, cols=c("Resp","Response","Group"), FUN=sum)
out.BoD <- unkeep(out.BoD, sources=TRUE, prevresults=TRUE)
out.BoD$Iter <- as.numeric(as.character(out.BoD$Iter))
out.BoD <- CollapseMarginal(out.BoD, cols = c("Country", "Ages","Gender"), fun = "sample")@output

write.csv(out.population, "population.csv", row.names=FALSE)
write.csv(out.salmonintake, "salmonintake.csv", row.names=FALSE)
write.csv(out.herringintake, "herringintake.csv", row.names=FALSE)
write.csv(out.omegaintake, "omegaintake.csv", row.names=FALSE)
write.csv(out.dioxinintake, "dioxinintake.csv", row.names=FALSE)
write.csv(out.omegatotal, "omegatotal.csv", row.names=FALSE)
write.csv(out.herringconc, "herringconc.csv", row.names=FALSE)
write.csv(out.salmonconc, "salmonconc.csv", row.names=FALSE)
write.csv(out.herringdioxinexp, "herringdioxinexp.csv", row.names=FALSE)
write.csv(out.salmondioxinexp, "salmondioxinexp.csv", row.names=FALSE)
write.csv(out.BoD, "BoD.csv", row.names=FALSE)

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

Health impact model (Monte Carlo)

Second part

Plot concentrations and survey

Interface for BBN model

References

Keywords

See also

Related files

Navigation menu