Silakan hyöty-riskiarvio

Edistymisluokitus

Opasnetissa lukuisat sivut ovat työn alla eri vaiheissa. Niiden tietosisältöön pitää siis suhtautua harkiten. Tämän sivun sisällön edistyminen on arvioitu:

Tämä sivu on täysluonnos: Sivu on kirjoitettu kertaalleen alusta loppuun eli kaikki olennaiset sisällöt ovat jo paikoillaan. Kuitenkaan sisältöjä ei ole vielä kunnolla tarkistettu, ja esimerkiksi tärkeitä viitteitä voi puuttua.

Tämä sivu on tutkimus. Sivutunniste: Op_fi3831
Moderaattori:Jouni (katso kaikki) Sivun edistyminen: Täysluonnos. Arvostuksen määrää ei ole arvioitu (ks. peer review).
Lisää dataa {{#opasnet_base_link:Op_fi3831}}

Aiheesta on kirjoitettu myös raportti Silakan hyöty-riskiarvio / raportti Eviran_tutkimuksia_1_2015_silakan_hyöty-riskiarviointi 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
Malli jossa silakkaodotusarvot korjattu, MeHg lisätty ja Iter = 5000 24.1.2015. Pysyvä ajo. Tuloksia
Sellaisia tuloksia 26.1.2015, joiden laskemiseen on tarvittu ei-julkista dataa.

Näytä tuloksia

Malliajo 15.10.2015 jossa oliot pysyvästi tallennettu. Tuloksia 19.10.2015, mukana muutama uusi englanninkielinen kuva
op_en:File:BRA of Baltic herring.pptx esitys englanniksi aiheesta

+ Näytä koodi - Piilota koodi

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

⇤--#: . Kun vastauksia tulkitaan, tulee huomata, että seuraavia tuoreita tietoja EI OLE KÄSITELTY.

Seychell-tutkimus: omega-3 estää elohopeavaikutusta

EFSAn hyöty-riskiarvio omega-3:sta ja elohopeasta, lehdistötiedote

EFSAn hyöty-riskiarvio omega-3:sta ja elohopeasta

Engeset et al. 2014: Fish consumption and mortality in the European Prospective Investigation into Cancer and Nutrition cohort --Jouni Tuomisto (keskustelu) 26. tammikuuta 2015 kello 11.48 (UTC) (type: truth; paradigms: science: attack)

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.

Silakkakyselyn vastaajat ikä vs sukupuoli
Silakkakyselyn vastaajat ikä vs koulutus
Silakkakyselyn vastaajien lasten iät
Silakansyönnin keskiarvo ikäryhmittäin ja sukupuolittain
Silakansyönnin liukuva keskiarvo iän mukaan sukupuolittain
Silakansyönnin kumulatiivinen jakauma ikäryhmittäin ja sukupuolittain
Ravinto- ja haitta-aineiden pitoisuuksia silakassa
Altsituminen silakassa oleville ravinto- ja haitta-aineille
Altistuminen silakan aineille sukupuolittain
Silakansyönti maakunnittain
Silakansyönnin tautitaakka ikäryhmittäin ja sukupuolittain
Silakansyönnin tautitaakka ikäryhmittäin ja sukupuolittain alle 50-vuotiailla
Yksilöllinen riski silakansyönnistä iän mukaan
Yksilöllinen riski silakansyönnistä syönnin mukaan
Toisiko silakan täyskielto terveyshyötyä?
Lisätiedon arvo silakkasuosituksen antamiseksi
Toisiko silakkarajoitus < 3 g/d terveyshyötyä
Silakkarajoituksen vaikutus tautitaakkaan alle 50-vuotiailla
Silakkarajoituksen vaikutus tautitaakkaan

Vanhoja kuvia (ajettu vanhalla malliversiolla eivätkä siis ole ajan tasalla:

Laske tulokset

Aiempia laskentaharjoituksia on 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.

Korjattu koodi 24.1.2015, 3000 iteraatiota Tulokset

+ Näytä koodi - Piilota koodi

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

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

----#: . Tällä hetkellä kuitenkin pitäydytään nykyisessä bootstrapissä, koska raporttiversio pitää saada matkaan 4.11.2014. --Jouni Tuomisto (keskustelu) 2. marraskuuta 2014 kello 10.15 (UTC) (type: truth; paradigms: science: comment)

Avainmuuttujat

Avainmuuttujat(Sarake#)
Obs	Muuttuja	Aikuinen	Lapsi
1	Ikä	3	123
2	Sukupuoli	2	2
3	Maakunta	80	80
4	Paino	125	125
5	Painokerroin	91	91
6	Kyselyjakso	126
7	Kokoaikaisuus	4	4
8	Ammatti	5	5
9	Siviilisääty	6	6
10	Tienaajaammatti	7	7
11	Tienaaja	8	8
12	Ruokapäättäjä	9	9
13	Elämäntilanne	10	10
14	Tuloluokka	30	30
15	Koulutus	40	40
16	Televisio	65	65
17	Asunto	75	75
18	Kuntakoko	76	76
19	Kokosilakka	92	126
20	Silakkaruoka	93	127
21	Silakkalisuke	94	128
22	Silakkalkm	95	129
23	Kalapuikot	109
24	Tonnikala	110
25	Kirjolohi	111
26	Itämerenlohi	112
27	Ahven	113
28	Hauki	114
29	Kuha	115
30	Siika	116
31	Purkkitonnikala	117
32	Äyriäiset	118
33	Muikku	119

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.^[1].

Saantioletukset(#,g /annos)
Obs	Muuttuja	Arvo	Tulkinta	Selite
1	V95	En syö silakkaa ollenkaan	0
2	V95	1-2	0.5-2
3	V95	3-5	3-5
4	V95	6-10	6-10
5	V95	Enemmän kuin 10	10-12
6	V95	EOS	0-4
7	Koko silakan paino		15-30
8	Raaka-ainesilakan paino		20-40
9	Raaka-ainesilakan paino, lapset		5-25
10	Lisukesilakan paino		5-10
11	TEQ-puoliintumisaika		2737.5	t0.5: 7.5 * 365, Dioxin half-life in adults (d)
12	TEQ-imeytymisosuus		0.8	f_ing: 0.80, Fraction of dioxin ingested that is actually absorbed
13	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
14	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]].

Tausta-altistus(-)
Obs	Tausta	Sukupuoli	Exposure_agent	Result	Yksikkö	Selite
1	Kyllä	Mies	Vitamin_D	11.7	µg /d	Finriski 12 - 0.3 silakasta
2	Kyllä	Nainen	Vitamin_D	8.5	µg /d	Finriski 8.7 - 0.2 silakasta
3	Kyllä	Mies	EPA	120	mg /d	Finriski 125 - 4.6 silakasta
4	Kyllä	Nainen	EPA	96	mg /d	Finriski 100 - 3.9 silakasta
5	Kyllä	Mies	DHA	118	mg /d	Finriski 125 - 6.7 silakasta
6	Kyllä	Nainen	DHA	94	mg /d	Finriski 100 - 5.4 silakasta
7	Kyllä		PCDDF	0	pg /d (TEQ)
8	Kyllä		PCB	0	pg /d (TEQ)
9	Kyllä		MeHg	0	µg /d
10	Kyllä		logTEQ	0	log(pg /g)
11	Ei			0

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

Laatupainotetut elinvuodet vaikutuksille(DALY /case)
Obs	Trait	Response_metric	DALY	Lisätietoja
1	CHD	Mortality	5 - 15	Oletetaan DW 1 ja D 10 U 50%
2	CHD2	Mortality	5 - 15	Oletetaan DW 1 ja D 10 U 50%
3	Stroke	Mortality	5 - 15	Oletetaan DW 1 ja D 10 U 50 %
4	Developmental dental defects incl. agenesis	Yes/No according to "Developmental Defects of Enamel Index"	0 - 0.12	DW 0.001 D 60 U 100 %
5	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 %
6	Vitamin D recommendation	1 if recommendation not met	0.0001 - 0.0101	DW 0.001 D 1 U 101x
7	Dioxin recommendation	1 if not met	0.0001 - 0.0101	DW 0.001 D 1 U 101x
8	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.
9	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.
10	Tooth defect	Developmental defect: caries or missing tooth	0.008 (0.003 - 0.017)	Periodontitis weight from IHME. D: 1. U from IHME
11	Dental defect	Developmental defect: caries or missing tooth	0.008 (0.003 - 0.017)	Periodontitis weight from IHME. D: 1. U from IHME

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) [3]. Assumed 50 years of condition and 30 points reduction to cause this disability weight.

Cohen ym. 2005^[2] käyttivät QALY-painokertoimia julkaisusta ^[3]. 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

Sarkoomatutkimus (Jouni)
Beneris: Salasanan takana olevia Beneris-tietoja (Projektitietoa)
Kalastajatutkimus (Anu)
DIPP (Suvi)
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
[4] 236 riviä PCDD/F- ja PCB-TEQ-pitoisuuksia kopioitu EU-kalat-hankkeen Exceleistä.
[5] 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)

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

Saantitietoa

Arviointeja

Kuluttaja-arvio metyylielohopeasta ja omega-3-rasvahapoista kalassa

Arviointimenetelmiä

op_en:Environmental burden of disease calculation

Tautiriskejä

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ää: [6] Ajettu 09.08.2013 06:45:44
Muuttujakuvaukset Ajettu 07.01.2014 14:13:57
Pitemmälle vietyjä laskelmia syyskuun ja joulukuun kyselytuloksista Ajettu 07.01.2014 16:37:39, mm. regressioanalyysi ja vastaajaprofiilit
Regressioanalyysi 28.3.2014 Päivitys 3.4.2014
Ajo sisältää terveysvaikutuslaskelmat ja tallennetun indivrisk-ovariablen Ajettu 30.08.2014 21:25:46
Malliversio sisältää myös epäonnistuneen syöpälaskennan
Laskennassa mm. hammasvasteet bootstrapillä Ajettu 02.11.2014 05:17:47
- Ajon kuvaajat Ajettu 02.11.2014 05:40:09
Malli 12.11.2014 varsin lopullisessa muodossaan mutta vain 100 iteraatiota.
- Ajon kuvaajat ajettu 13.11.2014 [7] 19.11.2014
Tulokset (28.12.2014), joissa oapplyssä ja totcasesssa olleet virheet on kaikki korjattu [8]
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ä kysymys
- Taulukko kysymysten perusteista ja niistä seuraavista ratkaisuista.
- [9] Kalankulutustilastoja Suomesta → siirretty sivulle Kalansyönti Suomessa
- Silakkakyselyn alustavat kysymykset. Lopussa on myös THL:n eri tutkimuksissa käytettyjä kalakysymyksiä.
- 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ä

Is the fear of dioxin cancer more harmful than dioxin? Jouko Tuomisto, Jouni T. Tuomisto: Toxicology Letters 2012.
Helldán et al. FINRAVINTO 2012 -tutkimus
Mozaffarian D1, Rimm EB. JAMA. 2006 Oct 18;296(15):1885-99. Fish intake, contaminants, and human health: evaluating the risks and the benefits.
EU kalat II - Itämeren tärkeimmät ympäristömyrkyt kalassa 2009-2011.
op_en:Benefit-risk assessment on farmed salmon
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
Nu vet fler att all fisk inte är nyttig Livsmedelverket 28.11.2013
Kirjolohin on parasta terveydelle mutta ahven on ympäristöystävällisin ruoka Helsingin Sanomat 8.3.2006
Hierarchical Bayesian model of Russian duma election 2011
- OpenBUGS model about election (does not work)
rjags test (password required)
rbugs test (password required) (rjags in CRAN)
Becker W, Darnerud PO, Petersson-Grawé K: Fiskkonsumtion – risk och nytta. Livsmedelsverket, Rapport 12/2007
Dietary Reference Values - A reference to whom?. Eufic, Food Today 04/2013

Viitteet

↑ Fineli
↑ 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. [1]
↑ 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 [2]

Aiheeseen liittyviä tiedostoja

[1] Fineli

[2] 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. [1]

[3] 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 [2]

[1]

[2]

[3]