Lähdekoodi sivulle Ruori

[[Luokka:Projekti]]
{{arviointi|moderaattori=Jouni}}
'''Ruori''' on VN-TEAS-hanke, jossa arvioidaan erilaisia ruokaan liittyviä riskitekijöitä, niiden vähentämispotentiaalia ja niiden terveys- ja talousvaikutuksia.

== Rajaus ==

=== Kysymys ===

Millaista tautitaakkaa Suomessa aiheuttavat Ruori-altisteet (tyydyttynyt rasva, vähäiset vihannekset, vähäiset hedelmät, liiallnen suola, dioksiinit, lyijy, toksoplasma, norovirus ja legionella?

=== Aikataulu ja käyttäjät ===

* Hanke alkoi 2018 ja loppuu 30.6.2019.
* Toteuttajina ovat Ruokavirasto, THL ja Helsingin yliopisto.
* Seuraavat skenaariot ovat tarkastelussa:
** 

== Vastaus ==

Vastaus kirjoitetaan tähän pian.

== Perustelut ==

=== Data ===

Mitä kaikkea kuuluu vähäiseen hedelmien tai vihannesten syöntiin?
* Vähähedelmäinen ruokavalio: hedelmien kulutus alle 3 annosta päivässä (310 g yhteensä) (sisältää tuoreet, pakastetut, keitetyt, säilötyt ja kuivatut hedelmät mutta ei sisällä hedelmämehuja tai suolaan tai etikkaan säilöttyjä hedelmiä) http://www.healthdata.org/terms-defined. Diet low in fruits: Consumption of less than 3 servings (310 g total) of fruits per day (includes fresh, frozen, cooked, canned, or dried fruit but excludes fruit juices and salted or pickled fruits). 
* Vähävihanneksinen ruokavalio: vihannesten kulutus alle 4 annosta  (400 g yhteensä) (sisältää tuoreet, pakastetut, keitetyt, säilötyt ja kuivatut vihannekset mukaan lukien palkokasvit mutta ei sisällä suolaan tai etikkaan säilöttyjä vihanneksia eikä pähkinöitä, siemeniä tai tärkkelyspitoisia vihanneksia kuten perunaa tai maissia). Diet low in vegetables: Consumption of less than 4 servings (400 g total) of vegetables per day (includes fresh, frozen, cooked, canned, or dried vegetables including legumes but excluding salted or pickled, juices, nuts and seeds, and starchy vegetables such as potatoes or corn).

Luken tilastoista löytyy tietoja kalansyönnistä Suomessa. Järvikalaa ei ole eritelty, mutta muut kuin viljellyt ja merilajit ovat yhteensä 2.6 kg/a henkeä kohti. https://stat.luke.fi/en/fish-consumption-2017_en

Voiko DALYt muuntaa euroiksi, ja miten se tehdään?
* Drake ehdottaa globaalin arvon päättämistä DALYn hinnaksi, samaan tapaan kuin 1,25 dollarin alittava päivätulo on määritelty absoluuttiseksi köyhyydeksi. Tällöin kaikki tuota hintaa kustannustehokkaammat toimet kannattaisi tehdä joko kansallisin, tai jos se ei jostain syystä onnistu, kansainvälisin toimin. Hän ei kuitenkaan ehdota suuruutta tälle hinnalle.<ref>Drake T. (2014) Priority setting in global health: towards a minimum DALY value. Health Economics Letter 23:2:248-252. https://doi.org/10.1002/hec.2925</ref>
* Brent on analysoinut implisiittisiä hintoja DALYlle Global Fund for AIDS, Tuberculosis, and Malaria -säätiön rahoituspäätöksistä. DALYn hinta näyttää olevan 6300 USD kaikille taudeille keskimäärin, ja 11900 USD erityisesti HIV/AIDSille<ref>Brent RJ. (2011) An implicit price of a DALY for use in a cost-benefit analysis of ARVs. Applied Economics 43:11:1413-1421. https://doi.org/10.1080/00036840802600475</ref>. Globaalit luvut ovat toki pienemmät kuin mitä rikkaissa länsimaissa katsottaisiin aiheelliseksi käyttää.
* Erilaisista arvioinneista löytyy vaihtelevia lukuja yhden DALYn rahalliselle hinnalle. Esimerkiksi IOMin Shecan-projekti käytti arvoa 50393 €/menetetty elinvuosi<ref>Minstry R. (2011) Methodology for valuing health impacts on the SHEcan project. IOM Research Project P937/96. http://ec.europa.eu/social/BlobServlet?docId=10178&langId=en</ref>, ja IGCB(N)-meluarviointiryhmä käyttää arvoa 60000 GBP/QALY (laatupainotettu elinvuosi) mutta samalla toteaa, että eri arvioinneissa arvot voivat vaihdella välillä 29000 - 130000 GBP/QALY<ref>The Interdepartmental Group on Costs and Benefits Noise Subject Group. (2014) Environmental Noise: Valuing impacts on sleep disturbance, annoyance, hypertension, productivity and quiet. https://assets.publishing.service.gov.uk/government/uploads/system/uploads/attachment_data/file/380852/environmental-noise-valuing-imapcts-PB14227.pdf</ref><ref>The Interdepartmental Group on Costs and Benefits Noise Subject Group (IGCB(N)). (2010) Noise & Health – Valuing the Human Health Impacts of Environmental Noise Exposure. https://khub.net/c/document_library/get_file?uuid=6a229977-e27a-43c5-a780-e224649bd2df&groupId=6197021</ref>.
* Berryn ja Flindellin mukaan Isossa-Britanniassa käytäntö on muodostunut sellaiseksi, että lääkkeet tai muut lääketieteelliset toimenpiteet saavat kansallisessa terveysjärjestelmässä helposti puollon, jos ne tuottavat yhden terveen elinvuoden alle 20000 GBP:n kustannuksilla. Tyypillisesti toteutetaan hankkeita tasolla 30000 GBP/QALY, mutta hankkeilta hinnaltaan yli 50000 GBP/QALY vaaditaan erityisiä perusteluja<ref>Berry BF, Flindell IH. (2009) Estimating  Dose-Response  Relationships between Noise Exposure and Human Health Impacts in the UK. BEL Technical Report 2009-002. https://webarchive.nationalarchives.gov.uk/20130123222353/http://archive.defra.gov.uk/environment/quality/noise/igcb/documents/tech-report.pdf</ref>.
* Toisaalta Hammitt todistelee, että hinta per tilastollinen elämä (value per statistical life, VSL) ja hinta per DALY muuttuvat epälineaarisesti suhteessa toisiinsa, eikä näin ollen olisi mahdollista käyttää hyvinvointimuutoksen mittarina vakiolla kerrottua DALY-arvoa, ainakaan taloudellisen hyvinvointiteorian (economic welfare theory) puitteissa.<ref>Hammitt, J.K. (2013) Admissible utility functions for health, longevity, and wealth: integrating monetary and life-year measures. J Risk Uncertain 47: 311. https://doi.org/10.1007/s11166-013-9178-4</ref>

Pitoisuusanalyysien kustannukset (poistettu).

Trikiinin valvontakustannukset: Tämän artikkelin mukaan trikiinin DALYt ovat vain luokkaa 100 DALY/miljardi ihmistä, joten valvonta ei ole mielekästä<ref>Brecht Devleesschauwer, Nicolas Praet, Niko Speybroeck, Paul R. Torgerson, Juanita A. Haagsma, KrisDe Smet, K. Darwin Murrell, Edoardo Pozio, Pierre Dorny. (2015) The low global burden of trichinellosis: evidence and implications. International Journal for Parasitology 45, 2–3, 95-99. [https://doi.org/10.1016/j.ijpara.2014.05.006] [https://www.sciencedirect.com/science/article/pii/S0020751914001374]</ref>.

==== Malliparametrit ====

<t2b name="Malliparametrit" index="Response,Exposure_agent,Type,Subgroup,Unit" obs="Result" desc="Description" unit="various">
CHD death||BoD|Age:Female 25-69|DALY|9876 (9103 - 10784)|Z:\Projects\RUORI\tautitaakka\Rasvat\IHD_data_IHME.csv
CHD death||BoD|Age:Female 70+|DALY|42750 (41007 - 45909)|Z:\Projects\RUORI\tautitaakka\Rasvat\IHD_data_IHME.csv
CHD death||BoD|Age:Male 25-69|DALY|48851 (54035 - 46123)|Z:\Projects\RUORI\tautitaakka\Rasvat\IHD_data_IHME.csv
CHD death||BoD|Age:Male 70+|DALY|48150 (46327 - 51255)|Z:\Projects\RUORI\tautitaakka\Rasvat\IHD_data_IHME.csv
Diet high in sodium||BoD|Age:Total population|DALY|27670 (1310 - 66420)|IHME GBD2017
Diet low in fruits||BoD|Age:Total population|DALY|36050 (20570 - 54800)|IHME GBD2017 http://ghdx.healthdata.org/gbd-results-tool?params=gbd-api-2017-permalink/da3bde44e863adb438c5fb47a89942fb
Diet low in vegetables||BoD|Age:Total population|DALY|28440 (14190 - 45960)|IHME GBD2017
Liver cancer||BoD|Age:Age 25-64|DALY|2499 (2018 - 3115)|From IHME (2014)
Liver cancer||BoD|Age:Age 65-74|DALY|2745 (2245 - 3311)|From IHME (2014)
dummy||case burden|Age:Female 25-69|DALY /case|0|Needed for case_burden to cover all Ages
dummy||case burden|Age:Female 70+|DALY /case|0|Needed for case_burden to cover all Ages
dummy||case burden|Age:Male 25-69|DALY /case|0|Needed for case_burden to cover all Ages
dummy||case burden|Age:Male 70+|DALY /case|0|Needed for case_burden to cover all Ages
dummy||case burden|Age:Age 25-64|DALY /case|0|Needed for case_burden to cover all Ages
dummy||case burden|Age:Age 65-74|DALY /case|0|Needed for case_burden to cover all Ages
IQ loss||case burden|Age:Age 1|DALY /IQ|0.11 (0.06 - 0.16)|Arja used 0.013 but here we use Goherr value instead
Listeria||case burden|Age:Total population|DALY/case|10 (5 - 13.3)|WHO 2015 report, European values (Table A8.2)
Noro virus||case burden|Age:Total population|DALY/case|0.0015 - 0.0025|WHO 2015 report, European values (Table A8.2)
Toxoplasma gondii||case burden|Age:Age 0 (congenital)|DALY/case|7|WHO 2015 report, European values (Table A8.2)
Toxoplasma gondii||case burden|Age:Age 1+ (acquired)|DALY/case|0.05|WHO 2015 report, European values (Table A8.2)
|Aflatoxin|exposure|Age:Age 25-64|ng /kg /d|0.85 - 1.14|Finravinto 2012
|Aflatoxin|exposure|Age:Age 65-74|ng /kg /d|0.5 - 0.67|Finravinto 2012
|Lead|exposure|Age:Age 1|ug /l|27.9|Measured as blood concentration. RASKURI, Z:\Projects\RUORI\tautitaakka\Lyijy\Lyijy_tautitaakkadata.xlsx
|Saturated fat|exposure|Age:Age 25-69|E%|13.1 (12.9 - 13.4)|Finland, 2010 situation from Wang et al. Supplementary
|Saturated fat|exposure|Age:Age 70+|E%|13.2 (12.8 - 13.6)|Finland, 2010 situation from Wang et al. Supplementary
|Aflatoxin|frexposed|Age:Age 25-64|fraction|1|dummy variable
|Aflatoxin|frexposed|Age:Age 65-74|fraction|1|dummy variable
|Lead|frexposed|Age:Age 1|fraction|0.0657|Population exposed to lead over threshold: 3126. RASKURI, Z:\Projects\RUORI\tautitaakka\Lyijy\Lyijy_tautitaakkadata.xlsx
IQ loss||incidence|Age:Age 1|IQ /100000py|596000|On average, a population has ca. 6 IQ points per person below 100: mean(abs(rnorm(10000, 100,15)-100))/2
Listeria||incidence|Age:Total population|# /100000py|1.22|Tartuntatautirekisteri 2016: 66 kpl. WHO 2015 report, European values (Table A8.2): 0.2 (0.2 - 0.3)
Liver cancer||incidence|Age:Age 25-64|# /100000py|4.06|Finnish Cancer Registry, average 2011-2015
Liver cancer||incidence|Age:Age 65-74|# /100000py|26.16|Finnish Cancer Registry, average 2011-2015
Noro virus||incidence|Age:Total population|# /100000py|1652 (630 - 3294)|WHO 2015 report, European values (Table A8.2)
Toxoplasma gondii||incidence|Age:Age 0 (congenital)|# /100000py|0.3 (0.2 - 0.7)|WHO 2015 report, European values (Table A8.2)
Toxoplasma gondii||incidence|Age:Age 1+ (acquired)|# /100000py|119 (77 - 188)|WHO 2015 report, European values (Table A8.2)
Diet low in fruits|Fruits|PAF|Age:Total population|fraction|1|dummy variable
Listeria|Listeria|PAF|Age:Total population|fraction|1|dummy variable
Noro virus|Noro virus|PAF|Age:Total population|fraction|1|dummy variable
CHD death|Saturated fat|PAF|Age:Age 25-69|fraction|0.064 (0.050 - 0.078)|Finland, 2010 situation from Wang et al. Supplementary
CHD death|Saturated fat|PAF|Age:Age 70+|fraction|0.048 (0.033 - 0.063)|Finland, 2010 situation from Wang et al. Supplementary
Diet high in sodium|Sodium|PAF|Age:Total population|fraction|1|dummy variable
Toxoplasma gondii|Toxoplasma gondii|PAF|Age:Total population|fraction|1|dummy variable
Diet low in vegetables|Vegetables|PAF|Age:Total population|fraction|1|dummy variable
||population|Age:Age 1||47577|Statistics Finland, 2018 http://stat.fi/til/synt/index.html
||population|Age:Age 25-64|#|2835089|Statistics Finland?
||population|Age:Age 65-74|#|615487|Statistics Finland?
||population|Age:Total population|#|5471753|Statistics Finland?
||population|Age:Age 0 (congenital)|#|47577|Statistics Finland?
||population|Age:Age 1+ (acquired)|#|5471753|Statistics Finland?
Hepatitis||prevalence|Hepatitis:Hepatitis B-|fraction|0.005|TerveSuomi
Hepatitis||prevalence|Hepatitis:Hepatitis B+|fraction|0.995|TerveSuomi
|Saturated fat|scenario exposure|Age:Male 18-24|E%|14.9|Finravinto 2017. Supplementary table 7.12. Average daily intake of saturated fats by gender and age.
|Saturated fat|scenario exposure|Age:Male 25-44|E%|15.2|Finravinto 2017. Supplementary table 7.12. Average daily intake of saturated fats by gender and age.
|Saturated fat|scenario exposure|Age:Male 45-64|E%|15.3|Finravinto 2017. Supplementary table 7.12. Average daily intake of saturated fats by gender and age.
|Saturated fat|scenario exposure|Age:Male 65-74|E%|14.4|Finravinto 2017. Supplementary table 7.12. Average daily intake of saturated fats by gender and age.
|Saturated fat|scenario exposure|Age:Female 18-24|E%|13.8|Finravinto 2017. Supplementary table 7.12. Average daily intake of saturated fats by gender and age.
|Saturated fat|scenario exposure|Age:Female 25-44|E%|14.8|Finravinto 2017. Supplementary table 7.12. Average daily intake of saturated fats by gender and age.
|Saturated fat|scenario exposure|Age:Female 45-64|E%|14.3|Finravinto 2017. Supplementary table 7.12. Average daily intake of saturated fats by gender and age.
|Saturated fat|scenario exposure|Age:Female 65-74|E%|14.0|Finravinto 2017. Supplementary table 7.12. Average daily intake of saturated fats by gender and age.
IQ loss|Lead|threshold|Age:Age 1|ug /l|24|Lanphear et al 2005 https://doi.org/10.1289/ehp.7688 CHECK THRESHOLD
Liver cancer|Aflatoxin|UR|Hepatitis:Hepatitis B-|# /(ng /kg /d /100000py)|0.01 (0.002 - 0.03)|WHO Is this per year or per lifetime?
Liver cancer|Aflatoxin|UR|Hepatitis:Hepatitis B+|# /(ng /kg /d /100000py)|0.3 (0.01 - 0.5)|WHO Is this per year or per lifetime?
IQ loss|Lead|UR|Age:Age 1|IQ l /ug|0.039|Lanphear et al 2005 https://doi.org/10.1289/ehp.7688 CHECK THRESHOLD
</t2b>

Annos-vasteet on tässä vain näytillä, ja oikeat käyttöön tulevat luvut ovat sivulla [[:op_en:ERFs of environmental pollutants]].

<t2b name="Decisions" index="Decision,Option,Variable,Cell,Change" obs="Result" desc="Description" unit="-">
Adjust|BAU|incidence||Multiply|0.00001|1/100000 py --> 1 py
Adjust|BAU|PAF||Identity|1|Scaled from Wang to Finravinto 2017: mean(15.2,15.3,14.8,14.0)/13.1
Adjust|BAU|PAF|Exposure_agent:Saturated fat|Multiply|1.132|Healthy people. Data from TerveSuomi
Hepatitis|Hepatitis B-|BoD||Multiply|0.995|Hepatitis B patients. Data from TerveSuomi
Hepatitis|Hepatitis B+|BoD||Multiply|0.005|1/100000py --> 1/py
Adjust|BAU|ERF||Identity|1|For completion
Adjust|BAU|ERF|Response:Liver cancer|Multiply|0.00001|For completion
</t2b>

<t2b name="CollapseMarginals" index="Variable,Index,Probs,Function" obs="Dummy" desc="Description" unit="-">
BoD|incidenceSource,disabilityweightSource,populationSource,BoDSource||sum|1|Remove redundant
PAF|Unit, Exposure, Scaling,Exposcen, ER_function, ERFchoiceSource, exposureSource, bgexposureSource, BWSource, doseSource, thresholdSource, ERFSource, RRSource, frexposedSource, incidenceSource, InpPAFSource||sum|1|Remove redundant
case_burden|case_burdenSource||sum|1|Fill missing Ages
</t2b>

=== Laskenta ===

==== Arviointimalli ====

* Malliajo 29.5.2019 [http://fi.opasnet.org/fi-opwiki/index.php?title=Toiminnot:RTools&id=dbscMpzNoaiytUi3]
* Malliajo 31.5.2019 [http://fi.opasnet.org/fi-opwiki/index.php?title=Toiminnot:RTools&id=JLSryIWqMiZEMhBh]

<rcode graphics=1>
library(OpasnetUtils)
library(ggplot2)

openv.setN(0)

DecisionTableParser(opbase.data("Op_fi5889",subset="Decisions"))
dat <- opbase.data("Op_fi5889", subset="Malliparametrit")[-1]

CTable <- opbase.data("Op_fi5889",subset="CollapseMarginals")
for(i in 1:ncol(CTable)) {CTable[[i]] <- as.character(CTable[[i]])}
CollapseTableParser(CTable)

ColBoD2$probs<-NULL # Fix CollapseTableParser! NA and "" in probs should work too.

#' create_subgroups adjusts the data table for ovariables
#' @param dat data.frame
#' @param type type of data that is used. Must match content in column Type
#' @param drop columns to remove
#' @return data.frame

#' Only works for one subgroup per cell. Should be generalised to match several subgroups separated by ;
#' @param data.frame that contains columns Subgroup and - if type is used - Type
#' @param type Selects given type from column Type
#' @param drop column names to drop

prepare <- function(dat, type=NULL, drop=NULL) {
  out <- dat
  if(!is.null(type)) out <- out[out$Type %in% type , ]
  if(!is.null(drop)) out <- out[!colnames(out) %in% drop]
  out$Subgroup <- as.character(out$Subgroup)
  tmp <- ifelse(out$Subgroup=="" | is.na(out$Subgroup), " : ", out$Subgroup)
  tmp <- lapply(strsplit(tmp,";"), function(x) {strsplit(x, ":")[[1]]})
  tmp <- data.frame(
    Sg_type = sapply(tmp, function(x) {x[[1]]}),
    Subgroup = sapply(tmp, function(x) {x[[2]]}),
    stringsAsFactors = FALSE
  )
  for(i in setdiff(unique(tmp$Sg_type)," ")) {
    out[[i]] <- ifelse(tmp$Sg_type==i, tmp$Subgroup,NA)
  }
  out$Subgroup <- NULL
  return(out)
}

ERFchoice <- Ovariable(
  "ERFchoice",
  dependencies = data.frame(Name="exposure"),
  formula = function(...) {
    out <- cbind(unique(exposure@output["Exposure_agent"]),Result=1)
    return(out)
  }
)

dec <- opbase.data("Op_fi5889", subset="Decisions")[-1]

population <- Ovariable("population", data = prepare(dat,"population",c("Type","Exposure_agent","Response")))

exposure <- Ovariable("exposure", data = prepare(dat, "exposure", c("Type","Response")))

frexposed <- Ovariable("frexposed", data=prepare(dat, "frexposed", c("Type","Response","Unit")))

incidence <- Ovariable("incidence", data = prepare(dat,"incidence",c("Type","Exposure_agent","Unit")))

disabilityweight <- Ovariable(
  "disabilityweight",
  ddata = "Op_en7748",
  subset = "DALYs of responses"
)
colnames(disabilityweight@data)[match(
  c("Resp","disabilityweightResult"),
  colnames(disabilityweight@data))] <- c("Response","Result")

disabilityweight@data <- orbind(
  disabilityweight@data,
  prepare(dat,"disability weight * duration",c("Type","Exposure_agent","Unit"))
)

PAF <- Ovariable(
  "PAF",
  data =  prepare(dat,"PAF","Type"),
  dependencies=data.frame(Name=c("PAF1","PAF2")),
  formula = function(...) {
    return(unkeep(PAF1,sources=TRUE,prevresults = TRUE))#OpasnetUtils::combine(PAF1, PAF2))# ,c("Subgroup","Adjust")
  }
)

duration <- 1

#objects.latest("Op_en7422", code_name="BoD") # [[Burden of disease]]

BoD2 <- Ovariable(
  "BoD2",
  dependencies=data.frame(Name=c("incidence","disabilityweight","population","duration")),
  formula=function(...) {
    out <- incidence * population * duration * disabilityweight
    return(out)
  }
)

BoDattr <- Ovariable(
  "BoDattr",
  dependencies = data.frame(Name=c("BoD2","PAF")),
  formula = function(...) {
    out <- BoD2 * PAF
    return(out)
  }
)

##########################################

# This is code Op_en2261/PAF1 on page [[Health impact assessment]].
library(OpasnetUtils)

### NOTE! The function does not work if the resultnames are not doseResult, ERFResult and thresholdResult.
PAF1 <- Ovariable(
  "PAF1", # This calculates the burden of disease for background-independent endpoints.
  dependencies = data.frame(
    Name = c(
      #      "population", # Population divided into subgroups as necessary
      "dose", # Exposure to the pollutants
      "ERF", # Other ERFs than those that are relative to background.
      "threshold", # exposure level below which the agent has no impact.
      "incidence", # background incidence of the disease. It cancels out if response not dependent on this
      "frexposed", # fraction of population that is exposed
      "sumExposcen", # function that calculates difference between exposure scenarios
      "mc2d",  # Function to run two-dimensional Monte Carlo
      "mc2dparam" # Parameters for mc2d function. Default: don't run.
    ), 
    Ident = c(
      #      "Op_en2261/population", # [[Health impact assessment]]
      "Op_en2261/dose",       # [[Health impact assessment]]
      "Op_en2031/initiate",   # [[Exposure-response function]]
      "Op_en2031/initiate",   # [[Exposure-response function]]
      "Op_en2261/incidence",  # [[Health impact assessment]]
      "Op_en2261/frexposed",  # [[Health impact assessment]]
      "Op_en2261/sumExposcen",# [[Health impact assessment]]
      "Op_en7805/mc2d",       # [[Two-dimensional Monte Carlo]]
      "Op_en7805/mc2dparam"   # [[Two-dimensional Monte Carlo]]
    )
  ),
  formula = function(...) {
    out <- Ovariable()
    exps <- unique(dose$Exposure_agent)
    
    tmp <- ERF[ERF$ER_function %in% c("UR", "CSF", "ERS") , ]
    if(!is.null(exps)) tmp <- tmp[tmp$Exposure_agent %in% exps , ]
    
    if(nrow(tmp@output)>0) {
      tmp2 <- dose - threshold
      result(tmp2) <- pmax(0,result(tmp2))
      out <- tmp2 * tmp * frexposed / incidence # Actual equation
    }
    
    # Step estimates: value is 1 below threshold and above ERF, and 0 in between.
    
    tmp <- ERF[ERF$ER_function %in% c("Step", "ADI", "TDI", "RDI", "NOAEL") , ]
    if(!is.null(exps)) tmp <- tmp[tmp$Exposure_agent %in% exps , ]
    if(nrow(tmp@output)>0) {
      tmp <- 1 - (dose >= threshold) * (dose <= tmp)
      tmp <- tmp * frexposed / incidence
      
      if(nrow(out@output)==0) out <- tmp else out <- OpasnetUtils::combine(out, tmp)
    }
    out <- oapply(out, NULL, sum, c("Exposure_agent", "Exposure", "ER_function", "Scaling"))
    out <- mc2d(out) # Run two-dimensional Monte Carlo if needed
    
    return(sumExposcen(out))
  }
)

#objects.store(PAF1)
cat("Ovariable casesabs stored. page: Op_en2261, code_name: casesabs.\n")

##################################################

# This is code Op_en2261/casesrr on page [[Health impact assessment]].
library(OpasnetUtils)

PAF2 <- Ovariable(
  "PAF2", # This calculates the burden of disease for endpoints using RR.
  dependencies = data.frame(
    Name = c(
      #      "population", # Population divided into subgroups as necessary
      "RR", # Relative risks for the given exposure
      "frexposed", # fraction of exposed population
      "sumExposcen" # function that calculates difference between exposure scenarios
    ), 
    Ident = c(
      #      "Op_en2261/population", # [[Health impact assessment]]
      "Op_en2261/RR",         # [[Health impact assessment]]
      "Op_en2261/frexposed",  # [[Health impact assessment]]
      "Op_en2261/sumExposcen" # [[Health impact assessment]]
    )
  ),
  formula = function(...) {
    out <- frexposed * (RR - 1) # RR diluted to get PAF instead of AF
    r <- result(out)
    result(out) <- (r > 0) * (r / (r + 1)) + (r <= 1) * r # Negative if RR < 1
    
    return(sumExposcen(out))
  }
)

#objects.store(PAF2)
cat("Ovariable PAF2 stored. page: Op_en2261, code_name: casesrr.\n")

###############################################

# This is code Op_en2031/initiate on page [[Exposure-response function]]

library(OpasnetUtils)

ERF <- Ovariable(
  "ERF",
  dependencies = data.frame(Name = c(
    "ERF_env", # [[ERF of environmental pollutants]]
    "ERF_omega3", # [[ERF_of_omega-3_fatty_acids]]
    "ERF_mehg", # [[ERF_of_methylmercury]]
    "ERF_diox", # [[ERF_of_dioxin]]
    "ERF_vit", # [[ERF_of_dioxin]]
    "ERFchoice" # an ovariable for choosing case-specific Responses, Exposure_units or ER_functions
  ), Ident = c(
    "Op_en5827/initiate",
    "Op_en5830/initiate",
    "Op_en5825/initiate",
    "Op_en5823/initiate",
    "Op_en6866/initiate",
    "Op_en2031/ERFchoice"
  )),
  formula = function(...) {
    
    out <- OpasnetUtils::combine(
      ERF_env, 
      ERF_omega3, 
      ERF_mehg, 
      ERF_diox,
      ERF_vit
    )
    #    out@output$Age[is.na(out@output$Age)] <- ""
    out@marginal[colnames(out@output) == "Exposure_unit"] <- FALSE
    if("Iter" %in% colnames(out@output)) {
      out@output <- fillna(out@output, marginals = "Iter")} # To make sure that deterministic ERFs work properly.
    out <- unkeep(out * ERFchoice, cols=c("ERFchoiceSource", "ERFchoiceResult")) #technical nuisance removed
    out@output <- prepare(out@output)
    
    return(out)                         
  }
)

threshold <- Ovariable(
  "threshold",
  dependencies = data.frame(Name = c(
    "threshold_env", # [[ERF of environmental pollutants]]
    "threshold_omega3", # [[ERF_of_omega-3_fatty_acids]]
    "threshold_mehg", # [[ERF_of_methylmercury]]
    "threshold_diox", # [[ERF_of_dioxin]]
    "threshold_vit" # [[ERFs of vitamins]]
  ), Ident = c(
    "Op_en5827/initiate",
    "Op_en5830/initiate",
    "Op_en5825/initiate",
    "Op_en5823/initiate",
    "Op_en6866/initiate"
  )),
  formula = function(...) {
    
    out <- OpasnetUtils::combine(
      threshold_env, 
      threshold_omega3, 
      threshold_mehg, 
      threshold_diox,
      threshold_vit
    )
    #    out@output$Age[is.na(out@output$Age)] <- ""
    out@output <- prepare(out@output)
    
    return(out)                         
  }
)

#objects.store(ERF, threshold)
cat("Ovariables ERF, threshold stored. Page: Op_en2031, code: initiate.\n")

##############################

# population <- EvalOutput(population, verbose=TRUE)

BoD2 <- EvalOutput(BoD2, forceEvaluation=TRUE, verbose=TRUE)

BoDattr <- EvalOutput(BoDattr, verbose=TRUE)

oprint(ls())

oprint(summary(exposure))
oprint(summary(dose))
oprint(summary(ERF))
oprint(summary(incidence))
oprint(summary(frexposed))
oprint(summary(PAF1))
oprint(summary(PAF2),digits=7)
oprint(summary(BoD2),digits=7)
oprint(summary(BoDattr),digits=7)

oprint(exposure@output)
oprint(dose@output)
oprint(ERF@output)
oprint(incidence@output)
oprint(frexposed@output)
oprint(PAF1@output)
oprint(PAF2@output)
oprint(PAF@output)
oprint(PAF@data)
oprint(BoD2@output,digits=7)
oprint(BoDattr@output,digits=7)

ggplot(ERF@output, aes(x=Response, weight=ERFResult))+geom_bar()+coord_flip()
ggplot(BoD2@output, aes(x=Response, weight=BoD2Result))+geom_bar()+coord_flip()
ggplot(BoDattr@output, aes(x=Response, weight=BoDattrResult))+geom_bar()+coord_flip()
</rcode>

==== Tautitaakkakuvia ====

<rcode label="Ruori ja tautitaakka -esityksen kuvat (aja omalla koneella)">
# Tämä on koodi Op_fi5889 sivulla [[Ruori]]

library(OpasnetUtils)
library(tidyverse)
library(treemap)
library(plotly)

#objects.latest("Op_en6007", code_name = "hnh2035") # [[OpasnetUtils/Drafts]] pushIndicatorGraph

##### Finnish translations
transl <- as_tibble(opbase.data("Op_fi3944", subset="Tautiluokittelu")) %>% # [[Tautitaakka Suomessa]]
  mutate(Id=as.integer(as.character(Id)))

BS <- 24

palet <- c(
  '#a6cee3','#1f78b4','#b2df8a','#33a02c','#fb9a99','#e31a1c','#fdbf6f',
  '#ff7f00','#cab2d6','#6a3d9a','#ffff99','#b15928', # First 12 colours from Colorbrewer Paired.
  'cyan2','cyan4','plum1','plum4', 'darkslategray4','darkslategray1','firebrick3'
)

###################### Graphs for the Tautitaakka auttaa hahmottamaan ... article

# DALYs by causes (risks not included)
dc <- 
  as_tibble(opasnet.csv("2/2d/IHME_Fin_Risks_by_Cause.zip", wiki="opasnet_en",
                        unzip="IHME-GBD_2017_DATA-8ce9adcf-1.csv",sep=",",dec=".",header=TRUE)) %>%
  left_join(transl[transl$Type=="Cause",], by=c("cause_id"="Id")) %>%
  mutate(cause_name = Result) %>% # transl combines some entries from cause_name to Name (in Finnish: Result)
  #  "Neglected tropical diseases and malaria",
  #  "HIV/AIDS and sexually transmitted infections",
  #  "Enteric infections" >> "Other infectious diseases"
  filter(measure_name=="DALYs (Disability-Adjusted Life Years)" & metric_name=="Number" & year == 2017) %>%
  group_by(cause_name) %>% 
  summarise(value=sum(val)) %>%
  mutate(cause_name = factor(cause_name, levels=cause_name[order(value)])) %>%
  mutate(valy = value/1000 + ifelse(value<250000, 20, -40)) # unit kDALY/a

plot_cause <- ggplot(dc2, aes(x=cause_name,weight=value/1000,fill=cause_name, label=round(value/1000)))+
  geom_bar(position="stack")+geom_text(aes(y=valy))+coord_flip()+
  guides(fill=FALSE)+
  scale_fill_manual(values = rev(palet))+ 
  labs(
    title="Tautitaakka Suomessa syittäin 2017",
    x="Tauti tai haitta",
    y="Tautitaakka (tuhatta DALYa vuodessa)")

plot_cause
# ggsave("Tautitaakka Suomessa 2017.svg", width=8, height=5) # Png conversion: 300 pixels/inch, font size 4 times larger

############## DALYs by causes and risks

dcr <- as_tibble(
  opasnet.csv("2/2d/IHME_Fin_Risks_by_Cause.zip", wiki="opasnet_en",
              unzip="IHME-GBD_2017_DATA-8c9ca17f-1.csv",sep=",",dec=".",header=TRUE)
) %>%
  left_join(transl[transl$Type=="Cause",], by=c("cause_id"="Id")) %>%
  filter(measure_name=="DALYs (Disability-Adjusted Life Years)" & metric_name=="Number" & year == 2017) %>%
  #  mutate(cause_name = factor(Result, levels=dc2$cause_name)) %>% # transl combines some entries from cause_name to Name (in Finnish: Result)
  left_join(transl[transl$Type=="Risk",], by=c("rei_id"="Id"))

dcr <- dcr %>%
  mutate(
    cause_name=factor(Result.x,levels=levels(dc$cause_name)),
    rei_name=factor(Result.y, levels = aggregate(dcr["val"],dcr["Result.y"],sum) %>% transmute(Result.y[order(val)]) %>% pull())) %>%
  group_by(rei_name, cause_name, Luokka.y) %>%
  summarise(val = sum(val))

#  mutate(rei_name = factor(rei_name, levels= dcr2$rei_name)) # transl combines some entries from cause_name to Name (in Finnish: Result)

# Strange discrepancy in two sources from IHME: especially in cardiovascular diseases. 
# Reason not known, use dcr anyway; because errors not found.
# ihme <- read.csv("C:/Users/jtue/AppData/Local/Temp/download(1).csv")
# ihme <- ihme[ihme$Measure=="DALYs",]
# sum(ihme$Value)
# [1] 911630.9
# sum(dcr$val)
# [1] 1205234
# View(
#   merge(
#     aggregate(ihme["Value"], by=ihme["Cause.of.death.or.injury"], FUN=sum),
#     aggregate(dcr["val"], by=dcr["cause_name"], FUN=sum),
#     by.x=c("Cause.of.death.or.injury"), by.y=c("cause_name")
#   )
# )


plot_risk <- ggplot(dcr, aes(x=rei_name, weight=val/1000, fill=cause_name))+geom_bar()+coord_flip()+
  scale_fill_manual(values = rev(palet)[-8])+ # Skin diseases missing, so remove that colour 
  guides(fill = guide_legend(reverse=TRUE, title=NULL, keyheight=0.8))+
  theme_gray(base_size=11)+ theme(legend.position=c(0.83,0.29)) + 
  labs(
    title = "Tautitaakka Suomessa 2017 tunnettujen riskitekijöiden mukaan",
    x="Riskitekijä", y="Tautitaakka (tuhatta DALYa vuodessa)")

plot_risk
# ggsave("Tautitaakka Suomessa 2017 riskitekijöittäin.svg", width=8, height=6) # Png conversion: 300 pixels/inch, font size 4 times larger

ruori <-   data.frame(
  Luokka.y = c(rep("Ympäristö",3),rep("Ravitsemus",4),rep("Ympäristö",3)),
  rei_name = c(
    "Aflatoksiini",
    "Lyijy",
    "Dioksiini",
    "Tyydyttyneet rasvat",
    "Suola",
    "Hedelmien liian vähäinen saanti",
    "Kasvisten liian vähäinen saanti",
    "Norovirus",
    "Listeria",
    "Toksoplasma"
  ),
  val = c(
    8,
    6,
    22,
    8190,
    27062,
    35314,
    27778,
    1,
    660,
    432
  ),
  lower = c(
    6,
    NA,
    NA,
    5857,
    1330,
    20197,
    13740,
    NA,
    330,
    270
  ),
  upper = c(
    13,
    NA,
    NA,
    11125,
    65542,
    53635,
    44716,
    NA,
    880,
    702
  )
)

dcr$rei_name <- as.character(dcr$rei_name)

# Subtract Ruori results from IHME background
dcr$val[dcr$rei_name=="Ravitsemusriskit" & dcr$cause_name=="Sydän- ja verisuonitaudit"] <-
  dcr$val[dcr$rei_name=="Ravitsemusriskit" & dcr$cause_name=="Sydän- ja verisuonitaudit"] - sum(ruori$val[ruori$Luokka.y=="Ravitsemus"])

dcr$val[dcr$rei_name=="Muut ympäristöriskit" & dcr$cause_name=="Syöpä"] <-
  dcr$val[dcr$rei_name=="Muut ympäristöriskit" & dcr$cause_name=="Syöpä"] - sum(ruori$val[ruori$Luokka.y=="Ympäristö"])

dcr <- rbind(dcr[colnames(dcr)!="cause_name"],ruori[!colnames(ruori) %in% c("upper","lower")])
dcr$rei_name[dcr$rei_name=="Ravitsemusriskit"] <- "Muut ravitsemusriskit"

pdf("Tautitaakka riskitekijöittäin.pdf", width=11.67, height=4.27)
treemap(dcr, index=c("Luokka.y","rei_name"),vSize="val",
        title=paste0("Tunnettujen riskitekijöiden tautitaakka Suomessa 2017 (",
                     round(sum(dcr$val),-3)," DALY)"),
        align.labels=list(c("center","top"),c("center","center")),
        aspRatio = 2, border.col = "gray")
dev.off()

tmp <- ruori %>%
  filter(rei_name %in% c("Aflatoksiini","Lyijy","Dioksiini","Norovirus","Listeria","Toksoplasma")) %>%
  rename(Altiste=rei_name) %>%
  mutate(Altiste = factor(Altiste, levels=ruori$rei_name[order(ruori$val)]))

thlBarPlot(tmp, xvar=Altiste, yvar=val)+
  coord_flip()+geom_errorbar(data=tmp, aes(ymin=lower, ymax=upper), width=0.25, size=0.7)+
  labs(
    title="Ruorissa tutkittujen riskien tautitaakka",
    subtitle = "Tautitaakka mikrobiologisille ja kemiallisille epäpuhtauksille (DALYa vuodessa)")

ggsave("Mikrobiologiset ja kemialliset riskit Ruorissa.svg", width=8, height=5)

tmp <- ruori %>%
  filter(!rei_name %in% c("Aflatoksiini","Lyijy","Dioksiini","Norovirus","Listeria","Toksoplasma")) %>%
  rename(Altiste=rei_name) %>%
  mutate(Altiste = factor(Altiste, levels=ruori$rei_name[order(ruori$val)]))

thlBarPlot(tmp, xvar=Altiste, yvar=val)+
  coord_flip()+geom_errorbar(data=tmp, aes(ymin=lower, ymax=upper), width=0.25, size=0.7)+
  labs(
    title="Ruorissa tutkittujen riskien tautitaakka",
    subtitle = "Tautitaakka ravitsemusriskeille (DALYa vuodessa)")

ggsave("Ravitsemusriskit Ruorissa.svg", width=8, height=5)

####################################

# Katso sivut [[:op_en:Goherr assessment]] ja [[:op_en:Health effects of Baltic herring and salmon: a benefit-risk assessment]]
objects.get('155401096341')

groups <- function(o) {
  o$Group <- paste(o$Gender, o$Ages)
  o$Group <- factor(o$Group, levels = c(
    "Female 18-45",
    "Male 18-45",
    "Female >45",
    "Male >45"
  ))
  return(o)
}

varit12 <- c(
  '#2f62ad', # THL-värit, koko väripaletti. Nro 2
  '#7cd0d8', # 4
  '#571259', # 9
  '#5faf2c', # 11
  '#bf4073', # 14
  '#3b007f', # 18
  '#16994a', # 20
  '#cccc72', # 22
  '#0e1e47', # 3
  '#25a5a2', # 5
  '#cc7acc', # 7
  '#244911' # 12
)

#### Burden of disease

# Figure 5. from Goherr assessment manuscript

if(FALSE) tmp <- groups(BoD * info) else tmp <- groups(BoDRaw * info)

tmp <- oapply(
  tmp[!grepl("TWI", tmp$Resp) & tmp$Background=="Yes",],
  INDEX=c("Resp","Group","Country","Cons.policy","Background"),
  mean
)@output

# levels(tmp$Group)
# [1] "Female 18-45" "Male 18-45"   "Female >45"   "Male >45"    
levels(tmp$Group) <- c("Naiset 18-45","Miehet 18-45","Naiset >45","Miehet >45")

#levels(tmp$Resp)
#[1] "Stroke"           "Heart (CHD)"      "Cancer"           "Child's IQ"       "Infertility"     
#[6] "Tooth defect"     "Dioxin TWI"       "TWI 2018"         "Vitamin D intake"
levels(tmp$Resp) <- c("Aivohalvaus","Sepelvaltimotauti","Syöpä","Lapsen ÄO","Hedelmättömyys",
                      "Hammasvaurio","","","D-vitamiini")

fig22 <- ggplot(tmp[tmp$Country=="FI" & tmp$Cons.policy=="BAU",],
                aes(x=Group, weight=Result, fill=Resp))+
  geom_bar()+facet_grid(Cons.policy ~ Country)+
  theme_grey(base_size=BS-8)+
  theme(
    legend.position = "bottom"
    #    axis.text.x = element_text(angle = -90)
  )+
  labs(
    title="Itämeren kalan aiheuttama tautitaakka Suomessa ryhmittäin",
    y = "Tautitaakka (DALY/a)",
    fill="",
    x=""
  )+
  scale_fill_manual(values=varit12[-(7:8)])

fig22
ggsave("Itämeren kalan tautitaakka Suomessa.svg", width=8, height=5)

######### Figure 6.

tmp <- groups(BoDRaw * info)

#> levels(tmp$Resp)
#[1] "Stroke"           "Heart (CHD)"      "Cancer"           "Child's IQ"       "Infertility"     
#[6] "Tooth defect"     "Dioxin TWI"       "TWI 2018"         "Vitamin D intake"
levels(tmp$Resp) <- c("Aivohalvaus","Sepelvaltimotauti","Syöpä","Lapsen ÄO","Hedelmättömyys",
                      "Hammasvaurio","TWI 2001","TWI 2018","D-vitamiini")

tmp@output <- rbind(
  cbind(tmp@output,Focusgroup = "Kaikki"),
  cbind(tmp@output[tmp$Group=="Female 18-45",],Focusgroup = "Nuoret naiset")
)

tmp <- oapply(
  tmp[tmp$Cons.policy=="BAU" & tmp$Country=="FI",],
  INDEX=c("Resp","Background","Focusgroup","Country"),
  mean
)@output

tmp <- rbind(
  cbind(
    Objective="Kokonaisterveys, oletus",
    tmp[tmp$Background=="Yes" & !grepl("TWI",tmp$Resp),]
  ),
  cbind(
    Objective="Kokonaisterveys, taustatta",
    tmp[tmp$Background=="No" & !grepl("TWI",tmp$Resp),]
  ),
  cbind(
    Objective=tmp$Resp,
    tmp
  )[grepl("TWI",tmp$Resp),]
)

sums <- aggregate(tmp["Result"], by=tmp[c("Objective","Focusgroup","Country")],sum)
fig22b <- ggplot(tmp, aes(x=Objective, weight=Result))+
  geom_bar(aes(fill=Resp))+facet_grid(Country ~ Focusgroup)+
  theme_grey(base_size=BS-8)+
  theme(
    legend.position = "bottom",
    axis.text.x = element_text(angle = -15, hjust=0.5)
  )+
  geom_text(data=sums, size=6, aes(label=round(Result,-2),y=pmax(200,Result+200)))+
  labs(
    title="Tautitaakka eri tavoitteiden näkökulmasta",
    y = "Tautitaakka (DALY/a)",
    fill="",
    x=""
  )+
  scale_fill_manual(values=varit12)
fig22b

ggsave("Itämeren kalan tautitaakka eri näkökulmista.svg", width=10, height=7)

############################ Health care costs

dat <- scrape.webtable(
  "https://en.wikipedia.org/wiki/List_of_countries_by_total_health_expenditure_per_capita",
  6) 
colnames(dat) <- c("Country","Y2000","Y2005","Y2010","Y2015")
dat$Y2015 <- as.numeric(gsub(",","",as.character(dat$Y2015)))
dat$Country <- as.character(dat$Country)
dat$Country[match(c(
  "Iran (Islamic Republic of)",
  "Lao People's Democratic Republic",
  "United Republic of Tanzania",
  "United States of America",
  "Venezuela (Bolivarian Republic of)",
  "Viet Nam"
),dat$Country)] <- c("Iran","Laos","Tanzania","United States","Venezuela","Vietnam")

dat2 <- read.csv("C:/Users/jtue/AppData/Local/Temp/IHME-GBD_2017_DATA-96ba4b83-1.csv")
levels(dat2$location_name)[match(c(
  "The Bahamas","The Gambia"),levels(dat2$location_name))] <- c("Bahamas","Gambia")

dat3 <- merge(dat,dat2, by.x="Country",by.y="location_name", all=TRUE)
dat3$label <- ifelse(dat3$Country %in% c("Finland","Sweden","United States","Poland"), dat3$Country, " ")

tmp <- dat3[dat3$metric_name=="Rate" & dat3$cause_name=="All causes",]
tmp$label[match(c("Finland","Sweden","United States","Poland"),tmp$label)] <- c("Suomi","Ruotsi", "USA","Puola")

ggplot(tmp, aes(x=Y2015,y=val/100000, label=label, colour=label))+
  geom_point()+geom_text()+
  labs(
    title="Tautitaakka terveysmenojen funktiona maittain 2015",
    y="Tautitaakka per henkilö (DALY)",
    x="Terveydenhuoltomenot per henkilö (USD)"
  )

ggsave("Tautitaakka terveysmenojen funktiona maittain.svg", width=9, height=5)

pl <- plot_ly(tmp, x=~Y2015, y=~val/100000,type="scatter", mode="markers", hovertext=~Country) %>% #, text=~label,color=~label) %>%
  layout(
    title="Tautitaakka terveysmenojen funktiona maittain 2015",
    xaxis=list(title="Terveydenhuoltomenot per henkilö (USD)"),
    yaxis=list(title="Tautitaakka per henkilö (DALY)")
  )

pl
# pushIndicatorGraph(pl, 103)#, API_KEY= apikey)
</rcode>

== Viitteet ==

<references/>

== Katso myös ==

* {{#l:Ruori-hankkeen tautitaakka-arviot.pptx}}