View source for Benefit-risk assessment of Baltic herring and salmon intake

[[Category:Goherr]] 
[[Category:Ongoing assessment]]
{{progression class
|progression = Full draft
|curator = THL
|date = 2020-08-12
}}
'''For the final version of this assessment, see the scientific article.<ref name="tuomisto2020"/>
{{summary box
|question=What are the current individual and 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 a) fish consumption recommendations or limitations, b) different management scenarios of Baltic sea fish stocks, or c) selection of fish sizes for human consumption?
|answer=BONUS GOHERR project (2015-2018) looked at this particular question. A health impact assessment was performed based on fish consumption survey done in the four target countries; EU Fish 2 study about dioxin and PCB concentrations; a dynamic growth model about Baltic herring stock sizes and dioxin accumulation; scientific literature about exposure-response functions of several compounds found in Baltic fish; and online models produced in Opasnet web-workspace. This page describes the benefit-risk assessment performed.<ref name="tuomisto2020">Tuomisto, J.T., Asikainen, A., Meriläinen, P., Haapasaari, P. Health effects of nutrients and environmental pollutants in Baltic herring and salmon: a quantitative benefit-risk assessment. BMC Public Health 20, 64 (2020). https://doi.org/10.1186/s12889-019-8094-1</ref>

Dioxin and PCB concentrations have been constantly decreasing in Baltic fish for 40 years, and now they are mostly below EU limits. Also Baltic herring consumption has been decreasing during the last decades and is now a few grams per day, varying between age groups (old people eat more), genders (males eat more) and countries (Estonians eat more and Danes less than others). People reported that better availability of easy products, recipes, and reduced pollutant levels would increase their Baltic herring consumption. In contrast, recommendations to reduce consumption would have little effect on average.

Health benefits of Baltic herring and salmon clearly outweigh health risks in age groups over 45 years. Benefits are higher even in the most sensitive subgroup, women at childbearing age. The balance is close to even, if exceedance of the tolerable daily intake is given weight in the consideration and if other omega-3 sources are given priority over fish. The analysis was robust in a sense that we did not find uncertainties that could remarkably change the conclusions and suggest postponing decisions in hope of new information.

Overall, main arguments from this health assessment and other disciplines studied in Goherr are in favour of getting rid of dioxin-based food restrictions related to Baltic herring and salmon, and promoting human consumption of Baltic fish.
}}
[[image:Goherr_WP5_HIA_structure.JPG|thumb|500px|Schematic picture of the health benefit-risk model for Baltic herring and salmon intake.]]
{{#l:Goherr health benefit-risk assessment.odp|Assessment presentation}} &middot;
[http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=kAei8Vhxm3cEHQ4j show ready-made model results]

== Scope ==
This assessment is part of the [[Goherr:_WP5_Linking_the_health_of_the_Baltic_Sea_with_health_of_humans:_Dioxin|WP5]] work in [[Goherr:_Project | 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 individual and 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 a) fish consumption recommendations or limitations, b) different management scenarios of Baltic sea fish stocks, or c) selection of fish sizes for human consumption?

=== 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_Building_a_decision_support_model_for_integrated_governance|Goherr WP6]].

=== Participants ===

* University of Helsinki: Sakari Kuikka, Päivi Haapasaari, Suvi Ignatius, Kirsi Hoviniemi, Inari Helle, Annukka Lehikoinen, Mika Rahikainen
* National Institute for Health and Welfare (THL): Jouni Tuomisto, Arja Asikainen, Päivi Meriläinen
* University of Oulu: Timo P. Karjalainen, Simo Sarkki, Mia Pihlajamäki
* Swedish University of Agricultural Sciences (SLU): Anna Gårdmark, Johan Östergren, Magnus Huss, Andreas Bryhn, Philip Jacobson
* University of Aalborg/Innovative Fisheries Management (IFM-AAU): Alyne Delaney, Jesper Raakjaer
* Stakeholders needed in the assessment: fisher's associations, agricultural/fisheries ministries in Finland, Sweden, Estonia, Denmark.

=== 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)
* Area considered: Sweden, Finland, Denmark, Estonia. For detailed herring/salmon stock modelling, only the Bothnian Sea and Gulf of Bothnia is considered.
* Policies considered: See the table below.
* This assessment looks at health only. Goherr as a project will consider wider objectives: threats to and state of the fish stocks; impacts and governance responses. See below.

=== Decisions and scenarios ===

In Goherr, the scope is wide and we are looking at scenarios about fragmented vs. integrated governance and high vs. low human impact at the Baltic Sea Region. These scenarios are described in more detail on a page about [[Goherr:_3.2._Online_description_of_the_scenarios_developed,_applicable_in_the_dioxin_model_(Task_5.1)_and_decision_support_model_of_WP6_(Month_16)._Responsible_partner:_UOULU. |management scenarios developed in Goherr WP3]].

In this assessment, will look at health only, but the decisions and options considered are based on the wider discussions about integrated governance. The main decisions considered are improvements about fish availability and recipes, food recommendations, and selection of fish sizes used for consumption. The changes caused by each option are estimated based on the [[Goherr: Fish consumption study|Goherr fish consumption study]].<ref name="pihlajamaki2019"/>
* Info.improvements with options
** BAU (Business as usual): current availability of and chemical levels in herring and salmon.
** Yes: Better availability of and less chemicals in herring and salmon.
* Recomm.herring with options
** BAU: authorities recommendations about Baltic herring do not change.
** Eat more: authorities recommend to eat more Baltic herring.
** Eat less: authorities recommend to eat less Baltic herring.
* Recomm.salmon with options
** BAU: authorities recommendations about Baltic salmon do not change.
** Eat more: authorities recommend to eat more Baltic salmon.
** Eat less: authorities recommend to eat less Baltic salmon.
* Coherent consumer policy (Cons.policy): This is a combination of the three decisions above, with options:
** BAU: BAU is chosen for all three decisions.
** Less: BAU for Info.improvements and Eat less for both herring and salmon recommendation
** More: Yes for Info.improvements and Eat more for both herring and salmon recommendation
** Inconsistent: All other combinations

There are also more technical scenarios: what if nobody eats fish more than 3 g/d on average, what time points will be considered, what pollutants will be considered, and will other sources of nutrients and pollutants be considered as primary or will fish be considered as primary source. The ordering in the latter scenario is important if exposure-response functions are non-linear, as is case with e.g. omega-3 fatty acids. The marginal omega-3 benefits are highest with low exposures, and marginal benefit becomes smaller as the exposure increases. Thus, the health impacts of a primary source are considered larger than those of a secondary source.

{{hidden|
<t2b name="Decisions" index="Num,Decision maker,Decision,Option,Variable,Cell,Change,Description" unit="-">
0||Adjust|BAU|exposure||Add|Add something to avoid zero|0.01
0||Adjust|BAU|incidence||Multiply|#/100000py --> #/py|0.00001
0|HELCOM|Info.improvements|BAU|effinfo||Replace|Current availability of and chemicals in herring and salmon|0
1|HELCOM|Info.improvements|Yes|effinfo||Identity|Better availability of and less chemicals in herring and salmon|
0|National food safety authority|Recomm.herring|BAU|effrecommRaw|Fish:Herring|Replace|Authorities recommendations about Baltic herring do not change|0
2|National food safety authority|Recomm.herring|Eat more|effrecommRaw|Fish:Herring;Recommendation:Eat less|Replace|Authorities recommend to eat more Baltic herring|0
2|National food safety authority|Recomm.herring|Eat less|effrecommRaw|Fish:Herring;Recommendation:Eat more|Replace|Authorities recommend to eat less Baltic herring|0
0|National food safety authority|Recomm.salmon|BAU|effrecommRaw|Fish:Salmon|Replace|Authorities recommendations about Baltic salmon do not change|0
3|National food safety authority|Recomm.salmon|Eat more|effrecommRaw|Fish:Salmon;Recommendation:Eat less|Replace|Authorities recommend to eat more Baltic salmon|0
3|National food safety authority|Recomm.salmon|Eat less|effrecommRaw|Fish:Salmon;Recommendation:Eat more|Replace|Authorities recommend to eat less Baltic salmon|0
0|HELCOM|Select.size|BAU|size|Fish:Baltic herring;Lower:70|Replace|Select herring size for consumption: only small, only large, or both|0
0|HELCOM|Select.size|BAU|size|Fish:Baltic herring;Lower:120|Replace|Eat only large herring (>17 cm)|0
4|HELCOM|Select.size|Ban large|size|Fish:Baltic herring;Lower:170|Replace|Eat only small herring (<17 cm)|0
4|HELCOM|Select.size|Ban large|size|Fish:Baltic herring;Lower:220|Replace|Eat only small herring|0
4|HELCOM|Select.size|Ban large|size|Fish:Baltic herring;Lower:270|Replace|Eat only small herring|0
4|HELCOM|Select.size|New products|size|Fish:Baltic herring|Identity|Eat both small and large herring|
0|National food safety authority|Background|Yes|expo_bg||Identity|Consider background from other sources|
5|National food safety authority|Background|No|expo_bg||Replace|Do not consider background from other sources|0
0|National food safety authority|Limit|No limit|amount||Identity|Do not limit fish concumption of heavy eaters|
6|National food safety authority|Limit|Max 3 g/d|amount||Identity|This limits consumption to max 3 g per day (requires also a function in code)|
0|Risk assessor|Time|2018|time||Replace|Short timeline (2018 only)|2018
7|Risk assessor|Time|2010|time||Replace|Long timeline (1990, 2000, 2010, 2018)|2010
7|Risk assessor|Time|2000|time||Replace|Long timeline|2000
7|Risk assessor|Time|1990|time||Replace|Long timeline|1990
0|Risk assessor|Mixtures|BAU|conc||Identity|All compounds as mixture in the model|
8|Risk assessor|Mixtures|MeHg|conc|Exposure_agent:MeHg|Replace|Consider situation without MeHg|0
8|Risk assessor|Mixtures|Dioxin|conc|Exposure_agent:PCDDF|Replace|Consider situation without PCDDF and PCB|0
8|Risk assessor|Mixtures|Dioxin|conc|Exposure_agent:PCB|Replace|Consider situation without PCDDF and PCB|0
8|Risk assessor|Mixtures|Dioxin|conc|Exposure_agent:TEQ|Replace|Consider situation without PCDDF and PCB|0
8|Risk assessor|Mixtures|Omega3|conc|Exposure_agent:DHA|Replace|Consider situation without DHA and EPA|0
8|Risk assessor|Mixtures|Omega3|conc|Exposure_agent:EPA|Replace|Consider situation without DHA and EPA|0
8|Risk assessor|Mixtures|Omega3|conc|Exposure_agent:Omega3|Replace|Consider situation without DHA and EPA|0
8|Risk assessor|Mixtures|Vitamin D|conc|Exposure_agent:Vitamin D|Replace|Consider situation without MeHg|0
9|Risk assessor|Time|2009|time||Replace|Same time as EU Fish2 study (2009)|2009
9|Risk assessor|Select.size|All sizes|size|Fish:Baltic herring|Identity|Size fractions to include EU Fish 2 study|
0|National food safety authority|Cons.policy|BAU|effinfo||Replace|Current availability of and chemicals in herring and salmon|0
10|National food safety authority|Cons.policy|Increase|effinfo||Identity|Better availability of and less chemicals in herring and salmon|
10|National food safety authority|Cons.policy|Reduce|effinfo||Replace|Current availability of and chemicals in herring and salmon|0
0|National food safety authority|Cons.policy|BAU|effrecommRaw|Fish:Herring|Replace|Authorities recommendations about Baltic herring do not change|0
10|National food safety authority|Cons.policy|Increase|effrecommRaw|Fish:Herring;Recommendation:Eat less|Replace|Authorities recommend to eat more Baltic herring|0
10|National food safety authority|Cons.policy|Reduce|effrecommRaw|Fish:Herring;Recommendation:Eat more|Replace|Authorities recommend to eat less Baltic herring|0
0|National food safety authority|Cons.policy|BAU|effrecommRaw|Fish:Salmon|Replace|Authorities recommendations about Baltic salmon do not change|0
10|National food safety authority|Cons.policy|Increase|effrecommRaw|Fish:Salmon;Recommendation:Eat less|Replace|Authorities recommend to eat more Baltic salmon|0
10|National food safety authority|Cons.policy|Reduce|effrecommRaw|Fish:Salmon;Recommendation:Eat more|Replace|Authorities recommend to eat less Baltic salmon|0
0|LUKE|Luke.scaling|BAU|amount||Identity|Use the Goherr fish consumption study data directly|
11|LUKE|Luke.scaling|Scale|amount|Fish:Salmon|Multiply|Average amount is scaled to actual Luke data (0.07 kg/a)|0.12
11|LUKE|Luke.scaling|Scale|amount|Fish:Herring|Multiply|Average amount is scaled to actual Luke data (0.31 kg/a)|0.22
12||OmegaERF|Cohen|ERF|Response:CHD2 mortality;Observation:ERF;ER_function:RR|Replace|ERF from Cohen et al 2005|1
0||OmegaERF|Cochrane|ERF|Response:CHD2 mortality;Observation:ERF;ER_function:Relative Hill|Replace|ERF from Cochrane 2018|0
12||OmegaERF|Cohen|ERF|Response:Stroke mortality;Observation:ERF;ER_function:Relative Hill|Identity|ERF from Cohen et al 2005|
0||OmegaERF|Cochrane|ERF|Response:Stroke mortality;Observation:ERF;ER_function:Relative Hill|Replace|ERF from Cochrane 2018|0
</t2b>

<t2b name="Marginals" index="Variable,Index,Probs,Function" obs="dummy" desc="Description" unit="-">
effrecomm|Recommendation||sum|1|Change the increase/reduce to increase/BAU/reduce index of Cons.policy
amount|Ages, Gender, oftenSource, muchSource, oftensideSource, muchsideSource, amountRawSource effinfoSource, effrecommSource, effrecommRawSource, infoSource||sum|1|Remove redundant. 
conc|infoSource, conc_vitSource, conc_pcddfResult||sum|1|Remove redundant columns. Somehow model thinks conc_pcddfResult is marginal
expo_dir|infoSource, expo_bgSource, concSource||sum|1|Remove Source columns to save memory
expo_indir|f_ingSource, t0.5Source,f_mtocSource, BFSource||sum|1|Remove redundant
PAF|Exposcen, Adjust, populationSource, infoSource, frexposedSource, populationSource, RRSource, exposureSource, BWSource, incidenceSource||sum|1|Remove redundant columns
BoD|BoDtSource, incidenceSource, populationSource, case_burdenSource, InpBoDSource||sum|1|Remove redundant columns
BoDattr|Response, Adjust, BoDSource, ERFSource, doseSource, PAFSource, infoSource, populationSource||sum|1|Remove redundant columns
</t2b>
}}

=== Timing ===
The assessment started in April 2015. The first stakeholder meeting was in February 2016 and the second in November 2017. The final results with the full assessment were finalised in June 2018.

== Answer ==
=== Results ===

The results of this assessment have been published in a peer-reviewed scientific journal.<ref name="tuomisto2020"/>

;NOTE!: Not all of the following graphs are from the final model version (run at September, 2019). Please check the date of the graph before making conclusions.

<gallery heights=300 widths=400>
Goherr benefit-risk assessment fig1.svg|Lengths of all fish species studied in EU Fish2 study. See Fig 2 for Baltic herring and salmon.
Goherr benefit-risk assessment fig2.svg|Fish lengths may vary a lot in different scenarios for Baltic herring. Goherr model is well in line with actual measurements in EU Fish2.
Goherr benefit-risk assessment fig3.svg|Dioxin concentrations have reduced a lot since 1990. The same trend have existed since the 1970's (data not shown).
Goherr benefit-risk assessment fig4.svg|Baltic herring has larger variation in concentration because it depends on the fish size. Nowadays half or even more of Baltic herring and salmon is below the EU limit value.
Goherr benefit-risk assessment fig5.svg|With active policy to not use large (>17 cm) Baltic herring, the concentrations on the plate would be even lower than nowadays.
Goherr benefit-risk assessment fig6.svg|Dioxin concentrations are strongly correlated with the size of Baltic herring, while the impact is much lower in salmon.
Goherr benefit-risk assessment fig7.svg|Limiting fish use to 3 g/d per species would theoretically reduce pollutant exposures a lot and have an impact on threshold endpoints such as dioxin TWI (tolerable weekly intake). However, this policy is not realistic (see Figs 8-11). It is an important value judgement whether fish is seen as a secondary or primary source of nutrients, i.e. whether to calculate background intake or not, respectively, due to non-linear dose-responses like with stroke and heart disease; for health impacts, see Fig 27.
Goherr benefit-risk assessment fig8.svg|Baltic herring consumption shows several trends: men eat more than women; the old eat more than the young; Estonians eat more than others. For salmon, Danes and Estonians report much more consumption than what is available on the market; clearly they don't know the origin of their salmon.
Goherr benefit-risk assessment fig9.svg|Consumption policies (improved information, recommendations) seem effective in increasing Baltic herring (and somewhat salmon) consumption. However, they are not effective in reducing consumption as only a small fraction obeys this recommendation. See also Figs 10 and 11.
Goherr benefit-risk assessment fig10.svg|Individual change in consumption depends on previous consumption but not much on population subgroup. Also, although some people obey the recommendation to reduce consumption, almost an equal amount does the opposite, and most do not change.
Goherr benefit-risk assessment fig11.svg|Fig 10 is about people consuming large amounts of fish. The exactly same phenomenon is visible in subgroups that eat moderate amounts (i.e. 5 g/d or less).
Goherr benefit-risk assessment fig12.svg|There is large individual variation (almost hundredfold) in fish consumption in most subgroups. There is also a large fraction of people who do not eat these fishes at all (25-90 % with Baltic herring, 70-90 % with Baltic salmon).
Goherr benefit-risk assessment fig13.svg|This is a detail from Fig 12 showing the subgroups that eat fish the most.
Goherr benefit-risk assessment fig14.svg|Most people are below the recommended (vitamin D) or tolerable (methylmercury, dioxin) daily intake. The differences between countries are clearly larger with dioxin (TEQ) than with methylmercury (MeHg).
Goherr benefit-risk assessment fig15.svg|Dioxin exposures from Fig 14 are shown by subgroup and country. Few young females exceed the tolerable daily intake any more, as both concentrations and consumption have reduced.
Goherr benefit-risk assessment fig16.svg|Average exposures do not go beyond tolerable limits of dioxin except with the most active increasing policy in Denmark and Estonia (which countries also over-reported salmon intake). Methylmercury exposures are far from the limit, but still it has significant health impacts (see Fig 26)
Goherr benefit-risk assessment fig17.svg|Exposure-response functions of the health impacts considered in this assessment. With yes/no impacts no (0) means benefit (or no harm) and yes (1) means harm (or no benefit).
Goherr benefit-risk assessment fig18.svg|Figures with numbers of cases must be interpreted with caution. The model unrealistically assumes that the whole population subgroup behaves like one individual, so only the average numbers are realistic. However, these figures are useful in demonstrating how the risks are distributed within subgroups: most impacts come to very few people.
Goherr benefit-risk assessment fig19.svg|Same as Fig 18 but the focus is on the upper corner of the figure, i.e. those with high impacts.
Goherr benefit-risk assessment fig20.svg|Heart disease and stroke are risks that actually affect large fractions of the older subgroups.
Goherr benefit-risk assessment fig21.svg|Burden of disease figures have a realistic x axis and are clearer than figures 18-20. Heart disease, stroke, and dioxin TWI (tolerable weekly intake) affect large fraction of population groups and are fairly large impacts, while other endpoints affect only a small fraction of the population. Dioxin TWI risk is very sensitive to the disability weight given.
Goherr benefit-risk assessment fig22.svg|When looking that the net health benefits, it is clear that old age groups benefit a lot from eating fish despite risks. The impacts overall are much smaller in young age groups, and in women the critical issue is effects of child's intelligence quotient (IQ) and tooth defects, not the health impacts to the woman herself.
Goherr benefit-risk assessment fig22b.svg|Outcome of interest using different objectives. The default objective (the main assessment of this article) focusses on total net health effect in the whole population. The second objective focusses on young women only. Tolerable weekly intakes from 2001 and 2018 are converted to DALYs based on the number of people exceeding the guidance value.
Goherr benefit-risk assessment fig23.svg|For the majority of the population subgroups (except in Denmark), the net health impact deviates from zero. In many people (35-75 % of the old age groups and less than 15 % in the young age groups), the impact of Baltic fish is beneficial but for some people (less than 10 % of the subgroup in practically all subgroups except young Estonian women) the risks are sligthly greater.
Goherr benefit-risk assessment fig24.svg|Size selection of Baltic herring clearly reduces dioxin intake and also risks to youg women, but the overall picture of net benefits changes little because it is dominated by omega3 effects, which do not change in these scenarios.
Goherr benefit-risk assessment fig25.svg|Time trends show that in each country, the risks have come down a lot and especially the tooth defect risk has reduced. Indeed, these were observed in Finland in children born in the 1980's<ref name="alaluusua1996">Satu Alaluusua, Pirjo-Liisa Lukinmaa, Terttu Vartiainen, Maija Partanen, Jorma Torppa, Jouko Tuomisto. (1996) Polychlorinated dibenzo-p-dioxins and dibenzofurans via mother's milk may cause developmental defects in the child's teeth. Environmental Toxicology and Pharmacology Volume 1, Issue 3, 15 May 1996, Pages 193-197. {{doi|10.1016/1382-6689(96)00007-5}}</ref> but not later.
Goherr benefit-risk assessment fig26.svg|When one exposure agent is removed at a time, it is possible to see the impacts on endpoints that are affected by several compounds, in this case child's IQ. Without methylmercury in fish, the health benefits of omega3 would be roughly double than what they are now. Note that the scales on the y axes vary.
Goherr benefit-risk assessment fig27.svg|If background is not considered, i.e. fish is seen as the primary source of omega3, the health benefits are clearly larger than in the business-as-usual scenario of this model. It is a value judgement to decide whether fish or other sources are considered primary. Limiting fish intake to 3 g/d would reduce risks a lot, but benefits would also clearly decrease.
Goherr benefit-risk assessment fig28.svg|In a bigger picture, Baltic fish and its health hazards are only one of the many environmental health risks. It is not even close to the largest ones, but it may be in the top 10 list.
Goherr benefit-risk assessment fig29.svg|This is a simplified version of Fig 22: it is only about Finland, inconsistent policies have been omitted, and health impacts have been categorised into only three groups: benefits, risks, and TWI dioxin (which is actually not a health endpoint even if we have given it a disability value.)
</gallery>

'''Value of information analyses

