+ Näytä koodi- Piilota koodi
library(sorvi)
library(OpasnetUtils)
library(ggplot2)
objects.latest('Op_en6007', code_name = 'answer') # [[OpasnetUtils/Drafts]] fetches fillna function.
N <- 2
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.
objects.latest('Op_fi2761', code_name = 'alusta') # [[Talotyypit Suomessa]], ovariable talot.
# Summaa taloista pois perhekoko ja muodosta uusi ovariable asukkaat.
levels(talot@output$Talotyyppi) <- ifelse(
levels(talot@output$Talotyyppi) == "Asuinkerrostalo",
"Kerrostalo",
"Pientalo"
)
talot@output <- dropall(talot@output)
dat <- as.data.frame(as.table(tapply(talot@output$Result, talot@output[c("Kunta", "Talotyyppi")], sum)))
colnames(dat)[colnames(dat) == "Freq"] <- "Result"
asukkaat <- EvalOutput(Ovariable("asukkaat", data = dat), N = N)
syopakuolleisuus@data <- syopakuolleisuus@data[syopakuolleisuus@data$Primaaripaikka == "Keuhkot, henkitorvi" , ] # ICD.10.koodi == "C33-34"
radonpit <- EvalOutput(radonpit, N = N)
kunnat <- EvalOutput(kunnat, N = N)
syopakuolleisuus <- EvalOutput(syopakuolleisuus, 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)
tupakkaaltistus <- EvalOutput(Ovariable("tupakkaaltistus",
data = data.frame(Tupakka = c("Kyllä", "Ei"), Result = c(1, 0))
), N = N)
dat <- tidy(opbase.data("Op_fi2760.maakunnat_ja_sairaanhoitopiirit"))
colnames(dat)[colnames(dat) == "Result"] <- "Sairaanhoitopiiri"
dat$Result <- 1
aluevastaavuus <- EvalOutput(Ovariable("aluevastaavuus", data = dat))
väestö <- kunnat * asukkaat * aluevastaavuus
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)
oprint(summary(radonalt, marginals = c("Maakunta", "Talotyyppi")))
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))
väestöosuus <- 1
PAF <- väestöosuus * (RR - 1) / (väestöosuus * (RR - 1) + 1)
PAF@output <- fillna(PAF@output, c("Tupakka", "Talotyyppi"))
#taustariski@output <- taustariski@output[taustariski@output$Tupakka == "Kyllä" , colnames(taustariski@output) != "Tupakka"]
vaikutus <- väestö * syopakuolleisuus * tupakointi * PAF /100000
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()
kokonaisvaikutus <- exp(log(rr.tupakka) * tupakkaaltistus + log(rr.radon) * radonpit)
radonvaikutus <- kokonaisvaikutus - exp(log(rr.tupakka) * tupakkaaltistus)
radonvaikutus <- radonvaikutus * syopakuolleisuus /100000 * väestö
kokonaisvaikutus <- kokonaisvaikutus * syopakuolleisuus /100000 * väestö
oprint(out[1:20, ])
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)
| |