+ Näytä koodi- Piilota koodi
library(OpasnetUtils)
library(ggplot2)
objects.latest('Op_en6007', code_name = 'answer') # [[OpasnetUtils/Drafts]] fetches fillna and collapsemarg function.
###### OSA 1: VÄESTÖ
### ASUNTOVÄESTÖ
objects.latest('Op_fi2761', code_name = 'alusta') # [[Talotyypit Suomessa]], ovariable talot.
# Yksinkertaista talotyyppiluokittelua ja summaa taloista pois perhekoko.
levels(asuntoväestö@data$Talotyyppi) <- ifelse(
levels(asuntoväestö@data$Talotyyppi) == "Asuinkerrostalo",
"Kerrostalo",
"Pientalo"
)
### KUNNAT
objects.latest('Op_fi3907', code_name = 'alusta') # [[Suomen kunnat]], ovariable kunnat.
### 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),
Tupakka = 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_fi2760.maakunnat_ja_sairaanhoitopiirit"))
colnames(dat)[colnames(dat) == "Result"] <- "Sairaanhoitopiiri"
dat$Result <- 1
dat
}
)
### VÄESTÖN LASKENTA
väestö <- Ovariable("väestö",
dependencies = data.frame(Name = c("kunnat", "asuntoväestö", "aluevastaavuus", "tupakointi")),
formula = function(...) {
asuntoväestö <- collapsemarg(asuntoväestö, c("Asuntokunta"), "sum")
out <- kunnat * asuntoväestö * aluevastaavuus * tupakointi
return(out)
}
)
###### OSA 2: RISKISUHTEET
### ALTISTUS
# RADON
objects.latest('Op_fi2759', code_name = 'alusta') # [[Radon sisäilmassa]], ovariable radonpit, data.frame radon.
radon$Altiste <- "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(Tupakka = c("Kyllä", "Ei"), Altiste = "Tupakka", Result = c(1, 0))
)
# ALTISTUS YHTEENSÄ
altistus <- Ovariable("altistus",
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("Tupakka", "Maakunta", "Talotyyppi")) # Tätä on muutettava, jos toiseen lisätään indeksejä.
return(out)
}
)
### SUHTEELLINEN RISKI
# Annosvasteiden perustiedot.
ERF <- Ovariable("ERF", data = data.frame(
Altiste = c("Radon", "Tupakka"),
Primaaripaikka = "Keuhkot, henkitorvi", # ICD.10.koodi == "C33-34"
Result = c(1.0016, 20)
))
tausta.altistus <- 0 # Radonille voisi olla 5.
RR <- Ovariable("RR",
dependencies = data.frame(Name = c("ERF", "altistus", "tausta.altistus")),
formula = function(...) {
out <- exp(log(ERF) * (altistus - tausta.altistus)) # Relative risk given the exposures
return(out)
}
)
### SYÖPÄKUOLLEISUUS
objects.latest('Op_fi3912', code_name = 'alusta') # [[Syöpäkuolleisuus Suomessa]], ovariable syopakuolleisuus.
# syopakuolleisuus@data <- syopakuolleisuus@data[syopakuolleisuus@data$Primaaripaikka == "Keuhkot, henkitorvi" , ] # ICD.10.koodi == "C33-34"
#syopakuolleisuus <- EvalOutput(syopakuolleisuus, N = N)
###### OSA 3: TAUTIKUORMA
### ARVIO PERUSTUEN POPULATION ATTRIBUTABLE FRACTIONIIN
# Tässä oletetaan, että riskisuhteiden avulla voidaan suoraan laskea eri teijöiden osuus tautikuormasta.
# Luku on yliarvio, koska jokainen altiste lasketaan ikään kuin muita altisteita ei olisi tautia aiheuttamassa.
# Etuna on, että jokaiselle altisteelle saadaan arvio erikseen.
väestöosuus <- 1
vaikutuspaf <- Ovariable("vaikutuspaf",
dependencies = data.frame(Name = c("väestöosuus", "RR", "väestö", "syopakuolleisuus")),
formula = function(...) {
PAF <- väestöosuus * (RR - 1) / (väestöosuus * (RR - 1) + 1)
out <- väestö * PAF * syopakuolleisuus
return(out)
}
)
### ARVIO PERUSTUEN TAUSTARISKIIN
# Tässä oletetaan, että on olemassa populaatiokohtainen taustariski, jota ei altistumisen takia havaita.
# Se voidaan kuitenkin laskea ottamalla havaittu tautikuorma ja poistamalla altistumisen vaikutus.
# Arvio lienee tarkempi kuin PAF-menetelmällä, mutta altistekohtainen tieto menetetään.
# Calculate subgroup-specific background risks by first calculating a temporary burden estimate.
# temp is a population-weighted average of relative risks. The actual disease risk is divided by this value and then multiplied
# by the subgroup-specific relative risk. In this way, the total burden of all subgroups equals to what is actually seen.
vaikutusrr <- Ovariable("vaikutusrr",
dependencies = data.frame(Name = c("väestö", "RR")),
formula = function(...) {
RRkok <- collapsemarg(RR, c("Altiste"), "prod")
temp1 <- RRkok * väestö # Population-weighted sum of the relative risk.
temp2 <- (RRkok * 0 + 1) * väestö # Population-weighted sum of ones with the same dimensions.
integrate <- colnames(temp1)[
colnames(temp1) %in% temp1@marginal & !
colnames(temp1) %in% c("Sairaanhoitopiiri", "Sukupuoli", "Iter")
]
#c("Tupakka", "Maakunta", "Talotyyppi", "Primaaripaikka", "ERFSource", "altistusSource", "RRSource", "Kunta",
# "kunnatSource", "talotSource", "väkiSource", "aluevastaavuusSource", "Sukupuoli", "tupakointiSource", "väestöSource"
#)
temp1 <- collapsemarg(temp1, integrate, "sum") # Aggregate to the same indices as kokonaisvaikutus.
temp2 <- collapsemarg(temp2, integrate, "sum")
tausta <- syopakuolleisuus / (temp1 / temp2) # syopakuolleisuus without any exposure
cat("väestö * tausta, väestö * RRkok, tausta * RRkok\n")
print(nrow((väestö * tausta)@output))
print(nrow((väestö * RRkok)@output))
print(nrow((tausta * RRkok)@output))
print(ograph(RRkok, x = "Maakunta", fill = "Tupakka"))
print(ograph(tausta, x = "Maakunta", fill = "Sukupuoli"))
out <- väestö * tausta * RRkok
# out <- collapsemarg(out, c("Iter", "Kunta", "Tupakka", "Talotyyppi"), "sum")
return(out)
}
)
#turhat <- c("Maakunta", "kunnatSource", "talotSource", "väkiSource", "Sairaanhoitopiiri", "aluevastaavuusSource",
# "tupakointiSource", "väestöSource", "ERFSource")
väestö <- EvalOutput(väestö, N = N)
#väestö <- collapsemarg(väestö, "Maakunta", fun = "pick", picks = "Itä-Uusimaa")
vaikutuspaf <- EvalOutput(vaikutuspaf, N = N)
vaikutusrr <- EvalOutput(vaikutusrr, N = N)
if(verbose) {
cat("Väestö\n")
print(nrow(väestö@output))
oprint(head(väestö@output))
cat("Kunnat\n")
print(nrow(kunnat@output))
oprint(head(kunnat@output))
cat("Asuntoväestö\n")
print(nrow(asuntoväestö@output))
oprint(head(asuntoväestö@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(altistus@output))
oprint(head(altistus@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(syopakuolleisuus@output))
oprint(head(syopakuolleisuus@output), digits = 6)
cat("Kokonaisvaikutus PAF-menetelmällä\n")
print(nrow(vaikutuspaf@output))
oprint(head(vaikutuspaf@output))
cat("Kokonaisvaikutus RR-menetelmällä\n")
print(nrow(vaikutusrr@output))
oprint(head(vaikutusrr@output))
# oprint(vaikutuspaf)
# oprint(vaikutusrr)
}
temp <- vaikutuspaf@output
temp$vaik <- temp$vaikutuspafResult / N
ggplot(temp, aes(weight = vaik, x = Tupakka, fill = Talotyyppi)) + geom_bar(position = 'dodge')
ograph(RR, x = erottelu1, fill = erottelu2, title = 'Suhteellinen yksilöriski altistumattomaan verrattuna') + coord_trans(y = 'log2')
ograph(RR, x = "Tupakka", fill = "Altiste") + coord_trans(y = 'log2')
ggplot(vaikutuspaf@output, aes(weight = vaikutuspafResult, x = Tupakka, fill = Talotyyppi)) + geom_bar(position = "dodge") +
theme_grey(base_size = 24) +
labs( # label names
title = "Vaikutus PAF-menetelmällä"
)
ggplot(vaikutusrr@output, aes(weight = vaikutusrrResult, x = Tupakka, fill = Talotyyppi)) + geom_bar(position = "dodge") +
theme_grey(base_size = 24) +
labs( # label names
title = "Vaikutus RR-menetelmällä"
)
ggplot(väestö@output, aes(x = Talotyyppi, weight = väestöResult)) + geom_bar() + labs(title = "Väestö")
ggplot(altistus@output[altistus@output$Altiste == "Radon" , ], aes(x = altistusResult, fill = Talotyyppi)) + geom_density(alpha = 0.2)
## POISTETTIIN KOKO KUNTAKARTTATOIMINNALLISUUS, KOSKA JOTKIN OSAT OVAT VANHENTUNEET. ON HELPOMPI RAKENTAA ALUSTA UUDESTAAN.
| |