For detailed results, see [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=vRYbR54ddT2UsBjg model run on 18.11.2018]. Value of information was looked at in three parts, where a bunch of similar decisions were considered together. In these VOI analyses, infertility was used as the outcome for the sperm concentration effect, while tolerable weeksly intakes (both the current "Dioxin TWI" and the new suggested "TWI 2018") were ignored.

Value of information was calculated for the total burden of disease in a random population subgroup in the four study countries, but using uncertainties for individual people. This approach ensures that value of information is not underestimated, because at population level many uncertainties are smaller than at individual level.

* Select herring size
** There is practically no expected value of perfect information (EVPI) (only 1.5 DALY/a) because ''Ban large'', i.e. switching to small herring is in most cases better than other alternatives. However, also other options are beneficial, and the expected value of including that option is 16 DALY/a.
** If that option is excluded, EVPI increases to 8 DALY/a. 
* Consider background and limit maximal fish intake to 3 g/d are evaluated at the same time.
** EVPI is slightly higher than with herring size, 51 DALY. This is because there is no obvious single decision option to choose.
** Dropping the option Background=No would cost 1880 DALY/a, demonstrating that that is clearly a good choice. However, whether background should be considered or not is not an actionable decision but rather a value judgement about how the situation should be seen. In practice, if you consider background intake (Background=Yes), you ignore a large amount of health benefits from omega3 fatty acids in fish. Some people may say that ignoring it is exactly what you should do because those omega3 fatty acids can easily be received from sources that do not have pollutants (the default in this assessment), while others say that fish and other natural foods are the primary source, and omega3 pills and other food supplements should only be used if undernourished.
** If you always consider background intake, then the model uncertainties decrease, and your EVPI is lower (10 DALY/a). The largest EVPI (42 DALY/a) is obtained when fish intake is not limited to 3 g/d; this is because there is more room for benefits leveling off and relative importance of risks increasing, thus increasing uncertainty to decision making.
* Improved information (including availability and usability of fish) and consumption recommendations.
** EVPI with these decisions is 40 DALY/a, so there is some uncertainty about what to do.
** The most important decision option is to increase information and fish availability (145 DALY/a), while any of the other options can be excluded without much change in expected value.

In a previous analysis we used tolerable weeksly intake instead of infertility ([http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=cPBeTD3VkgPcW0yV model run on 20.4.2018, data not shown]). The disability weight used for tolerable weekly dioxin intake is highly uncertain in the model (hundredfold uncertainty 0.0001 - 0.01 DALY/case of exceedance). Therefore, presumably it would be very important to know the actual value that the society wants to allocate to this impact. But actually it is not, as knowing the value has expected value of partial perfect information (EVPPI) of only 12 DALY/a. The reason for this seems to be that this disability weight rarely becomes so high that a decision maker would actually regret fish-promoting policies.

=== Conclusions ===

'''Overall, main arguments from this health assessment and other disciplines studied in Goherr are in favour of getting rid of dioxin-based food restrictions related to Baltic herring and salmon, and promoting human consumption of Baltic fish.

Fish is healthy food, and its use should be promoted. This applies to Baltic herring and salmon as well, even when considering the dioxin and methylmercury concentrations and vulnerable subgroups. This assessment has shown several reasons that support this conclusion:
*	Net health benefits are clear in older age groups with increased risk of cardiovascular diseases.
*	Even in the subgroup of young females, the risks are close to or smaller than benefits.
*	We also considered the tolerable weekly intake (TWI) of dioxin in the assessment. This is actually not a health risk ''per se'', but an indicator that the exposure is approaching levels where actual health harm may occur. However, we did give it a reasonable weight in the model, as we thought that it is something that people don’t want to exceed. Even when it is considered in the assessment, the previous conclusions prevail. If this outcome is ignored, the health risks become so small that there is little uncertainty about the conclusion.
* Dioxin concentrations have decreased dramatically since the 1990's, which was the starting point of our scenarios. In fact, the decrease started even earlier during the 1970's and 1980's when industrial emissions started to improve. The current levels are roughly one tenth of the worst levels. So, the remaining dioxin problem is way smaller than its historical reputation implies.
*	Young women are the risk group but they eat less Baltic herring and salmon than other population subgroups. Therefore, fish promotion policies are likely to increase the health benefits in other subgroups (especially the elderly) much more than health risks in young women.
*	According to the fish consumption survey (Task 5.3), policies promoting Baltic fish consumption seem to encourage people to increase their fish intake substantially. In contrast, discouraging policies were effective only in a minority, while some would paradoxically increase fish intake, and the average change would be negligible.
*	Recommendations to limit Baltic herring intake can be targeted to the risk group of young women only. However, it is possible that as a side effect, these recommendations affect also other age groups even if they would benefit from increasing herring consumption. The other groups could face risks that are larger than benefits in young women, and the overall health impact could be poor. (However, while people responded that such recommendations would not reduce their herring intake, herring consumption has been decreasing for decades in parallel with the discussion on the dioxin risks of herring. Thus, definitive conclusion on this point is not possible.)
*	There are still large uncertainties in both scientific and value-based issues. However, our results show that robust conclusions about decision options can be made with the current information. Also, it seems that there is no single source of additional information that would reveal crucial insights into this issue. In other words, there seems to be little or no value in postponing recommendations in the hope that we would, in the near future, learn something that would help us make wiser decisions.

The conclusions above are based on health considerations only. If we include economic, cultural, and food security aspects into this consideration, we can argue that:
* The price of Baltic herring for food is clearly higher than price of herring for feed. Therefore, for fishing industry it would be beneficial to start catching more fish for human consumption and therefore effects are synergistic with those of health.
* In many areas around the Baltic Sea, Baltic herring has high cultural value as a food and as a tradition. An increase in consumption would help maintaining the culture, and vice versa.
* Baltic Sea could be an important source of human food, but  currently it is rather a source of animal feed. If dioxin-based food regulations were abandoned, it would be easier to develop the Baltic Sea as a food reserve.

== Rationale ==

{{assessment|moderator=Arja}}
[[image:Goherr_WP5_HIA_structure.JPG|thumb|Schematic picture of the health benefit-risk model for Baltic herring and salmon intake.]]
[[File:BRA of Baltic herring and salmon.png|400px|thumb|Detailed modelling diagram of the health benefit-risk model. Green nodes are original data, red nodes are based on scientific literature, and blue nodes are computational nodes. Those with a number are generic nodes designed to be used in several assessments. the number refers to the page identifier in Op_en wiki (this Opasnet wiki). Nodes with red borders are places where it is possible to aggregate an assessment of random individuals to population-level examination. At those points, individual variation disappears and only population-level uncertainty remains in the model. This aggregation is done with mc2d function, but it is NOT used in the default model runs and therefore all results currently presented on this page are based on individual-level assessments.]]

=== Stakeholders ===

Who will be affected by the decisions?
* Professional fishers and their organisations
* Anglers and other recreational fishers and their organisations
* Land owner fishers and their organisation (they have inherited rights)
* Saami people also have inherited rights (in Sweden only?)
* Food and fish industry.
* Mink and fox farmers
* All people utilising recreational values of the Baltic Sea (relates to eutrophication)
* Farmers and agricultural sector
* Consumers of fish
* Producers of fish oil and fish meal
* Baltic Sea RAC (Regional Advisory Concil) represents fishers (located in Copenhagen)
* Hydropower plant owners

Who will affect the decisions?
* Food safety organisations (EVIRA, Livsmedelsverket, ...)
* Ministries (of agriculture, environment, and commerce)
* EU: DG Mare, EU Parliament?
* SWAM (Swedish Agency for Marine and Water Management)

Who are interested in the decisions?
* Scientists
* Environmental NGOs
* Bureaucrats 
* ICES (International Council for the Exploration of the Sea)
* Helcom

=== Dependencies ===

The assessment model is implemented in a modular way in Opasnet. In practice, this means that the data and code used for different parts of the model is located at different pages in Opasnet. In this section, we give the overview and links to the module pages, and all details can be found from there. The whole idea of Opasnet is that different modules can be used in different assessments simultaneously, and that updates in any module are fully reflected in all assessments when they are rerun. 

The assessment of '''[[Goherr:_Fish_consumption_study |consumption of fish]]''' was based on a survey conducted by Taloustutkimus Ltd in the four study countries. Around 500 adults were recruited from each country, and consumption of fish, especially that of Baltic salmon and herring, was asked. Also, we asked reasons for eating or not eating fish and factors that would increase or reduce fish consumption. Gender and age (18-45 or 45+) were separated in the model.

'''[[EU-kalat|Concentrations of dioxins and PCBs in fish]]''' were based on EU Fish 2 study conducted by THL in 2009. The results were based on pooled and individual fish samples (98 Baltic herring and 9 salmon samples) and  analysed for 17 dioxin and 37 PCB congeners. Toxic equivalent quantities (TEQ) were used by multiplying each congener concentration with its potential to induce dioxin-like effects ([[Toxic equivalency factor|toxic equivalency factor]], TEF) and summing up. Size-specific concentration distributions were estimated for salmon and herring, and dioxin and PCB TEQs using linear regression and hierarchical Bayesian modelling using R (version 3.4.3 and JAGS package, http://cran.r-project.org/). Herring sizes and dioxin concentrations in different scenarios came from the fish growth model by SLU applied in BONUS GOHERR project; those results are published elsewhere<ref name="pihlajamaki2019">Mia Pihlajamäki, Arja Asikainen, Suvi Ignatius, Päivi Haapasaari and Jouni T. Tuomisto. Forage Fish as Food: Consumer Perceptions on Baltic Herring. Sustainability 2019, 11(16), 4298; https://doi.org/10.3390/su11164298</ref>

'''Other concentrations'''. [[Concentrations of beneficial nutrients in fish]] were based on published data and a dataset obtained from the Finnish Food Safety Authority Evira. [[Mercury concentrations in fish in Finland]] were based on Kerty database produced by the Finnish Environment Institute.

'''Exposures to pollutants and nutrients''' were simply products of consumption amounts and concentrations in the consumed fish. An exception to this were the infant's exposures to [[Infant's dioxin exposure|dioxin]] and [[ERF of methylmercury#Data|methylmercury]] during pregnancy and breast-feeding, as they were derived from the mother's exposure using simple toxicokinetic models. 

'''[[Exposure-response function]]s''' were derived for all relevant pollutants and nutrients. [[ERF_of_dioxin | Exposure-response functions of dioxins]] were derived for several endpoints. Tooth defects were based on an epidemiological study in Finland<ref name="alaluusua1996"/>. Cancer morbidity was based on U.S.EPA dioxin risk assessment<ref>U.S.EPA. Guidance for Assessing Chemical Contaminant Data for Use in Fish Advisory. Volume 2: Risk Assessment and Fish Consumption Limits, 3rd Edition. 2000. Table 3-1.</ref>. Tolerable daily intake was based on EC Scientific Committee on Food recommendation<ref>EC Scientific Committee on Food. (2001) Opinion of the Scientific Committee on Food on the risk assessment of dioxins and dioxin-like PCBs in food. CS/CNTM/DIOXIN/20 final [https://ec.europa.eu/food/sites/food/files/safety/docs/cs_contaminants_catalogue_dioxins_out90_en.pdf]</ref>. [[ERF of omega-3 fatty acids|Exposure-response functions for omega-3 fatty acids]] on coronary heart disease and stroke mortalities were from a previous risk assessment<ref>Cohen, J.T., PhD, Bellinger, D.C, PhD, W.E., MD, Bennett A., and Shaywitz B.A. 2005b. A Quantitative Analysis of Prenatal Intake of n-3 Polyunsaturated Fatty Acids and Cognitive Development. American Journal of Preventive Medicine 2005;29(4):366–374).</ref>. [[ERF_of_methylmercury |Exposure-response function of methylmercury]] on child's intelligence quotient was based on a previous risk assessment<ref>Cohen JT, Bellinger DC, Shaywitz BA. A quantitative analysis of prenatal methyl mercury exposure and cognitive development. Am J Prev Med. 2005 Nov;29(4):353-65.</ref>. [[ERFs of vitamins|Exposure-response function for vitamin D]] was a step function based on the daily intake recommendations for adults in Finland<ref>Finnish Nutrition Recommendations 2014 [https://www.evira.fi/globalassets/vrn/pdf/ravitsemussuositukset_2014_fi_web.3_es-1.pdf]</ref>.

'''[[Burden of disease]]''' was estimated in two alternative ways: if the burden of a particular disease in the target population was known, the [[attributable risk|attributable fraction]] of a particular compound exposure was calculated. If it was not known, the excess number of cases due to the exposure was estimated using [[health impact assessment]], and this was multiplied by the years under disease per case and the disability weight of the disease. If the exposure-response function was relative to background risk of the disease, [[disease risk]]s from Finland were used for all countries. Population data was from Eurostat (http://ec.europa.eu/eurostat). Disability weights and durations of diseases were based on the estimates from the Institute for Health Metrics and Evaluation (https://healtdata.org), adjusted using author judgement when appropriate estimates were not available.

=== Model parameters ===

We assume that the background exposure is a uniform distribution between zero and the average Finnish intake for nutrients (according to the Finriski study). The typical nutrient intake from Baltic herring is subtracted from the average to avoid double counting. 

{| class="wikitable collapsible collapsed"
! Show details
|----
|

<t2b name="Background exposure" index="Country,Gender,Exposure_agent" desc="Unit,Description" unit="-">
FI|Male|Vitamin D|0 - 11.7|µg /d|Finriski 12 - 0.3 silakasta 
SE|Male|Vitamin D|0 - 11.7|µg /d|Finriski 12 - 0.3 silakasta 
EE|Male|Vitamin D|0 - 11.7|µg /d|Finriski 12 - 0.3 silakasta
DK|Male|Vitamin D|0 - 11.7|µg /d|Finriski 12 - 0.3 silakasta
|Female|Vitamin D|0 - 8.5|µg /d|Finriski 8.7 - 0.2 silakasta
|Male|EPA|0 - 120|mg /d|Finriski 125 - 4.6 silakasta
|Female|EPA|0 - 96|mg /d|Finriski 100 - 3.9 silakasta
|Male|DHA|0 - 118|mg /d|Finriski 125 - 6.7 silakasta
|Female|DHA|0 - 94|mg /d|Finriski 100 - 5.4 silakasta
||ALA|0|mg /d|
||PCDDF|0|pg /d (TEQ)|
||PCB|0|pg /d (TEQ)|
||MeHg|0|µg /d|
||logTEQ|0|log(pg /g)|
||Fish|0|g /d|
</t2b>

<t2b name="Exposure-response functions of interest" index="Exposure_agent,Resp,Response,ER_function,Scaling" obs="Dummy" unit="-">
TEQ|Tooth defect|Yes or no dental defect|ERS|None|1
TEQ|Cancer|Cancer morbidity|CSF|BW|1
TEQ|TWI 2001|Dioxin recommendation tolerable daily intake|TDI|BW|1
TEQ|TWI 2018|Dioxin recommendation tolerable daily intake 2018|TDI|BW|1
TEQ|Infertility|Sperm concentration|ERS|None|1
DHA|Child's IQ|Loss in child's IQ points|ERS|None|1
Omega3|Heart (CHD)|CHD2 mortality|Relative Hill|None|1
Omega3|Heart (CHD)|CHD2 mortality|RR|None|1
ALA|Heart (CHD)|CHD2 mortality|RR|None|1
Vitamin D|Vitamin D intake|Vitamin D recommendation|Step|None|1
MeHg|Child's IQ|Loss in child's IQ points|ERS|BW|1
Omega3|Cancer|Breast cancer|RR|None|1
Fish|Depression|Depression|RR|None|1
Fish|Mortality|All-cause mortality|RR|None|1
</t2b>

<t2b name="DALYs of responses" index="Response" obs="DALY" desc="Description" unit="DALY /case">
CHD2 mortality|5 - 15|Assumes disability weight 1 and duration 10 a (all durations are in years) with 50% uncertainty
Stroke mortality|5 - 15|Assumes disability weight 1 and duration 10 a with 50 % uncertainty
Yes or no dental defect|0 - 0.12|disability weight 0.001 and duration 60 a with 100 % uncertainty. 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?
Cancer morbidity|0.3937391 (0.3566650 - 0.4356150)|IHME estimate for breast cancer (the most important cancer in this analysis): 19.7 (95% CI 21.8 - 17.8) YLL/case. This comes from a lifetime exposure, so it is (linearly( assumed that 1/50 of this is caused by one-year exposure.
Vitamin D recommendation|0.0001 - 0.0101|disability weight 0.001 and duration 1 a with 101-fold uncertainty. 1 if not met
Dioxin recommendation tolerable daily intake|0.0001 - 0.0101|disability weight 0.001 and duration 1 a with 101-fold uncertainty. 1 if not met
Dioxin recommendation tolerable daily intake 2018|0.0001 - 0.0101|disability weight 0.001 and duration 1 a with 101-fold uncertainty. 1 if not met
Sperm concentration|0-5|disability weight 0.1 and duration 50 a in boys, with 2-fold uncertainty. See ERF of dioxin
Loss in child's IQ points|0.11 (0.06 - 0.16)|Mild intellectual disability (IQ<70) has disability weight 0.043 (95 % CI 0.026-0.064) based on IHME. This is scaled to one IQ point with duration 75 a. 
</t2b>

Disability weights of IHME<ref name="disabilityweight">Institute for Health Metrics and Evaluation. (2019) Disability weights for GBD2017 study. http://ghdx.healthdata.org/record/global-burden-disease-study-2017-gbd-2017-disability-weights. Accessed 27 March 2019.</ref>.

''Lower means the lower bound of the size group.

<t2b name="Size distribution of fish species" index="Fish,Lower" obs="Fraction" unit="mm,fraction">
Baltic herring|70|0.112540
Baltic herring|120|0.070343
Baltic herring|170|0.344871
Baltic herring|220|0.472215
Baltic herring|270|0.000031
Baltic herring|320|0
Salmon|700|1
Salmon|1000|0
</t2b>

<t2b name="Diox conc fishing areas" index="Area,Fish" desc="Concentrations used for countries" unit="TEQ">
24|Herring|3.08|DK
25-29|Herring|1.95|SE
28|Herring|2.84|
30-31|Herring|2.29|FI (reference: reflects EU-kalat study)
32|Herring|2.45|EE
24|Salmon|16.0|DK
25-29|Salmon|10.9|
28|Salmon|10.9|
30-31|Salmon|	9.3|FI, SE (reference: reflects EU-kalat study)
32|Salmon|10.7|EE
</t2b>

PCDD/F and PCB concentration distributions for each fish species and country are estimated in the following way:
* [[EU-kalat|EU Fish 2]] 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 Fish 2 study, except that the means are scaled based on the area of interest relative to the reference (see table above).

Population data from Eurostat database 
<t2b name="Population" index="Country,Group" unit="n">
DK|Male 18-45|1032266
DK|Male >45|1207260
DK|Female 18-45|1003587
DK|Female >45|1296678
EE|Male 18-45|252485
EE|Male >45|237413
EE|Female 18-45|239174
EE|Female >45|339881
FI|Male 18-45|970016
FI|Male >45|1182808
FI|Female 18-45|921887
FI|Female >45|1339537
SE|Male 18-45|1823042
SE|Male >45|2065334
SE|Female 18-45|1738408
SE|Female >45|2199156
</t2b>

{| {{prettytable}}
|+'''Landing of fish from different areas in Baltic Sea by country.
|----
!rowspan="2"| Fish species and area
!colspan="4"|Country
|----
! 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 %|| || || 
|----
|}

