+ Näytä koodi- Piilota koodi
library(sorvi)
library(OpasnetUtils)
library(ggplot2)
collapsemarg <- function(object, marginals, fun = sum) {
object@output <- dropall(object@output)
temp <- as.data.frame(as.table(tapply(result(object), object@output[marginals], fun)))
colnames(temp)[colnames(temp) == "Freq"] <- "Result"
out <- EvalOutput(Ovariable(object@name, data = temp))
return(out)
}
objects.latest('Op_en6007', code_name = 'answer') # [[OpasnetUtils/Drafts]] fetches fillna function.
N <- 2
###### VÄESTÖ
objects.latest('Op_fi3907', code_name = 'alusta') # [[Suomen kunnat]], ovariable kunnat.
objects.latest('Op_fi2761', code_name = 'alusta') # [[Talotyypit Suomessa]], ovariable talot.
# Summaa taloista pois perhekoko ja muodosta uusi ovariable asukkaat.
talot@output <- talot@output[talot@output$Kunta == "Akaa" , ] ########### POISTA TÄMÄ
levels(talot@output$Talotyyppi) <- ifelse(
levels(talot@output$Talotyyppi) == "Asuinkerrostalo",
"Kerrostalo",
"Pientalo"
)
talot@output <- dropall(talot@output)
talot <- collapsemarg(talot, c("Kunta", "Talotyyppi"), sum)
#colnames(dat)[colnames(dat) == "Freq"] <- "Result"
#asukkaat <- EvalOutput(Ovariable("asukkaat", data = dat), N = N)
kunnat <- EvalOutput(kunnat, N = N)
tupakointi <- EvalOutput(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)
)
), N = N)
aluevastaavuus <- EvalOutput(Ovariable("aluevastaavuus",
data = {
dat <- tidy(opbase.data("Op_fi2760.maakunnat_ja_sairaanhoitopiirit"))
colnames(dat)[colnames(dat) == "Result"] <- "Sairaanhoitopiiri"
dat$Result <- 1
dat
}
))
väestö <- kunnat * talot * aluevastaavuus * tupakointi
###### altistus
objects.latest('Op_fi2759', code_name = 'alusta') # [[Radon sisäilmassa]], ovariable radonpit, data.frame radon.
radonpit <- EvalOutput(radonpit, N = N)
radonpit@output <- radonpit@output[radonpit@output$Maakunta == "Pirkanmaa" , ] ############# POISTA TÄMÄ
tupakkaaltistus <- EvalOutput(Ovariable("tupakkaaltistus",
data = data.frame(Tupakka = c("Kyllä", "Ei"), Result = c(1, 0))
), N = N)
altistus <- orbind(
tupakkaaltistus * EvalOutput(Ovariable(data = data.frame(Altiste = "Tupakka", Result = 1)), N = N),
radonpit * EvalOutput(Ovariable(data = data.frame(Altiste = "Radon", Result = 1)), N = N)
)
altistus <- EvalOutput(Ovariable("altistus", data = altistus), N = N)
altistus@output <- fillna(altistus@output, c("Tupakka", "Maakunta", "Talotyyppi"))
#colnames(temp)[colnames(temp) == "altistusResult"] <- "Result"
#temp <- EvalOutput(Ovariable("temp", data = temp))
######## 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)
####### Relative risk
ERF <- EvalOutput(Ovariable("ERF", data = data.frame(Altiste = c("Radon", "Tupakka"), Result = c(1.006, 20))), N = N)
tausta.altistus <- 0 # Radonille voisi olla 5.
RR <- exp(log(ERF) * (altistus - tausta.altistus)) # Relative risk given the exposures
####### taustariski
# Calculate subgroup-specific background risks by first calculating a temporary burden estimate.
# temp is a population-weighted sum 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.
kokonaisvaikutus <- syopakuolleisuus * väestö
temp <- collapsemarg(RR, c("Maakunta", "Iter", "Tupakka", "Talotyyppi"), prod)
temp <- temp * väestö
temp <- temp / collapsemarg(temp, c("Sairaanhoitopiiri", "Sukupuoli", "Iter"), sum)
tausta <- kokonaisvaikutus * temp
########### PAF
väestöosuus <- 1
RR <- exp(log(ERF) * (altistus - tausta.altistus))
PAF <- väestöosuus * (RR - 1) / (väestöosuus * (RR - 1) + 1)
#PAF@output <- fillna(PAF@output, c("Tupakka", "Talotyyppi", "Maakunta"))
######### vaikutus
vaikutus <- väestö * syopakuolleisuus * PAF
vaikutus@output <- vaikutus@output[!is.na(vaikutus@output$Result) , ]
vaikutusrr <- väestö * tausta * RR
ggplot(asukkaat@output, aes(x = asukkaatResult)) + geom_density()
ggplot(altistus@output, aes(x = altistusResult)) + geom_density()
ggplot(altistus@output, aes(x = altistusResult, fill = Altiste)) + geom_density()
ggplot(vaikutus@output, aes(weight = vaikutusResult, fill = Altiste)) + geom_bar()
ggplot(vaikutus@output, aes(weight = Result, x = Altiste)) + geom_bar()
ggplot(vaikutusrr@output, aes(weight = Result, x = Altiste)) + geom_bar()
print("Yksikköriski: tapauksia tuhatta asukasta kohti vuodessa")
temp <- as.data.frame(as.table(tapply(out$Radonvaikutus/out$Asukkaita, list(out$Maakunta, out$Tupakka, out$Talo), sum)/n))
temp <- temp[!is.na(temp$Freq), ]
#head(temp)
temp$Freq <- temp$Freq*1000
colnames(temp) <- c("Maakunta", "Tupakka", "Talotyyppi", "Riski")
oprint(temp)
print("Radonin aiheuttama lisäriski, tapausta vuodessa")
oprint(as.data.frame(as.table(((tapply(out$Radonvaikutus, list(out$Talo, out$Tupakka), sum)/n)))))
ggplot(out, aes(x = Maakunta, weight = Kokonaisvaikutus, fill = Tupakka)) + geom_bar(position = "stack") # position = "dodge"
out$Radonvaikutus <- ifelse(out$Radonvaikutus < 0, 0, out$Radonvaikutus)
out$Radonsata <- round(as.numeric(as.character(out$Radon)), -2)
ggplot(out, aes(x = Radon, weight = Radonvaikutus)) + geom_density(alpha = 1)
ggplot(out, aes(x = Radonsata, weight = Radonvaikutus, fill = Talo)) + geom_bar(position = "stack")
###########################################################
# (C) 2011 Leo Lahti <leo.lahti[at]iki.fi> All rights reserved.
# License: FreeBSD, http://en.wikipedia.org/wiki/BSD_licenses
# Tama esimerkki on testattu sorvi-paketin versiolla 0.1.23
# Esimerkki Suomen kuntatason vaestonkasvutilastojen (Tilastokeskus)
# visualisoinnista Maanmittauslaitoksen karttadatalla (vuonna 2010)
###############################################
# Lue Suomen kuntarajat SpatialPolygon-muodossa
# (C) Maanmittauslaitos 2011
# http://www.maanmittauslaitos.fi/aineistot-palvelut/digitaaliset-tuotteet/ilmaiset-aineistot/hankinta
#LoadData(MML)
#sp <- MML[["1_milj_Shape_etrs_shape"]][["kunta1_p"]]
sp <- LoadMML(data.id = "kunta1_p", resolution = "1_milj_Shape_etrs_shape")
#################################################
# Lue kuntatason vaestonkasvutiedot tilastokeskuksen StatFin-tietokannasta
# http://www.stat.fi/tup/statfin/index.html
# PC Axis-muodossa
px <- statfin.px("vrm/synt/080_synt_tau_203.px")
# statfin.px-funktio käyttää soRvin read.px-funktiota hakemaan dataa URLista, joka alkaa "http://pxweb2.stat.fi/database/StatFin/"
# URLin loppuosa pitää antaa statfin.px:lle parametrina.
# Poimi taulukosta halutut tiedot
# pxs <- subset(as.data.frame(px), Väestönmuutos.ja.väkiluku == "Luonnollinen väestönlisäys" & Vuosi == year)
# pxs <- subset(as.data.frame(px), Väestönmuutos.ja.väkiluku == parameter & Vuosi == year)
# Putsaa data
# vaestonkasvu <- preprocess.px(pxs)
# head(vaestonkasvu)
################################################
vaestonkasvu <- out[out$Talo == talo & out$Tupakka == tupakka, ]
vaestonkasvu <- as.data.frame(as.table(tapply(vaestonkasvu[, vaste], list(vaestonkasvu$Kunta), sum)/n))
colnames(vaestonkasvu) <- c("Kunta", "Vaste")
head(vaestonkasvu)
vaestonkasvu[is.na(vaestonkasvu)] <- 0
##############################################
# Lisaa tiedot karttaobjektiin
sp@data$vaestonkasvu <- vaestonkasvu$Vaste[match(sp$Kunta.FI, vaestonkasvu$Kunta)]
# Korvaa puuttuvat arvot nollalla
sp[["vaestonkasvu"]][is.na(sp[["vaestonkasvu"]])] <- 0
#colnames(out)[vaste]
#vaestonkasvu$Vaste
################################################
# Maarittele varipaletti
my.palette <- colorRampPalette(c("white", "red"), space = "rgb") # , "blue"
ncol <- 40 # Number of colors
#################################################
# Piirra kuva
varname <- "vaestonkasvu"
int <- max(abs(vaestonkasvu$Vaste))
int
q <- spplot(sp, varname,
col.regions = my.palette(ncol),
main = paste(colnames(out)[vaste], talo, "Tupakointi: ", tupakka, sep = ", "),
colorkey = TRUE,
lwd = .4,
col = "black",
at = seq(0, 0 + int, length = ncol)
)
print(q)
| |