R-tools

> wiki_username <- "Jouni"
> server <- 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_fi3831", code_name = "compute") # Haetaan valmiiksi laskettu indivrisk, pop ja exposure.
> #objects.get("fLcPm8KyLhUuZC0c") # Haetaan tilapäinen versio, jossa on yli 3 g/d saannit typistetty.
> #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@output$Ikä <- as.numeric(as.character(amount@output$Ikä))
> amount@output$Ikäluokka <- cut(amount@output$Ikä, c(0, 20, 65, 80), include.lowest = TRUE)
> amount@output$Sukupuoli <- factor(amount@output$Sukupuoli, levels = c("Nainen", "Mies"))
> 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"))
> levels(amosummary$Ikäluokka) <- c("0-20", "21-65", "66-80")
> 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"
+ )
> 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))

Call:
lm(formula = Silakkamäärä ~ Ikä + Sukupuoli + Paino + Rannikko, 
    data = subset(silakka, Ikäryhmä == "Aikuinen"))
Residuals:
    Min      1Q  Median      3Q     Max 
 -2.898  -1.511  -0.942  -0.096 160.923 
Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)     -0.124397   0.617355  -0.202 0.840327    
Ikä              0.026126   0.006681   3.910 9.52e-05 ***
SukupuoliMies    0.172041   0.274870   0.626 0.531453    
Paino            0.007703   0.008162   0.944 0.345397    
RannikkoSisämaa -0.957231   0.246429  -3.884 0.000106 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.499 on 2037 degrees of freedom
Multiple R-squared:  0.01603,	Adjusted R-squared:  0.0141 
F-statistic: 8.297 on 4 and 2037 DF,  p-value: 1.23e-06

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

Call:
lm(formula = Silakkamäärä ~ Ikä + Sukupuoli + Paino + Maakunta, 
    data = subset(silakka, Ikäryhmä == "Aikuinen"))
Residuals:
    Min      1Q  Median      3Q     Max 
 -3.451  -1.567  -0.860   0.058 160.259 
Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)               -1.454100   0.839418  -1.732   0.0834 .  
Ikä                        0.026864   0.006726   3.994 6.73e-05 ***
SukupuoliMies              0.145004   0.275620   0.526   0.5989    
Paino                      0.008414   0.008185   1.028   0.3041    
MaakuntaEtelä-Pohjanmaa    1.072769   0.781956   1.372   0.1702    
MaakuntaEtelä-Savo        -0.263332   1.177478  -0.224   0.8231    
MaakuntaKainuu            -0.260139   0.849460  -0.306   0.7595    
MaakuntaKanta-Häme         0.307080   0.968020   0.317   0.7511    
MaakuntaKeski-Pohjanmaa    1.177495   1.008773   1.167   0.2432    
MaakuntaKeski-Suomi       -0.058500   0.788285  -0.074   0.9408    
MaakuntaKymenlaakso        1.043907   0.993823   1.050   0.2937    
MaakuntaLappi             -0.456959   1.164465  -0.392   0.6948    
MaakuntaPäijät-Häme       -0.068317   1.424061  -0.048   0.9617    
MaakuntaPirkanmaa          1.494552   0.700294   2.134   0.0329 *  
MaakuntaPohjanmaa          0.418393   1.363158   0.307   0.7589    
MaakuntaPohjois-Karjala    0.204066   0.803862   0.254   0.7996    
MaakuntaPohjois-Pohjanmaa  1.892090   0.739405   2.559   0.0106 *  
MaakuntaPohjois-Savo      -0.216787   0.728201  -0.298   0.7660    
MaakuntaSatakunta          1.394953   0.783358   1.781   0.0751 .  
MaakuntaUusimaa            1.220247   0.634451   1.923   0.0546 .  
MaakuntaVarsinais-Suomi    1.004783   0.723623   1.389   0.1651    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.498 on 2021 degrees of freedom
Multiple R-squared:  0.02413,	Adjusted R-squared:  0.01447 
F-statistic: 2.499 on 20 and 2021 DF,  p-value: 0.0002538

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

Call:
lm(formula = Silakkamäärä ~ Ikä + Sukupuoli + Paino + Rannikko, 
    data = subset(silakka, Ikäryhmä == "Lapsi"))
Residuals:
   Min     1Q Median     3Q    Max 
-1.509 -0.709 -0.391 -0.077 58.422 
Coefficients:
                Estimate Std. Error t value Pr(>|t|)  
(Intercept)      0.79276    0.30752   2.578   0.0101 *
Ikä              0.08800    0.06340   1.388   0.1655  
SukupuoliMies   -0.54956    0.25436  -2.161   0.0310 *
Paino           -0.01566    0.01744  -0.898   0.3693  
RannikkoSisämaa  0.35539    0.25537   1.392   0.1644  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.661 on 846 degrees of freedom
Multiple R-squared:  0.009633,	Adjusted R-squared:  0.004951 
F-statistic: 2.057 on 4 and 846 DF,  p-value: 0.08456

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

Call:
glm(formula = Suurkuluttaja ~ Ikä + Sukupuoli + Paino + Rannikko, 
    family = binomial, data = silakka)
Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-0.5374  -0.2545  -0.1721  -0.1159   3.5267  
Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
(Intercept)     -5.410008   0.476066 -11.364  < 2e-16 ***
Ikä              0.030484   0.006705   4.546 5.46e-06 ***
SukupuoliMies    0.471457   0.271016   1.740 0.081931 .  
Paino            0.006683   0.006660   1.003 0.315638    
RannikkoSisämaa -1.035546   0.289948  -3.571 0.000355 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Dispersion parameter for binomial family taken to be 1)
    Null deviance: 659.31  on 2892  degrees of freedom
Residual deviance: 600.75  on 2888  degrees of freedom
AIC: 610.75
Number of Fisher Scoring iterations: 7

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