{| {{prettytable}}
|+'''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
|----
|}

{| {{prettytable}}
|+'''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
|----
|}

<t2b name="Model parametres" index="Response,Exposure_agent,Type,Subgroup,Unit" obs="Result" desc="Description" unit="various">
Loss in child's IQ points||incidence|Group:Female 18-45|IQ /100000py|596000|On average, a population has ca. 6 IQ points per person below 100: mean(abs(rnorm(10000, 100,15)-100))/2
Cancer morbidity||incidence||# /100000py|657|Statistics Finland https://tilastot.syoparekisteri.fi/syovat/ applies to all subgroups because cancer is lifetime risk
Dioxin recommendation tolerable daily intake||incidence||# /100000py|11000|Average over the whole Goherr population is 11 %
Dioxin recommendation tolerable daily intake 2018||incidence||# /100000py|32000|Average over the whole Goherr population is 32 %
Yes or no dental defect||incidence|Group:Female 18-45|# /100000py|22400|Alaluusua et al 2004 found 11/49 cases in two lowest groups
Sperm concentration||incidence|Group:Female 18-45|# /100000py|7000|Male infertility rate is 7 % (Wikipedia)
Vitamin D recommendation||incidence||# /100000py|22000|Valsta et al.: 14 and 30 % of Finnish males and females, respectivey, are below average need (7.5 ug/d) http://urn.fi/URN:ISBN:978-952-343-238-3
</t2b>

