+ Näytä koodi- Piilota koodi
library(OpasnetUtils)
library(ggplot2)
library(reshape2)
openv.setN(200)
objects.latest("Op_en6007", code_name = "answer") # OpasnetUtils/Drafts # Tarvitaan vielä unkeep
##########################################################################
# KYSELYTIEDOT SILAKASTA
silakka <- opbase.data("Op_fi3831", subset = "Silakka") # [[:op_fi:Silakan hyöty-riskiarvio]]
silakka$Paino[silakka$Paino == 0 & silakka$Ikäryhmä == "Aikuinen"] <- 75 #Tehdään karkea inputointi
silakka$Paino[silakka$Paino == 0 & silakka$Ikäryhmä == "Lapsi"] <- 30
# Inputoinnin voisi tehdä myös iän perusteella, mutta pitäisi kaivaa lasten painotieto. (aikuiset datasta)
# Arvotaan henkilöitä N kpl jokaiseen sukupuoli-maakunta-ikäryhmään.
silakka$Age <- cut(silakka$Ikä, breaks = c(0,1,(1:20)*5), labels = c("0", "1-4", "5-9", "10-14", "15-19", "20-24",
"25-29", "30-34", "35-39", "40-44", "45-49", "50-54", "55-59", "60-64", "65-69", "70-74", "75-79", "80-84",
"85-89", "90-94", "95-"), right = FALSE
)
ehto <- unique(silakka[c("Age")]) #, "Sukupuoli", "Maakunta")]) Eiköhän painokerroin huolehdi sukupuolen ja maakunnan
out <- data.frame()
for(i in 1:nrow(ehto)) {
temp <- silakka[
silakka$Age == ehto$Age[i]
# & silakka$Sukupuoli == ehto$Sukupuoli[i]
# & silakka$Maakunta == ehto$Maakunta[i]
, ]
temp <- temp[sample(1:nrow(temp), get("N", envir = openv), replace = TRUE, prob = temp$Painokerroin) , ]
out <- rbind(out, data.frame(Iter = 1:nrow(temp), temp))
}
# Luo ovariablet arvotun datan perusteella
rannikko <- c("Uusimaa", "Pohjanmaa", "Kymenlaakso", "Etelä-Pohjanmaa", "Satakunta", "Keski-Pohjanmaa", "Pohjois-Pohjanmaa", "Varsinais-Suomi")
sisämaa <- c("Kanta-Häme", "Pirkanmaa", "Etelä-Karjala", "Pohjois-Savo", "Pohjois-Karjala", "Etelä-Savo", "Keski-Suomi", "Päijät-Häme", "Lappi", "Kainuu")
# pitkä = pitkä lista yksilökohtaisia määrittelyjä
pitkä <- out[c("Iter", "Age", "Sukupuoli", "Ikä", "Paino", "Maakunta")] # "Ikäryhmä", OTETAAN POIS TARPEETTOMANA.
pitkä$Rannikko <- ifelse(out$Maakunta %in% rannikko, "Rannikko", "Sisämaa")
# lyhyt = lyhyt lista yksilökohtaisia määrittelyjä eli vain välttämättömät Iter ja Ikäryhmä.
lyhyt <- pitkä[c("Iter", "Age", "Ikä", "Sukupuoli")] # Tähän on lisättävä Age ja Sukupuoli, jos/kun ERF tai muu tieto indeksoidaan niillä.
kokoaik <- Ovariable("kokoaik", data = data.frame(lyhyt, Arvo = out$Result, Result = out$Kokosilakka)) # Kokonaisten silakoiden syönti kertaa/3 kk
osanaik <- Ovariable("osanaik", data = data.frame(lyhyt, Result = out$Silakkaruoka)) # Silakkaruokien syönti kertaa/3 kk
lisuaik <- Ovariable("lisuaik", data = data.frame(lyhyt, Result = out$Silakkalisuke)) # Silakkalisukkeiden syönti kertaa/3 kk
BW <- Ovariable("BW", data = data.frame(lyhyt, Result = out$Paino)) # Ruumiinpaino
######################################################################################
# SILAKOIDEN ANNOSKOOT
solet <- opbase.data("Op_fi3831.kalansyonti") # Silakkaoletukset sivulta [[:op_fi:Silakan hyöty-riskiarvio]]
silakoita <- Ovariable("silakoita", data = solet[solet$Muuttuja == "V95" , c("Arvo", "Result")]) # Silakoita per silakka-annos
silakanpaino <- Ovariable("silakanpaino", data = solet[solet$Muuttuja == "Koko silakan paino" , ]["Result"])
raakaainepaino <- Ovariable("raakaainepaino", data = data.frame(
Ikä = 0:100,
Result = c(
rep(solet$Result[solet$Muuttuja == "Raaka-ainesilakan paino"], 15),
rep(solet$Result[solet$Muuttuja == "Raaka-ainesilakan paino, lapset"], 86)
)
))
lisukepaino <- Ovariable("lisukepaino", data = solet[solet$Muuttuja == "Lisukesilakan paino" , ]["Result"])
määrä <- (kokoaik * silakoita * silakanpaino + osanaik * raakaainepaino + lisuaik * lisukepaino ) / 91 # Per 3 kk -> per d
result(määrä)[result(määrä) == 0] <- 0.01 # Ei jätetä nollia saantiin
#määrä@marginal <- määrä@marginal & ! colnames(määrä@output) %in% c("Ikä", "Paino", "Arvo") # Tarvitaanko vai voiko poistaa?
# Luo muuttuja, jossa ovat kaikki kiinnostavat selittävät tekijät
tiedot <- määrä * Ovariable(data = data.frame(pitkä, Result = 1))
tiedot <- unkeep(tiedot, cols = "Arvo", prevresults = TRUE, sources = TRUE)
tietomarg <- colnames(tiedot@output)[tiedot@marginal]
colnames(tiedot@output)[colnames(tiedot@output) == "Result"] <- "Silakkamäärä"
tiedot@output$Result <- 1
tiedot@marginal <- colnames(tiedot@output) %in% tietomarg
#####Jointataan ikäryhmä, jotta saadaan johdonmukaiset terveysvasteet (tosin insidenssit pitäisi myös laskea ikäryhmittäin!)
# TÄTÄ EI SIIS TARVITA, KOSKA Age on MARGINAL!
#ikäjoint <- Ovariable(output = data.frame(määrä@output[c("Iter", "Ikäryhmä")], Result = 1))
#ikäjoint <- CollapseMarginal(ikäjoint, cols = "Ikäryhmä", fun = "sample", probs = c(
# 0.8358142231, # yli 15-vuotiaita (Aikuinen). Tiedot Tilastokeskus 2012.
# 0.1641857769 # 0-14-vuotiaita (Lapsi)
#))
#
#määrä <- unkeep(määrä * ikäjoint, cols = "Ikäryhmä")
#
#BW <- unkeep(BW * ikäjoint, cols = "Ikäryhmä")
#
#tiedot <- tiedot * ikäjoint
#################################################################################
# PITOISUUDET: PCDD/F, PCB, OMEGA3; JA MUUT TAPAUSKOHTAISET OVARIABLET
arvo <- function(x, measure.vars) {
out <- x[sample(1:nrow(x), get("N", envir = openv), replace = TRUE) , ]
out <- data.frame(Iter = 1:nrow(out), out)
out <- melt(
out,
measure.vars = measure.vars,
variable.name = "Exposure_agent",
value.name = "Result"
)
return(out)
}
# PCDD/F ja PCB-pitoisuudet # pg /g f.w. nd = LOQ
objects.get("x89WlgJlDLA02vnL") # PCDD/F-data table1
pcdd <- table1[2:nrow(table1), c(2,3)] # Aika tiputetaan pois toistaiseksi tarpeettomana (2002, 2009 dataa)
colnames(pcdd) <- c("PCDDF", "PCB")
pcdd$PCDDF <- as.numeric(levels(pcdd$PCDDF)[pcdd$PCDDF])
pcdd$PCB <- as.numeric(levels(pcdd$PCB)[pcdd$PCB])
pcdd$TEQ <- pcdd$PCDDF + pcdd$PCB
pcdd <- arvo(pcdd, measure.vars = c("PCDDF", "PCB", "TEQ"))
# Hae pitoisuustiedot: D-vitamiini, EPA, DHA, Omega-3. D-vitamiini: ug /100 g, rasvahapot: mg /100 g
objects.get("qhV9V72f9zLmXX50") # Ravintoainedata table
ravinteet <- table[2:nrow(table), 3:6]
colnames(ravinteet) <- c("Vitamin_D", "EPA", "DHA", "Omega3")
ravinteet <- arvo(ravinteet, measure.vars = c("Vitamin_D", "EPA", "DHA", "Omega3"))
ravinteet$Result <- as.numeric(ravinteet$Result) / 100 # per 100 g -> per g
concentration <- Ovariable("concentration", data = rbind(ravinteet, pcdd))
################ ALTISTUS
exposure <- unkeep(määrä, cols = "Arvo", prevresults = TRUE, sources = TRUE) * concentration
frexposed <- 1
bgexposure <- 0
######################################################### VÄESTÖ
# Tarkastellaan totcases-laskennassa aluksi yksilöriskiä, ja vasta lopussa kerrotaan yksilötulokset yksilöiden edustamien ryhmien koolla.
population <- Ovariable("population", data = data.frame(Result = 1))
# Väkimäärä. TÄTÄ KÄYTETÄÄN disincidence-ovariablen suhteuttamisessa CHD:n suhteen.
pop <- opbase.data("Op_en2949") # [[Population of Finland]]
pop$Obs <- NULL
pop$Observation <- NULL
age <- pop$Result[pop$Age == "0-4"]
pop$Result[pop$Age == "0-4"] <- 4/5 * age
pop <- orbind(pop, data.frame(Age = "0", Result = 1/5 * age))
levels(pop$Age)[levels(pop$Age) == "0-4"] <- "1-4"
pop <- Ovariable("pop", data = pop)
###############################################################################################
# ANNOSVASTEET JA TAUTIRISKIT
# dioksiinille,
# omega-3-rasvahapoille,
# vitamiineille (D-vitamiinille)
# ja metyylielohopealle ja niiden vaikutuksille
# Kuitenkaan metyylielohopeapitoisuuksia ei ole tässä malliversiossa, joten ne tippuvat pois.
# Annosvasteet
objects.latest("Op_en5823", code_name = "initiate") # [[ERF of dioxin]], ovariables ERF, threshold
temp1 <- ERF
temp2 <- threshold
objects.latest("Op_en5830", code_name = "initiate") # [[ERF of omega-3 fatty acids]], ovariables ERF, threshold
temp1@data <- orbind(ERF@data, temp1@data)
temp2@data <- orbind(threshold@data, temp2@data)
objects.latest("Op_en6866", code_name = "initiate") # [[ERFs of vitamins]], ovariables ERF, threshold
temp1@data <- orbind(ERF@data, temp1@data)
temp2@data <- orbind(threshold@data, temp2@data)
objects.latest("Op_en5825", code_name = "initiate") # [[ERF of methylmercury]], ovariables ERF, threshold
ERF@data <- orbind(ERF@data, temp1@data)
threshold@data <- orbind(threshold@data, temp2@data)
# Drop nuisance indices because they use a lot of memory in oapply.
dropp <- c("Response_metric", "Exposure_route", "Exposure_metric", "Exposure_unit")
ERF <- unkeep(EvalOutput(ERF), dropp, sources = TRUE)
threshold <- unkeep(EvalOutput(threshold), dropp, sources = TRUE)
######################################## Tautiriski
objects.latest("Op_en5917", code_name = "initiate") # [[:op_en:Disease risk]] ovariable disincidence
disincidence <- disincidence / 100000 # Change from 1/100000py to 1/py.
disincidence <- unkeep(disincidence, cols = c("Age", "Sex", "Population", "Unit", "Response_metric"))
# Näitä indeksejä ei tarvita tässä arvioinnissa.
disincidence@output <- disincidence@output[disincidence@output$Trait != "CHD" , ] # Tulee ikävakioidusta alta.
############################################################################
# Ikävakioidut kuolemansyyt
# Sairastuvuus (tarkemmin sanottuna [[Kuolemansyyt Suomessa]])
syy <- Ovariable("syy", data = opbase.data("Op_fi4558", subset = "2012", include = list(Kuolemansyy = c(
"04-21 Syövät (C00-C97)", "27 Iskeemiset sydäntaudit (I20-I25)"
))))
syy@data$Ikä <- gsub(" ", "", syy@data$Ikä)
colnames(syy@data)[colnames(syy@data) == "Ikä"] <- "Age"
# Syöpää ei haluta linkata ikäryhmittäin ERF:iin, vaan sitä käytetään muodostamaan diskontattu elinaikainen
# odotettu eliniän menetys jos kuolee syöpään tietyn ikäisenä. Ks. myöhemmin.
väli <- Ovariable(output = data.frame(
Kuolemansyy = c("27 Iskeemiset sydäntaudit (I20-I25)", "04-21 Syövät (C00-C97)"),
Trait = c("CHD", "Syövät"), # Syövät nimetään hassusti, koska sen EI haluta linkkautuvan ERF:iin.
Result = 1
))
syyt <- syy / (pop / 2) * väli # Pop ei ole sukupuolen mukaan jaoteltu mutta syy on, joten pistetään kahtia.
marginals <- union(colnames(syyt@output)[syyt@marginal], colnames(disincidence@output)[disincidence@marginal])
syyt@output <- orbind(disincidence, syyt)
syyt@marginal <- colnames(syyt@output) %in% marginals
syyt <- unkeep(syyt, cols = c("Kuolemansyy", "Observation"), prevresults = TRUE, sources = TRUE)
syyt@output$Sukupuoli[syyt@output$Sukupuoli == "Miehet"] <- "Mies" # Ei ole faktori vaan character
syyt@output$Sukupuoli[syyt@output$Sukupuoli == "Naiset"] <- "Nainen"
disincidence <- syyt # Nyt sisältää sekä alkuperäisen disincidencen että ikäluokittaisen CHD:n
###################################################################################
#### HIA-ovariablet
objects.latest("Op_en2261", code_name = "initiate") # [[:op_en:Health impact assessment]] ovariables totcases, AF
## Tilapäisesti totcases on täällä, koska se kaipaa korjaamista
totcases <- Ovariable("totcases", # This calculates the total number of cases in each population subgroup.
# The cases are calculated for specific (combinations of) causes. However, these causes are NOT visible in the result.
dependencies = data.frame(Name = c(
"population", # Population divided into subgroups as necessary
"disincidence", # Incidence of the disease of interest
"ERF", # Exposure-response function of the pollutants or agents (RR for unit exposure)
"exposure", # Exposure to the pollutants
"frexposed", # fraction of population that is exposed
"bgexposure", # Background exposure (a level below which you cannot get in practice)
"threshold", # exposure level below which the agent has no impact.
"BW" # body weight
)),
formula = function(...) {
#################################################################33
######### Body weight scaling: In some cases, exposure is given as per body weight and in some cases as absolute amounts.
# Here we add one index to define for this difference.
scaling <- Ovariable("scaling", data = data.frame(
ERF_parameter = unique(ERF@output$ERF_parameter),
Result = 1 * grepl("bw", unique(ERF@output$ERF_parameter)))
# if ERF parameter contains bw, scaling is TRUE, i.e. 1.
)
expo <- exposure / (((BW - 1) * scaling) + 1) # If scaling is 0, BW cancels out.
expo <- expo * Ovariable( # Create alternative scenario with zero exposure.
output = data.frame(Exposcen = c("BAU", "No exposure"), Result = c(1, 0)),
marginal = c(TRUE, FALSE)
)
####################################################################
####### This part is about risks relative to background.
# Calcualte the risk ratio to each subgroup based on the exposure in that subgroup.
# Combine pollutant-specific RRs by multiplying. For description, see [[Exposure-response function]].
test <- FALSE
#First take the relative risk estimates
ERFrr <- ERF
ERFrr@output <- ERF@output[ERF@output$ERF_parameter %in% c("RR", "RR bw") , ]
if(nrow(ERFrr@output) > 0) {
bigger <- (threshold > bgexposure) * threshold + (1 - (threshold > bgexposure)) * bgexposure # Choose bigger
RRrr <- exp(log(ERFrr) * (expo - bigger)) # Actual function
RR <- (RRrr - 1) * (expo > threshold) + 1 # RR is 1 below threshold
test <- TRUE
}
# Then take the relative Hill estimates
ED50 <- threshold
ED50@output <- ED50@output[ED50@output$ERF_parameter %in% c("Relative Hill", "Relative Hill bw") , ]
if(nrow(ED50@output) > 0) {
RRrhill <- 1 + (expo * ERF) / (expo + ED50) # ERF has parameter value for Imax. If Imax < 0, risk reduces.
if(!test) RR <- RRhill else {
marginals <- union(
colnames(RR@output)[RR@marginal],
colnames(RRrhill@output)[RRrhill@marginal]
)
RR@output <- orbind(RR * 1, RRrhill * 1)
RR@marginal <- colnames(RR@output) %in% marginals
}
test <- TRUE
}
# We need OR but not yet crucial, so let's postpone this. See [[:op_en:Converting between exposure-response parameters]]
# #Then take the odds ratio estimates
# OR <- ERF
# OR@output <- ERF@output[ERF@output$ERF_parameter %in% c("OR", "OR bw") , ]
# if(nrow(OR@output) > 0) {
# bigger <- (threshold > bgexposure) * threshold + (1 - (threshold > bgexposure)) * bgexposure # Choose bigger
# OR <- RR = OR/( 1-PX0+OR*PX0 ) # Actual function where PX0 is the background incidence. How to write a code?
# RR <- (RRrr - 1) * (expo > threshold) + 1 # RR is 1 below threshold
# test <- TRUE
# }
# Dilute the risk in the population if not all are exposed i.e. frexposed < 1.
RR <- frexposed * (RR - 1) + 1
RR <- unkeep(RR, prevresults = TRUE, sources = TRUE)
RR <- oapply(RR, cols = "Exposure_agent", FUN = prod)
# Remove redundant columns.
population@output <- dropall(population@output)
population <- unkeep(population, sources = TRUE, prevresults = TRUE)
disincidence <- unkeep(disincidence, sources = TRUE, prevresults = TRUE)
# takeout is a vector of column names of indices that are NOT in the disease incidence.
takeout <- setdiff(
colnames(population@output)[population@marginal],
colnames(disincidence@output)[disincidence@marginal]
)
# pci is the proportion of cases across different population subroups based on differential risks and
# population sizes. pci sums up to 1 for each larger subgroup found in disincidence.
# See [[Population attributable fraction]].
pci <- unkeep(population * RR, prevresults = TRUE, sources = TRUE)
# Divide pci by the values of the actually exposed group (discard nonexposed)
temp <- pci
temp@output <- temp@output[temp@output$Exposcen == "BAU" , ]
temp <- unkeep(temp, cols = "Exposcen")
pci <- pci / oapply(temp, cols = takeout, FUN = sum)
if(length(takeout > 0)) population <- oapply(population, cols = takeout, FUN = sum) # Aggregate to larger subgroups.
# population@output <- population@output[!is.na(result(population)) , ] # Not sure if this is necessary.
# The cases are divided into smaller subgroups based on weights in pci.
# Redundant columns are removed because totcases are used within attrcases which has same columns.
out <- disincidence * population * pci
##########################################################################
############# This part is about absolute risks (i.e., risk is not affected by background rates).
# Unit risk (UR) and cancer slope factor (CSF) estimates.
# Exposure-response slope (ERS)
UR <- ERF
UR@output <- UR@output[UR@output$ERF_parameter %in% c("UR", "CSF", "ERS", "UR bw", "CSF bw", "ERS bw") , ]
if(nrow(UR@output) > 0) {
UR <- threshold + UR * expo * frexposed # Actual equation
# threshold is here interpreted as the baseline response (intercept of the line). It should be 0 for
# UR and CSF but it may have meaningful values with ERS
UR <- unkeep(UR, prevresults = TRUE, sources = TRUE)
UR <- oapply(UR, cols = "Exposure_agent", FUN = sum)
UR@output <- UR@output[!is.na(result(UR)) , ] # Remove empty rows
out2 <- population * UR
} else {
out2 <- out
out2@output <- out2@output[FALSE , ]
}
# Step estimates: value is 1 below threshold and above ERF, and 0 in between.
Step <- ERF
Step@output <- Step@output[Step@output$ERF_parameter %in% c(
"Step", "Step bw", "ADI", "ADI bw", "TDI", "TDI bw", "RDI", "RDI bw", "NOAEL", "NOAEL bw") , ]
if(nrow(Step@output) > 0) {
Step <- 1 - (expo >= threshold) * (expo <= Step) # Actual equation
Step <- unkeep(Step, prevresults = TRUE, sources = TRUE)
Step <- oapply(Step, cols = "Exposure_agent", FUN = sum)
Step@output <- Step@output[!is.na(result(Step)) , ] # Remove empty rows
out3 <- (population * 0 + 1) * Step # This is to maintain the ovariable structure
} else {
out3 <- out
out3@output <- out3@output[FALSE , ]
}
#####################################################################
# Combining effects
marginals <- union(union(
colnames(out@output)[out@marginal],
colnames(out2@output)[out2@marginal]),
colnames(out3@output)[out3@marginal]
)
out@output <- orbind(out, out2)
out@output <- orbind(out, out3)
out@marginal <- colnames(out@output) %in% marginals
out <- unkeep(out, sources = TRUE, prevresults = TRUE)
out@output <- dropall(out@output[!is.na(result(out)) , ]) # Clean the output for missing locations
return(out)
}
)
totcases <- EvalOutput(totcases)
##################################################################################
# Tapauskohtaiset postprosessoinnit
indivrisk <- totcases
#indivrisk <- unkeep(totcases, cols = c("Population", "Unit", "Exposure_unit", "Exposure_metric", "Exposure_route", "Response_metric")) # "Age", "Sex",
## Poistetaan hammasvaste, koska se on matemaattisesti epäselvä ja ne annosvasteet jotka ovat kahdesti.
indivrisk@output <- indivrisk@output[ ! indivrisk@output$Trait %in% c(
"Developmental dental defects incl. agenesis",
"Cancer, total",
"Coronary heart disease",
"Child's intelligence"
), ]
indivrisk@output$Ikä <- as.numeric(levels(indivrisk@output$Ikä)[indivrisk@output$Ikä]) # orbind in totcases confuses.
indivrisk <- indivrisk * tiedot # Palautetaan kyselystä Ikä, Sukupuoli, Maakunta, Paino ja Silakkamäärä.
indivrisk@marginal[colnames(indivrisk@output) == "Age"] <- TRUE # Jostain syystä marginaali tippuu pois; palautetaan.
#Etsitään hedelmällisessä iässä olevat naiset ja lasten ÄO-vaste ja korjataan synnytyksen todennäköisyydellä.
mother <- indivrisk@output$Sukupuoli == "Nainen" & indivrisk@output$Ikä > 17 & indivrisk@output$Ikä < 39
iq <- indivrisk@output$Trait == "Child's IQ"
result(indivrisk) <- result(indivrisk) * ifelse(iq, ifelse(mother, 0.1, 0), 1) # Probability of birth during a year.
########################### Analyysit
cat("Yksilöriski: tapausmäärä / 1000 henkilövuotta. \n")
oprint(summary(indivrisk * 1000))
if(exists("server")) BS <- 24 else BS <- 18
ggplot(indivrisk@output, aes(x = Trait, y = Result, colour = Exposcen)) + geom_point()
ggplot(indivrisk@output, aes(x = Exposcen, weight = Result, fill = Trait)) + geom_bar()
ggplot(indivrisk@output, aes(x = Age, weight = Result, fill = Trait)) + geom_bar()
#############################
# DALYT
dweight <- Ovariable("dweight", ddata = "Op_fi3831.laatupainotetut_elinvuodet_vaikutuksille")
dweight@data$Response_metric <- NULL
#lkm <- Ovariable(output = data.frame(
# Trait = c("CHD", "Cancer", "Vitamin D recommendation", "Child's IQ"),
# Result = c(rep(5000000, 3), 50000) / get("N", envir = openv) # Other impacts apply to the population of Finland,
# # maternal impacts apply to 50000 infants per year.
#))
daly <- indivrisk * dweight * pop
cat("DALYjen määrä koko Suomen väestössä per vuosi\n")
oprint(summary(daly))
signer <- Ovariable(data = data.frame(Exposcen = c("BAU", "No exposure"), Result = c(1, -1)))
dalys <- unkeep(daly * signer, prevresults = TRUE, sources = TRUE)
dalys <- oapply(dalys, cols = "Exposcen", FUN = sum)
dalys@output <- dalys@output[!is.na(result(dalys)) , ]
ggplot(daly@output, aes(x = Exposcen, weight = Result, fill = Trait)) + geom_bar()
ggplot(dalys@output, aes(x = Age, weight = Result, fill = Trait)) + geom_bar()
ggplot(exposure@output, aes(x = Ikä, y = Result, colour = Sukupuoli)) + geom_point() +
facet_wrap(~ Exposure_agent, scales = "free_y")
#oprint(indivrisk)
#oprint(daly)
#oprint(dalys)
#dalyne@output <- dalyne@output[dalyne@output$Exposcen == "No exposure" ,]
#dalyne <- unkeep(dalyne, cols = "Exposcen")
#dalye <- daly - dalyne
#dalye@output <- dalye@output[dalye@output$Exposcen != "No exposure" ,]
#dalye <- unkeep(dalye, cols = "Exposcen")
#ggplot(daly@output, aes(x = Trait, weight = Result, fill = Response_metric)) + geom_bar()
#ggplot(dalys@output, aes(x = Trait, weight = Result, fill = Response_metric)) + geom_bar()
#ggplot(daly@output, aes(x = Silakkamäärä, y = Result, colour = Response_metric)) + geom_point()
#dalys <- oapply(daly, cols = c("Trait", "Response_metric", "Exposure_route", "Exposure_metric", "Exposure_unit", "totcasesSource", "dweightSource"), FUN = sum)
########### Diskontatut syövät EI KÄYTETÄ MIHINKÄÄN MUTTA TALLENNETAAN TOISTAISEKSI ARKISTOINTIMIELESSÄ.
# Päätelmä koko tästä jaksosta on, että en oikeastaan ymmärrä mitä syöville pitäisi tehdä. Niinpä taidan
# yksinkertaisesti skipata koko ikäkysymyksen, olettaa samansuuruisen syöpäriskin ja elinajan menetyksen kaikille.
# Syöpä on joka tapauksessa niin pieni, että isokaan virhe ei tuota virhettä pääselmiin.
#ikäero <- Ovariable("ikäero", data = data.frame(
# Ikä = c("0", "1-4", "5-9", "10-14", "15-19", "20-24", "25-29", "30-34", "35-39", "40-44", "45-49", "50-54",
# "55-59", "60-64", "65-69", "70-74", "75-79", "80-84", "85-89", "90-94", "95-"),
# Result = 1
#))
#syöpäsum <- oapply(syyt, cols = c("Age"), FUN = sum)
#syöpäsum@output <- syöpäsum@output[syöpäsum@output$Trait == "Syövät" , ]
#disk <- (ikäero * syyt)@output
#disk$Ikäv <- as.numeric(gsub("-[0-9]*", "", disk$Ikä))
#disk$Agev <- as.numeric(gsub("-[0-9]*", "", disk$Age))
#disk$Elinikä <- disk$Agev - disk$Ikäv # Elinikä nykyisestä iästä syöpäkuolemaan
#disk <- disk[disk$Trait == "Syövät" & disk$Elinikä >= 0 , ]
#disk$Diskontto <- 1/(1.03)^disk$Elinikä # Kuolinhetken diskonttokorko
#disk$Pd <- (80 - disk$Agev) * disk$Diskontto # Diskontatut meneteyt elinvuodet per kuolema olettaen 80 odotusikä
#disk$Yd <- disk$Pd * disk$Result / 0.03 # Elinvuosimenetys kertaa todennäköisyys kuolla tietyn ikäisenä
# Odotettavissa olevien menetettyjen elinvuosien määrä ko. ikävaiheessa.
#as.data.frame(as.table(tapply(disk$Pd, disk["Ikä"], sum))) # Ikävaiheet yhteenlaskettuna nykyhetkeen eli DALY-vaikutus
# nykyisillä syöpämäärillä eri ikäryhmissä. Kaikki syövät yhteenlaskettuna BoD on 1-5 vuotta ikäryhmästä riippuen.
# Kuitenkin tämä on hyödytön tieto, koska syöpävaikutusta ei lasketa BoD:stä PAF:lla. Sen sijaan tarvittaisiin menettyjen
# diskontattujeen vuosien määrä ehdolla että tulee syöpä. Tätä varten ei tarvita syöpäilmaantuvuutta ollenkaan vaan
# ainoastaan ikäerosta riippuva diskonttaus. Eli seuraavasti:
#disk <- ikäero
#colnames(disk@data)[colnames(disk@data) == "Ikä"] <- "Age"
#disk <- (ikäero * disk)@output
#disk$Ikäv <- as.numeric(gsub("-[0-9]*", "", disk$Ikä))
#disk$Agev <- as.numeric(gsub("-[0-9]*", "", disk$Age))
#disk$Elinikä <- disk$Agev - disk$Ikäv # Elinikä nykyisestä iästä syöpäkuolemaan
#disk <- disk[disk$Elinikä >= 0 , ]
#disk$Diskontto <- 1/(1.03)^disk$Elinikä # Kuolinhetken diskonttokorko
#disk$Pd <- (80 - disk$Agev) * disk$Diskontto # Diskontatut meneteyt elinvuodet per kuolema olettaen 80 odotusikä
#disk$Yd <- disk$Pd * 5 # Odotettavissa olevien menetettyjen elinvuosien määrä ko. ikävaiheessa.
#as.data.frame(as.table(tapply(disk$Pd, disk["Ikä"], sum))) # Ikävaiheet yhteenlaskettuna nykyhetkeen eli DALY-vaikutus
| |
