Lähdekoodi sivulle Silakan hyöty-riskiarvio

[[op_en:Benefit-risk assessment of Baltic herring]]
[[heande:Silakan riskiarvio]]
[[Category:Dioxin]]
[[Category:Fish]]
{{tutkimus|moderaattori=Jouni|edistyminen=Täysluonnos|edistymistaulu=Kyllä}}
:''Aiheesta on kirjoitettu myös raportti [[Silakan hyöty-riskiarvio / raportti]] {{#l:Eviran_tutkimuksia_1_2015_silakan_hyöty-riskiarviointi.pdf}} sekä luento [[:Tiedosto:Kala ja terveys.ppt]]. Osa työstä sekä lisäaineistoja löytyy sivulta [[:heande:Silakan riskiarvio]]. Englanniksi löytyy [[:op_en:Benefit-risk assessment of Baltic herring]] ja [[:op_en:File:BRA of Baltic herring.pptx]].

'''Silakan hyöty-riskiarvio''' on tutkimus, jossa selvitetään suomalaisten silakansyöntiä ja silakan sisältämien hyödyllisen ravintoaineiden hyötyjä haitallisten ympäristömyrkkyjen haittoja vastaan. Erityisen kiinnostuksen kohteena ovat lasten riskit. Tutkimus on kolmeosainen ja nyt ollaan menossa vaiheessa 3:
# Aiemman tiedon perusteella tehdään alustava arvio, jossa erityisesti yritetään tunnistaa kriittisiä kohderyhmiä kaulutuskyselyä varten.
# Tehdään kohdennettu kysely silakan syöntimäärien selvittämiseksi riskiryhmissä.
# Tehdään päivitetty riskiarvio aiempien vaikutustietojen ja altistustietojen sekä kyselystä saatujen uusien syöntitietojen perusteella.

==Kysymys==

# Millainen on suomalaisten, erityisesti riskiryhmien, altistuminen silakalle ja sen hyödyllisille ja haitallisille aineille?
# Kuinka suuret ovat silakan syönnin hyödyt ja haitat väestölle kokonaisuutena, eri väestöryhmille ja erityisesti eri riskiryhmille?

==Vastaus==

* Tuloksista on kirjoitettu raportti: [[Silakan hyöty-riskiarvio / raportti]]
* [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=OG0cNA6g7H9qmJp2 Malli jossa silakkaodotusarvot korjattu, MeHg lisätty ja Iter = 5000] 24.1.2015. Pysyvä ajo. [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=eElVa3RtPhXVEuZ2 Tuloksia]
* [http://fi.opasnet.org/fi_wiki/index.php?title=Special:RTools&id=ZqXmrRlobrfJhzny Sellaisia tuloksia 26.1.2015], joiden laskemiseen on tarvittu ei-julkista dataa.

=== Näytä tuloksia ===

* [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=PZQkBIvf8KJOoRok Malliajo 15.10.2015] jossa oliot pysyvästi tallennettu. [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=ByNEe1UbulZvkYWC Tuloksia 19.10.2015, mukana muutama uusi englanninkielinen kuva]
* [[:op_en:File:BRA of Baltic herring.pptx]] esitys englanniksi aiheesta

<rcode label="Näytä tuloksia" graphics=1 embed=0 variables="name:server|type:hidden|default:TRUE">
library(OpasnetUtils)
library(ggplot2)
library(reshape2)
library(plyr)

################## Suomentaja

suomenna <- function(ova) {
	d <- ova@output
	if("Trait" %in% colnames(d)) {
		levels(d$Trait)[levels(d$Trait) == "Cancer"] <- "Syöpä"
		levels(d$Trait)[levels(d$Trait) == "CHD2"] <- "Sydänkuolema"
		levels(d$Trait)[levels(d$Trait) == "Stroke"] <- "Aivohalvaus"
		levels(d$Trait)[levels(d$Trait) == "Child's IQ"] <- "Lapsen ÄO"
		levels(d$Trait)[levels(d$Trait) == "Tooth defect"] <- "Hammasvaurio"
		levels(d$Trait)[levels(d$Trait) == "Dental defect"] <- "Hammasvaurio"
		levels(d$Trait)[levels(d$Trait) == "Vitamin D recommendation"] <- "D-vitamiinin saantisuositus"
		levels(d$Trait)[levels(d$Trait) == "Dioxin recommendation"] <- "Dioksiinin saantisuositus"
	}
	if("Exposure_agent" %in% colnames(ova@output)) {
		levels(d$Exposure_agent)[levels(d$Exposure_agent) == "Vitamin_D"] <- "D-vitamiini"
		levels(d$Exposure_agent)[levels(d$Exposure_agent) == "EPA"] <- "EPA"
		levels(d$Exposure_agent)[levels(d$Exposure_agent) == "DHA"] <- "DHA"
		levels(d$Exposure_agent)[levels(d$Exposure_agent) == "Omega3"] <- "Omega3"
		levels(d$Exposure_agent)[levels(d$Exposure_agent) == "PCDDF"] <- "Dioksiini"
		levels(d$Exposure_agent)[levels(d$Exposure_agent) == "PCB"] <- "PCB"
		levels(d$Exposure_agent)[levels(d$Exposure_agent) == "TEQ"] <- "TEQ"
		levels(d$Exposure_agent)[levels(d$Exposure_agent) == "MeHg"] <- "Metyylielohopea"
	}
	if("Hedelm" %in% colnames(d)) {
		levels(d$Hedelm)[levels(d$Hedelm) == FALSE] <- "Ei"
		levels(d$Hedelm)[levels(d$Hedelm) == TRUE] <- "Kyllä"
	}
	
	colnames(d)[colnames(d) == "Age"] <- "Ikäryhmä"
	colnames(d)[colnames(d) == "Exposure_agent"] <- "Altiste"
	colnames(d)[colnames(d) == "Trait"] <- "Vaste"
	colnames(d)[colnames(d) == "IkÃ¤"] <- "Ikä"
	colnames(d)[colnames(d) == "SilakkamÃ¤Ã¤rÃ¤"] <- "Silakkamäärä"
	
	return(d)
}

##################### VOI-ANALYYSI

#objects.latest("Op_en2480", code_name = "Initiate") # Load VOI functions

VOI <- function(ova, dec) { # Use this simpler total VOI function instead
	ncpi <- oapply(oapply(ova, cols = dec, FUN = "min"), cols = "Iter", FUN = "mean")
	ncuu <- oapply(oapply(ova, cols = "Iter", FUN = "mean"), cols = dec, FUN = "min")
	out <- ncpi - ncuu
	return(out)
}


########### HAE RISKITULOKSET JA ALUSTA

#!!+++++++++++++++++++++++++++++++++++++++++++++++++
objects.latest("Op_en7392", code_name = "translate") # translate
objects.latest("Op_fi3831", code_name = "compute") # Haetaan valmiiksi laskettu indivrisk, pop ja exposure. 
objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]] oggplot

#objects.get("fLcPm8KyLhUuZC0c") # Haetaan tilapäinen versio, jossa on yli 3 g/d saannit typistetty.

korjaa <- function(o) {
	d <- o@output
	colnames(d)[colnames(d) == "IkÃ¤"] <- "Ikä"
	colnames(d)[colnames(d) == "SilakkamÃ¤Ã¤rÃ¤"] <- "Silakkamäärä"
	o@output <- d
	return(o)
}
amount <- korjaa(amount)
#daly <- korjaa(daly)
exposure <- korjaa(exposure)
indivrisk <- korjaa(indivrisk)

#ii+++++++++++++++++++++++++++++++++++++++++++++++++

if(exists("server")) BS <- 24 else BS <- 24

cat("Tässä ajossa on", length(unique(indivrisk@output$Iter)), "iteraatiota.\n")

################ VASTAAJAT 

silakka$Sukupuoli <- factor(silakka$Sukupuoli, levels = c("Nainen", "Mies"))

ggplot(subset(silakka, Ikäryhmä == "Aikuinen"), aes(x = Ikä, weight = 1, fill = Sukupuoli)) + geom_bar(position = "dodge") + 
	theme_gray(base_size = BS) + labs(# font etc. http://sape.inf.usi.ch/quick-reference/ggplot2/themes
		title = "Vastaajat iän ja sukupuolen mukaan", 
		x = "Ikä (a)", 
		y = "Lukumäärä"
)

ggplot(subset(silakka, Ikäryhmä == "Lapsi"), aes(x = Ikä, weight = 1, fill = "Red")) +
	geom_bar(position = "dodge") + theme_gray(base_size = BS) +
	labs(
		title = "Lapset iän mukaan", 
		x = "Ikä (a)", 
		y = "Lukumäärä"
	)

######################### Pitoisuudet

ggplot(suomenna(concentration), aes(x = Result, colour = Altiste))+geom_density(size = 1)+scale_x_log10() +
     theme_gray(base_size = BS) + labs(title = "Eri aineiden pitoisuudet silakassa", x = "Pitoisuus (g /g)", y = "Todennäköisyystiheys")
 
ggplot(suomenna(concentration), aes(y = Result, x = Altiste))+geom_boxplot()+scale_y_log10() +
     theme_gray(base_size = BS) + labs(title = "Eri aineiden pitoisuudet silakassa", y = "Pitoisuus (g /g)")

########################## SILAKANSYÖNTI

amount$Ikä <- as.numeric(as.character(amount$Ikä))
amount$Ikäluokka <- cut(amount$Ikä, c(0, 20, 65, 80), levels = c("0-20", "21-65", "66-80"), include.lowest = TRUE)
amount$Sukupuoli <- factor(amount$Sukupuoli, levels = c("Nainen", "Mies"))
rannikko <- c(
	"Uusimaa",
	"Pohjanmaa", 
	"Kymenlaakso",
	"Etelä-Pohjanmaa",
	"Satakunta",
	"Keski-Pohjanmaa",
	"Pohjois-Pohjanmaa", 
	"Varsinais-Suomi"
)

amount$Rannikko <- ifelse(amount$Maakunta %in% rannikko, TRUE, FALSE)

cat("Silakansyönti ikäryhmittäin ja sukupuolittain.\n")

amosummary <- summary(amount, marginals = c("Sukupuoli", "Ikäluokka", "Rajoita"), fun = c("mean", "sd", "Q0.025", "Q0.975", "length"))

oprint(amosummary)

ggplot(subset(amosummary, Rajoita == "Ei"), aes(x = Ikäluokka, weight = mean, fill = Sukupuoli))  +
geom_bar(position = "dodge") + theme_gray(base_size = BS) + 
     labs(
         title = "Silakansyönti ikäryhmittäin", 
         y = "Syönti g/d",
         x = "Ikäryhmä"
     )

amo <- suomenna(amount)

ggplot(subset(amo, Rajoita == "Ei"), aes(x = Result, colour = Ikäluokka)) + stat_ecdf(size = 1) + 
  facet_wrap( ~ Sukupuoli) + scale_x_log10() + theme_gray(base_size = BS) +
  labs(
	title = "Silakansyönti sukupuolen ja iän mukaan", 
	colour = "Ikäryhmä", 
	x = "Silakansyönti (g/d)", 
	y = "Kumulatiivinen todennäköisyysjakauma"
  )

amen <- translate(amount[amount$Rajoita == "Ei" , colnames(amount@output) != "Ikä"])
ggplot(amen@output, aes(x = Result, colour = Age_group)) + stat_ecdf(size = 2) + 
  facet_wrap( ~ Sex) + scale_x_log10() + theme_gray(base_size = BS) +
  labs(
	title = "Consumption of Baltic herring by sex and age", 
	colour = "Age group", 
	x = "Consumption (g/d)", 
	y = "Cumulative probability distribution"
  )

ggplot(amen@output, aes(x = Result, group = Maakunta, colour = Rannikko)) + stat_ecdf(size = 1) + 
coord_cartesian(ylim = c(0.45, 1.02), xlim = c(0.03, 100)) +
 scale_x_log10() + theme_gray(base_size = BS) +
  labs(
	title = "Consumption of Baltic herring by district", 
	colour = "Coastal area?", 
	x = "Consumption (g/d)", 
	y = "Cumulative probability distribution"
  )

ggplot(amo, aes(x = cut(Ikä, c(0, 20, 65, 80), include.lowest = TRUE), y = Result, fill = Sukupuoli)) + geom_boxplot() + 
  facet_wrap( ~ Sukupuoli, scale = "free", ncol = 2) + theme_gray(base_size = BS) +
  labs(
	title = "Silakansyönti sukupuolen ja iän mukaan", 
	colour = "Ikäryhmä", 
	y = "Silakansyönti (g/d)", 
	x = "Ikä"
  )

ggplot(subset(amo, as.integer(as.character(Iter)) <= 1000 & Rajoita == "Ei"), aes(x = Ikä, y = Result, colour = Sukupuoli)) + geom_smooth(size = 1.5) + 
	theme_gray(base_size = BS) + labs(title = "Silakansyönti iän mukaan", y = "Silakansyönti (g/d)")
	
ggplot(subset(amo, as.integer(as.character(Iter)) <= 1000 & Rajoita == "Ei"), aes(x = Maakunta, y = Result)) + geom_boxplot() + 
	theme_grey(base_size = BS) + # font etc. http://sape.inf.usi.ch/quick-reference/ggplot2/themes
	labs( # label names
		title 	= "Silakan kokonaissyönti maakunnittain"
	) +
	theme(axis.text.x = element_text(angle = 90, hjust = 1)) +# Turn text and adjust to right 
	coord_cartesian(ylim = c(0, 15)) 

ggplot(subset(amo, as.integer(as.character(Iter)) <= 1000 & Rajoita == "Ei"), aes(x = Maakunta, y = Result)) + geom_boxplot() + 
	theme_grey(base_size = BS) + # font etc. http://sape.inf.usi.ch/quick-reference/ggplot2/themes
	labs( # label names
		title 	= "Silakan kokonaissyönti maakunnittain"
	) +
	theme(axis.text.x = element_text(angle = 90, hjust = 1)) # Ei katkaista y-akselia

# Testataan todennäköisillä determinanteilla: sukupuoli, ikä, maakunta, paino.

cat("Regressioanalyysi silakansyönnin selittäjistä aikuisilla.\n")

fit <- lm(Silakkamäärä ~ Ikä + Sukupuoli + Paino + Rannikko, data = subset(silakka, Ikäryhmä == "Aikuinen"))

oprint(summary(fit))

fit <- lm(Silakkamäärä ~ Ikä + Sukupuoli + Paino + Maakunta, data = subset(silakka, Ikäryhmä == "Aikuinen"))

oprint(summary(fit))

fit <- lm(Silakkamäärä ~ Ikä + Sukupuoli + Paino + Rannikko, data = subset(silakka, Ikäryhmä == "Lapsi"))

cat("Regressioanalyysi silakansyönnin selittäjistä lapsilla.\n")

oprint(summary(fit))

#### Testataan suurkulutusta selittäviä tekijöitä

silakka$Suurkuluttaja <- silakka$Silakkamäärä >= 10

fit <- glm(Suurkuluttaja ~ Ikä + Sukupuoli + Paino + Rannikko, family = binomial, data = silakka)
oprint(summary(fit))

######################### ALTISTUKSET

expo <- suomenna(exposure)
expo$Sukupuoli <- factor(expo$Sukupuoli, levels = c("Nainen", "Mies"))

ggplot(subset(expo, Rajoita == "Ei" & as.integer(as.character(Iter)) <= 1000), aes(x = Ikä, y = Result, colour = Sukupuoli)) + geom_point() +
  facet_wrap(~ Altiste, scales = "free_y", ncol = 2) + theme_gray(base_size = BS) +
  labs(title = "Yksilötason altistuminen päivässä")

ggplot(subset(expo, Tausta == "Ei"), aes(x = Result, colour = Altiste)) + geom_density(size = 1) + scale_x_log10() + facet_grid(Rajoita ~ .) +
  labs(
    title = "Altistuminen silakassa oleville aineille", 
	y = "Todennäköisyystiheys",
    x = "Altistuminen, g /d") +
  theme_gray(base_size = BS)

temp <- translate(exposure@output[exposure$Tausta == "Ei" & exposure$Exposure_agent %in% c("TEQ", "Vitamin_D", "Omega3", "MeHg"), ])
ggplot(temp, 
aes(x = Result, colour = Exposure_agent)) + geom_density(size = 2) + scale_x_log10() + facet_grid(Rajoita ~ .) +
  labs(
    title = "Exposure to compounds in Baltic herring. Limit to 3 g/d?", 
	y = "Probability distribution",
    x = "Exposure, g /d",
colour = "Compound") +
  theme_gray(base_size = BS)


ggplot(subset(expo, Tausta == "Ei" & Rajoita == "Ei"), aes(x = Result, colour = Sukupuoli)) + geom_density(size = 1) + 
  facet_wrap( ~ Altiste, scale = "free", ncol = 2) + scale_x_log10() + theme_gray(base_size = BS) +
  labs(
	title = "Altistumisjakauma altisteittain ja sukupuolittain",
	y = "Todennäköisyystiheys",
	x = "Altistuminen (g/d)")

ggplot(subset(expo, Tausta == "Ei"), aes(x = Result, colour = Sukupuoli)) + geom_density(size = 1) + 
  facet_grid(Altiste ~ Rajoita , scale = "free") + scale_x_log10() + theme_gray(base_size = BS) +
  labs(title = "Altistumisjakauma sukupuolittain\nSuositellaanko rajoitusta?",
	x = "Altistuminen (g/d)",
	y = "Todennäköisyystiheys")

ggplot(subset(expo, Rajoita == "Ei" & Tausta == "Ei"), aes(x = Result, colour = Ikäryhmä)) + geom_density() + 
  facet_wrap( ~ Altiste, scale = "free", ncol = 2) + scale_x_log10() + theme_gray(base_size = BS) + 
  labs(title = "Altistumisen todennäköisyysjakauma ikäryhmittäin")

ggplot(subset(expo, Rajoita == "Ei" & Tausta == "Kyllä"), aes(x = Result, colour = Ikäryhmä)) + geom_density() + 
  facet_wrap( ~ Altiste, scale = "free", ncol = 2) + scale_x_log10() + theme_gray(base_size = BS) + 
  labs(title = "Altistumisen todennäköisyysjakauma ikäryhmittäin. Tausta mukana")

################################ VÄESTÖTIETO JA PAINOKERTOIMET

#!!++++++++++++++++++++++++++++++++++++++++++++++++
dweight <- Ovariable("dweight", ddata = "Op_fi3831.laatupainotetut_elinvuodet_vaikutuksille")
#ii++++++++++++++++++++++++++++++++++++++++++++++++

dweight@data$Response_metric <- NULL

daly <- unkeep(indivrisk * dweight, sources = TRUE, prevresults = TRUE)
daly@output <- subset(daly@output, Trait != "Dioxin recommendation")

# pophed is the population size for each Hedelm * Age subgroup. This can be summed over either Hedelm or Age.
# pophed is for a single sex, so it should be multiplied with an ovariable with index Sex to get values for the whole population.
pop <- EvalOutput(pop)
pophed <- pop / 2 * Ovariable(
	output = data.frame(unique(daly@output[c("Hedelm", "Age")]), Result = 1),
	marginal = c(TRUE, TRUE, FALSE)
)

############################# RYHMÄKOHTAISET DALYT

# daly: all DALYs (contains exposure and non-exposure scenarios)
## dalyp: DALYs on population level, i.e. daly multiplied by subgroup sizes.
## dalyi: individual DALYs summed over Traits

#Oikeastaan tässä pisäisi laskea joka ryhmälle indivriskistä bootsträppäyksellä keskiarvoja, joista sitten muodostuu
# populaatiolle keskiarvovaikutuksen jakauma. Tämä sitten kerrottaisiin dweightillä ja pophedillä.

dalyp <- daly * pophed 

cat("Silakansyönnin DALYt vasteen, sukupuolen ja politiikkojen mukaan.\n")

# Jokin vika eli ei käytetä J-unitseja
#library(scales)
#junit_trans = function() trans_new("junit", function(x) x^(1/3), function(x) x^3)
#ggplot(data.frame(Result = -50:50), aes(x = Result)) + geom_density() + coord_trans(x = "junit") + labs(title = "Kuva 7. Todennäköisyyden tiheysjakauma ikäryhmittäin")
#Error in if (zero_range(range)) { : missing value where TRUE/FALSE needed

############################# DALYT IKÄRYHMITTÄIN

d1 <- subset(suomenna(dalyp), Rajoita == "Ei")
Nd <- as.data.frame(as.table(tapply(d1$Iter, d1[c("Rajoita", "Vaste", "Sukupuoli", "Hedelm", "Tausta", "Ikäryhmä")], FUN = length)))
Nd <- Nd[!is.na(Nd$Freq) , ]
d1 <- merge(d1, Nd)

test1 <- !d1$Ikäryhmä %in% c("50-54", "55-59", "60-64", "65-69", "70-74", "75-79")
# Tämä varmistaa, että pylväät eivät pilkkoudu kirjaviksi
d1 <- d1[order(d1$Vaste) , ]

ggplot() +
  geom_bar(data = subset(d1, Result > 0), aes(x = Ikäryhmä, y = Result/Freq, fill = Vaste), stat = "identity") +
  geom_bar(data = subset(d1, Result <= 0), aes(x = Ikäryhmä, y = Result/Freq, fill = Vaste), stat = "identity") +
  facet_grid(Tausta ~ Sukupuoli) + theme_gray(base_size = BS) +
  labs(title = "Silakansyönnin tautitaakka Suomessa",
  y = "Menetettyjä elinvuosia (DALY/vuosi)")  + theme(axis.text.x = element_text(angle = 90))

temp <- dalyp[dalyp$Trait != "Vitamin D recommendation" , ]
levels(temp$Trait)[levels(temp$Trait) == "CHD2"] <- "Cardiac death"
temp@marginal[colnames(temp@output) == "Age"] <- TRUE
oggplot(temp[temp$Tausta == "Ei" , ], x = "Age", fill = "Trait", turnx = TRUE) + facet_grid(Rajoita ~ Sex)+
labs(y = "Life years lost (DALY/a",
title = "Health impacts of Baltic herring in Finland. Restrict to 3 g/d?",
colour = "Response"
)

d2 <- translate(d1)
#oggplot(temp, x = "Age", fill = "Response") +
#  facet_grid(Background ~ Sex) +
#  labs(title = "Burden of disease due to Baltic herring in Finland",
#  y = "Life years  lost (DALY/a)")

ggplot() +
  geom_bar(data = subset(d1, test1 & Result > 0), aes(x=Ikäryhmä, y = Result/Freq, fill = Vaste), stat = "identity") +
  geom_bar(data = subset(d1, test1 & Result <= 0), aes(x=Ikäryhmä, y = Result/Freq, fill = Vaste), stat = "identity") +
  facet_grid(Tausta ~ Sukupuoli) + theme_gray(base_size = BS) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  labs(title = "Silakansyönnin tautitaakka Suomessa alle 50-vuotiailla",
  y = "Menetettyjä elinvuosia (DALY/vuosi)")  + theme(axis.text.x = element_text(angle = 90))


############################## VÄESTÖTASON MUUTOKSET RAJOITUSPOLITIIKAN TAKIA

d2 <- dalyp * Ovariable(output = data.frame(Rajoita = c("Kyllä", "Ei"), Result = c(1, -1)), marginal = c(TRUE, FALSE))
d2 <- suomenna(oapply(d2, INDEX = c("Trait", "Sukupuoli", "Iter", "Tausta", "Age"), FUN = sum))
Nd <- as.data.frame(as.table(tapply(d2$Iter, d2[c("Vaste", "Sukupuoli", "Tausta", "Ikäryhmä")], FUN = length)))
Nd <- Nd[!is.na(Nd$Freq) , ]
d2 <- merge(d2, Nd)
d2 <- d2[order(d2$Vaste, d2$Ikäryhmä, d2$Sukupuoli) , ]

test2 <- !d2$Ikäryhmä %in% c("50-54", "55-59", "60-64", "65-69", "70-74", "75-79", "80-84", "85-89", "90-94", "95-")

ggplot() +
  geom_bar(data = subset(d2, Result > 0), aes(x = Ikäryhmä, y = Result/Freq, fill = Vaste), stat = "identity") +
  geom_bar(data = subset(d2, Result <= 0), aes(x = Ikäryhmä, y = Result/Freq, fill = Vaste), stat = "identity") +
  facet_grid(. ~ Sukupuoli) + theme_gray(base_size = BS) +
  labs(title = "Silakkarajoituksen 3 g/d vaikutus tautitaakkaan",
  y = "Menetettyjä elinvuosia (DALY/vuosi)")  + theme(axis.text.x = element_text(angle = 90))

ggplot() +
  geom_bar(data = subset(d2, test2 & Result > 0), aes(x=Ikäryhmä, y = Result/Freq, fill = Vaste), stat = "identity") +
  geom_bar(data = subset(d2, test2 & Result <= 0), aes(x=Ikäryhmä, y = Result/Freq, fill = Vaste), stat = "identity") +
  facet_grid(Tausta ~ Sukupuoli) + theme_gray(base_size = BS) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
  labs(title = "Silakkarajoituksen 3 g/d vaikutus alle 50-vuotiailla",
  y = "Menetettyjä elinvuosia (DALY/vuosi)")  + theme(axis.text.x = element_text(angle = 90))

######################### YKSILÖRISKIT ERI VASTEILLE

riski <- subset(suomenna(indivrisk), Tausta == "Ei" & Rajoita == "Ei" & Iter <= 1000)

ggplot(riski, aes(x = Ikä, y = Result)) + 
	geom_point() + geom_smooth() +
  facet_wrap(~ Vaste, scale = "free_y", ncol = 2)+  theme_gray(base_size = BS) +
  labs(title = "Yksilön silakkariski vasteittain iän mukaan",
	y = "Vaste (vasteilla eri yksiköt)")

ggplot(riski, aes(x = Silakkamäärä, y = Result)) + geom_point() + 
	geom_smooth() + theme_gray(base_size = BS) + scale_x_log10() +
	facet_wrap(~ Vaste, ncol = 2, scale = "free_y") +
	labs(title = "Yksilöriski vasteittain silakansyönnin mukaan",
		x = "Silakan syönti (g/d)",
		y = "Vaste (vasteilla eri yksiköt)")

ggplot(riski, aes(x = Ikä, y = Paino, colour = Sukupuoli)) + 
	geom_point() + theme_gray(base_size = BS) +
	labs(title = "Ruumiinpaino iän ja sukupuolen mukaan",
		x = "Ikä (a)",
		y = "Ruumiinpaino (kg)"
	)

######################## YKSILÖLLISET DALYT

dalyi <- oapply(daly, INDEX = c("Sukupuoli", "Iter", "Rajoita", "Tausta", "Age"), FUN = sum, use_plyr = TRUE)

ggplot(dalyi@output, aes(x = Age, y = Result, colour = Rajoita)) + geom_boxplot() + theme_gray(base_size = BS) +
	facet_grid(Tausta ~ Sukupuoli) +
  labs(title = "Yksilötason DALYt sukupuolittain ja iän mukaan")

ggplot(dalyi@output, aes(x = Age, y = Result, colour = Rajoita)) + geom_boxplot() + theme_gray(base_size = BS) +
	facet_grid(Tausta ~ Sukupuoli) + coord_cartesian(ylim = c(-0.005, 0.001)) +
  labs(title = "Yksilötason DALYt sukupuolittain ja iän mukaan", y = "DALY/a (Rajattu alue)")

# Look at individual risks and probability that the herring eating recommendation is correct.

dalyir <- oapply(daly, cols = c("Trait"), FUN = "sum")

dalyir3 <- dalyir * Ovariable(output = data.frame(Rajoita = c("Kyllä", "Ei"), Result = c(1, -1)), marginal = c(TRUE, FALSE))
dalyir3 <- oapply(dalyir3, INDEX = c("Sukupuoli", "Iter", "Hedelm"), FUN = sum)
dalyir3 <- oapply(dalyir3, INDEX = c("Sukupuoli", "Hedelm"), FUN = function(x) sum(x <= -0.00003) / length(x)) 
# Terveysvaikutuksen rajana 1 terve päivä elinaikana (1/365/80) enemmän.

ggplot(suomenna(dalyir3), aes(x = Hedelm, weight = Result, fill = Sukupuoli)) + 
	geom_bar(position = "dodge") + theme_gray(base_size = BS) + 
	coord_cartesian(ylim = c(0,1)) +
	labs(
		title = "Toisiko silakkarajoitus 3 g/d terveyshyötyä?", 
		y = "Todennäköisyys että toisi",
		x = "Kuuluuko hedelmälliseen ikäryhmään?"
	)

## Look at recommendation of total herring ban.

dalyir0 <- dalyir * Ovariable(output = data.frame(Rajoita = c("Kyllä", "Ei"), Result = c(0, -1)), marginal = c(TRUE, FALSE))
dalyir0 <- oapply(dalyir0, INDEX = c("Sukupuoli", "Iter", "Hedelm"), FUN = sum)
dalyir0 <- oapply(dalyir0, INDEX = c("Sukupuoli", "Hedelm"), FUN = function(x) sum(x <= -0.00003) / length(x)) 

ggplot(suomenna(dalyir0), aes(x = Hedelm, weight = Result, fill = Sukupuoli)) + 
	geom_bar(position = "dodge") + theme_gray(base_size = BS) + 
	coord_cartesian(ylim = c(0,1)) +
	labs(
		title = "Toisiko silakan täyskielto terveyshyötyä?", 
		y = "Todennäköisyys että toisi",
		x = "Kuuluuko hedelmälliseen ikäryhmään?"
	)

# Look at individual risks and probability that the herring eating recommendation is correct.
# WITHIN THE POPULATION THAT ATE >= 3 g/d

daly3n <- amount
daly3n@output <- subset(daly3n@output, Rajoita == "Ei")
daly3n <- dalyir3 / oapply(daly3n, INDEX = c("Sukupuoli", "Hedelm"), FUN = function(x) sum(x > 3) / length(x))

ggplot(suomenna(daly3n), aes(x = Hedelm, weight = Result, fill = Sukupuoli)) + 
	geom_bar(position = "dodge") + theme_gray(base_size = BS) + 
	coord_cartesian(ylim = c(0,1)) +
	labs(
		title = "Toisiko silakkarajoitus 3 g/d terveyshyötyä ylittäjille?", 
		y = "Todennäköisyys että toisi",
		x = "Kuuluuko hedelmälliseen ikäryhmään?"
	)

## Look at recommendation of total herring ban.
# WITHIN THE POPULATION THAT ATE > 0.01 g/d
## TÄTÄ VASTAAVAA KUVAA EI PIIRRETÄ, KOSKA JO NOLLA-ALTISTUKSELLA == 0.01 g/d TULEE TERVEYSVAIKUTUSTA 
## JA TÄMÄN SEURAUKSENA RAJOITUKSESTA HYÖTYVIÄ NÄYTTÄÄ OLEVAN ENEMMÄN KUIN KALAA SYÖVIÄ.

########################## DALYt suhteessa silakkamäärään

sil <- unique(suomenna(daly)[c("Iter", "Ikäryhmä", "Sukupuoli", "Rajoita", "Silakkamäärä")])
sil <- merge(sil, suomenna(dalyi))

ggplot(subset(sil, as.integer(as.character(Iter)) <= 1000) , aes(x = Silakkamäärä, y = Result, colour = Sukupuoli)) + geom_point() + geom_smooth() + 
	scale_x_log10() + facet_grid(. ~ Rajoita) + theme_gray(base_size = BS) + 
	labs(
		title = "Yksilökohtaiset DALYt silakansyönnin mukaan. Sarake: Rajoita", 
		x = "Silakansyönti (g/d)", 
		y = "DALY/hlö/a"
	)

############################# YKSILÖTASON VOI

dalyvoi <- oapply(daly, cols = "Trait", FUN = "sum")
dalyvoi <- VOI(dalyvoi, dec = "Rajoita") 

dalyvoi <- oapply(dalyvoi * pophed, cols = "Age", FUN = "sum") # Population per sex


ggplot(suomenna(dalyvoi), aes(x = Hedelm, weight = Result, fill = Sukupuoli))  +
	geom_bar(position = "dodge") + theme_gray(base_size = BS) + #theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
	labs(
		title = "Lisätiedon arvo silakkasuosituspäätökselle",
		y = "Lisätiedon arvo (DALY)",
		x = "Onko hedelmällisessä iässä?"
	)

######### Plot disability weights

dweight <- EvalOutput(dweight)
ggplot(subset(dweight@output, Trait %in% c(
	"Cancer",
	"CHD2",
	"Stroke",
	"Child's IQ",
	"Dental defect",
	"Vitamin D recommendation"
)), aes(x = dweightResult))+geom_density()+facet_wrap(~ Trait,  scale = "free")

######## Fill in the missing Age * Iter combinations by inputing
# This converts Age into marginal.
# This must be done, otherwise aggragation by oapply results in too small values.

dalypinp <- daly * pophed 

N <- length(unique(dalypinp@output$Iter))
rows <- unlist(tapply(1:nrow(dalypinp@output), 
	dalypinp@output[c("Trait", "Sukupuoli", "Hedelm", "Rajoita", "Age")],
	FUN = function(x) sample(x, N, replace = TRUE)
))
dalypinp@output <- dalypinp@output[rows , ]
dalypinp@output$Iter <- 1:N # From now on the Iter does not match that of any other variable
dalypinp@marginal[colnames(dalypinp@output) == "Age"] <- TRUE

oprint(summary(oapply(dalypinp, cols = c("Age", "Hedelm", "Trait"), FUN = sum), fun = c("mean", "sd", "Q0.025", "Q0.975")), digits = 1)

oprint(summary(oapply(dalypinp, cols = c("Age", "Hedelm", "Sukupuoli"), FUN = sum), fun = c("mean", "sd", "Q0.025", "Q0.975")), digits = 1)

oprint(summary(oapply(dalypinp, cols = c("Age", "Trait"), FUN = sum), fun = c("mean", "sd", "Q0.025", "Q0.975")), digits = 1)

oprint(summary(oapply(dalypinp, cols = c("Hedelm", "Trait"), FUN = sum), fun = c("mean", "sd", "Q0.025", "Q0.975")), digits = 1)

######## DALY distributions by trait

ggplot(subset(dalypinp@output, !Trait %in% c("Child's IQ", "Dental defect") | Sukupuoli == "Nainen") ,
	aes(x = Result, colour = Sukupuoli)) + geom_density() + facet_wrap(~ Trait , scale = "free") +
	labs(title = "Dalyjen jakaumat vasteittain aritmeettisella asteikolla")

ggplot(subset(dalypinp@output, !Trait %in% c("Child's IQ", "Dental defect") | Sukupuoli == "Nainen") ,
	aes(x = abs(Result), colour = Rajoita)) + geom_density() + facet_wrap(~ Trait , scale = "free") +
	scale_x_log10() + labs(title = "Absoluuttidalyt log-asteikolla")

cat("Yli 0.7 pg/kg/d dioksiineille altistuvat:\n")

nrow(subset(exposure@output, Sukupuoli == "Mies" & Exposure_agent == "TEQ" & Rajoita == "Ei" & Hedelm == FALSE & Age %in% c("0", "1-4", "5-9") & Result/Paino > 7E-13)) /
nrow(subset(exposure@output, Sukupuoli == "Mies" & Exposure_agent == "TEQ" & Rajoita == "Ei" & Hedelm == FALSE & Age %in% c("0", "1-4", "5-9")))

odag()

</rcode>

{{attack|1=# |2=Kun vastauksia tulkitaan, tulee huomata, että seuraavia tuoreita tietoja EI OLE KÄSITELTY.
* [http://www.urmc.rochester.edu/news/story/index.cfm?id=4238 Seychell-tutkimus: omega-3 estää elohopeavaikutusta]
* [http://www.efsa.europa.eu/en/press/news/150122.htm EFSAn hyöty-riskiarvio omega-3:sta ja elohopeasta, lehdistötiedote]
* [http://www.efsa.europa.eu/en/efsajournal/pub/3982.htm EFSAn hyöty-riskiarvio omega-3:sta ja elohopeasta]
* [http://link.springer.com/article/10.1007%2Fs10654-014-9966-4 Engeset et al. 2014]: Fish consumption and mortality in the European Prospective Investigation into Cancer and Nutrition cohort
|3=--[[Käyttäjä:Jouni|Jouni Tuomisto]] ([[Keskustelu käyttäjästä:Jouni|keskustelu]]) 26. tammikuuta 2015 kello 11.48 (UTC)}}

==Perustelut==

; Huomaa!
* Koulujen ja päiväkotien ruokailut vaikuttavat myös siihen, mitä lapset syövät. Näitä pitäisi selvittää.
* Miten selvitetään, mitä lapset oikeasti syövät pöydästä?
* Päiväkotien ruokailut?
* Mikä on silakan menekki koulussa, pitäisi selvittää jostain koulusta.

=== Tuloskuvat ===

:''26.1.2015 versiot kuvista ovat ne, jotka tulevat Eviran raporttiin. Huomaa, että kuvat saattavat päivittyä muiden hankkeiden yhteydessä senkin jälkeen, kun ''Silakan hyöty-riskiarvo'' -hanke on jo päättynyt.
<gallery widths=250px heights=200px>

File:Silakkakyselyvastaajat ikä vs sukupuoli.png|Silakkakyselyn vastaajat ikä vs sukupuoli
File:Silakkakyselyvastaajat ikä vs koulutus.png|Silakkakyselyn vastaajat ikä vs koulutus
File:Silakkakyselyvastaajat lasten iät.png|Silakkakyselyn vastaajien lasten iät
File:Silakansyönti ikäryhmittäin.png|Silakansyönnin keskiarvo ikäryhmittäin ja sukupuolittain
File:Silakansyönti Suomessa iän ja sukupuolen mukaan.png|Silakansyönnin liukuva keskiarvo iän mukaan sukupuolittain
File:Silakansyönnin jakauma Suomessa ikäryhmittäin.png|Silakansyönnin kumulatiivinen jakauma ikäryhmittäin ja sukupuolittain
File:Aineiden pitoisuuksia silakassa.png|Ravinto- ja haitta-aineiden pitoisuuksia silakassa
File:Altistuminen kalassa oleville aineille.png|Altsituminen silakassa oleville ravinto- ja haitta-aineille
File:Alistuminen silakan aineille sukupuolittain.png|Altistuminen silakan aineille sukupuolittain
File:Silakansyönti alueittain.png|Silakansyönti maakunnittain
File:Silakansyönnin DALYt ikäryhmittäin ja sukupuolittain.png|Silakansyönnin tautitaakka ikäryhmittäin ja sukupuolittain
File:Silakansyönnin DALYt nuorilla.png|Silakansyönnin tautitaakka ikäryhmittäin ja sukupuolittain alle 50-vuotiailla
File:Yksilöllinen riski silakansyönnistä iän_mukaan.png|Yksilöllinen riski silakansyönnistä iän mukaan
File:Yksilöllinen riski silakansyönnistä syönnin mukaan.png|Yksilöllinen riski silakansyönnistä syönnin mukaan
File:Olisiko silakansyönnin lopettamissuositus virhe.png|Toisiko silakan täyskielto terveyshyötyä?
File:Lisätiedon arvo silakkasuosituksen antamiseksi.png|Lisätiedon arvo silakkasuosituksen antamiseksi
File:Toisiko silakkarajoitus terveyshyötyä.png|Toisiko silakkarajoitus < 3 g/d terveyshyötyä
File:Silakkarajoituksen vaikutus tautitaakkaan nuorilla.png|Silakkarajoituksen vaikutus tautitaakkaan alle 50-vuotiailla
File:Silakkarajoituksen vaikutus tautitaakkaan.png|Silakkarajoituksen vaikutus tautitaakkaan

</gallery>

Vanhoja kuvia (ajettu vanhalla malliversiolla eivätkä siis ole ajan tasalla:
* [[:File:Dioksiinialtistuminen rannikolla ja sisämaassa.png]]
* [[:File:Dioksiinialtistuminen iän mukaan.png]]
* [[:File:PCDDF-altistus ikäryhmittäin.png]]
* [[:File:Kokonais-TEQ-altistus ikäryhmittäin.png]]
* [[:File:Omega3-rasvahappoaltistus ikäryhmittäin.png]]
* [[:File:D-vitamiinialtistus ikäryhmittäin.png]]
* [[:File:Silakansyönnin jakauma lapsilla.png]]
* [[:File:Silakansyönnin jakauma aikuisilla.png]]
* [[:File:Silakansyönnin DALYt ikäryhmittäin.png]]
* [[:File:SIlakansyönnin tautitaakka Suomessa nuorilla.png]]

=== Laske tulokset ===

:'' Aiempia laskentaharjoituksia on [http://fi.opasnet.org/fi_wiki/index.php?title=Silakan_hy%C3%B6ty-riskiarvio&oldid=28141#Laskenta arkistossa].

Syötemuuttujat käsitellään tällaisessa järjestyksessä (lopusta alkuun):
* out ulos ovariablesta
** out3 stepfunktiot
** out2 lineaariset funktiot
** out1 taustariskistä riippuvat funktiot
*** disincidence havaittu tautiriski
*** population
*** pci
**** RR yhdistelmä riskeistä
***** frexposed
***** RRor ei käytössä
****** ERF, threshold, exposure, bgexposure, BW
***** RRrhill relative Hill
****** ERF, threshold, exposure, bgexposure, BW
***** RRrr suhteellinen riski
****** ERF, threshold, exposure, bgexposure, BW
**** population

Laskee tuloksia lähtötiedoista alkaen ja tallentaa tulokset kantaan.  Tämän koodin ajaminen vie yli 5 min joten kannattaa käyttää ylempää koodia valmiille tuloksille. 

* [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=F1ih7JcxMh4bD7XJ Korjattu koodi 24.1.2015, 3000 iteraatiota] [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=IyWehjNgwxCYzAn3 Tulokset]

<rcode name="compute" label="Aja malli" embed=0 store=0>
library(OpasnetUtils)
library(ggplot2)
library(reshape2)
date()
openv.setN(300)

##########################################################################
# KYSELYTIEDOT SILAKASTA

#!!+++++++++++++++++++++++++++++++++++++++++++++++++++++
silakka <- opbase.data("Op_fi3831", subset = "Silakka") # [[:op_fi:Silakan hyöty-riskiarvio]]
#ii+++++++++++++++++++++++++++++++++++++++++++++++++++++

########## PREPROSESSOINTI

silakka$Paino[silakka$Paino < 4 & silakka$Ikäryhmä == "Aikuinen"] <- 75 #Tehdään karkea inputointi aikuisten painoon
silakka$Paino <- ifelse(silakka$Paino < 4 & silakka$Ikäryhmä == "Lapsi", 5+((60-5)/(15)*silakka$Ikä), silakka$Paino) 
# Lineaarinen ekstrapolaatio 5-60 kg

colnames(silakka)[colnames(silakka) == "Nro"] <- "Rivi"
levels(silakka$Result)[levels(silakka$Result) == "Ei syö silakkaa ollenkaan"] <- "En syö silakkaa ollenkaan"

rannikko <- c(
	"Uusimaa",
	"Pohjanmaa", 
	"Kymenlaakso",
	"Etelä-Pohjanmaa",
	"Satakunta",
	"Keski-Pohjanmaa",
	"Pohjois-Pohjanmaa", 
	"Varsinais-Suomi"
)

sisämaa <- c(
	"Kanta-Häme",
	"Pirkanmaa",
	"Etelä-Karjala",
	"Pohjois-Savo",
	"Pohjois-Karjala",
	"Etelä-Savo",
	"Keski-Suomi",
	"Päijät-Häme",
	"Lappi",
	"Kainuu"
)

silakka$Rannikko <- ifelse(silakka$Maakunta %in% rannikko, "Rannikko", "Sisämaa")

ages <- factor(c(
	"0", "1-4", "5-9", "10-14", "15-19", "20-24",
	"25-29", "30-34", "35-39", "40-44", "45-49", "50-54", 
	"55-59", "60-64", "65-69", "70-74", "75-79", "80-84", 
	"85-89", "90-94", "95-"), ordered = TRUE
)

silakka$Age <- cut(silakka$Ikä, breaks = c(0,1,(1:20)*5), labels = ages, right = FALSE)
silakka$Hedelm <- silakka$Ikä >= 20 & silakka$Ikä < 45 # Onko henkilö hedelmällisessä iässä?

# lyhyt = lyhyt lista yksilökohtaisia määrittelyjä eli vain välttämättömät.
lyhyt <- silakka[c("Age", "Ikä", "Hedelm", "Sukupuoli", "Rivi", "Maakunta")] 

kokoaik	<- Ovariable("kokoaik", data = data.frame(
	lyhyt, 
	Arvo = silakka$Result, 
	Result = silakka$Kokosilakka
)) # Kokonaisten silakoiden syönti kertaa/3 kk

osanaik	<- Ovariable("osanaik", data = data.frame(
	lyhyt,
	Result = silakka$Silakkaruoka
)) # Silakkaruokien syönti kertaa/3 kk

lisuaik	<- Ovariable("lisuaik", data = data.frame(
	lyhyt,
	Result = silakka$Silakkalisuke
)) # Silakkalisukkeiden syönti kertaa/3 kk

## Arvo yksi Iter jokaiselle kyselyn ihmiselle.

iterit <- Ovariable("iterit", data = data.frame(
	Iter = sample(get("N", envir = openv), nrow(silakka), replace = TRUE),
	Rivi = rownames(silakka),
	Result = 1
))

#### Arvo N riviä ehdot täyttävästä kyselyn osaryhmästä (tässä tapauksessa Hedelm+Sukupuoli-ryhmästä)

ehto <- unique(silakka[c("Sukupuoli", "Hedelm")]) # , "Age")]) Säästetään muistia #, "Maakunta")]) Eiköhän painokerroin huolehdi maakunnan

rivit <- data.frame()
for(i in 1:nrow(ehto)) {
	temp <- silakka[
		silakka$Sukupuoli == ehto$Sukupuoli[i] & silakka$Hedelm == ehto$Hedelm[i] , # Jätetään tästäkin Age pois ja Hedelm tilalle
		c("Sukupuoli", "Hedelm", "Painokerroin")
	] # Eikö tässä voisi yksinkertaistaa ja käyttää vain rivinumeroa eikä koko tempiä?
	temp <- temp[sample(1:nrow(temp), get("N", envir = openv), replace = TRUE, prob = temp$Painokerroin) , ]
	rivit <- rbind(rivit, data.frame(
		ehto[i , ],
		Iter = 1:get("N", envir = openv),
		Rivi = as.character(floor(as.numeric(rownames(temp)))),
		Result = 1
	))
}
rivit <- Ovariable(output = rivit, marginal = c(TRUE, TRUE, TRUE, TRUE, FALSE))

BW		<- rivit * Ovariable("BW", data = data.frame(lyhyt, Result = silakka$Paino)) # Ruumiinpaino
BW <- unkeep(BW, cols = c("Rivi", "Ikä"), prevresults = TRUE, sources = TRUE)

date()
######################################################################################
# SILAKOIDEN ANNOSKOOT

#!!+++++++++++++++++++++++++++++++++++++++++++++++++++
solet <- opbase.data("Op_fi3831.saantioletukset") # Silakkaoletukset sivulta [[:op_fi:Silakan hyöty-riskiarvio]]
#ii+++++++++++++++++++++++++++++++++++++++++++++++++++

silakoita	<- Ovariable("silakoita", data = solet[solet$Muuttuja == "V95" , c("Arvo", "Result")]) # Silakoita per silakka-annos

silakanpaino	<- Ovariable("silakanpaino", data = solet[solet$Muuttuja == "Koko silakan paino" , ]["Result"])

raakaainepaino	<- Ovariable("raakaainepaino", data = data.frame(
	Age = ages,
	Result = c(
		rep(solet$Result[solet$Muuttuja == "Raaka-ainesilakan paino, lapset"], 4), 
		rep(solet$Result[solet$Muuttuja == "Raaka-ainesilakan paino"], 17)
	)
))

lisukepaino	<- Ovariable("lisukepaino", data = solet[solet$Muuttuja == "Lisukesilakan paino" , "Result", drop = FALSE])

#!!++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Silakkaruokia (kpl/3 kk)
kouluruokailut <- Ovariable("kouluruokailut", data = data.frame(Result = "0;0;0.8;1.6"))

# Silakka-annoksen koko (g)
annos <- Ovariable("koulut", data = data.frame(Age = ages, Result = c(
	0,
	"20; 39; 40; 50; 60; 70; 80; 80; 100; 100; 100",
	"20; 39; 50; 65; 70; 70; 70; 80; 100; 120",
	"20; 39; 50; 65; 70; 70; 70; 80; 100; 120",
	"30; 39; 50; 80; 100; 110; 120; 120; 130; 160",
	rep(0, 16)))
)

# Syökö lapsi silakka-annoksensa vai ei?
into <- Ovariable("into", data = data.frame(Age = ages, Result = c(
	0,
	"1;1;1;1;1;1;1;1;1;1;0",
	"1;1;1;1;0",
	"1;1;1;1;0",
	"1;1;1;1;1;1;1;0;0;0",
	rep(0, 16)))
)
#ii++++++++++++++++++++++++++++++++++++++++++++++++++++

amount <- Ovariable("amount", 
	dependencies = data.frame(Name = c(
		"kokoaik",
		"silakoita",
		"silakanpaino",
		"osanaik",
		"raakaainepaino",
		"lisuaik",
		"lisukepaino",
		"kouluruokailut",
		"annos",
		"into"
	)),
	formula = function(...) {

		out <- (kokoaik * silakoita * silakanpaino + osanaik * raakaainepaino + lisuaik * lisukepaino) / 91 
		# Per 3 kk -> per d

		out <- out + kouluruokailut * annos * into / 91

		result(out)[result(out) == 0] <- 0.01 # Ei jätetä nollia saantiin
		return(out)
	}
)

########### Luodaan regressiomalleja varten data.frame, jossa joka yksilö mutta Iterin odotusarvo.

temp <- EvalOutput(amount, forceEval = TRUE, N = 200) # 200 iteraatiota riittää luotettavan odotusarvon laskemiseen.
temp <- oapply(temp, INDEX = temp@output["Rivi"], FUN = mean)@output # Silakansyönnin odotusarvo jokaiselle.
colnames(temp)[colnames(temp) == "amountResult"] <- "Silakkamäärä"
silakka <- merge(silakka, temp)

date()

########## Luodaan vaikutusarviointimallia varten ovariable, jossa arvotut yksilöt

riv <- rivit@output
riv$Result <- NULL

kokoaik@data <- merge(kokoaik@data, riv)
osanaik@data <- merge(osanaik@data, riv)
lisuaik@data <- merge(lisuaik@data, riv)

amount <- EvalOutput(amount, forceEval = TRUE)

# Muutetaan Age ja Maakunta epävarmaksi eli ei-marginaaliksi
amount@marginal[colnames(amount@output) %in% c("Age", "Maakunta")] <- FALSE 
# Sukupuoli ja Hedelm pidetään indekseinä koska niiden mukaan arvottiin

#### Luodaan Rajoita-päätös, joka toteutuessaan rajoittaa kaikkien silakansyönnin korkeintaan 3 g/d.

amount <- amount * Ovariable("rajoitus", data = data.frame(Rajoita = c("Kyllä", "Ei"), Result = 1))
amount <- unkeep(amount, cols = c("Rivi", "Ikä", "Arvo"), prevresults = TRUE, sources = TRUE)
result(amount) <- ifelse(result(amount) > 3 & amount@output$Rajoita == "Kyllä", 3, result(amount))

date()
#################################################################################
#  PITOISUUDET: PCDD/F, PCB, OMEGA3; JA MUUT TAPAUSKOHTAISET OVARIABLET

arvonta <- function(x, measure.vars) {
  out <- x[sample(1:nrow(x), get("N", envir = openv), replace = TRUE) , ]
  out <- data.frame(Iter = 1:nrow(out), out)
  out <- melt(
    out, 
    measure.vars = measure.vars, 
    variable.name = "Exposure_agent", 
    value.name = "Result"
  )
  return(out)
}

#!!++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# PCDD/F ja PCB-pitoisuudet # pg /g f.w. nd = LOQ
objects.get("x89WlgJlDLA02vnL") # PCDD/F-data table1
#ii++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

pcdd <- table1[2:nrow(table1), c(2,3)] # Aika tiputetaan pois toistaiseksi tarpeettomana (2002, 2009 dataa)
colnames(pcdd) <- c("PCDDF", "PCB")
pcdd$PCDDF <- as.numeric(as.character(pcdd$PCDDF))
pcdd$PCB <- as.numeric(as.character(pcdd$PCB))
pcdd$TEQ <- pcdd$PCDDF + pcdd$PCB
pcdd <- arvonta(pcdd, measure.vars = c("PCDDF", "PCB", "TEQ"))

#!!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
# Hae pitoisuustiedot: D-vitamiini, EPA, DHA, Omega-3. D-vitamiini: ug /100 g, rasvahapot: mg /100 g
objects.get("qhV9V72f9zLmXX50") # Ravintoainedata table
#ii+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

ravinteet <- table[2:nrow(table), 3:6]
colnames(ravinteet) <- c("Vitamin_D", "EPA", "DHA", "Omega3")
ravinteet <- arvonta(ravinteet, measure.vars = c("Vitamin_D", "EPA", "DHA", "Omega3"))
ravinteet$Result <- as.numeric(ravinteet$Result) / 100 # per 100 g -> per g

#!!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
#### Pitoisuustiedot metyylielohopeasta (tarkemmin sanottuna tämä on kokonaiselohopea).
objects.latest("Op_en4004", code_name = "initiate") # In mg /kg f.w. = ug /g f.w.
#ii+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

mehg <- concentration@data
colnames(mehg)[3] <- "MeHg"
mehg <- arvonta(mehg, "MeHg")
mehg <- mehg[! colnames(mehg) %in% c("Age", "Size")]

concentration <- Ovariable("concentration", data = rbind(ravinteet, pcdd, mehg))
date()
################ ALTISTUS

exposure <- amount * concentration

# temp <- oapply(exposure, INDEX = exposure@output[c("Sukupuoli", "Exposure_agent")], FUN = mean) # This is temporarily used
# to calculate the background intake from Baltic herring of vitamin D and omega3.

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

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

sumitem <- function(
	ova, #ovariable that has locations to sum
	cond, # index column that contains the locations to sum
	condvalue, # vector of locations to sum
	sumvalue # location to be given to the rows with the sums
) {
	d <- ova
	d@output <- d@output[d@output[[cond]] %in% condvalue , ]
	d <- oapply(d, cols = cond, FUN = sum)
	d@output[[cond]] <- sumvalue
	ova@output <- orbind(ova, d)
	return(ova)
}

addexposure <- sumitem(addexposure, "Exposure_agent", c("PCDDF", "PCB"), "TEQ")
addexposure <- sumitem(addexposure, "Exposure_agent", c("EPA", "DHA"), "Omega3")
addexposure <- unkeep(addexposure, prevresults = TRUE, sources = TRUE)

# Make the background exposure uncertain rather than an index.

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

addexposure <- addexposure * taustat

exposure <- exposure + addexposure

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

t0.5	<- Ovariable("t0.5", data = solet[solet$Muuttuja == "TEQ-puoliintumisaika" , ]["Result"])
f_ing	<- Ovariable("f_ing", data = solet[solet$Muuttuja == "TEQ-imeytymisosuus" , ]["Result"])
f_mtoc	<- Ovariable("f_mtoc", data = solet[solet$Muuttuja == "TEQ-osuus lapseen" , ]["Result"])
BF	<- Ovariable("BF", data = solet[solet$Muuttuja == "Imeväisen rasvamäärä" , ]["Result"])

temp <- exposure
temp@output <- temp@output[temp@output$Exposure_agent == "TEQ" , ]
temp <- log((temp * t0.5 * f_ing * f_mtoc / (log(2) * BF) ) + 1) # Actual conversion
temp@output$Exposure_agent <- "logTEQ"
temp <- unkeep(temp, prevresults = TRUE, sources = TRUE)
exposure@output <- orbind(exposure, temp)
exposure <- unkeep(exposure, prevresults = TRUE, sources = TRUE)

frexposed <- 1
bgexposure <- addexposure
date()
######################################################### VÄESTÖ

# Tarkastellaan totcases-laskennassa aluksi yksilöriskiä, ja vasta lopussa kerrotaan yksilötulokset yksilöiden edustamien ryhmien koolla.

population <- Ovariable("population", data = data.frame(Result = 1))

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

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

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

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

temp1 <- ERF
temp2 <- threshold

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

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

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

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

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

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

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

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

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


######################################## Tautiriski
#!!++++++++++++++++++++++++++++++++++++++++++++++++++
objects.latest("Op_en5917", code_name = "initiate") # [[:op_en:Disease risk]] ovariable disincidence
#ii++++++++++++++++++++++++++++++++++++++++++++++++++

disincidence <- disincidence / 100000 # Change from 1/100000py to 1/py.
disincidence <- unkeep(disincidence, cols = c("Age", "Sex", "Population", "Unit", "Response_metric")) 
# Näitä indeksejä ei tarvita tässä arvioinnissa.
disincidence@output <- disincidence@output[disincidence@output$Trait != "CHD" , ] # Tulee ikävakioidusta alta.

############################################################################
# Ikävakioidut kuolemansyyt

# Sairastuvuus (tarkemmin sanottuna [[Kuolemansyyt Suomessa]])
#!!++++++++++++++++++++++++++++++++++++++++++++++++++++
syy <- Ovariable("syy", data = opbase.data("Op_fi4558", subset = "2012", 
	include = list(Kuolemansyy = c(
		"04-21 Syövät (C00-C97)",
		"27 Iskeemiset sydäntaudit (I20-I25)",
		"29 Aivoverisuonien sairaudet (I60-I69)"
))))
#ii++++++++++++++++++++++++++++++++++++++++++++++++++++
syy@data$Ikä <- gsub(" ", "", syy@data$Ikä)
colnames(syy@data)[colnames(syy@data) == "Ikä"] <- "Age"

# Syöpää ei haluta linkata ikäryhmittäin ERF:iin, vaan sitä käytetään muodostamaan diskontattu elinaikainen
# odotettu eliniän menetys jos kuolee syöpään tietyn ikäisenä. Ks. myöhemmin.

väli <- Ovariable(output = data.frame(
  Kuolemansyy = c(
	rep("27 Iskeemiset sydäntaudit (I20-I25)", 2), 
	"04-21 Syövät (C00-C97)", 
	"29 Aivoverisuonien sairaudet (I60-I69)"
  ),
  Trait = c("CHD", "CHD2", "Syövät", "Stroke"), # Syövät nimetään hassusti, koska sen EI haluta linkkautuvan ERF:iin.
	# Lisäksi Tehdään kaksinkertainen sydäntaululistaus, koska on kaksi ERFiä: Mozaffarian ja Cohen.
  Result = 1
))

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

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

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

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

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

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

totcases <- EvalOutput(totcases)
date()
##################################################################################
# Tapauskohtaiset postprosessoinnit


### Muutetaan exposure yksikköön g /d kaikkien Exposure_agentien osalta.
skaala <- Ovariable(
	output = data.frame(
		Exposure_unit = "g /d", 
		Exposure_agent = c("Vitamin_D", "EPA", "DHA", "Omega3", "PCDDF", "PCB", "TEQ", "MeHg"), 
		Result = c(1E-6, 1E-3, 1E-3, 1E-3, 1E-12, 1E-12, 1E-12, 1E-6)),
	marginal = c(FALSE, TRUE, FALSE)
)

tiedot <- rivit * taustat * Ovariable(
	output = data.frame(
		silakka[c("Rivi", "Sukupuoli", "Age", "Hedelm", "Ikä", "Paino", "Silakkamäärä")],
		Result = 1
	),
	marginal = c(FALSE, TRUE, FALSE, TRUE, FALSE, FALSE, FALSE, FALSE)
)
tiedot <- unkeep(tiedot, cols = "Rivi", prevresults = TRUE, sources = TRUE)

date()

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

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

indivrisk <- totcases

indivrisk <- indivrisk * tiedot

result(indivrisk) <- result(indivrisk) * ifelse(
	indivrisk@output$Trait %in% c("Child's IQ", "Tooth defect", "Dental defect"), 
	ifelse(
		indivrisk@output$Sukupuoli == "Nainen" & indivrisk@output$Hedelm == "TRUE", 
		0.1, # Probability of birth during a year.
		0
	),
	1
)

#ages <- c(
#	"0", "1-4", "5-9", "10-14", "15-19", "20-24",
#	"25-29", "30-34", "35-39", "40-44", "45-49", "50-54", 
#	"55-59", "60-64", "65-69", "70-74", "75-79"
#)
indivrisk@output$Age <- factor(indivrisk@output$Age, levels = ages, ordered = TRUE)

indivrisk@output$Iter <- as.numeric(as.character(indivrisk@output$Iter))

#!!+++++++++++++++++++++++++++++++++++++++++++++++++
objects.store(indivrisk, pop, exposure, amount, ERF, threshold, silakka, disincidence, concentration, rivit)
#ii+++++++++++++++++++++++++++++++++++++++++++++++++
cat("indivrisk, pop, exposure, amount, ERF, threshold, silakka, disincidence, concentration, rivit stored. \n")
date()
</rcode>


Lasketaan jokaiselle kyselyn henkilölle oma silakansyönti, ja arvotaan vasta sen jälkeen. Tämä mahdollistaa myös sen, että tehdään LOESS-regressio datasta eikä bootstrap.

Koko kyselylomakkeen tietojen jäsennyksen voisi lohkaista omaksi mallikseen siten, että LOESS ajetaan N kertaa silakkamäärälle, jossa siinä on jo yhdistettynä kyselydata ja oletukset annosko'oista.  Näin saadaan parametrinen silakkamäärämalli, jota voidaan sämplätä haluttu määrä. Sen sijaan bootstrap ei sitten elegantisti toimisi, koska jo tätä preprosessointimallia ajettaessa pitäisi päättää iteraatioiden määrä ja tallentaa sämplätty määräovariable. Sen jälkeen voisi käyttää vain yhtä suuria tai pienempiä iteraatioita. Mutta tämä ei kyllä ole kummoinenkaan ongelma, jos alkuperäinen sample on riittävän iso.

Jos tehdään LOESS, pitää regressiossa olla kaikki ne parametrit, joita jatkossa käytetään. Näitä olisivat jatkuvina muuttujina Ikä, Sukupuoli, Paino ja Rannikko (neljä muuttujaa sallitaan). Sitten varsinaisessa mallissa käytetään predict.loess():ia. (Pitää opetella miten sellaisesta loess-objektista saa ovariablen.)

{{comment|# |Tällä hetkellä kuitenkin pitäydytään nykyisessä bootstrapissä, koska raporttiversio pitää saada matkaan 4.11.2014.|--[[Käyttäjä:Jouni|Jouni Tuomisto]] ([[Keskustelu käyttäjästä:Jouni|keskustelu]]) 2. marraskuuta 2014 kello 10.15 (UTC)}}

=== Avainmuuttujat ===

<t2b name="Avainmuuttujat" index="Muuttuja,Observation" locations="Aikuinen,Lapsi" unit="Sarake#">
Ikä|3|123
Sukupuoli|2|2
Maakunta|80|80
Paino|125|125
Painokerroin|91|91
Kyselyjakso|126|
Kokoaikaisuus|4|4
Ammatti|5|5
Siviilisääty|6|6
Tienaajaammatti|7|7
Tienaaja|8|8
Ruokapäättäjä|9|9
Elämäntilanne|10|10
Tuloluokka|30|30
Koulutus|40|40
Televisio|65|65
Asunto|75|75
Kuntakoko|76|76
Kokosilakka|92|126
Silakkaruoka|93|127
Silakkalisuke|94|128
Silakkalkm|95|129
Kalapuikot|109|
Tonnikala|110|
Kirjolohi|111|
Itämerenlohi|112|
Ahven|113|
Hauki|114|
Kuha|115|
Siika|116|
Purkkitonnikala|117|
Äyriäiset|118|
Muikku|119|
</t2b>

===Data===

; Malleissa käytettävät oletukset

* Muuttuja V95 (kysymystunnus SI2): Kuinka monta silakkaa syötte kerralla, kun syötte kokonaista silakkaa (kpl)?
* silakanpaino: yksittäisen silakkapihvin paino (kaksi silakkaa paistettuna vastakkain) (g).
* raakaainepaino: silakan paino annosta kohti, kun sitä käytetään raaka-aineena (esim. silakkalaatikossa) (g).
* lisukepaino: silakan paino annosta kohti, kun sitä käytetään lisukkeena (g).
* Silakkafileen keskipaino Finelin mukaan 15-30 g.<ref>[http://www.fineli.fi/foodunit.php?foodid=29207&lang=fi Fineli]</ref>. 

<t2b name='Saantioletukset' index='Muuttuja,Arvo' obs='Tulkinta' desc="Selite" unit='#,g /annos'>
V95|En syö silakkaa ollenkaan|0|
V95|1-2|0.5-2|
V95|3-5|3-5|
V95|6-10|6-10|
V95|Enemmän kuin 10|10-12|
V95|EOS|0-4|
Koko silakan paino||15-30|
Raaka-ainesilakan paino||20-40|
Raaka-ainesilakan paino, lapset||5-25|
Lisukesilakan paino||5-10|
TEQ-puoliintumisaika||2737.5|t0.5: 7.5 * 365, Dioxin half-life in adults (d)
TEQ-imeytymisosuus||0.8|f_ing: 0.80, Fraction of dioxin ingested that is actually absorbed
TEQ-osuus lapseen||0.1 - 0.25|f_mtoc: 0.25, Fraction of mother's dioxin burden that enters the child during breast-feeding. Alaluusua et Vartiainen
Imeväisen rasvamäärä||1580|BF: 0.2 * (7.692+8.141) / 2, body fat amount in the infant at six months of age (g) [[:op_en:Body weight of 0-24 months old Finnish children]].
</t2b>


<t2b name="Tausta-altistus" index="Tausta,Sukupuoli,Exposure_agent" desc="Yksikkö,Selite" unit="-">
Kyllä|Mies|Vitamin_D|11.7|µg /d|Finriski 12 - 0.3 silakasta
Kyllä|Nainen|Vitamin_D|8.5|µg /d|Finriski 8.7 - 0.2 silakasta
Kyllä|Mies|EPA|120|mg /d|Finriski 125 - 4.6 silakasta
Kyllä|Nainen|EPA|96|mg /d|Finriski 100 - 3.9 silakasta
Kyllä|Mies|DHA|118|mg /d|Finriski 125 - 6.7 silakasta
Kyllä|Nainen|DHA|94|mg /d|Finriski 100 - 5.4 silakasta
Kyllä||PCDDF|0|pg /d (TEQ)|
Kyllä||PCB|0|pg /d (TEQ)|
Kyllä||MeHg|0|µg /d|
Kyllä||logTEQ|0|log(pg /g)|
Ei|||0||
</t2b>

<pre>
Tämän silakkamallin tuottamat arviot keskimääräisestä silakan tuomasta altistuksesta. Tämä vähennetään Finriski-tutkimuksen
antamasta taustasaantiarviosta.
   Sukupuoli Exposure_agent      Result
1       Mies      Vitamin_D  0.25693803
2     Nainen      Vitamin_D  0.17339867
3       Mies            EPA  4.60868926
4     Nainen            EPA  3.87563627
5       Mies            DHA  6.74260515
6     Nainen            DHA  5.40637067
7       Mies         Omega3 15.66242617
8     Nainen         Omega3 12.95386085
9       Mies          PCDDF  9.33615627
10    Nainen          PCDDF  7.49031027
11      Mies            PCB  4.36167437
12    Nainen            PCB  3.50062804
13      Mies            TEQ 13.69783063
14    Nainen            TEQ 10.99093831
15      Mies           MeHg  0.04854054
16    Nainen           MeHg  0.03553385
</pre>

<t2b name="Laatupainotetut elinvuodet vaikutuksille" index="Trait,Response_metric" obs="DALY" desc="Lisätietoja" unit="DALY /case">
CHD|Mortality|5 - 15|Oletetaan DW 1 ja D 10 U 50%
CHD2|Mortality|5 - 15|Oletetaan DW 1 ja D 10 U 50%
Stroke|Mortality|5 - 15|Oletetaan DW 1 ja D 10 U 50 %
Developmental dental defects incl. agenesis|Yes/No according to "Developmental Defects of Enamel Index"|0 - 0.12|DW 0.001 D 60 U 100 %
Cancer|Lifetime probability|0 - 0.28|DW 0.1 D 20, lisäksi eliniän lyhenemä 5 a. Tämä syntyy elinikäisellä altistuksella joten oletetaan 1/50 syntyvän vuoden altistuksella. U 100 %
Vitamin D recommendation|1 if recommendation not met|0.0001 - 0.0101|DW 0.001 D 1 U 101x
Dioxin recommendation|1 if not met|0.0001 - 0.0101|DW 0.001 D 1 U 101x
Child's IQ|Change in IQ points|-0.0517 (-0.03 - -0.0817)|Intellectual disability, mild (IQ<70): 0.031 (0.018-0.049) From IHME. D 50 U from IHME.
Child´s IQ|Change in IQ points|-0.0517 (-0.03 - -0.0817)|Intellectual disability, mild (IQ<70): 0.031 (0.018-0.049) From IHME. D 50 U from IHME.
Tooth defect|Developmental defect: caries or missing tooth|0.008 (0.003 - 0.017)|Periodontitis weight from IHME. D: 1. U from IHME
Dental defect|Developmental defect: caries or missing tooth|0.008 (0.003 - 0.017)|Periodontitis weight from IHME. D: 1. U from IHME
</t2b>

*DW: disability weight (painokerroin: 0 ei haittaa, 1 kuolema)
*D: duration (vaikutuksen kesto vuosina)
* U: uncertainty (if > 100x from lower to upper, then assumes loguniform function)

Intellectual disability, mild (IQ<70): 0.031 (0.018-0.049) From GBD2010 (IHME) [http://ghdx.healthdata.org/record/global-burden-disease-study-2010-gbd-2010-disability-weights]. Assumed 50 years of condition and 30 points reduction to cause this disability weight.

Cohen ym. 2005<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. [http://www.ncbi.nlm.nih.gov/pubmed/16242602]</ref> käyttivät QALY-painokertoimia julkaisusta <ref>Nicole Mittmann, Kostas Trakas, Nancy Risebrough, Barbara A. Liu. Utility Scores for Chronic Conditions in a Community-Dwelling Population. PharmacoEconomics April 1999, Volume 15, Issue 4, pp 369-376 [http://link.springer.com/article/10.2165/00019053-199915040-00004]</ref>. Siinä oli kuitenkin vain 20 terveysvastetta ja kiinnostavia silakan kannalta syöpä, sydäntauti ja aivohalvaus. Näistä on käytettävissä myös DALY-painokertoimet.

; Käytettävissä olevat (?), hyödylliset datat

* [[:op_en:KTL sarcoma study|Sarkoomatutkimus]] (Jouni)
* Beneris: Salasanan takana olevia Beneris-tietoja ([http://www.opasnet.org/beneris/index.php/Beneris:Project_information Projektitietoa])
* Kalastajatutkimus (Anu)
* DIPP (Suvi)
* [http://www.thl.fi/fi_FI/web/fi/tutkimus/hankkeet/finriski Finriski]

==== Kouluruokailukysely ruokailuvastaaville ====

Tässä on kuvattu Taloustutkimuksen tekemän kouluruokakyselyn (21 ruokavastaavaa eri kunnista) tuloksia siltä osin, kuin ne ovat mallin kannalta tärkeitä. Luvut ovat vastaajien lukumääriä ja vastausten suuruuksia.

Kuinka usein silakkaa syödään?
* Ei ollenkaan 10
* Kerran kahdessa kuussa tai harvemmin 10
* Kerran kuussa 1

Mitä silakkaruokia on?
* Kokonaiset silakat (paistetut, silakkapihvit) 8
* Silakkalaatikko ym. 6
* Silakkaa lisukkeena 1

Keskimääräinen silakka-annos (g)?
* Päiväkodeissa 20, 39, 40, 50, 60, 70, 80, 80, 100, 100, 100
* Ala-asteella 20, 39, 50, 65, 70, 70, 70, 80, 100, 120
* Yläasteella 30, 39, 50, 80, 100, 110, 120, 120, 130, 160

Kuinka monta silakkaa per henkilö (kpl)?
* Päiväkodeissa 1, 1, 1, 2, 3
* Ala-asteella 1, 1, 2, 2
* Yläasteella 1, 2, 2, 2

Kuinka moni syö silakkaa?
* Päiväkodissa kaikki 6, suurin osa 5  = (6*1+5*0.8)/11 = 91 %
* Ala-asteella kaikki 4, suurin osa 4, puolet 1 alle puolet 1 = (4*1+4*0.8+1*0.5+1*0.3)/10 = 80 %
* Yläasteella kaikki 2, suurin osa 4, puolet 3, harvat 1 = (2*1+4*0.8+3*0.5+1*0.2)/10 = 69 %

Päätelmät mallia varten:
* kouluruoka: Oletetaan silakkaruokien runsaudeksi kolmessa kuukaudessa "0;0;0.8;1.6" eli jokainen listattu vaihtoehto on yhtä todennäköinen.
* annos: Oletetaan silakka-annoksen keskimääräiseksi kooksi nuo edellä mainitut seuraavin ikäryhmin:
** Age 1-4: päiväkodit
** Age 5-9: ala-aste
** Age 10-14: ala-aste
** Age 15-19: yläaste
* into: Oletetaan silakkaa syövien osuudeksi
** päiväkodeissa 1-4: "1;0.8"
** ala-asteella 5-14: "1;1;0.8;0.8;0.4"

=== Riippuvuudet ===

Vaikutusarviointimalli löytyy sivulta [[Silakan hyöty-riskiarvio#Laske tulokset]]. Se laskee ovariablet: indivrisk, pop, exposure, maara, ERF, threshold, tiedot, disincidence, totcases, maarakysely, concentration. Lähtötiedot tulevat useasta eri paikasta, ja ne on lueteltu tässä:
* [[Silakan hyöty-riskiarvio]] Silakkakyselydata
* [[Silakan hyöty-riskiarvio#Data]] Saantioletukset
* [[Silakan hyöty-riskiarvio#Data]] Taustasaannin määrä
* [[Silakan hyöty-riskiarvio#Data]] DALY-painokertoimet
* [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=x89WlgJlDLA02vnL] 236 riviä PCDD/F- ja PCB-TEQ-pitoisuuksia kopioitu EU-kalat-hankkeen Exceleistä.
* [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=qhV9V72f9zLmXX50] 61 riviä Evirasta saatua dataa EPA-, DHA-, omega3- ja D-vitamiinipitoisuuksista kopioitu Excelistä.
* [[:op_en:Mercury concentrations in fish in Finland]] 55 riviä elohopeadataa. Tiedot saatu EU-kalat II -hankkeesta 
* [[:op_en:Population of Finland]] Ikäryhmittäiset väkimäärät (miesten ja naisten oletetaan jakautuvan tasan joka ikäryhmässä)
* [[:op_en:ERF of dioxin]], ovariables ERF, threshold.  (Perustuu mm. tutkimukseen Alaluusua et al. 2004 (Seveso))
* [[:op_en:ERF of omega-3 fatty acids]], ovariables ERF, threshold.  (äidin raskaudenaikainen omega-3 rasvahappojen saanti vs. lapsen ÄO, Omega-3 rasvahapot vs. neurologiset häiriöt aikuisilla, etenkin aivoinfarkti, Omega-3 rasvahapot vs. sydäntautiriski aikuisilla)
* [[:op_en:ERFs of vitamins]], ovariables ERF, threshold
* [[:op_en:ERF of methylmercury]], ovariables ERF, threshold
* [[:op_en:Disease risk]] Tautiriski joillekin vasteille 
* [[Kuolemansyyt Suomessa]] Mutta sydäntautien ja aivohalvauksen ikäryhmittäiset kuolleisuusluvut tulevat täältä
* [[:op_en:Health impact assessment]] laskentakaavat tapausmäärille
* Mallissa myös oletetaan, että nainen saa 18-38 ikävuoden välillä keskimäärin 2.1 lasta. (Data upotettu koodiin.)

Vaikutusarviointimalli on hidas, ja siksi tulokset lasketaan ensin valmiiksi, ja kaikki DALY-laskennat ja kuvat lasketaan erikseen tällä koodilla: [[Silakan hyöty-riskiarvio#Näytä tuloksia]]. Tämä koodi tarvitsee ainoastaan nämä lähtötiedot:
* Vaikutusarviointimallin laskentatulokset komennolla objects.latest("Op_fi3831", code_name = "compute")
* [[Silakan hyöty-riskiarvio#Laatupainotetut elinvuodet vaikutuksille]] DALY-kertoimet

==== Katso myös (ei ole käytetty tässä arvioinnissa) ====

;Annosvasteita (englanniksi):
* [[:heande:Category:Review| Henna Karvosen Beneris-reviewt]]
* [http://www.efsa.europa.eu/en/efsajournal/pub/2813.htm Scientific Opinion on the Tolerable Upper Intake Level of vitamin D. EFSA Journal 2012;10(7):2813]
* [http://www.plosone.org/article/info%3Adoi%2F10.1371%2Fjournal.pone.0010530;jsessionid=49A1B8540750D6DC328960ACDB481914 Suominen-Taipale AL, Partonen T, Turunen AW, Männistö S, Jula A, Verkasalo PK: Fish Consumption and Omega-3 Polyunsaturated Fatty Acids in Relation to Depressive Episodes: A Cross-Sectional Analysis, PlosOne 2010, 5(5), e10530]
* [[:heande:All-cause mortality dose-response of fish consumption]]
* [http://www.ncbi.nlm.nih.gov/pmc/?term=10.1289/ehp.10399 Dioxin Exposure, from Infancy through Puberty, Produces Endocrine Disruption and Affects Human Semen Quality]

; Ravitsemustietoa yleisesti, mutta erityisesti kalasta: [[Kalansyönti Suomessa]]
; Pitoisuustietoa: [[Silakan ravintoainepitoisuudet]]

;Saantitietoa
* [[:heande:Estimated intake levels for Finnish children of methylmercury from fish]]
* [[:op_en:Fish oil intake in Beneris]]
* [[:beneris:Dioxin intakes in Finnish children]]

;Arviointeja
* [[:heande:Category:Review| Henna Karvosen Beneris-reviewt]]
* [http://www.slv.se/upload/dokument/rapporter/mat_naring/2007_12_fiskkonsumtion_risk_och_nytta.pdf Becker W, Darnerud PO, Petersson-Grawé K: Fiskkonsumtion – risk och nytta. Livsmedelsverket, Rapport 12/2007]
* [http://www.eufic.org/article/en/page/FTARCHIVE/artid/Dietary-Reference-Values-reference-for-whom/ Dietary Reference Values - A reference to whom?. Eufic, Food Today 04/2013]
* [http://www.fao.org/docrep/014/ba0136e/ba0136e00.pdf Joint FAO/WHO expert consultation on the risks and benefits of fish consumption, Rome, 25–29 January 2010, FAO Fisheries and Aquaculture Report No. 978, 2011]
* [http://www.vkm.no/dav/83cf7bd765.pdf Fish and seafood consumption in Norway – Benefits and risks, Norwegian scientific committee for food safety, March 2006]
* [http://www.norden.org/en/publications/publikationer/2008-568 Sand S, Becker W, Darnerud PO: Aspects of Risk-Benefit Assessment of Food Consumption,Directions for the future. TemaNord 2008:568]
* [http://www.sacn.gov.uk/pdfs/fics_sacn_advice_fish.pdf the Scientific Advisory Committee on Nutrition (SACN): Advice on fish consumption: benefits & risks. 2004]
* [http://www.efsa.europa.eu/en/supporting/pub/116e.htm EFSA Scientific Colloquium 6, 13-14 July 2006 - Tabiano, Italy: Risk-benefit analysis of foods: Methods and approaches.]
* [[:op_en:Benefit-risk assessment of fish consumption for Beneris]]

* [[Kuluttaja-arvio metyylielohopeasta ja omega-3-rasvahapoista kalassa]]

; Arviointimenetelmiä
* [[:op_en:Environmental burden of disease calculation]]

; Tautiriskejä
* [[:op_en:Risk of dental aberrations in children in Beneris]]
* [[:op_en:Total cancer risk in adults in Beneris]]
* [[:op_en:Intelligence quotient (IQ) in children in Beneris]]
* [[:op_en:Myocardial infarction risk in adults in Beneris]]
* [[:op_en:Remaining lifetime risk of CVD death in adults]]

==Katso myös==

Aiheesta on kirjoitettu myös raportti [[Silakan hyöty-riskiarvio / raportti]] sekä luento [[:Tiedosto:Kala ja terveys.ppt]]. Osa työstä sekä lisäaineistoja löytyy sivulta [[:heande:Silakan riskiarvio]].

; Aiempia malliajoja ja tuloksia
* Alustavia laskelmia ja hahmotelmia siitä, millaisessa muodossa tuloksia voisi esittää: [http://fi.opasnet.org/fi_wiki/index.php/Special:R-tools?id=gtSrX4eE9D12Of9H] Ajettu 09.08.2013 06:45:44
* [http://fi.opasnet.org/fi_wiki/index.php/Special:R-tools?id=RngqJaptYwZEVl2B Muuttujakuvaukset] Ajettu 07.01.2014 14:13:57
* [http://fi.opasnet.org/fi_wiki/index.php/Special:R-tools?id=DP2JmgvbY70DXNdJ Pitemmälle vietyjä laskelmia syyskuun ja joulukuun kyselytuloksista] Ajettu 07.01.2014 16:37:39, mm. regressioanalyysi ja vastaajaprofiilit
* [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=GyRdpSUBWnwkg0fo Regressioanalyysi 28.3.2014] [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=8kObRzacY4IGfBiC Päivitys 3.4.2014]
* [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=cSKlfueRfeTXOxGu Ajo sisältää terveysvaikutuslaskelmat ja tallennetun indivrisk-ovariablen] Ajettu 30.08.2014 21:25:46
* [http://fi.opasnet.org/fi_wiki/index.php?title=Silakan_hy%C3%B6ty-riskiarvio&oldid=26061 Malliversio sisältää myös epäonnistuneen syöpälaskennan]
* [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=LatLfvRmjrkVwAif Laskennassa mm. hammasvasteet bootstrapillä] Ajettu 02.11.2014 05:17:47
** [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=wPQ1pLFGLEI21gnG Ajon kuvaajat] Ajettu 02.11.2014 05:40:09
* [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=qe4YvVMcQCx2esgj Malli 12.11.2014] varsin lopullisessa muodossaan mutta vain 100 iteraatiota.
** [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=6D9ih0TQFx1YvWcJ Ajon kuvaajat] ajettu 13.11.2014 [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=VRgyoVfgpjS2CU7L] 19.11.2014
* Tulokset (28.12.2014), joissa oapplyssä ja totcasesssa olleet virheet on kaikki korjattu [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=XsQBGwOLgxuascO3]
* [http://fi.opasnet.org/fi_wiki/index.php?title=Toiminnot:RTools&id=8GSelUYVqZoUxnSl Tuloste, jossa taulukoissa ollut aggragointivirhe on korjattu] 7.1.2015.



* Aiemmassa versiossa oli seuraavanlaista tietoa, mutta se arkistoitiin koska sillä ei ollut käyttöä lopullisessa arvioinnissa:
** Kyselyn rakentamiseen liittyvä [http://fi.opasnet.org/fi_wiki/index.php?title=Silakan_hy%C3%B6ty-riskiarvio&oldid=25310#Kysymys kysymys]
** Taulukko kysymysten perusteista ja niistä seuraavista ratkaisuista.
** [http://fi.opasnet.org/fi_wiki/index.php?title=Silakan_hy%C3%B6ty-riskiarvio&oldid=25310#Data] Kalankulutustilastoja Suomesta &rarr; siirretty sivulle [[Kalansyönti Suomessa]]
** [http://fi.opasnet.org/fi_wiki/index.php?title=Silakan_hy%C3%B6ty-riskiarvio&oldid=25310#Kysymykset Silakkakyselyn alustavat kysymykset]. Lopussa on myös THL:n eri tutkimuksissa käytettyjä kalakysymyksiä.
** [http://fi.opasnet.org/fi_wiki/index.php?title=Silakan_hy%C3%B6ty-riskiarvio&oldid=25310#Muiden tutkimusten tuloksia Muiden tutkimusten tuloksia]. Sisältää lähinnä luettelon mahdollisesti relevanteista tutkimuksista ja avainhenkilöistä. Myös joitakin julkaisuja on linkattuna.

;Muuta hyödyllistä
* [http://www.sciencedirect.com/science/article/pii/S037842741200063X Is the fear of dioxin cancer more harmful than dioxin?] Jouko Tuomisto, Jouni T. Tuomisto: Toxicology Letters 2012.
* [https://www.julkari.fi/bitstream/handle/10024/110839/THL_RAP2013_016_%26sliitteet.pdf?sequence=1 Helldán et al. FINRAVINTO 2012 -tutkimus]
* [http://www.ncbi.nlm.nih.gov/pubmed/17047219 Mozaffarian D1, Rimm EB. JAMA. 2006 Oct 18;296(15):1885-99. Fish intake, contaminants, and human health: evaluating the risks and the benefits.]
* [http://www.evira.fi/portal/fi/tietoa+evirasta/asiakokonaisuudet/vierasaineet/tutkimukset+ja+projektit/eu-kalat+ii+-hanke/ EU kalat II - Itämeren tärkeimmät ympäristömyrkyt kalassa] 2009-2011.
* [[:op_en:Benefit-risk assessment on farmed salmon]]
* [http://www.ncbi.nlm.nih.gov/pmc/articles/PMC1567218/?page=4 Human Milk as a Source of Methylmercury Exposure in Infants.] Grandjean et al 1994.
* [[Hyöty-riskiarviointi metyylielohopeasta ja omega-3 rasvahapoista kalassa]]
* [http://www.efsa.europa.eu/en/search/doc/1673.pdf EFSA Journal 2010; 8(7):1673
Suggested citation: EFSA Scientific Committee; Guidance on human health risk-benefit assessment of food.] EFSA Journal 2010; 8(7)1673. [41 pp.]. doi:10.2093/j.efsa.2010.1673. Available online:  European Food Safety Authority, 2010 1. SCIENTIFIC OPINION. Guidance on human health risk-benefit assessment of foods. EFSA Scientific Committee2, 3. European Food Safety Authority (EFSA), Parma, Italy.
* [[:heande:Category:Fish]]
* [[:op_en:Beneris deliverables]]
* Food web Estonia: Mikä tämä on, pitäisi löytää kunnon linkki
* [http://www.slv.se/sv/grupp3/Pressrum/Nyheter/Pressmeddelanden/Nu-vet-fler-att-all-fisk-inte-ar-nyttig/ Nu vet fler att all fisk inte är nyttig] Livsmedelverket 28.11.2013
* [http://www.hs.fi/omaelama/artikkeli/Kirjolohi+on+parasta+terveydelle+mutta+ahven+onymp%C3%A4rist%C3%B6yst%C3%A4v%C3%A4llisin+ruoka/HS20060308SI1AU02gqo Kirjolohin on parasta terveydelle mutta ahven on ympäristöystävällisin ruoka] Helsingin Sanomat 8.3.2006
* [[Ven%C3%A4j%C3%A4n_vaalit_2011#Hierarkkinen_Bayes-malli_-_Hierarchial_Bayesian_model|Hierarchical Bayesian model of Russian duma election 2011]]
** [[:Tiedosto:RussianElections2011.ods|OpenBUGS model about election]] (does not work)
* [http://www.opasnet.org/eractest/index.php/Sandbox#JAGS_test rjags test] (password required)
* [http://www.opasnet.org/eractest/index.php/Sandbox#rbugs_test rbugs test] (password required) ([http://cran.r-project.org/web/packages/rjags/index.html rjags in CRAN])
* [http://www.slv.se/upload/dokument/rapporter/mat_naring/2007_12_fiskkonsumtion_risk_och_nytta.pdf Becker W, Darnerud PO, Petersson-Grawé K: Fiskkonsumtion – risk och nytta. Livsmedelsverket, Rapport 12/2007]
* [http://www.eufic.org/article/en/page/FTARCHIVE/artid/Dietary-Reference-Values-reference-for-whom/ Dietary Reference Values - A reference to whom?. Eufic, Food Today 04/2013]

==Viitteet==

<references/>

==Aiheeseen liittyviä tiedostoja==

* {{#l:Eviran_tutkimuksia_1_2015_silakan_hyöty-riskiarviointi.pdf}}
* {{#l:Kala ja terveys.ppt}}
* {{#l:Sopimus.pdf}}
* {{#l:Tietoisuus suosituksista2013.pdf}}
* {{#l:Linda Birnbaum EUssa2012.pdf}}
* {{#l:Linda Birnbaum Dioxin EDC 2012.pptx}}
* {{#l:Eviran_tutkimuksia_2_2015_raskasmetallit.pdf}}