* Tooth defect disease burden comes also from [[Disease risk#Calculations]]. But also incidence is needed because the ERF is not related to background: Alaluusua et al 2004 found 11/49 cases in two lowest groups

<t2b name="Model dependencies" index="Name,Ident,label,type,Description,Child" obs="Page" unit="-">
addexposure|Op_en7748/hia|additional exposure||Exposure to relevant agents from other sources than Baltic fish||
amount|Op_en7748/hia|consumption amount|key ovariable|Amount of fish consumed|exposure|
amountRaw|Op_en7748/hia|||Amount of fish consumed without background exposure||
assump|Op_en7748/hia|||Assumptions about fish meal sizes and other things||
bgexposure|Op_en7748/hia|background exposure||Subgroup-specific background intakes||
BoD|Op_en7748/hia|disase burden|objective|Total burden of disease (the objective of the assessment)||
BoDcase|Op_en7748/hia|disease burden based on cases||Burden of disease for responses that are estimated using numbers of cases||
BoDpaf|Op_en7748/hia|disase burden based on PAF||Burden of disease for responses that are estimated using population attributable fraction ||
BoDRaw|Op_en7748/hia|||Burden of disease calculated as if people in each subgroup are identical||
BoDt|Op_en7748/hia|observerd disase burden||Disease-specific burden of disease due to all causes||
BW|Op_en7748/hia|body weight||Body weight of an individual (nominal 70 kg)||
casesabs|Op_en7748/hia|cases of absolute ERFs||Number of cases from responses that have background-independent ERFs||
casesrr|Op_en7748/hia|cases of relative ERFs||Number of cases from responses that have background-dependent ERFs||
conc|Op_en7748/hia|concentrations in fish|key ovariable|Concentrations of relevant exposure agents|exposure|
conc_mehg|Op_en7748/hia|||Concentration of methylmercury||
conc_pcddf|Op_en7748/hia|||Concentration of dioxin and PCB||
conc_vit|Op_en7748/hia|||Concentration of vitamin D||
Decamount|Op_en7748/hia|||Decision about limiting intake to 10 g/d or not||
Decbgexposure|Op_en7748/hia|||Decision about considering background exposures from other sources or not||
Deceffinfo|Op_en7748/hia|||Decision about better availability and less chemicals in herring and salmon||
Deceffrecomm|Op_en7748/hia|||Decision about recommending more, less or BAU amount of consumption||
Decsizes|Op_en7748/hia|||Decision about the size range of consumed fish||
disabilityweight|Op_en7748/hia|disability weights||Disability weights for relevant diseases||
dose|Op_en7748/hia|dose|key ovariable|Dose of relevant exposure agents (scaled to body weight if necessary)|BoD|
duration|Op_en7748/hia|duration of disease||Technical ovariable for nominal disease duration (always 1 because burden is in disabilityweight)||
dx.expo.child|Op_en7748/hia|dioxin exposure in child|key ovariable|Child exposure to dioxin during pregnancy and breast feeding||
effinfo|Op_en7748/hia|effectiveness of information||Effectiveness of information on fish consumption||
effrecomm|Op_en7748/hia|effectiveness of recommendations||Effectiveness of recommendations on fish consumption||
ERF|Op_en7748/hia|exposure-response function|key ovariable|Exposure-response functions for relevant agents and responses|BoD|
ERF_diox|Op_en7748/hia|||ERF for dioxin (Ahr-mediated)||
ERF_env|Op_en7748/hia|||ERF for several environmental pollutants (not used in this model)||
ERF_mehg|Op_en7748/hia|||ERF of methyl mercury on child IQ||
ERF_omega3|Op_en7748/hia|||ERF of omega-3 fatty acids on child IQ, coronary heart disease and stroke||
ERF_vit|Op_en7748/hia|||ERF of vitamin D on intake recommendation compliance||
ERFchoice|Op_en7748/hia|||Technical ovariable slicing the relevant ERFs for this model||
expoRaw|Op_en7748/hia|||Exposure to relevant agents, does not include child's exposure||
exposure|Op_en7748/hia|exposure|key ovariable|Exposure to relevant agents, includes child's exposure||
f_ing|Op_en7748/hia|||Fraction of dioxins actually absorbed from food||
f_mtoc|Op_en7748/hia|||Fraction of mother's dioxins entering the child during pregnancy and breast feeding||
frexposed|Op_en7748/hia|fraction of exposed population||Fraction of the exposed population||
Hg|Op_en7748/hia|mercury data||Data of mercury concentrations from Kerty database||Op_en4004
incidence|Op_en7748/hia|incidence||Incidences of relevant diseases (CHD and stroke in this assessment)||
info|Op_en7748/hia|||Technical ovariable with non-marginal information that may drop during calculations (age, gender, country)||
jsp|Op_en7748/hia||key ovariable|Results of the key questions in the Goherr consumer survey|amount|
lengt|Op_en7748/hia|length||Length of the fish consumed (h missing to avoid confusion with function length)||
mc2d|Op_en7748/hia|||Function for 2D Monte Carlo||
mc2dparam|Op_en7748/hia|||Model parameters for 2D Monte Carlo||
much|Op_en7748/hia|how much fish||How much fish is eaten during a meal (Goherr consumer survey)||
muchside|Op_en7748/hia|how much side dishes||How much side meal fish is eaten during a meal (Goherr consumer survey)||
often|Op_en7748/hia|how often fish||How often fish meals are eaten fish (Goherr consumer survey)||
oftenside|Op_en7748/hia|how often side dishes||How often fish is eaten as side dish (Goherr consumer survey)||
population|Op_en7748/hia|population size||Population size in subgroups (country, age, gender)||
RR|Op_en7748/hia|relative risk||Relative risk of relevant responses given the modelled dose||
samps.j|Op_en7748/hia|||Parameters from Bayesian model for EU Fish study||
sizes|Op_en7748/hia|||Fractions of fish populations in each size group||
sumExposcen|Op_en7748/hia|||Function for exposure scenarios||
survey1|Op_en7748/hia|consumer survey||Data of the Goherr consumer survey||Op_en7749
t0.5|Op_en7748/hia|dioxin half-life||Half-life of dioxin in mother during breast-feeding (elimination+excretion to milk)||
threshold|Op_en7748/hia|threshold||Threshold for ERF of all relevant agents||
threshold_diox|Op_en7748/hia|||Threshold for dioxin ERF||
threshold_env|Op_en7748/hia|||Threshold for ERF of several pollutants (not used in this model)||
threshold_mehg|Op_en7748/hia|||Threshold for mrthylmercury ERF||
threshold_omega3|Op_en7748/hia|||Threshold for omega-3 fatty acid ERF||
threshold_vit|Op_en7748/hia|||Threshold for ERF of vitamin D||
time|Op_en7748/hia|time||Time points considered in the assessment||
vit.param|Op_en7748/hia|||Parameters from Bayesian model for vitamin concentration data||
fig1|Op_en7748/plots|||Fish length in the EU Fish 2 study of Goherr MC model|lengt|Op_en7748
fig2|Op_en7748/plots|||Fish length in the EU Fish 2 study of Goherr MC model, salmon and herring|lengt|Op_en7748
fig3|Op_en7748/plots|||Concentrations of relevant compounds in Baltic fish|conc|Op_en7748
fig4|Op_en7748/plots|||Dioxin in Baltic fish based on EU Fish 2 study or models|conc_pcddf|Op_en7748
fig5|Op_en7748/plots|||PCDD/F in Baltic fish based on EU Fish 2 study of models, salmon and herring|conc|Op_en7748
fig6|Op_en7748/plots|||Fish dioxins vs length by EU Fish 2 data or model|conc_pcddf|Op_en7748
fig7|Op_en7748/plots|||Exposure to compounds in fish by background and limit decisions|exposure|Op_en7748
fig8|Op_en7748/plots|||Consumption of Baltic fish by subgroup and policy|amount|Op_en7748
fig9|Op_en7748/plots|||Fish intake change due to consumption policies|amount|Op_en7748
fig10|Op_en7748/plots|||Change in individual fish intake after all consumption policies|amount|Op_en7748
fig11|Op_en7748/plots|||Change in individual fish intake after all consumption policies|amount|Op_en7748
fig12|Op_en7748/plots|||Consumption of Baltic fish by country and group|amount|Op_en7748
fig13|Op_en7748/plots|||Consumption of Baltic fish by country and group, detail|amount|Op_en7748
fig14|Op_en7748/plots|||Exposure to compounds from Baltic fish and background|exposure|Op_en7748
fig15|Op_en7748/plots|||Exposure to dioxin from Baltic fish|exposure|Op_en7748
fig16|Op_en7748/plots|||Exposure to compounds from fish and background|exposure|Op_en7748
fig17|Op_en7748/plots|||Exposure-response functions for health impacts|ERF|Op_en7748
fig18|Op_en7748/plots|||Numbers of cases for non-relative responses. All individuals are identical in population|casesabs|Op_en7748
fig19|Op_en7748/plots|||Numbers of cases for non-relative responses, detail. All individuals are identical in population|casesabs|Op_en7748
fig20|Op_en7748/plots|||Numbers of cases for relative responses. All individuals are identical in population|casesrr|Op_en7748
fig21|Op_en7748/plots|||Burden of disease by response and group|BoD|Op_en7748
fig22|Op_en7748/plots|||Burden of disease by country, group and policy|BoD|Op_en7748
fig23|Op_en7748/plots|||Net disease burden by country and population group|BoD|Op_en7748
fig24|Op_en7748/plots|||Burden of disease by country, group and herring size selection|BoD|Op_en7748
fig25|Op_en7748/plots|||Burden of disease in young females by country and time|BoD|Op_en7748
fig26|Op_en7748/plots|||Burden of disease by group and one compound removed at a time|BoD|Op_en7748
fig27|Op_en7748/plots|||Burden of disease for subgroups by background and limit|BoD|Op_en7748
fig28|Op_en7748/plots|||Environmental burden of disease in Finland|BoDRaw|Op_en7748
fig29|Op_en7748/plots|||Burden of disease by group and policy in Finland|BoD|Op_en7748
</t2b>
|}

=== Analyses ===

=== Indices ===
* Country (Denmark, Estonia, Finland, Sweden)
* Year (current [concentration data from 2009 but projected to 2018], future)
* Gender (female, male)
* Age (18-45 years, >45 years)
* Fish species (Baltic herring, Baltic salmon)
* Health end-point (coronary heart disease and stroke mortality, tooth defect caused by dioxins, intelligence quotient change in child of a pregnant or nursing woman, exceedance of tolerable daily intake of dioxin, cancer morbidity, and compliance with vitamin D recommendation)
* 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 Dependencies section above. Model results are presented as tables and figures, and the most important ones are shown on this page.

* 20.8.2019: [http://en.opasnet.org/en-opwiki/index.php?title=Benefit-risk_assessment_of_Baltic_herring_and_salmon_intake&oldid=43215#Calculations Archived previous models] that used old HIA model with casesabs and casesrr ovariables. The following codes were removed from the current page:
** 
* 18.5.2017: Archived exposure model ''Op7748/exposure'' by Arja (used separate ovariables for salmon and herring) [http://en.opasnet.org/en-opwiki/index.php?title=Benefit-risk_assessment_of_Baltic_herring_and_salmon_intake&oldid=40413#Exposure_model]
* Sketches about modelling determinants of eating (spring 2018) [http://en.opasnet.org/en-opwiki/index.php?title=Benefit-risk_assessment_of_Baltic_herring_and_salmon_intake&oldid=41947#Sketches_about_modelling_determinants_of_eating]

===== Manuscript graphs =====

* [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=G76irKppw7B3Py5T]
* Model run 29.8.2019 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=d2SV6NU4EOcqa0E6]
* Model run 1.9.2019 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=iwV2cemKKnFBl1uG]
* Model run 10.9.2019 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=pTwPumzOaotRjlRN] [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=A9pSQtOn7PQsAnGI]
* Model run 11.9.2019 official [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=2816VWIKZPtMKBAn]
* Model run 19.9.2019 with agent-specific disease burdens [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=JaizK4lnwxBw631r]
* Model run 24.9.2019 with EPA 2004 CSF. New official [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=HwNhImyTS3iauapg]

<rcode name="pregraph" graphics=1>
# This is code Op_en7748/pregraph on page [[Benefit-risk assessment of Baltic herring and salmon intake]]
library(OpasnetUtils)
library(ggplot2)
library(reshape2)

##### Parameters from user interface

BS <- 30

# This is a fork from code Op_en7919/ on page [[Health effects of Baltic herring and salmon: a benefit-risk assessment]]

############# VOI functions

objects.latest("Op_en2480",code_name="VOI") # [[Value of information]] VOI, binning

################ Fetch all the model ovariables

objects.latest("Op_en7748", code_name="model") # All objects of the model

cat("CollapseMarginals for all runs.\n")
oprint(margi)

##############  Draw graphs

######## Check first that everything is in order

if(FALSE) {
  oempty(all=TRUE) # Empty all ovariable outputs, marginals, and decision marks from openv
  DecisionTableParser(initiate(decsel=NULL))
  InpBoD <- EvalOutput(InpBoD, verbose=TRUE)
  BoDattrPer1 <- EvalOutput(BoDattrPer1,verbose=TRUE)
  
  ### Plot and print each ovariable 
  ovas <- list( # ovariables
    addexposure, amount, amountRaw, assump, BF, BoD, BoDattr, BoDattrPer1, BoDt, BW,
    case_burden, conc, conc.param, conc_mehg, conc_pcddf, conc_vit, dose,
    effinfo, effrecomm, ERF, ERF_diox, ERF_env, ERF_mehg, ERF_micr,
    ERF_omega3, ERF_vit, ERFchoice, expo_bg, expo_dir, expo_indir, exposure,
    f_ing, f_mtoc, frexposed, incidence, info, InpBoD, jsp, lengt, much,
    muchside, often, oftenside, PAF, population, RR, size, t0.5, time
  )
  for(i in 1:length(ovas)) {
    tmp <- ovas[[i]]
    cat("Ovariable", tmp@name, "\n")
    print(tryCatch(oggplot(tmp)+theme_gray(base_size=16), error=function(e1) cat("Cannot plot ",tmp@name)))
    oprint(tryCatch(summary(tmp), error=function(e1) cat("Cannot print summary for ",tmp@name)))
  }
}

########## Concentration

# Long timeline 
# Initiate the model: update decisions and set N.

oempty(all=TRUE)
DecisionTableParser(initiate(decsel=c(7)))
conc <- CheckCollapse(EvalOutput(conc))

### Figure 3. 

tmp <- conc
#tmp@output <- fillna(tmp@output, "Select.size")
tmp <- (BAUt * tmp[
  tmp$Exposure_agent %in% c("Vitamin D", "MeHg","TEQ","Omega3") &
    tmp$Iter %in% 1:1000,])@output

fig3 <- ggplot(tmp, aes(x=Result, colour=Fish, linetype=Time))+
  stat_ecdf(size=1.5)+theme_grey(base_size=BS)+
  facet_wrap(~Exposure_agent, scales="free_x")+
  scale_x_log10(labels=remzero)+
  labs(
    #    title="Concentrations of compounds in Baltic fish",
    x="Concentration (TEQ: pg/g; MeHg, vit D: µg/g; omega3: mg/g [f.w.])",
    y="Cumulative probability"
  )+scale_colour_manual(values=varit12)+
  theme(axis.text.x=(element_text(hjust=1)))
fig3
## Different sizes
# Initiate the model: update decisions.

oempty(all=TRUE)
DecisionTableParser(initiate(decsel=c(4)))
InpBoD <- EvalOutput(InpBoD)
BoDattrPer1 <- EvalOutput(BoDattrPer1)

cat("VOI for BoD and Select.size.\n")
tmp <- oapply(
  BoDattrPer1[!grepl("TWI", BoDattrPer1$Resp) , ],
  INDEX = c("Iter","Select.size","Group","Country"),
  FUN = sum
)

oprint(VOI(
  tmp,
  decision = c("Select.size")
))

##### Amount of consumption

# Initiate the model: update decisions and set N.

#openv.setN(3000)

oempty(all=TRUE)

if(TRUE) { # TRUE: direct study results, FALSE: average amount scaled to match Luke statistics
  tmp <- initiate(decsel=c(10)) # Info.improvement, Recommendations, i.e. Cons.policy #, OmegaERF, Background
} else {
  tmp <- initiate(decsel=c(5,10,11,12)) # Info.improvement, Recommendations, i.e. Cons.policy, OmegaERF, Background, Luke scaling
}
DecisionTableParser(tmp)
InpBoD <- EvalOutput(InpBoD, verbose=FALSE)
BoDattrPer1 <- EvalOutput(BoDattrPer1,verbose=FALSE)

cat("Indices used in BoDattrPer1\n")
tmp <- setdiff(colnames(BoDattrPer1@output)[BoDattrPer1@marginal],"Iter")
oprint(data.frame(
  Index = tmp,
  Locations = sapply(tmp, function(i) paste(unique(BoDattrPer1@output[[i]]),collapse=", "))
))

cat("Disease burden of infertility at individual level (mDALY/a per person)\n")

oprint(oapply( 
  BoDattrPer1[BoDattrPer1$Resp=="Infertility"&BoDattrPer1$Time=="2018",],
  #           info[info$Group=="Female 18-45",],
  c("Resp","Country","Group","Time"),
  mean
)@output)

cat("Disease burden of infertility at population level (DALY/a)\n")

oprint(oapply(
  BoDattr*info[info$Group=="Female 18-45",],
  c("Resp","Country","Group","OmegaERF"),
  mean
)@output)

cat("Background disease burden\n")

oprint(summary(oapply(BoD, NULL, sum, c("Group","BoDSource"))))#, marginals=c("Country","Response")))

amo <- info * amount * 365.25 / 1000 # g/d --> kg/a

cat("Fish consumption average (kg/a).\n")
oprint(oapply(amo*BAU,c("Country","Fish"),mean)@output)
oprint(oapply(amo*BAU,c("Country","Group","Fish"),mean)@output)

if(FALSE) {
  ova <- amo
  out <- aggregate(
    ova$Row,
    by = ova@output[c("Group")],
    FUN = function(x) {
      sample(x, 3000, replace=TRUE)
    }
  )
  #  FUN = function(x) {
  #    apply(array(as.numeric(sample(
  #      as.character(x),
  #      3000,
  #      replace = TRUE
  #    )), dim = c(1, 3000)), MARGIN = 2, FUN = mean)
  #    })
  temp <- melt(out[[length(out)]])
  out[[length(out)]] <- 1:nrow(out)
  colnames(temp) <- c("Nrow", "Iter", "Row")
  out <- merge(out, temp, by.x = "x", by.y = "Nrow")
  out$x <- NULL
  tmp2 <- amo@output[order(amo$Iter),]
  tst <- merge(
    amo@output[
      !duplicated(tmp2[c("Fish","Group","Row","Cons.policy")]) , # There may be several identical Rows on both sides
      colnames(tmp2)!="Iter"],
    out
  )
  out <- Ovariable(output=tst, marginal=colnames(tst) %in% c("Group", colnames(amo@output)[amo@marginal]))
  out <- oapply(out, c("Group","Fish","Cons.policy","Iter"),FUN=sum)
} else { # Endif
  out <- amo
}
amo.diff <- oapply(out*BAU, FUN=sum, cols="Cons.policy")
amo.diff$Amorig <- amo.diff$Result
amo.diff <- out - amo.diff
amo.diff <- amo.diff[amo.diff$Cons.policy %in% c("Increase","Reduce"),]
amo.diff <- amo.diff[amo.diff$Amorig < 5 , ] # Large values not interestng

if(nrow(amo.diff@output)>0) {
  
  # Figure 5.
  
  fig10 <- ggplot(data=amo.diff@output[amo.diff$Iter %in% 1:3000,],#amo.diff@output, # Changes >10 g/d=3.6 kg/a are uninteresting.
                  aes(x=Amorig, y=Result, colour=Group))+
    geom_point()+#data=amo.diff@output[amo.diff$Iter %in% 1:3000,])+
    geom_smooth(method="lm", se=FALSE)+
    theme_gray(base_size=BS)+facet_grid(Fish~Cons.policy)+
    theme(legend.position="bottom")+
    coord_cartesian(ylim=c(-4,8),xlim=c(0,4))+
    labs(
      #      title="Individuals' fish intake after all consumption policies",
      x="Current fish intake (kg/a)",
      y="Change in fish intake (kg/a)"
    )+scale_colour_manual(values=varit12)
  
  # Figure 4. 
  
  fig12 <- ggplot((BAU*amo+0.05)@output, aes(x=Result, colour=Country))+
    stat_ecdf(size=1.5)+
    geom_point(
      data = oapply(BAU*amo, c("Fish","Group","Country"),mean)@output,
      aes(x=Result, y=0.03, colour=Country), shape=1, size=5, stroke=3
    )+
    theme_grey(base_size=BS)+
    theme(legend.position = "bottom")+
    facet_grid(Fish~Group, scales="free_x")+scale_x_log10(labels=remzero)+
    labs(
      #    title="Consumption of Baltic fish by country and subgroup",
      x="Fish consumption of a random individual (kg/a)",
      y="Cumulative probability"
    )+scale_colour_manual(values=varit12)

}

######## Exposure

tmp <- info * unkeep(exposure[
  paste(exposure$Exposure_agent, exposure$Exposure) %in% c(
    "MeHg To child",
    "Vitamin D To eater",
    "Omega3 To eater",
    "TEQ To eater"
  ),],sources=TRUE)

# Figure 6. 

tmp <- BAU*tmp[tmp$Exposure_agent=="TEQ",]*7/70+0.1

fig15 <- ggplot(tmp@output, aes(x = Result, colour=Country))+
  stat_ecdf(size=1.5) + scale_x_log10(labels=remzero) +
  geom_point(
    data = oapply(tmp, c("Group","Country"),mean)@output,
    aes(x=Result, y=0.03, colour=Country), shape=1, size=5, stroke=3
  )+
  theme_gray(base_size = BS)+
  theme(legend.position = "bottom")+
  geom_vline(data=recom[recom$Unit=="pg/kg/week",], aes(xintercept=Result, linetype=Line), colour="red", size=1)+
  facet_wrap(~ Group)+
  labs(
    #    title="Exposure to dioxin from Baltic fish",
    y="Cumulative probability",
    x = "Exposure (pg/kg/week)",
    colour="",
    linetype=""
  )+scale_colour_manual(values=varit12)

#### Burden of disease

varit12[3] <- "#DDC0D3FF" # Make the pink brighter to separate depression from infertility

# Figure 7.

tmp <- oapply(
  BoDattrPer1[!grepl("TWI", BoDattrPer1$Resp) & BoDattrPer1$Background=="Yes" & BoDattrPer1$OmegaERF=="Cochrane",],
  INDEX=c("Resp","Group","Country","Cons.policy"),
  mean
)@output

fig22 <- ggplot(tmp,
                aes(x=Group, weight=BoDattrPer1Result, fill=Resp))+
  geom_bar()+facet_grid(Cons.policy ~ Country)+
  theme_grey(base_size=BS-8)+
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = -90)
  )+
  labs(
    #    title="Disease burden of Baltic fish by country, group, and Cons.policy",
    y = "Disease burden (mDALY/a per person)",
    fill=""
  )+
  scale_fill_manual(values=varit12)

######### Figure 8.

tmp <- oapply(BoDattrPer1 * BAU, INDEX=c("Resp","Country","OmegaERF","Group"),mean)@output

tmp <- rbind(
  cbind(
    Objective="Net health, default",
    tmp[!grepl("TWI",tmp$Resp),]
  ),
  cbind(
    Objective="Net health, young women",
    tmp[!grepl("TWI",tmp$Resp) & tmp$Group=="Female 18-45",]
  ),
  cbind(
    Objective=tmp$Resp,
    tmp
  )[grepl("TWI",tmp$Resp),]
)

sums <- aggregate(tmp["Result"], by=tmp[c("Objective","Country","OmegaERF")],sum)

ggplot(tmp, aes(x=Objective, weight=Result))+
  geom_bar(aes(fill=Resp))+facet_grid(OmegaERF ~ Country)+
  theme_grey(base_size=BS-8)+
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = -90, hjust=0)
  )+
  geom_text(data=sums, size=6, aes(label=round(Result,1),y=pmax(1,Result+1)))+
  labs(
    title="Disease burden using different objectives for all",
    y = "Disease burden (mDALY/a per person)",
    fill=""
  )+
  scale_fill_manual(values=varit12)

fig22b <- ggplot(tmp[tmp$OmegaERF=="Cochrane",], aes(x=Objective, weight=Result))+
  geom_bar(aes(fill=Resp))+facet_grid(. ~ Country)+
  theme_grey(base_size=BS-8)+
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = -90, hjust=0)
  )+
  geom_text(data=sums[sums$OmegaERF=="Cochrane",], size=6, aes(label=round(Result,1),y=pmax(1,Result+1)))+
  labs(
    #    title="Disease burden using different objectives",
    y = "Disease burden (mDALY/a per person)",
    fill=""
  )+
  scale_fill_manual(values=varit12)

if(FALSE) {
  # Why are risks so high in Sweden?
  
  ggplot((info * BoDattr * BAU)@output, aes(x=Result, color=Country))+
    stat_ecdf()+
    facet_wrap(~Resp, scales="free_x")
  
  ggplot((info * BoDattrPer1 * BAU)@output, aes(x=Result, color=Country))+
    stat_ecdf()+
    facet_wrap(~Resp, scales="free_x")
  
  ggplot((info * exposure * BAU)@output, aes(x=Result, color=Country))+
    stat_ecdf()+
    facet_wrap(~Exposure_agent, scales="free_x")
  
  #  oapply(info * exposure * BAU,c("Country","Group","Exposure_agent"),max)@output
  #  oapply(info * amount * BAU,c("Country","Group","Fish"),max)@output
  ## There is young woman in Sweden who reports a lot of both Baltic herring and salmon eating, and her TEQ exposure
  ## is 2215 pg/d while the respective maxima in other countries are <660 pg/d.
  ## This explains the high risks specifically in Sweden: young women don't show much benefit, just risk.
  ## sort((-1*info*exposure[exposure$Exposure_agent=="TEQ",]*BAU)$Result)
  
  ggplot((info*much*BAU)@output,aes(x=Result, color=Country))+stat_ecdf()+facet_grid(Fish~Group)
  oprint(paramtable(info*much*BAU,c("Country","Group","Fish")))
  
  ggplot((info*muchside*BAU)@output,aes(x=Result, color=Country))+stat_ecdf()+facet_grid(Fish~Group)
  oprint(paramtable(info*muchside*BAU,c("Country","Group","Fish")))
}

############### Value of information for Cons.policy

cat("VOI for BoDattr and Cons.policy.\n")
tmp <- oapply(
  info*BoDattr[!grepl("TWI", BoDattr$Resp) & BoDattr$Background=="Yes" & BoDattr$OmegaERF=="Cochrane" , ],
  INDEX = c("Iter","Cons.policy","BoDattrSource","Group","Country"),
  FUN = sum
) 

oprint(VOI(
  tmp,
  decision = c("Cons.policy"),
  scenario = "BoDattrSource"
))

cat("VOI for BoDattr and Cons.policy for separate Groups.\n")
oprint(VOI(
  tmp,
  decision = c("Cons.policy"),
  scenario = "Group"
))

################ Comparison to other environmental health risks

# Disease burden statistics for Finland

a<-opbase.data("Op_fi3944", subset="Tautitaakka Suomessa THL:n tutkimuksessa") # [[Tautitaakka Suomessa]]
a$Result <- as.numeric(as.character(a$Result))
a <- a[order(-a$Result) , ]

lab <- c(
  "Fine particles, death",
  "Noise, annoyance+sleep disorder",
  "Indoor radon, cancer",
  "Environmental tobacco smoke, CHD",
  "UV radiation, cancer",
  "Fine particles, chr bronchitis",
  "Fine particles, symptom days",
  "BONUS GOHERR estimate",
  "MeHg in fish, IQ in children",
  "Moisture in homes, asthma",
  "Moisture in homes, respir symptoms",
  "Lead, IQ in children",
  "Environmental tobacco smoke, cancer"
)

