Ajettu 19.10.2015 21:22:50+ Näytä viestit ja virheet- Piilota viestit ja virheet
stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
stat_bin: binwidth defaulted to range/30. Use 'binwidth = x' to adjust this.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
Warning message:
In grid.Call.graphics(L_polygon, x$x, x$y, index) :
semi-transparency is not supported on this device: reported only once per page
Warning messages:
1: Stacking not well defined when ymin != 0
2: Stacking not well defined when ymin != 0
3: Stacking not well defined when ymin != 0
4: Stacking not well defined when ymin != 0
Warning messages:
1: Stacking not well defined when ymin != 0
2: Stacking not well defined when ymin != 0
3: Stacking not well defined when ymin != 0
4: Stacking not well defined when ymin != 0
Warning messages:
1: Stacking not well defined when ymin != 0
2: Stacking not well defined when ymin != 0
3: Stacking not well defined when ymin != 0
4: Stacking not well defined when ymin != 0
Warning messages:
1: Stacking not well defined when ymin != 0
2: Stacking not well defined when ymin != 0
3: Stacking not well defined when ymin != 0
4: Stacking not well defined when ymin != 0
Warning messages:
1: Stacking not well defined when ymin != 0
2: Stacking not well defined when ymin != 0
Warning messages:
1: Stacking not well defined when ymin != 0
2: Stacking not well defined when ymin != 0
3: Stacking not well defined when ymin != 0
4: Stacking not well defined when ymin != 0
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
Warning message:
In grid.Call.graphics(L_polygon, x$x, x$y, index) :
semi-transparency is not supported on this device: reported only once per page
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
Warning message:
In grid.Call.graphics(L_polygon, x$x, x$y, index) :
semi-transparency is not supported on this device: reported only once per page
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
geom_smooth: method="auto" and size of largest group is >=1000, so using gam with formula: y ~ s(x, bs = "cs"). Use 'method = x' to change the smoothing method.
Warning message:
In grid.Call.graphics(L_polygon, x$x, x$y, index) :
semi-transparency is not supported on this device: reported only once per page
Warning messages:
1: Removed 13000 rows containing non-finite values (stat_density).
2: Removed 13000 rows containing non-finite values (stat_density).
3: Removed 12000 rows containing non-finite values (stat_density).
4: Removed 12000 rows containing non-finite values (stat_density).
5: Removed 12000 rows containing non-finite values (stat_density).
6: Removed 12000 rows containing non-finite values (stat_density).
7: Removed 12000 rows containing non-finite values (stat_density).
8: Removed 12000 rows containing non-finite values (stat_density).
9: Removed 33936 rows containing non-finite values (stat_density).
10: Removed 34000 rows containing non-finite values (stat_density).
Warning messages:
1: In if (class(get(i)) == "ovariable") { :
the condition has length > 1 and only the first element will be used
2: In if (class(get(i)) == "odecision") { :
the condition has length > 1 and only the first element will be used
3: In odag() : No edges found!
+ Näytä koodi- Piilota koodi
> 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_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")
Tässä ajossa on 5000 iteraatiota.
+ Näytä koodi- Piilota koodi
> ################ 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")
Silakansyönti ikäryhmittäin ja sukupuolittain.
+ Näytä koodi- Piilota koodi
> amosummary <- summary(amount, marginals = c("Sukupuoli", "Ikäluokka", "Rajoita"), fun = c("mean", "sd", "Q0.025", "Q0.975", "length"))
> oprint(amosummary)
+ Näytä koodi- Piilota koodi
> 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ä"])
dictionary fetched successfully!
+ Näytä koodi- Piilota koodi
> 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")
Regressioanalyysi silakansyönnin selittäjistä aikuisilla.
+ Näytä koodi- Piilota koodi
> fit <- lm(Silakkamäärä ~ Ikä + Sukupuoli + Paino + Rannikko, data = subset(silakka, Ikäryhmä == "Aikuinen"))
> oprint(summary(fit))
+ Näytä koodi- Piilota koodi
> fit <- lm(Silakkamäärä ~ Ikä + Sukupuoli + Paino + Maakunta, data = subset(silakka, Ikäryhmä == "Aikuinen"))
> oprint(summary(fit))
+ Näytä koodi- Piilota koodi
> fit <- lm(Silakkamäärä ~ Ikä + Sukupuoli + Paino + Rannikko, data = subset(silakka, Ikäryhmä == "Lapsi"))
> cat("Regressioanalyysi silakansyönnin selittäjistä lapsilla.\n")
Regressioanalyysi silakansyönnin selittäjistä lapsilla.
+ Näytä koodi- Piilota koodi
> oprint(summary(fit))
+ Näytä koodi- Piilota koodi
> #### 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))
+ Näytä koodi- Piilota koodi
> ######################### 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")
Silakansyönnin DALYt vasteen, sukupuolen ja politiikkojen mukaan.
+ Näytä koodi- Piilota koodi
> # 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, binwidth = 5) + 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)
> ggplot() +
+ geom_bar(data = subset(d2, Result > 0), aes(x = Age, y = Result/Freq, fill = Response), stat = "identity") +
+ geom_bar(data = subset(d2, Result <= 0), aes(x = Age, y = Result/Freq, fill = Response), stat = "identity") +
+ facet_grid(Background ~ Sex) + theme_gray(base_size = BS) +
+ labs(title = "Burden of disease due to Baltic herring in Finland",
+ y = "Life years lost (DALY/a)") + 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)
+ Näytä koodi- Piilota koodi
> oprint(summary(oapply(dalypinp, cols = c("Age", "Hedelm", "Sukupuoli"), FUN = sum), fun = c("mean", "sd", "Q0.025", "Q0.975")), digits = 1)
+ Näytä koodi- Piilota koodi
> oprint(summary(oapply(dalypinp, cols = c("Age", "Trait"), FUN = sum), fun = c("mean", "sd", "Q0.025", "Q0.975")), digits = 1)
+ Näytä koodi- Piilota koodi
> oprint(summary(oapply(dalypinp, cols = c("Hedelm", "Trait"), FUN = sum), fun = c("mean", "sd", "Q0.025", "Q0.975")), digits = 1)
+ Näytä koodi- Piilota koodi
> ######## 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")
Yli 0.7 pg/kg/d dioksiineille altistuvat:
+ Näytä koodi- Piilota koodi
> 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")))
[1] 0.08622909
+ Näytä koodi- Piilota koodi
> odag()
NULL
