+ Näytä koodi- Piilota koodi
# This is code Op_fi5923/ on page [[Ruori]]
library(OpasnetUtils)
library(ggplot2)
library(thlGraphs)
library(plotly)
openv.setN(10)
objects.latest("Op_fi5923", code_name="model")
# First empty all objects for a fresh start. Otherwise may be problems with CheckDecisions.
oempty(all=TRUE)
InpBoD <- EvalOutput(InpBoD)
InpPAF <- EvalOutput(InpPAF)
utility <- EvalOutput(utility, verbose=TRUE)
# Sample from default and sensitivity scenario about lead threshold.
cat("Elintarvikeperäisen lyijyn vaikutukset herkkyystarkastelussa.\n")
oprint(summary(BoDattr[BoDattr$Exposure_agent=="Lyijy",]))
BoDattr <- CollapseMarginal(BoDattr,"Threshold","sample")
utility <- CollapseMarginal(utility,"Threshold","sample")
levels(BoDattr$Exposure_agent)[levels(BoDattr$Exposure_agent)=="Vihannesvaje"] <- "Kasvisvaje"
#levels(BoDattr$Response)
#[1] "Cancer morbidity" "IQ loss" "Listeriosis" "Liver cancer"
#[5] "Noro infection" "Sperm concentration" "Toxoplasmosis" "Yes or no dental defect"
#[9] "CHD death" "Diet high in sodium" "Diet low in fruits" "Diet low in vegetables"
levels(BoDattr$Response) <- c(
"Syöpä",
"Älykkyysosamäärän lasku",
"Listerioosi",
"Maksasyöpä",
"Noroinfektio",
"Miehen hedelmättömyys",
"Toksoplasmoosi",
"Hammasvaurio",
"Sydäntauti",
"Liika suola",
"Hedelmävaje",
"Kasvisvaje"
)
cat("Elintarvikeperäisiä tautitaakkoja Suomessa arpoen lyijylle oletetun tai matalamman kynnysarvon.\n")
tmp <- summary(oapply(BoDattr[BoDattr$Scenario=="BAU",],NULL,sum,c("Age","Response")))
oprint(data.frame(
Altiste = tmp$Exposure_agent,
Keskiarvo = signif(tmp$mean,2),
"95 luottamusväli" = paste0(signif(tmp$Q0.025,2)," - ", signif(tmp$Q0.975,2)),
Keskihajonta = signif(tmp$sd,2)
)[rev(match(lev, tmp$Exposure_agent)),])
cat("Elintarvikeperäisiä tautitaakkoja Suomessa ikä- ja tautiryhmittäin arpoen lyijylle oletetun tai matalamman kynnysarvon.\n")
tmp <- summary(BoDattr[BoDattr$Scenario=="BAU",])
oprint(data.frame(
Altiste = tmp$Exposure_agent,
Ikä = tmp$Age,
Vaste = tmp$Response,
Keskiarvo = signif(tmp$mean,2),
Mediaani = signif(tmp$median,2),
"95 luottamusväli" = paste0(signif(tmp$Q0.025,2)," - ", signif(tmp$Q0.975,2)),
Keskihajonta = signif(tmp$sd,2)
))
cat("Ruori-skenaarioiden vaikutus tautitaakkaan\n")
tmp <- summary(utility)
oprint(data.frame(
Altiste = tmp$Exposure_agent,
Keskiarvo = signif(tmp$mean,2),
"95 luottamusväli" = paste0(signif(tmp$Q0.025,2)," - ", signif(tmp$Q0.975,2)),
Keskihajonta = signif(tmp$sd,2)
)[rev(match(lev, tmp$Exposure_agent)),])
dodge <- position_dodge(width=0.7)
if(FALSE) {
gg <- ggplot(summary(oapply(BoDattr[BoDattr$Scenario=="BAU",],NULL,sum,"Age")),
aes(x=Exposure_agent, weight=unlist(mean), fill=Response))+geom_bar()+
theme(legend.position = "bottom")+
labs(
title="Elintarvikkeiden tautitaakkoja Suomessa",
subtitle="Haittapainotettua elinvuotta vuodessa (DALY/a)")+
coord_flip()
print(gg)
gg <- ggplot(summary(oapply(BoDattr, NULL, sum,c("Age","Response"))),
aes(x=Exposure_agent, weight=unlist(mean), fill=Scenario))+geom_bar(position="dodge")+
coord_flip(ylim=c(0,70000))+
labs(
title="Elintarvikeperäisiä tautitaakkoja Suomessa",
subtitle="Haittapainotettua elinvuotta vuodessa (DALY/a)")+
geom_errorbar(aes(ymin=unlist(Q0.025),ymax=unlist(Q0.975),group=Scenario),position=dodge, width=0.3)+
geom_text(aes(label=signif(unlist(mean),2), y=unlist(Q0.975)+5000, group=Scenario), position=dodge)
print(gg)
# Utility of actions
gg <- ggplot(summary(utility),aes(x=Exposure_agent, weight=unlist(mean)))+geom_bar(fill="lightblue")+
coord_flip(ylim=c(-9000,0))+
labs(
title="Ruori-skenaarioiden vaikutus tautitaakkaan",
subtitle="Haittapainotettua elinvuotta vuodessa (DALY/a)")+
geom_errorbar(aes(ymin=unlist(Q0.025),ymax=unlist(Q0.975)), width=0.3)+
geom_text(aes(label=signif(unlist(mean),2), y=unlist(Q0.025)-600))
print(gg)
} else {
###### RUN THESE ON OWN COMPUTER WITH thlGraphs PACKAGE
# levels(BoDattr$Exposure_agent)
# [1] "Aflatoksiini" "Dioksiini" "Norovirus" "Toksoplasma"
# [5] "Lyijy" "Listeria" "Tyydyttynyt rasva" "Vihannesvaje"
# [9] "Suola" "Hedelmävaje"
levels(BoDattr$Exposure_agent) <- c(
"Aflatoxin","Dioxin","Noro virus","Toxoplasma", "Lead","Listeria",
"Saturated fat","Lack of vegetables","Sodium","Lack of fruits")
thlBarPlot(summary(oapply(BoDattr[BoDattr$Scenario=="BAU",],NULL,sum,"Age")),
xvar=Exposure_agent, yvar=unlist(mean), groupvar=Response, legend.position="bottom",
colors=thlColors(n=12),
title="Elintarvikkeiden tautitaakkoja Suomessa", subtitle="Haittapainotettua elinvuotta vuodessa (DALY/a)")+coord_flip()
thlBarPlot(summary(oapply(BoDattr, NULL, sum,c("Age","Response"))),xvar=Exposure_agent, yvar=unlist(mean),
groupvar=Scenario,stacked=FALSE, ylimits=c(0,70000),
# title="Elintarvikeperäisiä tautitaakkoja Suomessa",
# subtitle="Haittapainotettua elinvuotta vuodessa (DALY/a)")+
title="Burden of disease of selected food-mediated risk factors in Finland",
subtitle="Disability-adjusted life years per year (DALY/a)")+
coord_flip()+
geom_errorbar(aes(ymin=unlist(Q0.025),ymax=unlist(Q0.975),group=Scenario),position=dodge, width=0.3)+
geom_text(aes(label=signif(unlist(mean),2), y=unlist(Q0.975)+5000, group=Scenario), position=dodge)
ggsave("Ruori Burden of disease.png",width=10,height=6)
# ggsave("Ruori-tautitaakka.png",width=10,height=6)
# Error bars are not used in this plotly because it is unclear what it means in a stacked bar.
p <-plot_ly(
summary(oapply(BoDattr[BoDattr$Scenario=="BAU",], NULL, sum,c("Age"))),
y=~Exposure_agent,
x=~signif(mean,3),
text=~Response,
name=~Response,
type="bar",
orientation = "h"
) %>%
layout(
barmode="stack",
title="Elintarvikeperäisiä tautitaakkoja Suomessa",
xaxis=list(title="Haittapainotettua elinvuotta vuodessa (DALY/a)"),
yaxis=list(title="")
)
# pushIndicatorGraph(p, 117)
# ggsave("Ruori-tautitaakat.png",width=10/1.2,height=6/1.2)
# Utility of actions
thlBarPlot(summary(utility),xvar=Exposure_agent, yvar=unlist(mean), ylimits=c(-9000,0),
title="Ruori-skenaarioiden vaikutus tautitaakkaan",
subtitle="Haittapainotettua elinvuotta vuodessa (DALY/a)")+coord_flip()+
geom_errorbar(aes(ymin=unlist(Q0.025),ymax=unlist(Q0.975)), width=0.3)+
geom_text(aes(label=signif(unlist(mean),2), y=unlist(Q0.025)-600))
# ggsave("Ruori-toimenpideiden vaikutus.png",width=10/1.2,height=6/1.2)
}
| |