levels(a$Impact)[1:12] <- lab[-8]
levels(a$Source)[levels(a$Source)=="Food"] <- "Food source"

# Figure 8.

tmp <- info * BoDattr * BAU

tmp <- oapply(
  tmp[tmp$Country=="FI" & !grepl("TWI", tmp$Resp),],
  cols="Iter",
  FUN=mean
)@output
levels(tmp$Resp)[levels(tmp$Resp) %in% c("Infertility","Tooth defect")] <- "Dioxin risk"
colnames(tmp)[colnames(tmp)=="Resp"] <- "Source"
tmp$Impact <- "BONUS GOHERR estimate"
tmp <- orbind(tmp,a[a$Result >= 380 , ])
tmp$Impact <- factor(tmp$Impact, lab)
tmp$Source <- factor(tmp$Source, sort(levels(tmp$Source)))

fig28 <- ggplot(data=tmp, aes(x=Impact, weight=Result,fill=Source))+
  geom_bar() + 
  coord_flip() + 
  theme_gray(base_size = BS) +
  labs(
    #    title = "Environmental disease burden in Finland", 
    x = "",
    y = "Disease burden (DALY/a in whole population)"
  ) + 
  scale_fill_manual(values = varit12[c(5,6,3,8,14,2,11,1,7,4)]) + #4,5,8,1,6,7,2,3,9)]) + 
  theme(legend.position = c(0.8, 0.7)) +
  guides(fill=guide_legend(title=NULL))+
  theme(axis.title.x = element_text(hjust=1))

################### Calculate ratio of benefits and risks at individual level
# NOTE! THIS IS NOT A RATIO OF MARGINAL BENEFIT THAT WE SHOULD LOOK AT. RATHER, THIS IS ABSOLUTE AVERAGE RISK
# DIFFERENCE BETWEEN TWO GROUPS: WINNERS AND LOSERS.
if(FALSE) {
  ratio <- BoDattr[!grepl("TWI",BoDattr$Resp), ] * BAU * info
  
  ratio$Type <- ifelse(result(ratio)>0,"risk","benefit")
  ratio <- oapply(ratio, c("Type","Country","Iter","Group"),sum)
  
  net <- ratio[ratio$Type=="benefit" , colnames(ratio@output)!="Type"] +
    ratio[ratio$Type=="risk" , colnames(ratio@output)!="Type"]
  net@name <- "net"
  colnames(net@output)[colnames(net@output)=="Result"] <- "netResult"
  
  ratio <- -1 * ratio[ratio$Type=="benefit" , colnames(ratio@output)!="Type"] /
    ratio[ratio$Type=="risk" , colnames(ratio@output)!="Type"]
  ratio@name <- "ratio"
  colnames(ratio@output)[colnames(ratio@output)=="Result"] <- "ratioResult"
  
  ratio <- ratio + net
  ratio$Rank <- rank(ratio$ratioResult)
  
  cat("Benefit/risk ratio calculated at individual level and averaged for Groups and Countries.\n")
  
  oprint(oapply(ratio, c("Group","Country"),mean)@output)
  
  ggplot(ratio@output,aes(x=ratioResult, colour=Group))+stat_ecdf()+
    scale_x_log10()+geom_vline(xintercept=1)
  
  ggplot(ratio@output,aes(x=netResult, colour=Group))+stat_ecdf()
  
  ggplot(ratio@output,aes(x=netResult,y=ratioResult,colour=Group))+
    geom_point()+
    geom_density_2d()+
    scale_y_log10()+geom_hline(yintercept=1)
  #  coord_cartesian(xlim=c(-2,2))
  #  facet_wrap(~Country)
  
  ggplot(ratio@output[ratio$Group=="Female 18-45",][1:100,],
         aes(x=ratioResult,y=Rank,colour=abs(netResult)))+
    geom_point(size=3,shape=21)+
    scale_x_log10()+geom_vline(xintercept=1)
}

#############################################

cat("Concentrations in fish: mean (95 % CI)\n")
oprint(paramtable(info * conc,c("Country","Fish","Exposure_agent")))

cat("Exposure to pollutants and nutrients in fish: mean (95 % CI)\n")
oprint(paramtable(info * exposure,c("Country","Exposure_agent")))

cat("Exposure response functions: mean (95 % CI)\n")
oprint(paramtable(ERF[ERF$Observation=="ERF",],c("Exposure_agent","Resp")))

cat("Total burden of disease of selected causes from all risk factors (kDALY/a): mean (95 % CI)\n")
oprint(
  reshape(
    paramtable(
      oapply(BoD/1000,NULL,sum,c("Group","BoDSource")) * ERFchoice,
      c("Country","Resp")
    ),
    timevar="Country", idvar="Resp",v.names="Estimate", direction="wide"
  )
)

cat("Individual attributable burden of disease of fish (mDALY/a per person): mean (95 % CI)\n")
oprint(paramtable(BoDattrPer1*BAU,c("Resp","Group","Country")))

cat("Attributable burden of disease of fish by country (DALY/a): mean (95 % CI)\n")
oprint(
  reshape(
    paramtable(
      info*BoDattr*BAU,
      c("Resp","Exposure_agent","Country")
    ),
    timevar="Country", idvar=c("Resp","Exposure_agent"),v.names="Estimate", direction="wide"
  )
)


cat("Burden of disease per case (DALY): mean (95 % CI)\n")
oprint(paramtable(case_burden*ERFchoice,c("Resp")))

################ Store results

images <- c("expotable", "varit12", ls()[grepl("fig",ls())])

if(FALSE) {
  objects.store(list=images)
  cat("Objects", images, "stored.\n")
}

pre <- "" # "Goherr benefit-risk assessment "
post <- ".pdf"
do <- FALSE

fig3
if(do) ggsave(paste0(pre,"Figure 3",post), width=15, height=10.5)
fig12
if(do) ggsave(paste0(pre,"Figure 4",post), width=14, height=10.5)
fig10
if(do) ggsave(paste0(pre,"Figure 5",post), width=11.5, height=10.5)
fig15
if(do) ggsave(paste0(pre,"Figure 6",post), width=11.5, height=10.5)
fig22
if(do) ggsave(paste0(pre,"Figure 7",post), width=11.5, height=10.5)
fig22b
if(do) ggsave(paste0(pre,"Figure 8",post), width=14, height=10.5)
fig28
if(do) ggsave(paste0(pre,"Figure 9",post), width=14, height=10.5)

################ Insight network

if(do) {
  gr <- scrape(type="assessment")
  objects.latest("Op_en3861", "makeGraph") # [[Insight network]]
  gr <- makeGraph(gr)
  export_graph(gr, "ruori.svg")
  render_graph(gr)
}
</rcode>

==== All result graphs ====

<rcode label="Run the model with preselected decisions" graphics=1>
# This is code Op_en7748/ on page [[Benefit-risk assessment of Baltic herring and salmon intake]]
library(OpasnetUtils)
library(ggplot2)
library(reshape2)

##### Parameters from user interface
decsel <- c(5,10)
N <- 1000
groups <- function(o) {return(o)} # Not needed any more but function still shows up in the code

objects.latest("Op_en7748", code_name="model") # All objects of the model

rm(list=ls()[grepl("^Dec", ls())])

# Decisions and marginal collapses
#decs <- deci[deci$Num %in% c("0",decsel) , ] # Choose decisions
#if("10" %in% decsel) {decs <- decs[!decs$Decision %in% c("Info.improvements","Recomm.herring","Recomm.salmon"),]} # Cons.policy #covers these
#DecisionTableParser(decs)

initiate(c(5,10))

cat("Model iterations run:",N,"\n")
cat("Decisions used in model:\n")
oprint(deci)

InpBoD <- EvalOutput(InpBoD)

BoDattr <- EvalOutput(BoDattr)

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

# Disease burden statistics for Finland

a<-opbase.data("Op_fi3944", subset="Tautitaakka Suomessa THL:n tutkimuksessa")
a$Result <- as.numeric(as.character(a$Result))
a <- a[order(-a$Result) , ]

############### Limits for fish concentrations (pg/g fresh weight)

recomconc <- data.frame(
  Exposure_agent=c("PCDDF","TEQ"),
  Compound=c("PCDDF","TEQ"),
  Result=c(3.5,6.5) # http://eur-lex.europa.eu/legal-content/EN/TXT/?uri=CELEX:02006R1881-20140701
)

############ Limits for dioxin intake (ug or pg/d per person)

recomintake <- data.frame(
  Exposure_agent=c("MeHg","Vitamin D","TEQ"),
  Result = c(
    13, # EFSA 2012: 1.3 µg/kg/wk, assuming 70 kg
    7.5, # Recommendation for adults in Finland
    140 # TWI 14 pg /kg /week, assuming 70 kg
  ))

eu3 <- eu2[
  eu2$Fish %in% c("Baltic herring","Salmon") &
    eu2$Compound %in% c("PCDDF","PCB","TEQ") ,
  ]

########### Print basic output

basicoutput <- function() return(oprint(data.frame(
  Parameter = c(
    "Number of iterations: ",
    "Info.improvements: ",
    "Recomm.salmon: ",
    "Recomm.herring: ",
    "Limit: ",
    "Select.size: ",
    "Time: ",
    "Background: ",
    "Mixtures: "
  ),
  Values = c(
    max(as.numeric(as.character(BoD$Iter))),
    paste(unique(BoD$Info.improvements), collapse=", "),
    paste(unique(BoD$Recomm.salmon), collapse=", "),
    paste(unique(BoD$Recomm.herring), collapse=", "),
    paste(unique(BoD$Limit), collapse=", "),
    paste(unique(BoD$Select.size), collapse=", "),
    paste(unique(BoD$Time), collapse=", "),
    paste(unique(BoD$Background), collapse=", "),
    paste(unique(BoD$Mixtures), collapse=", ")
  )
)
))

############# VOI functions

objects.latest("Op_en2480",code_name="VOI") # [[Value of information]] VOI, binning

##############  Draw graphs

############## Fish length

objects.get(runs[1]) # Long timeline 
if(sum(colnames(BoD@output) %in% "BoDResult")==2) { # This is needed for model runs before 20.4.2018
  colnames(BoDRaw@output)[colnames(BoDRaw@output)=="BoDResult"] <- "BoDRawResult"
  BoDRaw@name <- "BoDRaw"
  colnames(BoD@output)[colnames(BoD@output)=="BoDResult"] <- c("BoDRawResult","BoDResult")
}

cat("Basic output for figures 3 - 4, 25\n")
basicoutput()

########## Concentration

fig3 <- ggplot((BAUt*conc)@output, aes(x=concResult, colour=Fish, linetype=Time))+
  stat_ecdf()+theme_grey(base_size=BS)+
  facet_wrap(~Exposure_agent, scales="free_x")+scale_x_log10()+
  labs(
    title="Fig 3. Concentrations of relevant compounds in Baltic fish",
    x="Concentration (dioxin: pg/g; MeHg and vit D: µg/g; omega3: mg/g)",
    y="Cumulative probability"
  )

fig4 <- ggplot()+
  stat_ecdf(data=eu3@output, aes(x=eu2Result, colour="Data 2009"))+
  stat_ecdf(data=melt(exp(samps.j$pred[,,,1])), aes(x=value, colour="Bayes estimate"))+
  stat_ecdf(data=conc_pcddf@output, aes(x=conc_pcddfResult, colour="Goherr estimate", linetype=Time))+
  geom_vline(data=recomconc,aes(xintercept=Result), colour="pink", size=2)+
  theme_grey(base_size=BS)+
  facet_grid(Compound ~ Fish)+scale_x_log10()+
  labs(
    title="Fig 4. Dioxin in Baltic fish based on EU Fish2 study or models",
    y="Cumulative probability",
    x="Concentration (pg/g fw)"
  )

############### Burden of disease with time trend and size selection

BoD <- groups(BoD)

fig25 <- ggplot(oapply(
  BAUt*BoD[BoD$Group =="Female 18-45",],
  c("Group","Time","Resp","Country"),
  mean
)@output,
aes(x=Time, weight=Result, fill=Resp))+
  geom_bar()+facet_grid(.~Country)+
  theme_grey(base_size=BS)+
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = -90)
  )+
  labs(
    title="Fig 25. Burden of disease in young females by country and time",
    y = "Disease burden (mDALY/a per person)"
  )+
  scale_fill_manual(values=varit12)

######## Size selections

objects.get(runs[2]) # Herring size selections
if(sum(colnames(BoD@output) %in% "BoDResult")==2) { # This is needed for model runs before 20.4.2018
  colnames(BoDRaw@output)[colnames(BoDRaw@output)=="BoDResult"] <- "BoDRawResult"
  BoDRaw@name <- "BoDRaw"
  colnames(BoD@output)[colnames(BoD@output)=="BoDResult"] <- c("BoDRawResult","BoDResult")
}
cat("Basic output for figures 1 - 2, 5, 24\n")
basicoutput()

cat("VOI of BoDRaw using infertility rather than TWI by decision Select.size.\n")
tmp <- oapply(BoDRaw[!BoDRaw$Resp %in% c("Dioxin TWI", "TWI 2018"),],c("Iter","Select.size","BoDSource"),sum)
oprint(VOI(tmp, "Select.size",scenario="BoDSource"))

fig1 <- ggplot()+
  stat_ecdf(data=eu2@output, aes(x=Length, colour="EU fish data"))+
  stat_ecdf(data=lengt@output, aes(x=lengtResult, colour=Select.size))+#, colour="Goherr model"))+
  facet_wrap(~Fish, scales="free_x")+theme_gray(base_size=BS)+
  labs(
    title="Fig 1. Fish length in the EU Fish 2 study or Goherr MC model",
    y="Cumulative probability",
    x="Fish length (mm)"
  )

fig2 <- ggplot()+
  stat_ecdf(data=eu2@output[eu2$Fish %in% c("Baltic herring","Salmon"),],
            aes(x=Length, colour="EU fish data"))+
  stat_ecdf(data=lengt@output, aes(x=lengtResult, colour=Select.size))+#, colour="Goherr model"))+
  facet_wrap(~Fish, scales="free_x")+theme_gray(base_size=BS)+
  labs(
    title="Fig 2. Fish length in the EU Fish 2 study or Goherr MC model",
    y="Cumulative probability",
    x="Fish length (mm)"
  )

fig5 <- ggplot()+
  stat_ecdf(data=eu3@output[eu3$Compound=="PCDDF",],
            aes(x=eu2Result, colour="Data 2009"))+
  stat_ecdf(data=melt(exp(samps.j$pred[,1,,1])), 
            aes(x=value, colour="Bayes estimate"))+
  stat_ecdf(data=conc@output[conc$Exposure_agent=="PCDDF",],
            aes(x=concResult, colour="Goherr estimate"))+
  geom_vline(data=recomconc[recomconc$Exposure_agent=="PCDDF",],
             aes(xintercept=Result), colour="pink", size=2)+
  theme_grey(base_size=BS)+
  facet_grid(Select.size ~ Fish)+scale_x_log10()+
  labs(
    title="Fig 5. PCDD/F in Baltic fish based on EU Fish2 study or models",
    y="Cumulative probability",
    x="Concentration (pg/g fw)"
  )

BoD <- groups(BoD)

fig24 <- ggplot(oapply(
  BoD,#[BoD$Time == "2018",],
  INDEX=c("Resp","Group","Country","Select.size"),
  mean)@output,
  aes(x=Select.size, weight=BoDResult, fill=Resp))+
  geom_bar()+facet_grid(Country ~ Group)+
  theme_grey(base_size=BS)+
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = -90)
  )+
  labs(
    title="Fig 24. Burden of disease by country, group and herring size selection",
    y = "Disease burden (mDALY/a per person)"
  )+
  scale_fill_manual(values=varit12)

######### All herring sizes in 2009

objects.get(runs[3]) # Herring size selections and Time==2009
cat("Basic output for figure 6\n")
basicoutput()

temp <- (lengt+conc_pcddf[
  conc_pcddf$Time=="2009" & # Same as EU Fish2 study
    conc_pcddf$Compound!="TEQ" &
    conc_pcddf$Select.size=="All sizes",])@output
if(nrow(temp)>0) {
  fig6 <- ggplot()+
    geom_point(data=temp, cex=0.2, aes(x=lengtResult, y=conc_pcddfResult, colour="Goherr model"))+
    geom_point(data=eu3@output, aes(x=Length, y=eu2Result, colour="EU Fish2 data"),cex=2)+
    geom_smooth(data=eu3@output, aes(x=Length, y=eu2Result, colour="EU Fish2 data"),colour="blue")+
    facet_grid(Compound~Fish, scales="free_x")+scale_y_log10()+
    theme_gray(base_size=BS)+coord_cartesian(ylim=c(0.3,12))+
    labs(
      title="Fig 6. Fish dioxins vs length by EU Fish2 data or model",
      y="Dioxin concentration (pg/g fw)",
      x="Fish length (mm)"
    )
}

######### Exposure with varying background and max intake limit

objects.get(runs[4]) # Background, Limit
if(sum(colnames(BoD@output) %in% "BoDResult")==2) { # This is needed for model runs before 20.4.2018
  colnames(BoDRaw@output)[colnames(BoDRaw@output)=="BoDResult"] <- "BoDRawResult"
  BoDRaw@name <- "BoDRaw"
  colnames(BoD@output)[colnames(BoD@output)=="BoDResult"] <- c("BoDRawResult","BoDResult")
}
cat("Basic output for figure 7, 27.\n")
basicoutput()

cat("VOI for BoDRaw using infertility rather than TWI by decisions Background and Limit.\n")
tmp <- oapply(BoDRaw[!BoDRaw$Resp %in% c("Dioxin TWI", "TWI 2018"),],c("Iter","Background","Limit","BoDSource"),sum)
oprint(VOI(tmp, c("Background","Limit"),scenario="BoDSource"))

