+ Näytä koodi- Piilota koodi
library(OpasnetBaseUtils)
# Input korjailua basesta, jos niin halutaan
if(!(i.city == "Custom")) {
city.default.inputs <- op_baseGetData("opasnet_base", i.city)
if(!is.na(city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.dose"])) {
Cl.dose <- city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.dose"]}
if(!is.na(city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.delay"])) {
Cl.delay <- city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.delay"]}
if(!is.na(city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.byprod.t"])) {
Cl.byprod.t <- city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.byprod.t"]}
if(!is.na(city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.pH"])) {
Cl.pH <- city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.pH"]}
if(!is.na(city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.temp"])) {
Cl.temp <- city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.temp"]}
if(!is.na(city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.Br"])) {
Cl.Br <- city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.Br"]}
if(!is.na(city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.DOC"])) {
Cl.DOC <- city.default.inputs$Result[city.default.inputs$Havainto == "arvo" & city.default.inputs$Parametri == "Cl.DOC"]}
}
# Vesiopas -malli
# Niputetaan raakaveden ja sen käsittelyn inputit siistiksi
i.raw.pat.conc.val <- list(i.raw.pat.conc.Kamp, i.raw.pat.conc.Ecoli, i.raw.pat.conc.Rota, i.raw.pat.conc.Noro, i.raw.pat.conc.Crpt,
i.raw.pat.conc.Giardia)
# Annetaan raakaveden konsentraatiolle joku arvo inputtien perusteella, 2 mahdollista metodia: estimaatti, kirjallisuus
## Luokitukset
raw.pat.conc.lit <- op_baseGetData("opasnet_base", "Op_fi2655")
##
Pathogen <- c("Kampylobakteeri","E.coli O157:H7","Rotavirus","Norovirus","Cryptosporidium","Giardia")
raw.pat.conc <- data.frame(Pathogen, Raw.pat.conc = NA)
for(i in 1:length(i.raw.pat.conc.val)) {
if(is.character(i.raw.pat.conc.val[[i]])) raw.pat.conc[i,2] <- raw.pat.conc.lit[raw.pat.conc.lit$Luokitus == i.raw.class &
raw.pat.conc.lit$Patogeeni == Pathogen[i], "Result"] else raw.pat.conc[i,2] <- i.raw.pat.conc.val[[i]]
}
# Rakennetaan objekti, joka kertoo mitkä puhdistusmenetelmät ovat käytössä ja miten tehokkaita ne ovat
treat.used <- data.frame(Treatment = c("Perinteinen puhdistus","Hyvin toimiva puhdistus","Tehostettu puhdistus","Hidas hiekkasuodatus",
"Kalkkikivisuodatus","Aktiivihiilisuodatus","UV"), Active = 0)
treat.used$Active[i.treat.used] <- 1
if(1 %in% i.treat.used) {if(2 %in% i.treat.used) {if(3 %in% i.treat.used){treat.used$Active[c(1,2)] <- 0} else {treat.used$Active[c(1,3)] <- 0}}
else {treat.used$Active[c(2,3)] <- 0}} else {treat.used$Active[c(1,2,3)] <- 0}
treat.pat.eff <- op_baseGetData("opasnet_base", "Op_fi2656")
treat.pat.eff <- treat.pat.eff[,c("Patogeeni", "Vedenpuhdistusmenetelmä","Result")]
colnames(treat.pat.eff)[1:2] <- c("Pathogen", "Treatment")
treat.pat.reduct <- merge(treat.pat.eff, treat.used)
treat.pat.reduct$active.reduct <- treat.pat.reduct$Result * treat.pat.reduct$Active
treat.pat.reduct.tot <- as.data.frame(as.table(tapply(treat.pat.reduct$active.reduct, treat.pat.reduct$Pathogen, sum))) # lasketaan käytössä olevat
# log vähenemät yhteen
colnames(treat.pat.reduct.tot) <- c("Pathogen", "Treat.pat.reduct")
# Desinfiointi, klooraus
Cl.test.sens <- 0
f.Cl.pat.red <- function(ct, sens) { # Free chlorine disinfection log reduct estimation function
out <- NA
if (sum(sens) == 0) out <- 0
if (ct < sens[1] & sens[1] > 0) {out <- ct * (1 / sens[1])}
if (ct >= sens[1] & sens[1] > 0) {out <- 1 + (ct - sens[1]) * (1 / (sens[2] - sens[1]))}
if (ct >= sens[2] & sens[2] > 0) {out <- 2 + (ct - sens[2]) * (1 / (sens[3] - sens[2]))}
if (ct >= sens[3] & sens[3] > 0 & sens[4] > 0) {out <- 3 + (ct - sens[3]) * (1 / (sens[4] - sens[3]))}
if (ct >= sens[3] & sens[3] > 0 & sens[4] == 0) {out <- 3 + (ct - sens[3]) * (1 / (sens[3] - sens[2]))}
if (ct >= sens[4] & sens[4] > 0 & sens[5] > 0) {out <- 4 + (ct - sens[4]) * (1 / (sens[5] - sens[4]))}
if (ct >= sens[4] & sens[4] > 0 & sens[5] == 0) {out <- 4 + (ct - sens[4]) * (1 / (sens[4] - sens[3]))}
if (ct >= sens[5] & sens[5] > 0) {out <- 5}
out
}
f.tap.tthm.conc <- function(Method, ph, temp, doc, react.time, free.cl.dose, br.conc) { # TTHM concentration estimation function
if(Method == "Amy") return(0.0412 * (doc ^ 1.098) * (free.cl.dose ^ 0.152) * (br.conc ^ 0.068) * (temp ^ 0.609) * (ph ^ 1.601) * (react.time ^ 0.263))
if(Method == "Rodriguez") return(0.044 * (doc ^ 1.030) * (react.time ^ 0.262) * (ph ^ 1.149) * (free.cl.dose ^ 0.277) * (temp ^ 0.968))
}
#if(Cl.used == 1) { # parameter to be removed
Mrt <- 12
N.cstr <- 6
Ttimes <- c(0.001,0.01,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4,1.5,1.6,2,2.5,4)
Times <- Ttimes * Mrt
Times.inc <- Times - c(0,Times[1:20])
Probability <- (Times^(N.cstr - 1) * exp(-N.cstr * Times / Mrt))/(factorial(N.cstr - 1) * (Mrt / N.cstr) ^ N.cstr)
if(Cl.test.sens == 1) Cl.conc <- 0 else Cl.conc <- Cl.dose
Cl.t.inact <- 0.13
Cl.conc.t <- Cl.conc * exp(-Cl.t.inact * Times)
Cl.conc.t.inc <- Cl.conc.t * Times.inc
Cl.conc.t.cum <- cumsum(Cl.conc.t.inc)
Cl.conc.t.dist <- mean(sample(Cl.conc.t.cum, 10000, TRUE, Probability))
Cl.sens.pat <- data.frame(Cl.credits = rep(1:5, each = 6), Pathogen, Cl.sens = c(0.152,0.17,0.12,0.09,0,75,0.294,0.34,0.16,0.18,0,150,0.436,
0.52,0.2,0.245,0,216,0,1.06,0.3,0.314,0,0,0,0,0,0,0,0)) # log reductions for different "credits"
Cl.pat.reduct <- NA
for (i in 1:length(levels(Cl.sens.pat$Pathogen))) {
Cl.pat.reduct[i] <- f.Cl.pat.red(Cl.conc.t.dist, Cl.sens.pat[Cl.sens.pat$Pathogen == Pathogen[i], "Cl.sens"]) # epärobustia koodia
}
Cl.pat.reduct <- data.frame(Pathogen, Cl.pat.reduct)
tthm.tap.conc <- f.tap.tthm.conc("Amy", Cl.pH, Cl.temp, Cl.DOC, Cl.byprod.t, Cl.dose, Cl.Br)
tthm.tap.conc <- data.frame(Tthm.tap.con = tthm.tap.conc)
#}
# Vedenkulutus
k.v.kulutus <- 0.67 * v.kulutus / 1000 # keittämättömän veden kulutus (l)
# Patogeeneille altistuminen, pääketju
Exposure <- merge(raw.pat.conc, treat.pat.reduct.tot)
Exposure <- merge(Exposure, Cl.pat.reduct)
Exposure <- merge(Exposure, tthm.tap.conc)
Exposure$Tap.pat.conc <- Exposure$Raw.pat.conc * 10 ^ (-(Exposure$Treat.pat.reduct + Exposure$Cl.pat.reduct))
Exposure$Exp.pat <- Exposure$Tap.pat.conc * k.v.kulutus
# Terveysvaikutukset
# vaeston rakenne, yht 99.17678 %, problem?
vaesto <- op_baseGetData("opasnet_base", "Op_fi2652")[,c("Ikä","Result")]
colnames(vaesto) <- c("Age", "Osuus")
vaesto$Populaatio <- vaesto$Osuus * vaeston.koko
odotettu.elinika <- 81
f.exact.beta.poisson <- function(Param1, Param2, Dose) 1 - exp(-(Param1 / (Param1 + Param2)) * Dose)
f.beta.poisson.approx <- function(Param1, Param2, Dose) 1 - (1 + Dose / Param2)^(-Param1)
f.exponential <- function(Param1, Dose) 1 - exp(-Param1 * Dose)
f.alt.exponential <- function(Param1, Param2, Dose) 1 - (1 + (Dose / Param1))^(Param2)
dose.response <- op_baseGetData("opasnet_base", "Op_fi2653")
dose.response <- data.frame(Pathogen = unique(dose.response$Pathogen), Method = dose.response$Result.Text[dose.response$Parameter ==
"Function type"], Parameter1 = dose.response$Result[dose.response$Parameter == "Param1"], Parameter2 = dose.response$Result[
dose.response$Parameter == "Param2"])
dose.response <- merge(dose.response, Exposure[, c("Pathogen", "Exp.pat")])
for (i in 1:nrow(dose.response)) {
if(dose.response$Method[i] == "Exact beta Poisson") dose.response$P.inf[i] <- f.exact.beta.poisson(dose.response$Parameter1[i],
dose.response$Parameter2[i], dose.response$Exp.pat[i])
if(dose.response$Method[i] == "Beta Poisson approximation") dose.response$P.inf[i] <- f.beta.poisson.approx(dose.response$Parameter1[i],
dose.response$Parameter2[i], dose.response$Exp.pat[i])
if(dose.response$Method[i] == "Exponential") dose.response$P.inf[i] <- f.exponential(dose.response$Parameter1[i],
dose.response$Exp.pat[i])
if(dose.response$Method[i] == "Alternative exponential") dose.response$P.inf[i] <- f.alt.exponential(dose.response$Parameter1[i],
dose.response$Parameter2[i], dose.response$Exp.pat[i])
}
# infektioiden terveysvaikutukset
P.ill.g.inf <- data.frame(Pathogen, P.ill.g.inf = c(0.33, 1 - (270 / 1540), 0.9, 0.7, 0.71, 1)) # todennäköisyys sairastua kun saa infektion
# Kampylobakteeri, DALYt per infektio
P.treat.g.ill.Kamp.Gastr <- data.frame(Pathogen = Pathogen[c(1,1,1)], Outcome = "Gastroenteritis", ill.treat = c("Untreated",
"General practitioner", "Hospitalised", "Unspecified")[c(1,2,3)], P.treat.g.ill = c(0.7627, 0.2373, 0.0097))
P.treat.ill.g.inf.Kamp.Gastr <- merge(P.treat.g.ill.Kamp.Gastr, P.ill.g.inf)
P.treat.ill.g.inf.Kamp.Gastr$P.treat.ill.g.inf <- P.treat.ill.g.inf.Kamp.Gastr$P.ill.g.inf *
P.treat.ill.g.inf.Kamp.Gastr$P.treat.g.ill
duration.ill.treat.Kamp.Gastr <- data.frame(Outcome = c("Gastroenteritis"), ill.treat = c("Untreated", "General practitioner",
"Hospitalised", "Unspecified")[c(1,2,3)], dur.ill = c(5.1 / 365, 8.4 / 365, 14.39 / 365))
severity.ill.treat.Kamp.Gastr <- data.frame(Outcome = c("Gastroenteritis"), ill.treat = c("Untreated", "General practitioner",
"Hospitalised", "Unspecified")[c(1,2,3)], sev.ill = c(0.067, 0.393, 0.393))
daly.ill.treat.Kamp.Gastr <- merge(P.treat.ill.g.inf.Kamp.Gastr, duration.ill.treat.Kamp.Gastr)
daly.ill.treat.Kamp.Gastr <- merge(daly.ill.treat.Kamp.Gastr, severity.ill.treat.Kamp.Gastr)
daly.ill.treat.Kamp.Gastr$dalys <- daly.ill.treat.Kamp.Gastr$P.treat.ill.g.inf * daly.ill.treat.Kamp.Gastr$dur.ill *
daly.ill.treat.Kamp.Gastr$sev.ill
P.death.g.ill.Gastr <- 0.0004
P.death.g.inf.Gastr <- P.death.g.ill.Gastr * P.ill.g.inf$P.ill.g.inf[P.ill.g.inf$Pathogen == "Kampylobakteeri"]
death.Gastr.life.lost <- 13.2
daly.death.Kamp.Gastr <- P.death.g.inf.Gastr * death.Gastr.life.lost
## GBS Kamp.
P.gbs.g.ill <- 2e-004
P.gbs.g.inf <- P.gbs.g.ill * P.ill.g.inf$P.ill.g.inf[P.ill.g.inf$Pathogen == "Kampylobakteeri"]
dur.sev.factor.gbs <- data.frame(Outcome = c("Clinical GBS", "Residual GBS"), dur.sev.factor = c(0.29, 5.8)) # duration * severity * fraction?
daly.Kamp.gbs <- data.frame(dur.sev.factor.gbs$Outcome, dalys = dur.sev.factor.gbs$dur.sev.factor * P.gbs.g.inf)
P.death.g.gbs <- 0.08 / 3 # triangular 0.01, 0.02, 0.05
P.death.g.inf.gbs <- P.death.g.gbs * P.gbs.g.inf
death.gbs.life.lost <- 18.7
daly.death.Kamp.gbs <- P.death.g.inf.gbs * death.gbs.life.lost
## reactive arthritis Kamp.
P.arth.g.ill <- 0.02 # triangluar 0.01, 0.02, 0.03
P.arth.g.inf <- P.arth.g.ill * P.ill.g.inf$P.ill.g.inf[P.ill.g.inf$Pathogen == "Kampylobakteeri"]
duration.arth <- 6 / 52
severity.arth <- 0.21
daly.Kamp.arth <- P.arth.g.inf * duration.arth * severity.arth
# E.coli
P.wd.g.ill <- 0.53 # watery diarrhea
P.wd.g.inf <- P.wd.g.ill * P.ill.g.inf$P.ill.g.inf[P.ill.g.inf$Pathogen == "E.coli O157:H7"]
severity.wd <- 0.067
duration.wd <- 3.4 / 365
daly.wd.Ecoli <- P.wd.g.inf * severity.wd * duration.wd
P.hc.g.ill <- 0.47
P.hc.g.inf <- P.hc.g.ill * P.ill.g.inf$P.ill.g.inf[P.ill.g.inf$Pathogen == "E.coli O157:H7"]
severity.hc <- 0.39
duration.hc <- 5.6 / 365
daly.hc.Ecoli <- P.hc.g.inf * severity.hc * duration.hc
P.death.g.ill.Ecoli <- 0.00027
P.death.g.inf.Ecoli <- P.death.g.ill.Ecoli * P.ill.g.inf$P.ill.g.inf[P.ill.g.inf$Pathogen == "E.coli O157:H7"]
age.death.Ecoli <- 81 - 13.2
daly.death.Ecoli <- P.death.g.inf.Ecoli * (odotettu.elinika - age.death.Ecoli)
## Haemolytic uraemic syndrome (HUS)
P.hus.g.ill <- 0.01
P.hus.g.inf <- P.hus.g.ill * P.ill.g.inf$P.ill.g.inf[P.ill.g.inf$Pathogen == "E.coli O157:H7"]
severity.hus <- 0.93
duration.hus <- 21 / 365
daly.hus.Ecoli <- P.hus.g.inf * severity.hus * duration.hus
P.death.g.hus <- 0.04
P.death.hus.g.inf <- P.death.g.hus * P.hus.g.inf
age.death.hus.Ecoli <- 81 - 26.2
daly.death.hus.Ecoli <- P.death.hus.g.inf * (odotettu.elinika - age.death.hus.Ecoli)
## End Stage Renal Disease (ESRD)
P.esrd.g.hus <- 0.118
P.esrd.g.inf <- P.hus.g.inf * P.esrd.g.hus
severity.duration.hus <- 8.7 # severity * duration
daly.esrd.Ecoli <- P.esrd.g.inf * severity.duration.hus
P.death.g.esrd <- 0.0252
P.death.esrd.g.inf <- P.esrd.g.inf * P.death.g.esrd
age.death.esrd.Ecoli <- 81 - 34
daly.death.esrd.Ecoli <- P.death.esrd.g.inf * (odotettu.elinika - age.death.esrd.Ecoli)
# Rotavirus
P.treat.g.ill.Rotavirus <- data.frame(Pathogen = "Rotavirus", ill.treat = c("Untreated",
"General practitioner", "Hospitalised")[rep(1:3, each = 82)], Age = rep(0:81, 3), P.treat.g.ill = c(rep(0.82,5),
rep(0.95, 10), rep(0.99, 50), rep(0.97, 17), rep(0.137, 5), rep(0.0244, 5), rep(0.0511, 5), rep(0.0127, 50),
rep(0.0299, 17), rep(0.0416, 5), rep(0.0213, 5), rep(0, 72)))
P.treat.ill.g.inf.Rotavirus <- merge(P.treat.g.ill.Rotavirus, P.ill.g.inf)
P.treat.ill.g.inf.Rotavirus$P.treat.g.inf <- P.treat.ill.g.inf.Rotavirus$P.ill.g.inf * P.treat.ill.g.inf.Rotavirus$P.treat.g.ill
duration.ill.treat.Rotavirus <- data.frame(ill.treat = c("Untreated", "General practitioner","Hospitalised"), dur.ill = c(4.9 / 365,
7.1 / 365, 7.7 / 365))
severity.ill.treat.Rotavirus <- data.frame(ill.treat = c("Untreated", "General practitioner", "Hospitalised"), sev.ill = c(0.067,
0.393, 0.393))
daly.ill.treat.Rotavirus <- merge(P.treat.ill.g.inf.Rotavirus, duration.ill.treat.Rotavirus)
daly.ill.treat.Rotavirus <- merge(daly.ill.treat.Rotavirus, severity.ill.treat.Rotavirus)
daly.ill.treat.Rotavirus$dalys <- daly.ill.treat.Rotavirus$P.treat.g.inf * daly.ill.treat.Rotavirus$dur.ill *
daly.ill.treat.Rotavirus$sev.ill
P.death.Rotavirus <- data.frame(Age = 0:81, P.death.g.ill = c(rep(2.13e-005, 5), rep(0, 77)))
P.death.Rotavirus$P.death.g.inf <- P.death.Rotavirus$P.death.g.ill * P.ill.g.inf$P.ill.g.inf[P.ill.g.inf$Pathogen == "Rotavirus"]
P.death.Rotavirus$Life.lost <- odotettu.elinika - P.death.Rotavirus$Age
daly.death.Rotavirus <- data.frame(Age = P.death.Rotavirus$Age, dalys = P.death.Rotavirus$P.death.g.inf * P.death.Rotavirus$Life.lost)
# Norovirus
P.treat.g.ill.Norovirus <- data.frame(Pathogen = "Norovirus", ill.treat = c("Untreated",
"General practitioner", "Hospitalised")[rep(1:3, each = 82)], Age = rep(0:81, 3), P.treat.g.ill = c(rep(0.94876706,5),
rep(0.9902, 5), rep(0.98239, 5), rep(0.98434, 51), rep(0.992741, 16), rep(0.0448,5), rep(8.6e-003, 5), rep(0.0154, 5),
rep(0.0137, 51), rep(6.17e-003, 16), rep(6.43e-003,5), rep(1.2e-003, 5), rep(2.21e-003, 5), rep(1.96e-003, 51),
rep(8.85e-004, 16)))
P.treat.ill.g.inf.Norovirus <- merge(P.treat.g.ill.Norovirus, P.ill.g.inf)
P.treat.ill.g.inf.Norovirus$P.treat.g.inf <- P.treat.ill.g.inf.Norovirus$P.ill.g.inf * P.treat.ill.g.inf.Norovirus$P.treat.g.ill
duration.ill.treat.Norovirus <- data.frame(ill.treat = c("Untreated", "General practitioner","Hospitalised"), dur.ill = c(3.8 / 365,
5.73 / 365, 7.23 / 365))
severity.ill.treat.Norovirus <- data.frame(ill.treat = c("Untreated", "General practitioner", "Hospitalised"), sev.ill = c(0.067,
0.393, 0.393))
daly.ill.treat.Norovirus <- merge(P.treat.ill.g.inf.Norovirus, duration.ill.treat.Norovirus)
daly.ill.treat.Norovirus <- merge(daly.ill.treat.Norovirus, severity.ill.treat.Norovirus)
daly.ill.treat.Norovirus$dalys <- daly.ill.treat.Norovirus$P.treat.g.inf * daly.ill.treat.Norovirus$dur.ill *
daly.ill.treat.Norovirus$sev.ill
P.death.Norovirus <- data.frame(Age = 0:81, P.death.g.ill = c(rep(2.94e-006, 5), rep(0, 61), rep(2.04e-004, 16)))
P.death.Norovirus$P.death.g.inf <- P.death.Norovirus$P.death.g.ill * P.ill.g.inf$P.ill.g.inf[P.ill.g.inf$Pathogen == "Norovirus"]
P.death.Norovirus$Life.lost <- odotettu.elinika - P.death.Norovirus$Age
daly.death.Norovirus <- data.frame(Age = P.death.Norovirus$Age, dalys = P.death.Norovirus$P.death.g.inf * P.death.Norovirus$Life.lost)
# Cryptosporidium
P.treat.g.ill.Crypt <- data.frame(Pathogen = "Cryptosporidium", ill.treat = c("Untreated",
"General practitioner", "Hospitalised")[rep(1:3, each = 82)], Age = rep(0:81, 3), P.treat.g.ill = c(rep(0.9175730049999999,5),
rep(0.80937, 5), rep(0.6810499999999999, 5), rep(0.9774191, 50), rep(0.94706, 17), rep(0.082,5), rep(0.188, 5), rep(0.316, 5),
rep(0.0209, 50), rep(0.0367, 17), rep(4.26e-004,5), rep(2.63e-003, 5), rep(2.95e-003, 5), rep(1.66e-003, 50), rep(0.0146, 17)))
P.treat.ill.g.inf.Crypt <- merge(P.treat.g.ill.Crypt, P.ill.g.inf)
P.treat.ill.g.inf.Crypt$P.treat.g.inf <- P.treat.ill.g.inf.Crypt$P.ill.g.inf * P.treat.ill.g.inf.Crypt$P.treat.g.ill
duration.ill.treat.Crypt <- data.frame(ill.treat = c("Untreated", "General practitioner","Hospitalised"), dur.ill = c(3.5 / 365,
7 /365, 18.4 / 365))
severity.ill.treat.Crypt <- data.frame(ill.treat = c("Untreated", "General practitioner", "Hospitalised"), sev.ill = c(0.067,
0.393, 0.393))
daly.ill.treat.Crypt <- merge(P.treat.ill.g.inf.Crypt, duration.ill.treat.Crypt)
daly.ill.treat.Crypt <- merge(daly.ill.treat.Crypt, severity.ill.treat.Crypt)
daly.ill.treat.Crypt$dalys <- daly.ill.treat.Crypt$P.treat.g.inf * daly.ill.treat.Crypt$dur.ill *
daly.ill.treat.Crypt$sev.ill
P.death.Crypt <- data.frame(Age = 0:81, P.death.g.ill = c(rep(9.95e-007, 5), rep(0, 10), rep(2.09e-005, 50), rep(1.64e-003, 17)))
P.death.Crypt$P.death.g.inf <- P.death.Crypt$P.death.g.ill * P.ill.g.inf$P.ill.g.inf[P.ill.g.inf$Pathogen == "Cryptosporidium"]
P.death.Crypt$Life.lost <- odotettu.elinika - P.death.Crypt$Age
daly.death.Crypt <- data.frame(Age = P.death.Crypt$Age, dalys = P.death.Crypt$P.death.g.inf * P.death.Crypt$Life.lost)
# Giardia
P.treat.g.ill.Giardia <- data.frame(Pathogen = "Giardia", ill.treat = c("Untreated",
"General practitioner", "Hospitalised")[rep(1:3, each = 82)], Age = rep(0:81, 3), P.treat.g.ill = c(rep(0.9376,5),
rep(0.91034, 5), rep(0.72642, 5), rep(0.92486, 50), 0.54596, rep(0.5365, 16), rep(0.0609,5), rep(0.0852, 5), rep(0.272, 5),
rep(0.0721, 50), rep(0.451, 17), rep(1.5e-003,5), rep(4.46e-003, 5), rep(1.58e-003, 5), rep(3.04e-003, 51), rep(0.0125, 16)))
P.treat.ill.g.inf.Giardia <- merge(P.treat.g.ill.Giardia, P.ill.g.inf)
P.treat.ill.g.inf.Giardia$P.treat.g.inf <- P.treat.ill.g.inf.Giardia$P.ill.g.inf * P.treat.ill.g.inf.Giardia$P.treat.g.ill
duration.ill.treat.Giardia <- data.frame(ill.treat = c("Untreated", "General practitioner","Hospitalised"), dur.ill = c(10 / 365,
10 /365, 30 / 365))
severity.ill.treat.Giardia <- data.frame(ill.treat = c("Untreated", "General practitioner", "Hospitalised"), sev.ill = c(0.067,
0.393, 0.393))
daly.ill.treat.Giardia <- merge(P.treat.ill.g.inf.Giardia, duration.ill.treat.Giardia)
daly.ill.treat.Giardia <- merge(daly.ill.treat.Giardia, severity.ill.treat.Giardia)
daly.ill.treat.Giardia$dalys <- daly.ill.treat.Giardia$P.treat.g.inf * daly.ill.treat.Giardia$dur.ill *
daly.ill.treat.Giardia$sev.ill
# yhteenveto DALYistä
Health.effects <- vaesto[,c("Age","Populaatio")]
Health.effects$Untreated.Gastr.Kamp <- daly.ill.treat.Kamp.Gastr[daly.ill.treat.Kamp.Gastr$ill.treat == "Untreated", c("dalys")]
Health.effects$GP.Gastr.Kamp <- daly.ill.treat.Kamp.Gastr[daly.ill.treat.Kamp.Gastr$ill.treat == "General practitioner", c("dalys")]
Health.effects$Hospitalised.Gastr.Kamp <- daly.ill.treat.Kamp.Gastr[daly.ill.treat.Kamp.Gastr$ill.treat == "Hospitalised", c("dalys")]
Health.effects$Death.Gastr.Kamp <- daly.death.Kamp.Gastr
Health.effects$Clinical.GBS.Kamp <- daly.Kamp.gbs$dalys[1]
Health.effects$Residual.GBS.Kamp <- daly.Kamp.gbs$dalys[2]
Health.effects$Death.GBS.Kamp <- daly.death.Kamp.gbs
Health.effects$Arth.Kamp <- daly.Kamp.arth
Health.effects$WD.Ecoli <- daly.wd.Ecoli
Health.effects$HC.Ecoli <- daly.hc.Ecoli
Health.effects$Death.Ecoli <- daly.death.Ecoli
Health.effects$HUS.Ecoli <- daly.hus.Ecoli
Health.effects$Death.HUS.Ecoli <- daly.death.hus.Ecoli
Health.effects$ESRD.Ecoli <- daly.esrd.Ecoli
Health.effects$Death.ESRD.Ecoli <- daly.death.esrd.Ecoli
Health.effects <- merge(Health.effects, daly.ill.treat.Rotavirus[daly.ill.treat.Rotavirus$ill.treat == "Untreated", c("Age", "dalys")])
colnames(Health.effects)[ncol(Health.effects)] <- "Untreated.Rotavirus"
Health.effects <- merge(Health.effects, daly.ill.treat.Rotavirus[daly.ill.treat.Rotavirus$ill.treat == "General practitioner", c("Age", "dalys")])
colnames(Health.effects)[ncol(Health.effects)] <- "GP.Rotavirus"
Health.effects <- merge(Health.effects, daly.ill.treat.Rotavirus[daly.ill.treat.Rotavirus$ill.treat == "Hospitalised", c("Age", "dalys")])
colnames(Health.effects)[ncol(Health.effects)] <- "Hospitalised.Rotavirus"
Health.effects <- merge(Health.effects, daly.death.Rotavirus)
colnames(Health.effects)[ncol(Health.effects)] <- "Death.Rotavirus"
Health.effects <- merge(Health.effects, daly.ill.treat.Norovirus[daly.ill.treat.Norovirus$ill.treat == "Untreated", c("Age", "dalys")])
colnames(Health.effects)[ncol(Health.effects)] <- "Untreated.Norovirus"
Health.effects <- merge(Health.effects, daly.ill.treat.Norovirus[daly.ill.treat.Norovirus$ill.treat == "General practitioner", c("Age", "dalys")])
colnames(Health.effects)[ncol(Health.effects)] <- "GP.Norovirus"
Health.effects <- merge(Health.effects, daly.ill.treat.Norovirus[daly.ill.treat.Norovirus$ill.treat == "Hospitalised", c("Age", "dalys")])
colnames(Health.effects)[ncol(Health.effects)] <- "Hospitalised.Norovirus"
Health.effects <- merge(Health.effects, daly.death.Norovirus)
colnames(Health.effects)[ncol(Health.effects)] <- "Death.Norovirus"
Health.effects <- merge(Health.effects, daly.ill.treat.Crypt[daly.ill.treat.Crypt$ill.treat == "Untreated", c("Age", "dalys")])
colnames(Health.effects)[ncol(Health.effects)] <- "Untreated.Crypt"
Health.effects <- merge(Health.effects, daly.ill.treat.Crypt[daly.ill.treat.Crypt$ill.treat == "General practitioner", c("Age", "dalys")])
colnames(Health.effects)[ncol(Health.effects)] <- "GP.Crypt"
Health.effects <- merge(Health.effects, daly.ill.treat.Crypt[daly.ill.treat.Crypt$ill.treat == "Hospitalised", c("Age", "dalys")])
colnames(Health.effects)[ncol(Health.effects)] <- "Hospitalised.Crypt"
Health.effects <- merge(Health.effects, daly.death.Crypt)
colnames(Health.effects)[ncol(Health.effects)] <- "Death.Crypt"
Health.effects <- merge(Health.effects, daly.ill.treat.Giardia[daly.ill.treat.Giardia$ill.treat == "Untreated", c("Age", "dalys")])
colnames(Health.effects)[ncol(Health.effects)] <- "Untreated.Giardia"
Health.effects <- merge(Health.effects, daly.ill.treat.Giardia[daly.ill.treat.Giardia$ill.treat == "General practitioner", c("Age", "dalys")])
colnames(Health.effects)[ncol(Health.effects)] <- "GP.Giardia"
Health.effects <- merge(Health.effects, daly.ill.treat.Giardia[daly.ill.treat.Giardia$ill.treat == "Hospitalised", c("Age", "dalys")])
colnames(Health.effects)[ncol(Health.effects)] <- "Hospitalised.Giardia"
Health.effects <- reshape(Health.effects, idvar = c("Age"), times = colnames(Health.effects)[-c(1,2)], timevar = "Outcome",
varying = list(colnames(Health.effects)[-c(1,2)]), direction = "long")
colnames(Health.effects)[4] <- "P.daly.g.inf"
Health.effects$Pathogen <- NA
Health.effects$Pathogen[grep(".Kamp", Health.effects$Outcome)] <- Pathogen[1]
Health.effects$Pathogen[grep(".Ecoli", Health.effects$Outcome)] <- Pathogen[2]
Health.effects$Pathogen[grep(".Rotavirus", Health.effects$Outcome)] <- Pathogen[3]
Health.effects$Pathogen[grep(".Norovirus", Health.effects$Outcome)] <- Pathogen[4]
Health.effects$Pathogen[grep(".Crypt", Health.effects$Outcome)] <- Pathogen[5]
Health.effects$Pathogen[grep(".Giardia", Health.effects$Outcome)] <- Pathogen[6]
Health.effects <- merge(Health.effects, dose.response[,c("Pathogen", "P.inf")])
Health.effects$DALYs <- (1 - (1 - Health.effects$P.inf * Health.effects$P.daly.g.inf)^365) * Health.effects$Populaatio
# TTHM terveysvaikutukset (ylimääräiset virtsarakon syövät)
bladder.cancer.incidence <- data.frame(Age = 0:81, Incidence = c(rep(0, 5), rep(1.7, 5), rep(0, 5), rep(1.3, 5), rep(1, 5), rep(0, 5),
rep(1.1, 5), rep(5.5, 5), rep(11.7, 5), rep(30.6, 5), rep(53, 5), rep(85.2, 5), rep(173.7, 5), rep(245.3, 5), rep(304.6, 5),
rep(389.8, 5), rep(441.3, 2)))
bladder.cancer.slope <- log(1.006134)
bladder.cancer.OR <- data.frame(OR = exp(bladder.cancer.slope * tthm.tap.conc[1,1] * k.v.kulutus))
bladder.cancer.excess <- merge(bladder.cancer.incidence, vaesto)
bladder.cancer.excess <- merge(bladder.cancer.excess, bladder.cancer.OR)
bladder.cancer.excess$Cases <- bladder.cancer.excess$Populaatio * bladder.cancer.excess$Incidence * (bladder.cancer.excess$OR - 1) /
bladder.cancer.excess$OR / 100000
bladder.cancer.dalys <- bladder.cancer.excess[, c("Age", "Cases")]
bladder.cancer.dalys$Survival.rate <- c(rep(0.83, 55), rep(0.7, 10), rep(0.65, 10), rep(0.51, 7))
bladder.cancer.dalys$Survivors <- bladder.cancer.dalys$Cases * bladder.cancer.dalys$Survival.rate
bladder.cancer.dalys$Terminal.cases <- bladder.cancer.dalys$Cases * (1 - bladder.cancer.dalys$Survival.rate)
bladder.cancer.dalys <- reshape(bladder.cancer.dalys[,c("Age", "Survivors", "Terminal.cases")], idvar = "Age", times = c("Survivors",
"Terminal.cases"), timevar = "Cases", varying = list(c("Survivors", "Terminal.cases")), direction = "long")
colnames(bladder.cancer.dalys)[3] <- "Count"
cancer.clinical.free.dur <- data.frame(Cases = c("Survivors", "Terminal.cases"), Clinical.free.state = c(4.88, 2.55))
cancer.initial.treatment.dur <- data.frame(Initial.treatment = 0.12)
cancer.pre.terminal.dur <- data.frame(Cases = "Terminal.cases", Pre.terminal = 0.75)
cancer.terminal.dur <- data.frame(Cases = "Terminal.cases", Terminal = 0.08)
cancer.death.dur <- data.frame(Cases = "Terminal.cases", Age = 0:81, Death = 81:0)
cancer.dur <- merge(bladder.cancer.dalys[,c("Age", "Cases")], cancer.clinical.free.dur, all = TRUE)
cancer.dur <- merge(cancer.dur, cancer.initial.treatment.dur, all = TRUE)
cancer.dur <- merge(cancer.dur, cancer.pre.terminal.dur, all = TRUE)
cancer.dur <- merge(cancer.dur, cancer.terminal.dur, all = TRUE)
cancer.dur <- merge(cancer.dur, cancer.death.dur, all = TRUE)
cancer.dur <- reshape(cancer.dur, idvar = c("Age", "Cases"), times = colnames(cancer.dur)[-c(1,2)], timevar = "Outcome",
varying = list(colnames(cancer.dur)[-c(1,2)]), direction = "long")
colnames(cancer.dur)[4] <- "Duration"
cancer.clinical.free.sev <- data.frame(Clinical.free.state = 0.18)
cancer.initial.treatment.sev <- data.frame(Initial.treatment = 0.27)
cancer.pre.terminal.sev <- data.frame(Cases = "Terminal.cases", Pre.terminal = 0.64)
cancer.terminal.sev <- data.frame(Cases = "Terminal.cases", Terminal = 0.93)
cancer.death.sev <- data.frame(Cases = "Terminal.cases", Age = 0:81, Death = 1)
cancer.sev <- merge(bladder.cancer.dalys[,c("Age", "Cases")], cancer.clinical.free.sev, all = TRUE)
cancer.sev <- merge(cancer.sev, cancer.initial.treatment.sev, all = TRUE)
cancer.sev <- merge(cancer.sev, cancer.pre.terminal.sev, all = TRUE)
cancer.sev <- merge(cancer.sev, cancer.terminal.sev, all = TRUE)
cancer.sev <- merge(cancer.sev, cancer.death.sev, all = TRUE)
cancer.sev <- reshape(cancer.sev, idvar = c("Age", "Cases"), times = colnames(cancer.sev)[-c(1,2)], timevar = "Outcome",
varying = list(colnames(cancer.sev)[-c(1,2)]), direction = "long")
colnames(cancer.sev)[4] <- "Severity"
bladder.cancer.dalys <- merge(bladder.cancer.dalys, cancer.dur)
bladder.cancer.dalys <- merge(bladder.cancer.dalys, cancer.sev)
bladder.cancer.dalys$DALYs <- bladder.cancer.dalys$Count * bladder.cancer.dalys$Duration * bladder.cancer.dalys$Severity
# Matala syntymäpaino
syntymia.vuodessa <- vaeston.koko * 1170 / 1e5
sga.esiintyvyys <- syntymia.vuodessa * 0.1
sga.slope <- 0.0009955
sga.OR <- exp(sga.slope * tthm.tap.conc[1,1] * k.v.kulutus)
sga.excess <- (sga.OR - 1) * sga.OR * sga.esiintyvyys
#
temp <- merge(dose.response, P.ill.g.inf)
############# TULOKSET #########################################################################################################
raw.pat.conc # Patogeenien konsentraatio raakavedessä
treat.pat.reduct.tot # Patogeenien log vähenemä puhdistuksessa
dose.response[,c("Pathogen", "Exp.pat", "P.inf")] # Patogeeneille altistuminen ja infektion todennäköisyys
cat(sum((1 - (1 - temp$P.ill.g.inf * temp$P.inf)^365) * vaeston.koko), " vatsatautia vuodessa \n")
cat(sum(Health.effects$DALYs), " DALY:ä vatsataudeista \n")
#cat(sum(bladder.cancer.dalys$DALYs, na.rm = TRUE), " DALY:ä virtsarakon syövistä vuodessa johtuen TTHM pitoisuudesta \n")
#cat(sga.excess, " vastasyntynyttä vuodessa kärsii matalasta syntymäpainosta johtuen TTHM pitoisuudesta \n")
| |
