+ Näytä koodi- Piilota koodi
N <- 100
erottelu1 <- 'Maakunta'
erottelu2 <- 'Causes'
verbose <- FALSE
rm(radonpit, disincidence, exposurenew, envir = openv)
library(OpasnetUtils)
library(ggplot2)
################################## YLEISET OSIOT
openv.setN(N)
objects.latest('Op_en6007', code_name = 'answer') # [[OpasnetUtils/Drafts]] fetches fillna function.
objects.latest('Op_en2261', code_name = 'initiate') # [[op_en:Health impact assessment]] ovariables totcases, attrcases.
objects.latest('Op_fi2761', code_name = 'alusta') # [[Talotyypit Suomessa]], ovariable asunnot, asuntovaesto.
objects.latest('Op_fi3907', code_name = 'alusta') # [[Suomen kunnat]], ovariable kunnat.
objects.latest('Op_fi2759', code_name = 'alusta') # [[Radon sisäilmassa]], ovariable radonpit, data.frame radon.
objects.latest('Op_fi3912', code_name = 'alusta') # [[Syöpäkuolleisuus Suomessa]], ovariable syopakuolleisuus.
######################## TAPAUSKOHTAISET SÄÄDÖT
####### PÄÄTÖKSET
#DecisionTableParser(opbase.data("Op_fi2760.paatokset"))
###### OSA 1: VÄESTÖ
### ASUNTOVÄESTÖ
# Yksinkertaista talotyyppiluokittelua ja summaa taloista pois perhekoko.
levels(asuntovaesto@data$Talotyyppi) <- ifelse(
levels(asuntovaesto@data$Talotyyppi) == "Asuinkerrostalo",
"Kerrostalo",
"Pientalo"
)
asuntovaesto <- oapply(EvalOutput(asuntovaesto), cols = "Asuntokunta", FUN = sum)
### TUPAKOINTI JA ALUEVASTAAVUUS
tupakointi <- Ovariable("tupakointi", # Tupakoivien osuus väestöstä. Result-sarake summautuu ykköseen.
data = data.frame(
Sukupuoli = rep(c("Miehet", "Naiset"), 2),
Tupakoija = rep(c("Kyllä", "Ei"), each = 2),
Result = c(0.11, 0.08, 0.39, 0.42)
)
)
aluevastaavuus <- Ovariable("aluevastaavuus",
data = {
dat <- tidy(opbase.data("Op_fi3907.maakunnat_ja_sairaanhoitopiirit"))
colnames(dat)[colnames(dat) == "Result"] <- "Sairaanhoitopiiri"
dat$Result <- 1
dat
}
)
### VÄESTÖN LASKENTA
population <- Ovariable("population",
dependencies = data.frame(Name = c("kunnat", "asuntovaesto", "aluevastaavuus", "tupakointi")),
formula = function(...) {
out <- kunnat * asuntovaesto * aluevastaavuus * tupakointi
return(out)
}
)
########## OSA 2: ALTISTUS
bkgrexposure <- 0 # Radonille voisi olla 5.
# RADON
radon$Pollutant <- "Radon" # Uusi sarake kuvaamaan sitä, minkä altisteen altistumisesta tai riskistä on kyse ko. rivillä.
# TUPAKKA. Väestö on jo jaettu tupakoiviin ja tupakoimattomiin.
tupakkaaltistus <- Ovariable("tupakkaaltistus",
data = data.frame(Tupakoija = c("Kyllä", "Ei"), Pollutant = "Tupakka", Result = c(1, 0))
)
# ALTISTUS YHTEENSÄ
radonpit <- EvalOutput(radonpit, N = N)
exposure <- Ovariable("exposure",
dependencies = data.frame(Name = c("tupakkaaltistus", "radonpit")),
formula = function(...) {
out <- orbind(tupakkaaltistus * 1, radonpit * 1)
out <- out[!colnames(out) %in% c("tupakkaaltistusSource", "radonpitSource")]
out <- fillna(out, c("Tupakoija", "Maakunta", "Talotyyppi")) # Tätä on muutettava, jos toiseen lisätään indeksejä.
return(out)
}
)
###### OSA 3: RISKISUHTEET
### SUHTEELLINEN RISKI
# Annosvasteiden perustiedot.
ERF <- Ovariable("ERF", data = data.frame(
Primaaripaikka = "Keuhkot, henkitorvi", # ICD.10.koodi == "C33-34"
Pollutant = c("Radon", "Tupakka"),
Result = c(1.0016, 20)
))
### SYÖPÄKUOLLEISUUS
disincidence <- syopakuolleisuus
disincidence@name <- "disincidence"
disincidence@data <- disincidence@data[disincidence@data$Primaaripaikka == "Keuhkot, henkitorvi" , ]
# ICD.10.koodi == "C33-34"
###################### ACTUAL MODEL
attrcases <- EvalOutput(attrcases)
casescen <- Ovariable("casescen",
dependencies = data.frame(Name = c(
"population",
"attrcases",
"ERF",
"exposurenew",
"bkgrexposure"
)),
formula = function(...) {
population2 <- unkeep(population, sources = TRUE, prevresults = TRUE)
disbkgr <- attrcases / population2
disbkgr@output <- disbkgr@output[!is.nan(result(disbkgr)) & ! grepl("\\+", disbkgr@output$Causes) , ]
diskbkgr <- unkeep(disbkgr, cols = c("Tupakoija", "Causes"))
RR <- unkeep(exp(log(ERF) * (exposurenew - bkgrexposure)), sources = TRUE, prevresults = TRUE)
RR <- oapply(RR, cols = "Pollutant", FUN = prod)
casescen <- disbkgr * population2 * RR
return(casescen)
}
)
exposurenew <- exposure * EvalOutput(Ovariable("Policy", data = data.frame(
Radonleikkaus = c("BAU", "Yli 400 leikataan", "Yli 200 leikataan"),
Result = 1
)))
result(exposurenew) <- ifelse(
exposurenew@output$Radonleikkaus == "Yli 400 leikataan" &
exposurenew@output$Pollutant == "Radon" &
result(exposurenew) > 400,
400,
result(exposurenew)
)
result(exposurenew) <- ifelse(
exposurenew@output$Radonleikkaus == "Yli 200 leikataan" &
exposurenew@output$Pollutant == "Radon" &
result(exposurenew) > 200,
200,
result(exposurenew)
)
exposurenew <- exposurenew * EvalOutput(Ovariable("Policy2", data = data.frame(
Tupakkarajoitus = c("BAU", "5 % vähenemä tupakoinnissa"),
Result = 1
)))
result(exposurenew) <- ifelse(
exposurenew@output$Tupakkarajoitus == "5 % vähenemä tupakoinnissa" &
exposurenew@output$Pollutant == "Tupakka",
result(exposurenew) * 0.95,
result(exposurenew)
)
ggplot(exposurenew@output, aes(x = Radonleikkaus, y = Result)) + geom_boxplot() + labs(y = "Radonpitoisuus (Bq/m3)")
casescen <- EvalOutput(casescen)
#ggplot(disbkgr@output, aes(x = Result, colour = Sukupuoli)) + geom_density()
#ggplot(disbkgr@output, aes(x = Result, colour = Talotyyppi)) + geom_density()
ggplot(casescen@output, aes(x = Tupakkarajoitus, weight = casescenResult, fill = Radonleikkaus)) + geom_bar(position = "dodge") + labs(y = "Keuhkosyöpäkuolemia (Lung ca deaths) (#/a)")
#ggplot(casescen@output, aes(x = Tupakoija, weight = Result, fill = Radonleikkaus)) + geom_bar()
########### TÄMÄ POIS TAI KORJATTAVA
#result(exposurenew) <- ifelse(
# exposurenew@output$Tupakkarajoitus == "5 % tupakkavähenemä radonpientaloissa" &
# exposurenew@output$Pollutant == "Tupakka" &
# exposurenew@output$Talotyyppi == "Pientalo" &
# exposurenew@output$radonpit > 100,
# result(exposurenew) * 0.95,
# result(exposurenew)
#)
####################### OUTPUT GRAPHS AND TABLES
oprint(head(attrcases@output))
ggplot(attrcases@output, aes(weight = attrcasesResult / N, x = Causes, fill = Talotyyppi)) + geom_bar(position = 'stack') + labs(y = "Keuhkosyöpäkuolemia (lung cancer deaths) (#/a)")
ggplot(attrcases@output, aes(weight = attrcasesResult / N, x = Causes, fill = Radontoimi)) + geom_bar(position = 'dodge') + labs(y = "Keuhkosyöpäkuolemia (lung cancer deaths) (#/a)")
ggplot(attrcases@output, aes(weight = attrcasesResult / N, x = Causes, fill = Tupakoija)) +
geom_bar(position = 'stack') + labs(y = "Keuhkosyöpäkuolemia (lung cancer deaths) (#/a)")
ggplot(attrcases@output, aes(weight = attrcasesResult / N, x = Causes, fill = Maakunta)) +
geom_bar(position = 'stack') + labs(y = "Keuhkosyöpäkuolemia (lung cancer deaths) (#/a)")
ggplot(attrcases@output, aes(weight = attrcasesResult / N, x = Tupakoija, fill = Talotyyppi)) + geom_bar(position = "dodge") +
theme_grey(base_size = 24) + labs(title = "Tapauksia vuodessa")
ggplot(population@output, aes(x = Talotyyppi, weight = populationResult, fill = Tupakoija)) + geom_bar() +
theme_grey(base_size = 24) + labs(title = "Väestö")
ggplot(exposure@output[exposure@output$Pollutant == "Radon" , ], aes(x = exposureResult, colour = Talotyyppi)) + geom_density(adjust = 4, size = 1) + labs(title = "Radonpitoisuuden jakauma talotyypeittäin")
ggplot(radonpit@output, aes(x = Maakunta, y = radonpitResult, fill = Talotyyppi)) + geom_boxplot() + coord_flip(ylim = c(0, 1000)) + labs(y = "Radonpitoisuus (Radon concentration) (Bq/m3)")
if(verbose) {
cat("Väestö\n")
print(nrow(population@output))
oprint(head(population@output))
cat("Kunnat\n")
print(nrow(kunnat@output))
oprint(head(kunnat@output))
cat("Asuntoväestö\n")
print(nrow(asuntovaesto@output))
oprint(head(asuntovaesto@output))
cat("Aluevastaavuus\n")
print(nrow(aluevastaavuus@output))
oprint(head(aluevastaavuus@output))
cat("Tupakointi\n")
print(nrow(tupakointi@output))
oprint(head(tupakointi@output))
cat("Altistus\n")
print(nrow(exposure@output))
oprint(head(exposure@output))
cat("ERF\n")
print(nrow(ERF@output))
oprint(head(ERF@output), digits = 4)
cat("RR\n")
print(nrow(RR@output))
oprint(head(RR@output))
cat("Syöpäkuolleisuus\n")
print(nrow(disincidence@output))
oprint(head(disincidence@output), digits = 6)
cat("Kokonaisvaikutus PAF-menetelmällä\n")
print(nrow(attrcases@output))
oprint(head(attrcases@output))
}
## POISTETTIIN KOKO KUNTAKARTTATOIMINNALLISUUS, KOSKA JOTKIN OSAT OVAT VANHENTUNEET. ON HELPOMPI RAKENTAA ALUSTA UUDESTAAN.
####################################################3
attrisk <- Ovariable("attrcases", # Probabilities of disease attributed to specific (combinations of) causal exposures.
dependencies = data.frame(Name = c(
"ERF", # Exposure-response function
"exposure", # Total exposure to an agent or pollutant
"bkgrexposure" # Background exposure to an agent (a level below which you cannot get in practice)
)),
formula = function(...) {
# First calculate risk ratio and remove redundant columns because they cause harm when operated with itself.
RR <- unkeep(exp(log(ERF) * (exposure - bkgrexposure)), prevresults = TRUE, sources = TRUE)
PAF <- (RR - 1) / RR
# pollutants is a vector of pollutants considered.
pollutants <- unique(exposure@output$Pollutant)
pollutants <- levels(pollutants)[pollutants]
out <- 1
for(i in 1:length(pollutants)) {
# Attributable fraction of a particular pollutant is combined with all pollutant AFs.
# The combination has 2^n rows (n = number of pollutants). Pollutant is either + or - depending on
# whether it caused the disease or not.
temp <- Ovariable(data = data.frame(
Pollutant = pollutants[i],
Temp1 = c(paste(pollutants[i], "-", sep = ""), paste(pollutants[i], "+", sep = "")),
Result = c(-1, 1) # Non-causes are temporarily marked with negative numbers.
))
temp <- temp * PAF
# Non-causes are given the remainder (1-AF) of attributable fraction AF.
result(temp) <- ifelse(result(temp) > 0, result(temp), 1 + result(temp))
# Causes with 0 AF are marked 1. This must be corrected.
result(temp) <- ifelse(result(temp) == 1 & grepl("\\+", temp@output$Temp1), 0, result(temp))
out <- out * temp
out <- unkeep(out, cols = "Pollutant", sources = TRUE, prevresults = TRUE)
# Combine and rename columns.
if(i == 1) {
colnames(out@output)[colnames(out@output) == "Temp1"] <- "Causes"
} else {
out@output$Causes <- paste(out@output$Causes, out@output$Temp1)
out@output$Temp1 <- NULL
}
}
return(out)
}
)
############################
totrisk <- Ovariable("totrisk", # This calculates the annual risk of the disease in each population subgroup.
# The cases are NOT attributed to specific (combinations of) causes.
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
"exposure", # Exposure to the pollutants
"bkgrexposure"
)), # Background exposure (a level below which you cannot get in practice)
formula = function(...) {
# Calcualte the risk ratio to each subgroup based on the exposure in that subgroup.
# Combine pollutant-specific RRs by multiplying.
RR <- unkeep(exp(log(ERF) * (exposure - bkgrexposure)), prevresults = TRUE, sources = TRUE)
RR <- oapply(RR, cols = "Pollutant", 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 <- population * RR # unkeep(population * RR, prevresults = TRUE, sources = TRUE)
pci <- pci / oapply(unkeep(pci, prevresults = TRUE, sources = TRUE), cols = takeout, FUN = sum)
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 * pci
# out <- unkeep(out, sources = TRUE, prevresults = TRUE)
return(out)
}
)
rad <- exposure * EvalOutput(Ovariable(data = data.frame(Pollutant = "Radon" , Result = 1)))
tup <- exposure * EvalOutput(Ovariable(data = data.frame(Pollutant = "Tupakka", Result = 1)))
tup <- unkeep(tup, cols = c("Pollutant", "exposureResult", "radonpitResult"))
rad <- unkeep(rad, cols = c("Pollutant", "exposureResult", "tupakkaaltistusResult"))
temp <- totrisk * tup * rad * 0
temp@output$tupakkaaltistusResult <- as.factor(temp@output$tupakkaaltistusResult)
ggplot(temp@output, aes(x = radonpitResult, y = totriskResult, colour = tupakkaaltistusResult)) + geom_point() + scale_x_log10() + scale_y_log10()
#### hirisk on alaryhmä, jolla henkilökohtainen riski on erityisen suuri. Huom1 Tämä EI OLE väestöpainotettu otos,
# vaan perusyksikkönä on väestöryhmä. Ilmeisesti tämän takia painottuvat pienet alueet, jossa satunnaistaminen
# voi saada aikaan isoja riskejä.
hirisk <- totrisk
hirisk@output <- hirisk@output[result(hirisk) > 1E-5 , ]
nrow(totrisk@output)
nrow(hirisk@output)
summary(hirisk@output)
ggplot(totrisk@output, aes(x = Maakunta, y = totriskResult * 1E6, fill = Tupakoija)) + geom_boxplot(point = 1) + labs(title = "Keuhkosyöpäriski maakunnittain", y = "Kuoleman riski (Risk of death) (micromort/a)") + coord_flip(ylim = c(-2, 50))
| |
Näkökulmia radonhaitan torjuntaan