tmp <- groups(unkeep(exposure[exposure$Exposure=="To eater",],sources=TRUE) * info)

fig7 <- ggplot(
  oapply(tmp,INDEX=c("Limit","Exposure_agent","Background"),mean)@output,
  aes(x = Limit, weight=Result, fill=Background))+
  geom_bar() + theme_gray(base_size = BS)+
  facet_grid(Exposure_agent ~ Background, scales="free_y")+
  labs(
    title="Fig 7. Exposure to compounds in fish by background and limit",
    y="Exposure (dioxin: pg/d; vit D & MeHg: µg/d; omega3: mg/d)"
  )

BoD <- groups(BoD)

fig27 <- ggplot(
  oapply(BoD,INDEX=c("Limit","Resp","Group","Background"),mean)@output,
  aes(x = Limit, weight=BoDResult, fill=Resp))+
  geom_bar() + theme_gray(base_size = BS)+
  facet_grid(Group ~ Background, scales="free_y")+
  scale_fill_manual(values=varit12)+
  labs(
    title="Fig 27. Burden of disease for subgroups by background and limit",
    y="Disease burden (mDALY/a per person)"
  )

##### Amount of consumption

objects.get(runs[5]) # Info.improvement, Recommendations, i.e. Policy
if(sum(colnames(BoD@output) %in% "BoDResult")==2) { # This is needed for model runs before 20.4.2018
  colnames(BoDRaw@output)[colnames(BoDRaw@output)=="BoDResult"] <- "BoDRawResult"
  BoDRaw@name <- "BoDRaw"
  colnames(BoD@output)[colnames(BoD@output)=="BoDResult"] <- c("BoDRawResult","BoDResult")
}
cat("Basic output for figures 8 - 23, 28, 29.\n")
basicoutput()

cat("VOI for BoDRaw and Policy.\n")
tmp <- oapply(groups(info*BoDRaw[!BoDRaw$Resp %in% c("Dioxin TWI", "TWI 2018"),]),c("Iter","Info.improvements","Recomm.herring",
                                                                                    "Recomm.salmon","Policy","BoDSource"), sum) +
  disabilityweight[disabilityweight$Resp=="Dioxin TWI",]
tmp@name <- "BoDRaw"
oprint(VOI(
  unkeep(tmp,c("Resp","disabilityweightSource")),
  c("Info.improvements","Recomm.herring","Recomm.salmon","Policy"),
  "disabilityweightResult",
  "BoDSource"
))

amount <- groups(amount)

amount.diff <- unkeep(amount*BAU, prevresults=TRUE,sources=TRUE, cols=policies)
amount.diff$Amorig <- amount.diff$Result
amount.diff <- amount - amount.diff
amount.diff <- amount.diff[
  amount.diff$Policy %in% c("Increase","Reduce"),]

#cat("87.5th percentile (top 1/8) of fish eating in different population groups. 
#    There is a large subgroup of young Swedish females who say that they eat a lot of Baltic salmon.
#    (Probably they don't because it is not that much available.)\n")

# 85 % of young female swedes say that they eat more than 13 g/d of Baltic salmon.
# This is unlikely but it is the main reason why net health loss is positive
# in this group.

#oprint(oapply(
#  unkeep(amount*BAU, sources=TRUE),
#  INDEX=c("Country","Group","Fish"),
#  FUN=function(x) quantile(x, 0.875)
#)@output)

fig8 <- ggplot(oapply(amount, INDEX=c("Country","Group","Fish","Policy"),mean)@output,
               aes(x=Policy, weight=amountResult, fill=Fish))+geom_bar()+
  facet_grid(Country~Group)+
  theme_gray(base_size=BS)+
  theme(axis.text.x = element_text(angle=90, hjust=1, vjust=0.5))+
  labs(
    title="Fig 8. Consumption of Baltic fish by subgroup and policy",
    x="Population subgroup",
    y="Average fish consumption (g/d)"
  )

if(nrow(amount.diff@output)>0) {
  fig9 <- ggplot(amount.diff@output, aes(x=Result, colour=Group))+stat_ecdf()+
    theme_grey(base_size=BS)+
    theme(legend.position = "bottom")+
    coord_cartesian(xlim=c(-5,10))+facet_grid(Fish~Policy)+
    labs(
      title="Fig 9. Fish intake change due to consumption policies",
      x="Change in amount consumed by a random individual (g/d)",
      y="Cumulative probability"
    )
}

if(nrow(amount.diff@output)>0) {
  fig10 <- ggplot(amount.diff@output[amount.diff$Amorig<50,], # Changes >50 are unrealistic because intake always <100.
                  aes(x=Amorig, y=Result, colour=Group))+geom_point()+geom_smooth()+
    theme_grey(base_size=BS)+facet_grid(Fish~Policy)+
    coord_cartesian(ylim=c(-40,60))+
    labs(
      title="Fig 10. Change in individual fish intake after all consumption policies",
      x="Current fish intake (g/d)",
      y="Change in fish intake (g/d)"
    )
  
  fig11 <- ggplot(amount.diff@output[amount.diff$Amorig<5,], # Changes >50 are unrealistic because intake always <100.
                  aes(x=Amorig, y=Result, colour=Group))+geom_point()+geom_smooth()+
    theme_grey(base_size=BS)+facet_grid(Fish~Policy)+
    coord_cartesian(ylim=c(-5,6))+
    labs(
      title="Fig 11. Change in individual fish intake after all consumption policies",
      x="Current fish intake (g/d)",
      y="Change in fish intake (g/d)"
    )
}

fig12 <- ggplot((BAU*amount+0.05)@output, aes(x=Result, colour=Country))+stat_ecdf()+
  theme_grey(base_size=BS)+
  facet_grid(Fish~Group, scales="free_x")+scale_x_log10()+
  labs(
    title="Fig 12. Consumption of Baltic fish by country and group",
    x="Fish consumption of a random individual in BAU (g/day)",
    y="Cumulative probability"
  )

fig13 <- ggplot((BAU*amount+0.05)@output, aes(x=Result, colour=Country))+stat_ecdf()+
  theme_grey(base_size=BS)+
  facet_grid(Fish~Group, scales="free_x")+scale_x_log10()+
  labs(
    title="Fig 13. Consumption of Baltic fish by country and group",
    x="Fish consumption of a random individual in BAU (g/day)",
    y="Cumulative probability"
  )+coord_cartesian(xlim=c(10,100),ylim=c(0.8,1))

######## Exposure

tmp <- groups(unkeep(exposure[
  paste(exposure$Exposure_agent, exposure$Exposure) %in% c(
    "MeHg To child","Vitamin D To eater","Omega3 To eater","TEQ To eater"
  ),],sources=TRUE) * info)

fig14 <- ggplot((BAU*tmp+0.1)@output,
                aes(x = Result, colour=Country))+
  stat_ecdf() + scale_x_log10() + theme_gray(base_size = BS)+
  geom_vline(data=recomintake, aes(xintercept=Result), colour="pink", size=2)+
  facet_wrap(~ Exposure_agent, scales="free_x")+
  labs(
    title="Fig 14. Exposure to compounds from Baltic fish and background",
    y="Cumulative probability",
    x = "Exposure (dioxin: pg/d; vit D & MeHg: µg/d; omega3: mg/d)"
  )

fig15 <- ggplot((BAU*tmp[tmp$Exposure_agent=="TEQ",]+0.1)@output,
                aes(x = Result, colour=Country))+
  stat_ecdf() + scale_x_log10() + theme_gray(base_size = BS)+
  geom_vline(data=recomintake[recomintake$Exposure_agent=="TEQ",], aes(xintercept=Result), colour="pink",size=2)+
  facet_wrap(~ Group)+
  labs(
    title="Fig 15. Exposure to dioxin from Baltic fish",
    y="Cumulative probability",
    x = "Exposure (pg/d)"
  )

fig16 <- ggplot(oapply(tmp,INDEX=c("Country","Exposure_agent","Policy"),mean)@output,
                aes(x = Policy, weight=Result, fill=Policy))+
  geom_bar() + theme_gray(base_size = BS)+
  geom_hline(data=recomintake, aes(yintercept=Result), colour="pink",size=2)+
  facet_grid(Exposure_agent ~ Country, scales="free_y")+
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = -90)
  )+
  labs(
    title="Fig 16. Exposure to compounds from fish and background",
    y="Exposure (dioxin: pg/d; vit D & MeHg: µg/d; omega3: mg/d)"
  )

######### Table with exposures

if(TRUE) {
  BAU <- Ovariable(data=data.frame(
    Info.improvements = "BAU",
    Recomm.herring = "BAU",
    Recomm.salmon = "BAU",
    Limit = "No limit",
    Select.size = "BAU",
    Time = 2018,
    Mixtures = "BAU",
    Background = "Yes",
    Result = 1
  ))
} else {
  BAU <- 1
}
mc2dparam$run2d<-TRUE
mc2dparam$newmarginals<-"Country"
mc2dparam$N2 <- 1000

expotable <- mc2d(exposure[exposure$Exposure=="To eater"&exposure$Exposure_agent=="TEQ",]*BAU*info)
#summary(tmp) # Error with colnames?!
expotable <- rbind(
  data.frame(
    oapply(
      expotable,
      c("Exposure_agent","Country","Exposure"),
      FUN=function(x) quantile(x, probs=0.05)
    )@output,
    Prob = 0.05
  ),
  data.frame(
    oapply(
      expotable,
      c("Exposure_agent","Country","Exposure"),
      FUN=function(x) quantile(x, probs=0.25)
    )@output,
    Prob = 0.25
  ),
  data.frame(
    oapply(
      expotable,
      c("Exposure_agent","Country","Exposure"),
      FUN=function(x) quantile(x, probs=0.75)
    )@output,
    Prob=0.75
  ),
  data.frame(
    oapply(
      expotable,
      c("Exposure_agent","Country","Exposure"),
      FUN=function(x) quantile(x, probs=0.95)
    )@output,
    Prob = 0.95
  )
)
expotable$Result <- expotable$Result/70
expotable <- expotable[order(expotable$Country),]
cat("Exposure to dioxin (pg/kg/day in TEQ)")
oprint(expotable)

####### Numbers of cases

fig18 <- ggplot(groups(BAU*casesabs*info)@output, aes(x=Result, colour=Group))+
  stat_ecdf()+theme_grey(base_size=BS-2)+
  facet_grid(Country~ Resp, scales="free_x")+
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = -90)
  )+
  labs(
    title="Fig 18. Numbers of cases for non-relative responses (individual=population)",
    y="Cumulative probability",
    x="# of cases/a (for IQ: each change by 1 point is a case)"
  )

fig19 <- ggplot(groups(BAU*casesabs*info)@output, aes(x=Result, colour=Group))+
  stat_ecdf()+theme_grey(base_size=BS-2)+
  facet_grid(Country~ Resp, scales="free_x")+
  coord_cartesian(ylim=c(0.8,1))+
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = -90)
  )+
  labs(
    title="Fig 19. Numbers of cases for non-relative responses (individual=population)",
    y="Cumulative probability",
    x="# of cases/a (for IQ: each change by 1 point is a case)"
  )

fig20 <- ggplot(groups(BAU*casesrr*info)@output, aes(x=Result, colour=Group))+
  stat_ecdf()+theme_grey(base_size=BS)+
  facet_grid(Country~ Resp)+
  coord_cartesian(xlim=c(-12000,0))+
  theme(
    legend.position = "bottom"
  )+
  labs(
    title="Fig 20. Numbers of cases for relative responses (individual=population)",
    y="Cumulative probability",
    x="# of cases/a "
  )

#### Burden of disease

BoD <- groups(BoD)

fig21 <- ggplot((BAU*BoD)@output,
                aes(x=BoDResult, colour=Group))+stat_ecdf()+
  theme_grey(base_size=BS)+
  facet_wrap(~ Resp, scales = "free_x")+
  labs(
    title="Fig 21. Burden of disease by response and group",
    x="Disease burden (mDALY/a per a random individual)",
    y="Cumulative probability")

fig22 <- ggplot((oapply(BoD, INDEX=c("Resp","Group","Country","Policy"), mean))@output,
                aes(x=Group, weight=BoDResult, fill=Resp))+
  geom_bar()+facet_grid(Country ~ Policy)+
  theme_grey(base_size=BS)+
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = -90)
  )+
  labs(
    title="Fig 22. Burden of disease by country, group and policy",
    y = "Disease burden (mDALY/a per person)"
  )+
  scale_fill_manual(values=varit12)

fig23 <- ggplot((BAU * info * oapply(BoD, cols=c("Resp","Response"), FUN = sum))@output,
                aes(x=Result, colour=Group))+stat_ecdf()+
  theme_grey(base_size=BS)+
  coord_cartesian(xlim=c(-10,10))+facet_wrap(~Country)+
  labs(
    title="Fig 23. Net disease burden by country and population group",
    x="Disease burden (mDALY/a per a random individual)",
    y="Cumulative probability"
  )

tmp <- BoD[BoD$Policy!="Inconsistent" & BoD$Country=="FI",]
tmp$Impact <- ifelse(tmp$Resp %in% c("Tooth defect","Cancer","Dioxin TWI","TWI 2018"),"Risk","Benefit") # Does IQ go right in the code as it may be risk or benefit?
tmp$Impact[tmp$Resp %in% c("Dioxin TWI", "TWI 2018")] <- "TWI dioxin"
fig29 <- ggplot((oapply(tmp, INDEX=c("Resp","Group","Country","Policy","Impact"), mean))@output,
                aes(x=Group, weight=BoDResult, fill=Impact))+
  geom_bar()+facet_wrap( ~ Policy)+
  theme_grey(base_size=BS)+
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = -90)
  )+
  labs(
    title="Fig 29. Burden of disease by group and policy in Finland",
    y = "Disease burden (mDALY/a per person)"
  )+
  scale_fill_manual(values=varit12[c(1,4,3)])

################ Comparison to other environemental health risks

fig28 <- ggplot()+
  geom_bar(data=oapply(BoDRaw[BoDRaw$Country=="FI",]*BAU, cols="Iter",FUN=mean)@output,
           aes(x="Goherr MC model estimate", weight=Result,fill=Resp))+
  geom_bar(data=a[a$Result >= 380 , ], 
           aes(x = Impact, weight = Result, fill = Source)) + 
  coord_flip() + 
  theme_gray(base_size = BS-4) +
  labs(
    title = "Fig 28. Environmental burden of disease in Finland", 
    x = "",
    y = "Disease burden (DALY/a in whole population)"
  ) + 
  scale_fill_manual(values = varit12[c(3,4,6,10,2,9,8,1,5,7,11,12)]) + 
  theme(legend.position = c(0.8, 0.7)) 

####################### Comparison of intake to EFSA tolerable weekly intakes

tmp <- expoRaw[expoRaw$Exposure_agent=="TEQ",] + amount[amount$Fish=="Baltic herring",]
fig30 <- ggplot(tmp@output, aes(x=amountResult*30/100, y=expoRawResult/70*7))+geom_smooth()+
  coord_cartesian(xlim=c(0,7),ylim=c(0,20))+
  geom_hline(yintercept=c(2,14))+
  labs(title="How much can you eat herring and still stay below TWI?",
       x="Consumption of herring (100 g portions/mo)",
       y="Exposure to TEQ (pg/kg/week)")

###################### Mixtures 

objects.get(runs[6]) # Mixtures
if(sum(colnames(BoD@output) %in% "BoDResult")==2) { # This is needed for model runs before 20.4.2018
  colnames(BoDRaw@output)[colnames(BoDRaw@output)=="BoDResult"] <- "BoDRawResult"
  BoDRaw@name <- "BoDRaw"
  colnames(BoD@output)[colnames(BoD@output)=="BoDResult"] <- c("BoDRawResult","BoDResult")
}
cat("Basic output for figure 26.\n")
basicoutput()

BoD <- groups(BoD)

fig26 <- ggplot(oapply(BoD, c("Group","Resp","Mixtures"),mean)@output,
                aes(x=Mixtures, weight=BoDResult, fill=Resp))+
  geom_bar()+facet_wrap( ~ Group, scales="free_y")+
  theme_grey(base_size=BS-2)+
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = -90)
  )+
  labs(
    title="Fig 26. Burden of disease by group and one compound removed at a time",
    y = "Disease burden (mDALY/a per person)",
    x="Which compound is removed from the model?"
  )+
  scale_fill_manual(values=varit12)

######## Exposure-response functions

objects.get(runs[7]) # Exposure-response functions

fig17 <- ggplot(erfs, aes(x=Dose,y=Result,colour=Resp,group=Group,linetype=Expo))+
  geom_line(size=2)+facet_wrap(~Group, scales="free_y")+
  scale_x_log10()+theme_gray(base_size=BS-2)+
  labs(
    title="Fig 17. Exposure-response functions for health impacts",
    y="Relative risk (CHD, stroke) or number of cases",
    x="Dose (Omega3: mg/d, TEQ: pg/kg/d except tooth pg/d, MeHg: µg/kg/d, vitamin D: µg/d)"
  )

#################### VOI for mcd2=TRUE and Policy

objects.get(runs[8]) # Info.improvement, Recommendations, i.e. Policy with mc2d=TRUE

cat("VOI for BodRaw and Policy separately for Gender, Age, Country subgroups.\n")
tmp <- oapply(groups(info*BoDRaw[!BoDRaw$Resp %in% c("Dioxin TWI", "TWI 2018"),]),c("Iter","Info.improvements","Recomm.herring",
                                                                                    "Recomm.salmon","Policy","BoDSource"), sum) +
  disabilityweight[disabilityweight$Resp=="Dioxin TWI",]
tmp@name <- "BoDRaw"
oprint(VOI(
  unkeep(tmp,c("Resp","disabilityweightSource")),
  c("Info.improvements","Recomm.herring","Recomm.salmon","Policy"),
  "disabilityweightResult",
  c("Gender","Age","Country")
))
</rcode>

==== Initiate model ====

* Model run 27.7.2019 with oempty [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=RxDWFi4cGd9Abeo2]
* Model run 21.8.2019 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=ReixLHfE7cAve8uV]
* Model run 29.8.2019 without groups function [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=h1rfkXYfbx11iwhG]
* Model run 4.9.2019 with paramtable [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=UTOZiw5xXcUZMNLk]
* Model run 9.9.2019 with new InpBoD from IHME [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=MsnFNoLDlJpmCJmS]
* Model run 11.9.2019 official [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=1KgSFeuC2e9bYLQ9]
* Model run 23.9.2019 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=IHWhmO0ZDoahnyVP]

<rcode name="model" label="Initiate the whole model" graphics=1>
# This is code Op_en7748/model on page [[Benefit-risk assessment of Baltic herring and salmon intake]]
library(OpasnetUtils)

##### Parameters from user interface
decsel <- NULL # c(5,10)
N <- 10

############### Limits and recommendations

recom <- data.frame(
  Line=c(rep("Concentration limit",2),"TWI 2001","TWI 2018","MeHg TWI","recommended vitamin D"),
  Exposure_agent=c("PCDDF",rep("TEQ",3),"MeHg","Vitamin D"),
  Compound=c("PCDDF",rep("TEQ",3),"MeHg","Vitamin D"),
  Unit=c(rep("pg/g fresh weight",2),rep("pg/kg/week",2),"µg/d","µg/d"),
  Result=c(
    3.5, # http://eur-lex.europa.eu/legal-content/EN/TXT/?uri=CELEX:02006R1881-20140701
    6.5,
    14, # SCF 2001 TWI 14 pg /kg /week
    2,  # EFSA 2018 TWI 2 pg/kg/week
    13, # EFSA 2012: 1.3 µg/kg/wk, assuming 70 kg
    7.5 # Recommendation for adults in Finland
  ) 
)

###### Preprocess (these are typically needed when graphs are plotted)

initiate <- function(decsel =NULL) {
  
  # Decisions and marginal collapses
  rm(list=ls(envir=.GlobalEnv)[grepl("^Dec", ls(envir=.GlobalEnv))], envir=.GlobalEnv)
  decs <- deci[deci$Num %in% c("0",decsel) , ] # Choose decisions
  if("10" %in% decsel) {decs <- decs[!decs$Decision %in% c("Info.improvements","Recomm.herring","Recomm.salmon"),]} # Cons.policy covers these
  
  cat("Decisions used in model:\n")
  oprint(decs)
  
  return(decs)
}

remzero <- function(n){
  gsub("\\.0+$","", format(n, scientific = FALSE))
}

paramtable <- function(ova, INDEX) {
  tmp <- cbind(
    oapply(
      ova,
      INDEX=INDEX,
      FUN=function(x) quantile(x, c(0.025, 0.5, 0.975))
    )@output,
    oapply(
      ova,
      INDEX=INDEX,
      FUN= mean
    )@output
  )
  tmp1 <- gsub("\\.$","",gsub("(0*)$", "", sprintf("%.6f", signif(tmp[[(length(INDEX)+1)*2]], 2))))
  tmp2 <- gsub("\\.$","",gsub("(0*)$", "", sprintf("%.6f", signif(tmp[[length(INDEX)+1]][,1], 2))))
  tmp3 <- gsub("\\.$","",gsub("(0*)$", "", sprintf("%.6f", signif(tmp[[length(INDEX)+1]][,3], 2))))
  tmp$Estimate <- ifelse(tmp1==tmp2 & tmp2==tmp3,tmp1, paste0(tmp1, " (",tmp2, ", ",tmp3,")"))
  return(tmp[c(INDEX,"Estimate")])
}

policies <- c(
  "Info.improvements",
  "Recomm.herring",
  "Recomm.salmon",
  "Cons.policy",
  "Limit",
  "Select.size",
  "Time",
  "Background",
  "Mixtures"
)

BAU <- Ovariable(
  "BAU",
  data = data.frame(
    Cons.policy = "BAU",
    Info.improvements = "BAU",
    Recomm.herring = "BAU",
    Recomm.salmon = "BAU",
    Limit = "No limit",
    Select.size="BAU",
    Time="2018",
    Background="Yes",
    Mixtures="BAU",
    OmegaERF="Cochrane",
    BAUResult = 1
  )
  #  marginal = c(rep(TRUE,8),FALSE)
)

BAUt <- Ovariable(
  "BAU",
  data = data.frame(
    Cons.policy = "BAU",
    Info.improvements = "BAU",
    Recomm.herring = "BAU",
    Recomm.salmon = "BAU",
    Limit = "No limit",
    Select.size="BAU",
    Background="Yes",
    Mixtures="BAU",
    BAUResult = 1
  )
  #  marginal = c(rep(TRUE,8),FALSE)
)

resplev <- c(
  "Mortality",
  "Heart (CHD)",
  "Depression",
  "Vitamin D intake",
  "Cancer",
  "Child's IQ",
  "Infertility",
  "Tooth defect",
  "TWI 2001",
  "TWI 2018"
)

varit12 <- c(
  "#74AF59FF", # green Colours 1-6 from [[Goherr: Fish consumption study]]
  "#2F62ADFF", # blue
  "#D692BDFF", # pink
  "#29A0C1FF", # light blue
  "#BE3F72FF", # red
  "#FBB848FF", # yellow
  "#cc77ac", # roosa # Colours 7-12 from THL old palette?
  "#606060", # tumma harmaa
  "#cc77ac", # roosa 
  "#cc77ac", # roosa 
  "#c3c2c6", # keskiharmaa
  "#dcdfe2", # vaalea harmaa
  "#cc7acc", # 7
  "#244911" # 12
) 

################3# Decisions and marginal collapses
deci <- opbase.data("Op_en7748", subset = "Decisions")
# Decisions and marginal collapses
decs <- deci[deci$Num %in% c("0",decsel) , ] # Choose decisions
if("10" %in% decsel) {decs <- decs[!decs$Decision %in% c("Info.improvements","Recomm.herring","Recomm.salmon"),]} # Cons.policy covers these
DecisionTableParser(decs)

cat("Model iterations run:",N,"\n")
cat("Decisions used in model:\n")
oprint(deci)

margi <- as.data.frame(sapply(opbase.data("Op_en7748", subset = "Marginals"), as.character),stringsAsFactors=FALSE)
CollapseTableParser(margi)

cat("Marginals summed up in model:\n")
oprint(margi)

# For development and testing, use smaller ovariables
openv.setN(N) 

size <- Ovariable("size", ddata="Op_en7748", subset="Size distribution of fish species")

time <- Ovariable(
  "time",
  data = data.frame(Result=c(2018))
)

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]]
      "info"
    ), 
    Ident = c(
      "Op_en3104/initiate",
      "Op_en1838/initiate",
      "Op_en4004/conc_mehg",
      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("Baltic herring", "Salmon"), each = 4),
      Country = c("DK", "EE", "FI", "SE"),
      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=="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("Baltic herring", "Salmon"),
      Result = c(
        d$Result[d$Fish=="Herring"&d$Area=="30-31"],
        d$Result[d$Fish=="Salmon"&d$Area=="30-31"]
      )
    ))
    
    pcddf <- conc_pcddf * coeff
    colnames(pcddf@output)[colnames(pcddf@output)=="Compound"] <- "Exposure_agent"
    
    pcddf <- unkeep(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")
    levels(out$Fish)[levels(out$Fish)=="Herring"] <- "Baltic herring"
    
    out <- combine(out, pcddf)
    if("Time" %in% colnames(out@output)) out@output <- fillna(out@output, "Time")
    result(out) <- pmax(0.01,result(out))
    return(out * info)
  },
  unit = "PCDDF, PCB, TEQ: pg /g; Vitamin D, MeHg: µg /g; DHA, EPA, Omega3: mg /g fresh weight"
)

######## Standard amount but adjusted for country-specific herring size

objects.latest("Op_en7749", code_name = "initiate") # [[Goherr: Fish consumption study]]
## Variables assump, often, much, oftenside, muchside, amountRaw, effinfo, effrecomm, amount

amountOrigFormula <- amount@formula
amount@formula <- function(...) {
  amount <- amountOrigFormula()
  result(amount)[result(amount)>3 & amount$Limit=="Max 3 g/d"] <- 3
  levels(amount$Fish)[levels(amount$Fish)=="Herring"] <- "Baltic herring"
  amount <- info * amount
  amount@marginal[colnames(amount@output)=="Recommendation"]<- FALSE
  return(amount)
}
amount@meta$unit <- "g /d"

## Exposure with background exposure but without mother's exposure to child

expo_dir <- Ovariable(
  "expo_dir",
  dependencies=data.frame(Name=c("amount","conc","expo_bg")),
  formula = function(...) {
    out <- conc[conc$Exposure_agent=="TEQ",] * 0 + 1
    out$Exposure_agent <- "Fish"
    out <- combine(conc, out, name="conc")
    out <- oapply(amount * out, cols="Fish", FUN=sum)
    out <- Ovariable(output = data.frame(
      Exposcen = c("BAU", "No exposure"),
      Result = c(1, 0)
    ), marginal=c(TRUE,FALSE)) * out + expo_bg
    out$Exposure <- as.factor(
      ifelse(
        out$Exposure_agent %in% c("DHA", "MeHg"),
        "To child",
        "To eater"
      )
    )
    return(out)
  },
  unit = "PCDDF, PCB, TEQ: pg /d; Vitamin D, MeHg: µg /d; DHA, EPA, Omega3: mg /d; Fish: g /d"
)

## 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",
  unit = "PCDDF, PCB, TEQ: pg /d; Vitamin D, MeHg: µg /d; DHA, EPA, Omega3: mg /d"
)

# Should the background be specific for gender and country? At the moment it is.
expo_bg <- Ovariable(
  "expo_bg",
  dependencies = data.frame(Name="addexposure","info"),
  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 * info, prevresults = TRUE, sources = TRUE)
    
    return(out)
  },
  unit = "PCDDF, PCB, TEQ: pg /d; Vitamin D, MeHg: µg /d; DHA, EPA, Omega3: mg /d"
)

# Stores non-marginal columns for further use.
info <- Ovariable(
  "info",
  dependencies = data.frame(Name = c("jsp")),
  formula = function(...) {
    out <- jsp
    out$Group <- factor(
      paste(out$Gender, out$Ages),
      levels = c("Female 18-45", "Male 18-45", "Female >45", "Male >45")
    )
    out$Country <- factor(out$Country, ordered=FALSE)
    out <- unique(out@output[c("Iter","Country","Group","Gender","Row")])
    out$Result <- 1
    return(out)
  }
)

exposure <- Ovariable(
  "exposure",
  dependencies = data.frame(
    Name = c(
      "expo_dir", # direct exposure, i.e. the person eats or breaths the exposure agent themself
      "expo_indir", # indirect exposure, i.e. the person (typically fetus or infant) is exposed via someone else (mother)
      "mc2d" # 2D Monte Carlo function
    ),
    Ident = c(
      NA,
      "Op_en7797/expo_indir", # [[Infant's dioxin exposure]] # expo_indir
      "Op_en7805/mc2d") # [[Two-dimensional Monte Carlo]]
  ),
  formula = function(...) {
    out <- combine(expo_dir, expo_indir)
    out <- unkeep(out, "Source.1", sources=TRUE)
    out <- mc2d(out)
    return(out)
  },
  unit = "PCDDF, PCB, TEQ: (To eater: pg /day; to child: pg /g fat); Vitamin D, MeHg: µg /day; DHA, EPA, Omega3: mg /day"
)

# Parameters for 2-dimensional Monte Carlo (not used by default)

mc2dparam<- list(
  N2 = 1000, # Number of iterations in the new Iter
  strength = 50, # Sample size to which the fun is to be applied. Resembles number of observations
  run2d = FALSE, # Should the mc2d function be used or not?
  info = info, # Ovariable that contains additional indices, e.g. newmarginals.
  newmarginals = c("Group","Exposure"), # 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.
)

# Pick only the interesting dose-responses (saves time in calculations)
ERFchoice <- Ovariable(
  "ERFchoice",
  ddata = "Op_en7748", # [[Benefit-risk assessment of Baltic herring and salmon intake]]
  subset = "Exposure-response functions of interest"
)

# 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))
  },
  unit="fraction"
)

###### ERF

objects.latest("Op_en2031", code_name = "ERF2") # [[Exposure-response function]] ERF, subgrouping

ERFOrigFormula <- ERF@formula
ERF@formula <- function(...) {
  out <- ERFOrigFormula()
  out$Exposure <- as.factor(ifelse(
    out$Exposure %in% c(
      "Maternal intake through placenta",
      "Maternal ingested intake",
      "Intake through placenta and mother's milk"
    ),
    "To child",
    "To eater"
  ))
  out$Resp <- out$Resp <- factor(out$Resp, levels=resplev)
  out@marginal[colnames(out@output) %in% c("ER_function","Exposure_unit","Scaling")] <- FALSE
  return(out)
}
ERF@meta$unit <- "several"

dat <- opbase.data("Op_en7748", subset="Model parametres")[-1]

incidence <- Ovariable(
  "incidence",
  data = subgrouping(dat[dat$Type=="incidence" , c("Response","Subgroup","Result")], "Subgroup"),
  unit="# /person-year"
)

#### POPULATION
population <- Ovariable(
  "population", 
  ddata = "Op_en7748",
  subset = "Population",
  unit="#"
)

case_burden <- Ovariable(
  "case_burden",
  ddata = "Op_en7748",
  subset = "DALYs of responses",
  unit="DALY /case"
)

objects.latest("Op_en5917", code_name="InpBoD") # [[Disease risk]] InpBoD

InpBoD <- EvalOutput(InpBoD) # Evaluated because is not a dependency but an Input

BoDattrPer1 <- Ovariable(
  "BoDattrPer1",
  dependencies = data.frame(Name="BoDattr", Ident="Op_en2261/BoDattr2"), # [[Health impact assessment]] BoDattr
  formula =  function(...) {
    out <- info * BoDattr / population * 1000
    return(out)
  },
  unit="mDALY per person/year with continuous exposure"
)

BoDattrPer1 <- EvalOutput(BoDattrPer1, verbose=TRUE)

BoDattrOrigFormula <- BoDattr@formula
BoDattr@formula <- function(...) {
  out <- BoDattrOrigFormula()
  out$Resp <- as.character(out$Resp) # Subtracts mortalities from total mortality.
  tmp <- -1 * out[out$Resp %in% c("Heart (CHD)", "Cancer"),]
  tmp$Resp <- "Mortality"
  out <- combine(out, tmp)
  out$Resp <- factor(out$Resp, levels=resplev)
  out <- info * out # Make Group non-marginal as it is
  return(out)
}
BoDattr@meta$unit="DALY /a per subpopulation"

### Plot and print each ovariable 
ovas <- list( # ovariables
  addexposure, amount, amountRaw, assump, BF, BoD, BoDattr, BW,
  case_burden, conc, conc.param, conc_mehg, conc_pcddf, conc_vit, dose,
  effinfo, effrecomm, ERF, ERF_diox, ERF_env, ERF_mehg, ERF_micr,
  ERF_omega3, ERF_vit, ERFchoice, expo_bg, expo_dir, expo_indir, exposure,
  f_ing, f_mtoc, frexposed, incidence, info, InpBoD, jsp, lengt, much,
  muchside, often, oftenside, PAF, population, RR, size, t0.5, time
)
for(i in 1:length(ovas)) {
  tmp <- ovas[[i]]
  cat("Ovariable", tmp@name, "\n")
  print(tryCatch(oggplot(tmp)+theme_gray(base_size=16), error=function(e1) cat("Cannot plot ",tmp@name)))
  oprint(tryCatch(summary(tmp), error=function(e1) cat("Cannot print summary for ",tmp@name)))
}

################ Insight network

#gr <- scrape(type="assessment")
#objects.latest("Op_en3861", "makeGraph") # [[Insight network]]
#gr <- makeGraph(gr)
#export_graph(gr, "ruori.svg")
#render_graph(gr)

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

### Empty each ovariable for storage and store

oempty(all=TRUE)

objects.store(list=setdiff(ls(), "wiki_username"))
cat("Objects stored without output:", setdiff(ls(), "wiki_username"), ".\n")
</rcode>

==== Other codes ====
{|class="wikitable collapsible collapsed"
!See details
|----
|
===== Select and run =====

Test run 21.7.2019 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=uQxNDeVX2EIBADjQ]

<rcode label="Run model" variables="
name:decsel|description:Choose which decisions to add to model (BAU of each is always considered)|type:checkbox|options:
  '1';Information improvements: BAU (no), yes (size:2x);
  '2';Recommendations about herring: BAU (current), eat more, eat less (size:3x);
  '3';Recommendations about salmon: BAU (current), eat more, eat less (size:3x);
  '4';Select herring size for human consumption: BAU (large only), ban large, new products (both) (size:3x);
  '5';Should background be considered? BAU (yes), no (size:2x);
  '6';Is fish consumption limited to 3 g/day? BAU (no), yes (size:2x);
  '7';Long timeline: BAU (2018), 1990, 2000, 2010 (size:4x);
  '8';Compound mixture: BAU (all), MeHg, Dioxin, Omega3 (size:4x);
  '9';Scenario for length vs TEQ (all sizes in 2009);
  '10';Coherent consumer policy: BAU, increase, reduce (size:3x);
  '11';Scale average amount to Luke data: BAU (size:1x)|default:'5';'10'|
category:Model parameters|
name:N|description:How many iterations to run?|default:10
" graphics=1>
# This is code Op_en7748/ on page [[Benefit-risk assessment of Baltic herring and salmon intake]]
library(OpasnetUtils)

objects.latest("Op_en7748", code_name="model") # All objects of the model

rm(list=ls()[grepl("^Dec", ls())])

BAU <- Ovariable(
  "BAU",
  output = data.frame(
    Info.improvements = "BAU",
    Recomm.herring = "BAU",
    Recomm.salmon = "BAU",
    Limit = "No limit",
    Select.size="BAU",
    Time="2018",
    Background="Yes",
    Mixtures="BAU",
    BAUResult = 1
  ),
  marginal = c(rep(TRUE,8),FALSE)
)

BAUt <- BAU
BAUt$Time <- NULL

# Decisions and marginal collapses
decs <- deci[deci$Num %in% c("0",decsel) , ] # Choose decisions
if("10" %in% decsel) {decs <- decs[!decs$Decision %in% c("Info.improvements","Recomm.herring","Recomm.salmon"),]} # Cons.policy covers these
DecisionTableParser(decs)

cat("Model iterations run:",N,"\n")
cat("Decisions used in model:\n")
oprint(deci)

InpBoD <- EvalOutput(InpBoD)

BoDattr <- EvalOutput(BoDattr)

BS <- 24

oggplot(BoDattr)

fig21 <- ggplot((BAU*BoDattr)@output,
                aes(x=BoDattrResult, colour=Group))+stat_ecdf()+
  theme_grey(base_size=BS)+
  facet_wrap(~ Resp, scales = "free_x")+
  labs(
    title="Fig 21. Burden of disease by response and group",
    x="Disease burden (mDALY/a per a random individual)",
    y="Cumulative probability")

fig22 <- ggplot((oapply(BoDattr, INDEX=c("Resp","Group","Country","Policy"), mean))@output,
                aes(x=Group, weight=BoDattrResult, fill=Resp))+
  geom_bar()+facet_grid(Country ~ Policy)+
  theme_grey(base_size=BS)+
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = -90)
  )+
  labs(
    title="Fig 22. Burden of disease by country, group and policy",
    y = "Disease burden (mDALY/a per person)"
  )#+
#  scale_fill_manual(values=varit12)

fig23 <- ggplot((BAU * info * oapply(BoDattr, cols=c("Resp","Response"), FUN = sum))@output,
                aes(x=Result, colour=Group))+stat_ecdf()+
  theme_grey(base_size=BS)+
  coord_cartesian(xlim=c(-10,10))+facet_wrap(~Country)+
  labs(
    title="Fig 23. Net disease burden by country and population group",
    x="Disease burden (mDALY/a per a random individual)",
    y="Cumulative probability"
  )

</rcode>

===== Exposure-response functions =====

* Model run 12.4.2018 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=1JmEzjmNBvoXEo14] 152354803523
* Model run 13.10.2018 with infertility [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=PGei5FcnuUPwcTcz]

<rcode label="Plot Fig 17" graphics=1>
# This is code Op_en7748/erf on page [[Goherr assessment]]

library(OpasnetUtils)
library(ggplot2)

openv.setN(100)

# Fetch objects to show
objects.latest("Op_en2261",code_name="RR") # [[Health impact assessment]] RR
objects.latest("Op_en2261",code_name="casesabs") # [[Health impact assessment]] casesabs

# Pick only the interesting dose-responses (saves time in calculations)
ERFchoice <- Ovariable(
  "ERFchoice",
  ddata = "Op_en7748", # [[Benefit-risk assessment of Baltic herring and salmon intake]]
  subset = "Exposure-response functions of interest"
)

dose <- Ovariable(
  "dose",
  data=data.frame(
    Exposure_agent=rep(c("DHA","Omega3","MeHg","TEQ","Vitamin D"),each=18),
    Expo=rep(c("DHA","Omega3","MeHg","TEQ","Vit D"),each=18),
    Dose = c(0.03, 0.1, 0.3, 1, 1.1, 1.9, 2, 2.1, 3, 7, 8, 10, 30, 100, 110, 300, 1000, 3000),
    Result = c(0.03, 0.1, 0.3, 1, 1.1, 1.9, 2, 2.1, 3, 7, 8, 10, 30, 100, 110, 300, 1000, 3000)
  )
)

mc2dparam<- list(
  run2d = FALSE # Should the mc2d function be used or not?
)

RR <- EvalOutput(RR)
casesabs <- EvalOutput(casesabs)
erfs <- unkeep(combine(RR,casesabs), sources=TRUE)
erfs$Sex <- NULL
erfs <- oapply(erfs, cols="Iter",FUN=mean)@output
erfs$Group <- paste(erfs$Resp, "by", erfs$Expo)

ggplot(erfs, aes(x=Dose,y=Result,colour=Resp,group=Group,linetype=Expo))+
  geom_line(size=2)+facet_wrap(~Group, scales="free_y")+
  scale_x_log10()+theme_gray(base_size=18)+
  labs(
    title="Fig 17. Exposure-response functions for health impacts",
    y="Relative risk (CHD, stroke) or number of cases",
    x="Dose (Omega3: mg/d, TEQ: pg/kg/d except tooth pg/d, MeHg: µg/kg/d, vitamin D: µg/d)"
  )

objects.store(erfs)
cat("Data.frame erfs stored.\n")
</rcode>

===== Initiate goherr assessment ovariable =====
<span id="goherr"> </span>

<rcode name="goherr2" label="Initiate goherr (for developers only)" embed=1>
# This is code Op_en7748/goherr on page [[Benefit-risk assessment of Baltic herring and salmon intake]]
library(OpasnetUtils)

goherr <- Ovariable(
  "goherr",
  dependencies = data.frame(opbase.data("Op_en7748", subset="Model dependencies"))[-1],
  formula = function(...) {
    return(BoD)
  }
)
colnames(goherr@dependencies)[colnames(goherr@dependencies)=="Result"] <- "Page"

objects.store(goherr)
cat("Assessment ovariable goherr stored.\n")
</rcode>

===== Plot concentrations and survey =====

* Requires codes Op_en7748/bayes and indirectly Op_en7748/preprocess.
* Model run 1.3.2017 [http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=sgkkXkgFEuvaXT9R]

<rcode label="Plot concentrations and answers" graphics=1>
#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
</rcode>

===== 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. The following ovariables are used to produce the data for BBN: bgexposure, amount, expoRaw, conc, lengt, BoD. There are merged and then the following information are selected:
*  Information improvements
*  Food recommendations herring amount
*  Food recommendations salmon amount
*  Food recommendations for herring size
*  Consumer country
*  Consumer gender
*  Consumer age group
*  Size of an indiv. herring in human diet
*  Effect of improved information herring
*  Human consumption of herring
*  Human consumption of salmon
*  Effect of improved information salmon
*  Other intake of Omega3
*  Omega3 intake total
*  Dioxin conc herring
*  Dioxin conc salmon
*  Human intake of dioxin from herring
*  Human intake of dioxin from salmon
*  Dioxin intake total
*  Net burden of disease
*  Omega3 from herring
*  Omega3 from salmon

As a result a data.frame with 22 columns is produced. This is used to produce 13 csv files, one for each BBN node.

<rcode label="Produce BBN interface data (works on own computer only)">
#This is code Op_en7748/ on page [[Benefit-risk assessment of Baltic herring and salmon intake]]

library(OpasnetUtils)
library(reshape2)
objects.latest("Op_en7748", code_name="hia")
# http://en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=9aoTj5dnW9vGL0wX with 5000 iterations

if(TRUE) {
  BAU <- Ovariable(data=data.frame(
    Info.improvements = "BAU",
    Recomm.herring = "BAU",
    Recomm.salmon = "BAU",
    Limit = "No limit",
    Select.size = "BAU",
    Time = 2018,
    Mixtures = "BAU",
    Background = "Yes",
    Result = 1
  ))
} else {
  BAU <- 1
}
# Collect all ovariables together that contain input for BBN.
# Without Iter, it contains 72 unique rows: Exposure_agent 2, Fish 3 (effectively 2),
# Info.improvements 2, Recomm.herring 3, Recomm.salmon 3
d3 <-  BAU + bgexposure[bgexposure$Exposure_agent %in% c("TEQ", "Omega3"),] +
  amount + expoRaw + conc + lengt +
  oapply(BoD, cols=c("Resp","Response"), FUN=sum)

d3 <- d3@output[c(
  "Iter",
  "Country",
  "Gender",
  "Ages",
  "Exposure_agent",
  "Fish",
  "bgexposureResult",
  "amountResult",
  "expoRawResult",
  "concResult",
  "BoDResult",
  "Info.improvements",
  "Recomm.herring",
  "Recomm.salmon",
  "effinfoResult",
  "Select.size",
  "lengtResult"
)
]

d3$Group <- paste(d3$Fish, d3$Exposure_agent, sep="_")
d4 <- reshape(
  d3,
  v.names=c(
    "concResult",
    "amountResult",
    "expoRawResult",
    "bgexposureResult",
    "effinfoResult",
    "lengtResult"
  ),
  idvar=c("Iter","Info.improvements","Recomm.herring","Recomm.salmon","Select.size"),
  timevar="Group",
  drop=c("Exposure_agent","Fish"),
  direction="wide"
)

# Match the columns above to the final columns needed (below)
d5 <- d4[c(
  "Info.improvements", # 1
  "Recomm.herring",
  "Recomm.salmon",
  "Select.size",
  "Country",         # 5
  "Gender",
  "Ages",
  "lengtResult.Baltic herring_TEQ",
  "effinfoResult.Baltic herring_TEQ", # 9a
  "amountResult.Baltic herring_TEQ",   # 10
  "amountResult.Salmon_TEQ",
  "effinfoResult.Salmon_TEQ", # 9b
  "bgexposureResult.Baltic herring_Omega3",
  "Iter",
  "concResult.Baltic herring_TEQ",
  "concResult.Salmon_TEQ",
  "expoRawResult.Baltic herring_TEQ",
  "expoRawResult.Salmon_TEQ",
  "Iter",
  "BoDResult",
  "expoRawResult.Baltic herring_Omega3",
  "expoRawResult.Salmon_Omega3"
)]

colnames(d5) <- c(
  "Information_improvements",
  "Food_recommendations_herring_amount",
  "Food_recommendations_salmon_amount",
  "Food_recommendations_for_herring_size",
  "Consumer_country",                   # 5
  "Consumer_gender",
  "Consumer_age_group",
  "Size_of_an_indiv._herring_in_human_diet",
  "Effect_of_improved_information_herring", # 9a
  "Human_consumption_of_herring",     # 10
  "Human_consumption_of_salmon",
  "Effect_of_improved_information_salmon", # 9b
  "Other_intake_of_Omega3",
  "Omega3_intake_total",  
  "Dioxin_conc_herring",            # 15
  "Dioxin_conc_salmon",
  "Human_intake_of_dioxin_from_herring",
  "Human_intake_of_dioxin_from_salmon",
  "Dioxin_intake_total",
  "Net_burden_of_disease",       # 20
  "Omega3_from_herring",
  "Omega3_from_salmon"
)
d5$Omega3_intake_total <- d5$Omega3_from_herring + d5$Omega3_from_salmon
d5$Dioxin_intake_total <- d5$Human_intake_of_dioxin_from_herring + d5$Human_intake_of_dioxin_from_salmon

write.csv(d5[c(4,5,8)],        "Size_of_an_indiv._herring_in_human_diet.csv", row.names=FALSE)
write.csv(d5[c(1,5,6,7,9)],    "Effect_of_improved_information_herring.csv", row.names=FALSE)
write.csv(d5[c(1,5,6,7,12)],   "Effect_of_improved_information_salmon.csv", row.names=FALSE)
write.csv(d5[c(2,5,6,7,9,10)], "Human_consumption_of_herring.csv", row.names=FALSE)
write.csv(d5[c(3,5,6,7,12,11)],"Human_consumption_of_salmon.csv", row.names=FALSE)
write.csv(d5[c(13)],           "Other_intake_of_Omega3.csv", row.names=FALSE)
write.csv(d5[c(10,11,13,14)],  "Omega3_intake_total.csv", row.names=FALSE)
write.csv(d5[c(5,8,15)],       "Dioxin_conc_in_herring.csv", row.names=FALSE)
write.csv(d5[c(5,16)],         "Dioxin_conc_in_salmon.csv", row.names=FALSE)
write.csv(d5[c(10,15,17)],     "Human_intake_of_dioxin_from_herring.csv", row.names=FALSE)
write.csv(d5[c(11,16,18)],     "Human_intake_of_dioxin_from_salmon.csv", row.names=FALSE)
write.csv(d5[c(17,18,19)],     "Dioxin_intake_total.csv", row.names=FALSE)
write.csv(d5[c(6,7,14,19,20)], "Net_burden_of_disease.csv", row.names=FALSE)
</rcode>
|}

== See also ==

* Sara M. Pires, Géraldine Boué, Alan Boobis, Hanna Eneroth, Jeljer Hoekstra, Jeanne-Marie Membré, Inez Maria Persson, Morten Poulsen, Juliana Ruzante, Jacobvan Klaverena, Sofie T. Thomsen, Maarten J. Nauta. Risk Benefit Assessment of foods: Key findings from an international workshop. Food Research International, Available online 10 September 2018. {{doi|10.1016/j.foodres.2018.09.021}}
* [[Goherr:5.4 Benefit-risk assessment of previous, current and future fish intake]]
* [[Omega-3 content in salmon]]
* [http://www.livsmedelsverket.se/globalassets/rapporter/2013/2013_livsmedelsverket_21_risk_benefit_herring_salmonid_fish_ver2.pdf?id=4013 Risk and Benefit Assessment of Herring and Salmonid Fish from the Baltic Sea Area]: National Food Agency in Sweden performed a risk assessment about the Baltic herring exemption in 2011 and concluded that "In conclusion, a cessation of the exemption from maximum limits would be more beneficial from a public health point-of view than a continued exemption. In the case of no exemption there would be a decreased exposure of the population to dioxins and dl-PCB without any limitation of the intake of beneficial nutrients." The conclusion was based on two criteria: comparison of TWI and intake in women in childbearing age and in children; and a view that large herring with high dioxin levels can be banned and replaced with small herring with low dioxin levels without any other change.<ref>Anders Glynn, Salomon Sand and Wulf Becker. Risk and Benefit Assessment 
of Herring and Salmonid Fish from the Baltic Sea Area. Report 21/2013. Livsmedelsverket, Sweden.[http://www.livsmedelsverket.se/globalassets/rapporter/2013/2013_livsmedelsverket_21_risk_benefit_herring_salmonid_fish_ver2.pdf?id=4013]</ref>
* [http://www.livsmedelsverket.se/globalassets/rapporter/2010/2012_livsmedelsverket_7_market_basket_2010.pdf Swedish Market Basket 2010]: Swedish Food Agency performed a market basket study and concluded that the per capita dioxin+PSB intake was 39 (28-49) pg/d TEQ in 2010. This is about one fourth of the tolerable weekly intake (140 pg/d TEQ in a 70-kg person), and therefore it was concluded that such levels do not markedly increase the risk of health effects<ref>National Food Agency. Market Basket 2010 - chemical analysis, exposure estimation and health-related assessment of nutrients and toxic compounds in Swedish food baskets. Report 7/2012. Livsmedelsverket, Sweden. [http://www.livsmedelsverket.se/globalassets/rapporter/2010/2012_livsmedelsverket_7_market_basket_2010.pdf]</ref>.
* [http://www.livsmedelsverket.se/globalassets/matvanor-halsa-miljo/kostrad-matvanor/matvaneundersokningar/riksmaten_2010_20111.pdf?id=3588 Riksmaten 2010]
* [http://www.food.dtu.dk/-/media/institutter/foedevareinstituttet/publikationer/pub-2015/rapport_danskernes-kostvaner-2011-2013.ashx Danskernes kostvaner 2011-2013] 

* Assmuth Timo and Jalonen Pauliina 2005: Risks and management of dioxin-like compounds in Baltic Sea fish: An integrated assessment. Nordic Council of Ministers, Copenhagen. {{#l:Assmuth Jalonen Dioxin risk assessment 2005.pdf}} [http://www.diva-portal.org/smash/get/diva2:700892/FULLTEXT01.pdf]
* EFSA 2012. Update of the monitoring of levels if dioxins and PCBs in food and feed. EFSA Journal 10(7):2832. doi: 10.2903/j.efsa.2012.2832

== References ==
<references/>

== Keywords ==

Terminology
; Normative scenarios: paths you need to take to reach a defined goal
; Expolorative scenarios: identify key uncertainties and dependencies to describe coherent paths into the future.
; Governance types: How things are managed (e.g. top down command or co-management).
; Management action: Actions to be taken based on decision-maker's decision (i.e. decisions)

== See also ==
*[http://www.livsmedelsverket.se/globalassets/rapporter/2013/2013_livsmedelsverket_21_risk_benefit_herring_salmonid_fish_ver2.pdf?id=4013 Risk and Benefit Assessment of Herring and Salmonid Fish from the Baltic Sea Area]
* [http://www.livsmedelsverket.se/globalassets/rapporter/2010/2012_livsmedelsverket_7_market_basket_2010.pdf Swedish Market Basket 2010]
* [http://www.livsmedelsverket.se/globalassets/matvanor-halsa-miljo/kostrad-matvanor/matvaneundersokningar/riksmaten_2010_20111.pdf?id=3588 Riksmaten 2010]
* [http://www.google.fi/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0CB4QFjAA&url=http%3A%2F%2Fwww.food.dtu.dk%2F~%2Fmedia%2FInstitutter%2FFoedevareinstituttet%2FPublikationer%2FPub-2015%2FRapport%2520-%2520Danskernes%2520Kostvaner%25202011-2013.ashx%3Fla%3Dda&ei=jr9mVa2xJMKgsgH2hoDABw&usg=AFQjCNHQj6ft9PB2fcs-fVz_XKGYKJF-lg&sig2=rZ-d9u4I58xUqBKg3KjuOQ&bvm=bv.93990622,d.bGg Danskernes kostvaner 2011-2013]
* Shared understanding of Goherr policy situation using RDF approach (a detail about Baltic herring) [https://docs.google.com/drawings/d/1rQWHnMBBuYFUEZCyCok_DsM14t_MOZoF2FYuMLer5MI/edit]

{{Goherr}}

== Related files ==