+ Show code- Hide code
############## TILAPÄISESTI HIA-FUNKTIOITA TÄHÄN
library(OpasnetUtils)
dose <- Ovariable("dose", # This calculates the body-weight-scaled exposure or "dose" to be used with ERFs.
dependencies = data.frame(Name = c(
"ERF", # Exposure-response function of the pollutants or agents (RR for unit exposure)
"exposure", # Exposure to the pollutants
"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 <- unique(ERF@output[c("ERF_parameter", "Exposure_agent")])
scaling <- Ovariable("scaling", data = data.frame(
scaling,
Result = 1 * grepl("bw", scaling$ERF_parameter) # if ERF parameter contains bw, scaling is TRUE, i.e. 1.
))
out <- exposure / (((BW - 1) * scaling) + 1) # If scaling is 0, BW cancels out.
out <- out * Ovariable( # Create alternative scenario with zero exposure.
output = data.frame(Exposcen = c("BAU", "No exposure"), Result = c(1, 0)),
marginal = c(TRUE, FALSE)
)
return(out)
}
)
RR <- Ovariable("RR", # 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(
"ERF", # Exposure-response function of the pollutants or agents (RR for unit exposure)
"dose", # 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.
)),
formula = function(...) {
####################################################################
####### 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 <- list()
marginals <- character()
#First take the relative risk estimates
ERFrr <- ERF
ERFrr@output <- ERF@output[ERF@output$ERF_parameter %in% c("RR", "RR bw") , ]
if(nrow((ERFrr * dose)@output) > 0) {
bigger <- (threshold > bgexposure) * threshold + (1 - (threshold > bgexposure)) * bgexposure # Choose bigger
RRrr <- exp(log(ERFrr) * (dose - bigger)) # Actual function
RRrr <- (RRrr - 1) * (dose > threshold) + 1 # RR is 1 below threshold
test <- c(test, RRrr)
marginals <- c(marginals, colnames(RRrr@output)[RRrr@marginal])
}
# 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 * dose)@output) > 0 & nrow(ED50@output) > 0) {
RRrhill <- 1 + (dose * ERF) / (dose + ED50) # ERF has parameter value for Imax. If Imax < 0, risk reduces.
test <- c(test, RRrhill)
marginals <- c(marginals, colnames(RRrhill@output)[RRrhill@marginal])
}
# 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 * dose)@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?
# RRor <- (RRrr - 1) * (dose > threshold) + 1 # RR is 1 below threshold
# test <- c(test, RRor)
# marginals <- c(marginals, colnames(OR@output)[OR@marginal]
# }
if(length(test) == 0) return(data.frame())
if(length(test) == 1) out <- test[[1]]@output
if(length(test) == 2) out <- orbind(test[[1]], test[[2]])
if(length(test) == 3) out <- orbind(orbind(test[[1]]), test[[2]], test[3])
out <- Ovariable(output = out, marginal = colnames(out) %in% marginals)
# Dilute the risk in the population if not all are exposed i.e. frexposed < 1.
out <- frexposed * (out - 1) + 1
out <- unkeep(out, prevresults = TRUE, sources = TRUE)
out <- oapply(out, cols = "Exposure_agent", FUN = prod)
return(out)
}
)
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
"RR", # Relative risks for the given exposure
"ERF", # Other ERFs than those that are relative to background.
"dose", # Exposure to the pollutants
"frexposed", # fraction of population that is exposed
"threshold" # exposure level below which the agent has no impact.
)),
formula = function(...) {
test <- list()
marginals <- character()
############### First look at the relative risks based on RR
if(nrow((RR * dose)@output) > 0) {
# if(class(population) != "ovariable") population <- EvalOutput(Ovariable("population", data = data.frame(Result = population)))
# # Remove redundant columns and locations.
# 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 in population but NOT in the disease incidence.
# However, populationSource is kept because oapply does not run if there are no indices.
if(class(population) == "ovariable") {
takeout <- setdiff(colnames(population@output), c(colnames(disincidence@output), "populationSource"))
pop <- oapply(population, cols = takeout, FUN = sum) # Aggregate to larger subgroups.
} else {
takeout <- character()
pop <- population
}
# 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
# 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", prevresults = TRUE, sources = TRUE)
pci <- pci / oapply(temp, cols = takeout, FUN = sum)
# The cases are divided into smaller subgroups based on weights in pci.
# This is why the larger groups of population are used (pop instead of population).
out1 <- disincidence * pop * unkeep(pci, prevresults = TRUE, sources = TRUE)
out1 <- unkeep(out1, cols = "populationResult") # populationResult comes from pop and not from pci that actually contains
# the population weighting for takeout indices. Therefore it would be confusing to leave it there.
test <- c(test, out1)
marginals <- c(marginals, colnames(out1@output)[out1@marginal])
}
##########################################################################
############# This part is about absolute risks (i.e., risk is not affected by background rates).
# Unit risk (UR), cancer slope factor (CSF), and Exposure-response slope (ERS) estimates.
UR <- ERF
UR@output <- UR@output[UR@output$ERF_parameter %in% c("UR", "CSF", "ERS", "UR bw", "CSF bw", "ERS bw") , ]
if(nrow((UR * dose)@output) > 0) {
UR <- threshold + UR * dose * 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
UR <- population * UR
test <- c(test, UR)
marginals <- c(marginals, colnames(UR)[UR@marginal])
}
# Step estimates: value is 1 below threshold and above ERF, and 0 in between.
# frexposed cannot be used with Step because this may be used at individual and maybe at population level.
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 * dose)@output) > 0) {
Step <- 1 - (dose >= threshold) * (dose <= Step) # Actual equation
Step <- unkeep(Step, prevresults = TRUE, sources = TRUE)
Step <- oapply(Step, cols = "Exposure_agent", FUN = sum)
# out3 <- (population * 0 + 1) * Step # This is to maintain the ovariable structure # Does not work because population may be numeric.
test <- c(test, Step)
marginals <- c(marginals, colnames(Step@output)[Step@marginal])
}
#####################################################################
# Combining effects
if(length(test) == 0) return(data.frame())
if(length(test) == 1) out <- test[[1]]
if(length(test) == 2) out <- orbind(test[[1]], test[[2]])
if(length(test) == 3) out <- orbind(orbind(test[[1]], test[[2]]), test[[3]])
out <- Ovariable(output = out, marginal = colnames(out) %in% marginals)
return(out)
}
)
################################## Varsinainen malli alkaa tästä
library(OpasnetUtils)
library(ggplot2)
library(rgdal)
library(maptools)
library(RColorBrewer)
library(classInt)
library(OpasnetUtilsExt)
library(RgoogleMaps)
language <- "EN"
openv.setN(1)
objects.latest("Op_en6007", code_name = "answer") # [[OpasnetUtils/Drafts]] findrest
objects.latest("Op_en6289", code_name = "initiate") # [[Building model]] buildings, heatingEnergy. # Generic building model.
objects.latest("Op_en5417", code_name = "initiate") # [[Population of Kuopio]] population
objects.latest("Op_en5932", code_name = "initiate") # [[Building stock in Kuopio]] Building ovariables: buildingStock, renovation,
#renovationShares, construction, constructionAreas, buildingTypes, heatingShares, heatingSharesNew, eventyear, efficiencies
objects.latest("Op_en3435", code_name = "disperse") # [[Exposure to PM2.5 in Finland]] iF, emissionLocations
#Summarised Piltti matrix, another copy of the code on a more reasonable page
# Default run: en.opasnet.org/en-opwiki/index.php?title=Special:RTools&id=aXDIVDboftr1bTEd
objects.latest("Op_en2791", code_name = "initiate") # [[Emission factors for burning processes]] emissionFactors, fuelTypes
objects.latest("Op_en5488", code_name = "initiate") # [[Energy use of buildings]] energyUse, efficienciesNew, savingPotential,
#objects.latest('Op_en2261', code_name = 'initiate') # [[Health impact assessment]] ovariables totcases, AF, totcases2.
objects.latest('Op_en5917', code_name = 'initiate') # [[Disease risk]] disincidence
objects.latest('Op_en5827', code_name = 'initiate') # [[ERFs of environmental pollutants]] ERF, threshold
#objects.latest('Op_en5453', code_name = 'initiate') # [[Burden of disease in Finland]] BoD
frexposed <- 1 # fraction of population that is exposed
bgexposure <- 0 # Background exposure to an agent (a level below which you cannot get in practice)
BW <- 70 # Body weight (is needed for RR calculations although it is irrelevant for PM2.5)
directs <- opbase.data("Op_en5461", subset = "Direct inputs")
colnames(iF@output)[colnames(iF@output) == "City.area"] <- "Emission_site"
colnames(iF@output)[colnames(iF@output) == "Emissionheight"] <- "Emission_height"
colnames(emissionLocations@data) <- gsub(" ", "_", colnames(emissionLocations@data))
# Calculate intake fractions for Kuopio districts
districts <- tidy(opbase.data("Op_en3435.kuopio_city_districts"), widecol = "Location") # [[Exposure to PM2.5 in Finland]]
colnames(districts) <- gsub("\\.", "_", colnames(districts))
districts <- Ovariable("districts", data = data.frame(districts, Result = 1))
if(language == "EN") deci <- "Decisions" else deci <- "Päätökset"
decisions <- opbase.data('Op_en5461', subset = deci) # [[Climate change policies and health in Kuopio]]
DecisionTableParser(decisions)
# Remove previous decisions, if any.
rm(
"buildings",
"buildingStock",
"buildingTypes",
"construction",
"efficiencies",
"efficienciesNew",
"energyUse",
"heatingShares",
"heatingSharesNew",
"renovation",
"renovationShares",
"fuelTypes",
"year",
envir = openv
)
fuelTypes <- CheckDecisions(EvalOutput(fuelTypes, verbose = TRUE))
#################### Manage the data before calculating
# The building stock is measured as m^2 floor area.
if(!exists("language")) language <- "FI"
obsyear <- (192:205) * 10
####### Remove columns and rows that are not needed from data. These must be done before EvalOutput.
# This is the population of Kuopio (i.e. population living in the building stock)
population <- EvalOutput(population, verbose = TRUE)
pop <- oapply(population, cols = "City_area", FUN = "sum") # Sum across city areas.
pop <- unkeep(pop * 1, prevresults = TRUE) # Remove populationResult column because it would cause trouble in merge.
pop <- population / pop
buildingStock@data <- buildingStock@data[
buildingStock@data$Observation == "AreaHR" ,
colnames(buildingStock@data) != "Observation"
]
buildingStock <- buildingStock * pop
construction@data <- construction@data[
construction@data$Observation == "Area" ,
colnames(construction@data) != "Observation"
]
construction <- construction * constructionAreas / 100 / 3 # Statistics are for three years (2010-2012)
energyUse@data <- energyUse@data[
energyUse@data[["Energy use"]] == "Heat" ,
colnames(energyUse@data) != "Energy use"
]
savingPotential@data <- savingPotential@data[
savingPotential@data$Observation == "Relative" ,
colnames(savingPotential@data) != "Observation"
]
savingPotential <- 1 - savingPotential / 100
# Note: If you add decisions to savingPotential, this formula must be updated.
# The data are not yet specific to construction year, so remove index:
heatingSharesNew@data <- heatingSharesNew@data[colnames(heatingSharesNew@data) != "Constructed"]
# Fill in Heating types and convert from % to fraction.
heatingSharesNew <- findrest(heatingSharesNew, cols = "Heating", total = 100)
result(heatingSharesNew) <- result(heatingSharesNew) / 100
heatingShares <- EvalOutput(heatingShares, verbose = TRUE)
result(heatingShares) <- result(heatingShares) / 100 # From % to fraction
# Fill in the rest of the data (the last emission class was omitted because it is determined by the total).
efficienciesNew <- findrest(efficienciesNew, cols = "Efficiency", total = 100)
result(efficienciesNew) <- result(efficienciesNew) / 100
# When renovations are done, which type are they?
renovationShares <- findrest(renovationShares, cols = "Renovation", total = 100)
result(renovationShares) <- result(renovationShares) / 100
# What fraction of buildings is renovated per year?
renovation <- EvalOutput(renovation, verbose = TRUE)
result(renovation) <- result(renovation) / 100 # From % to fraction
disincidence <- disincidence / 100000 # Change units from 1/100000py to 1/py.
###################### Actual model
buildings <- EvalOutput(buildings, verbose = TRUE)
####### Sort factors
buildings@output <- buildings@output[
buildings@output$EfficiencyPolicy == "BAU" ,
colnames(buildings@output) != "EfficiencyPolicy"
]
buildings@output$Renovation <- factor(
buildings@output$Renovation,
levels = c("None", "General", "Windows", "Technical systems", "Sheath reform"),
ordered = TRUE
)
buildings@output$Efficiency <- factor(
buildings@output$Efficiency,
levels = c("Old", "New", "Low-energy", "Passive"),
ordered = TRUE
)
buildings@output$RenovationPolicy <- factor(
buildings@output$RenovationPolicy,
levels = c("BAU", "Active renovation", "Effective renovation"),
ordered = TRUE
)
heatingEnergy <- EvalOutput(heatingEnergy, verbose = TRUE)
heatingEnergy <- unkeep(heatingEnergy, cols = "Building2", sources = TRUE, prevresults = TRUE)
if(exists("server")) # The following code only works from Opasnet server.
{
temp <- buildings * districts
temp@output <- temp@output[temp@output$Year == 2040 , ]
temp <- unkeep(temp, sources = TRUE, prevresults = TRUE)
temp@output <- dropall(temp@output)
temp <- oapply(temp, cols = c("Building", "Heating", "Efficiency", "Renovation"), FUN = "sum", na.rm = TRUE)
if(language == "EN") lege <- "Floor area" else lege <- "Rakennuspinta-ala"
MyRmap(
ova2spat(temp, coord = c("E", "N"), proj4string = "+init=epsg:3067"), # National Land Survey uses EPSG:3067 (ETRS-TM35FIN)
plotvar = "Result", legend_title = lege, numbins = 8, pch = 19, cex = 2
)
}
####### Calculate emissions
emis <- heatingEnergy
emis@output <- emis@output[emis@output$Year >= 1980 , ]
emis <- oapply(emis, cols = c("Efficiency", "Renovation", "Building"), FUN = "sum")
emis@output$Year <- as.numeric(as.character(emis@output$Year))
emis <- emis * fuelTypes * emissionFactors * 3.6 * 1E-9 # convert from kWh /a to MJ /a and mg to ton
emis <- unkeep(emis, sources = TRUE, prevresults = TRUE)
emis <- emis * emissionLocations
emis@output$Emission_site <- ifelse(
emis@output$Emission_site == "At site of consumption",
as.character(emis@output$City_area),
as.character(emis@output$Emission_site)
)
emis2 <- oapply(
emis,
cols = c("Heating", "Burner", "Fuel", "City_area", "Emission_site", "Emission_height", "emissionLocationsSource"),
FUN = sum
)
ggplot(emis@output, aes(x = Year, weight = Result, fill = Emission_site)) + geom_bar() + facet_grid( Pollutant ~ ., scales = "free_y")
### Use these population and iF values in health impact assessment. Why?
pop <- 5E+5
iF <- 1E-7
exposure <- Ovariable("exposure",
# emis is in ton /a
# iF = conc (g /m3) * pop (#) * BR (m3 /s) / emis (g /s) <=> conc = emis * iF / BR / pop # conc is the exposure concentration
dependencies = data.frame(Name = c("emis2", "iF", "pop")),
formula = function(...) {
BR <- 20 # Nominal breathing rate (m^3 /d)
BR <- BR / 24 / 3600 # m^3 /s
out <- 1E+12 / 365 / 24 / 3600 # Emission scaling from ton /a to ug /s.
out <- (emis2 * out) * iF / BR / pop # the actual equation
out <- unkeep(out, prevresults = TRUE, sources = TRUE)
return(out)
}
)
exposure <- EvalOutput(exposure, verbose = TRUE)
colnames(exposure@output)[colnames(exposure@output) == "Pollutant"] <- "Exposure_agent"
ggplot(exposure@output, aes(x = Year, weight = exposureResult)) + geom_bar() + facet_grid(Exposure_agent ~ FuelPolicy, scales = "free_y")
#totcases <- totcases2 # This is temporary, until totcases2 is established and replaces totcases.
totcases <- EvalOutput(totcases)
DW <- Ovariable("DW", data = data.frame(directs[c("Exposure.agent", "Disease")], Result = directs$DW))
L <- Ovariable("L", data = data.frame(directs[c("Exposure.agent", "Disease")], Result = directs$L))
#ggplot(cases@output[cases@output$Causes == "PM2.5+" , ], aes(x = Year, weight = Result/10, fill = Disease)) + geom_bar() + facet_grid(FuelPolicy ~ RenovationPolicy, scale = "free_y")
DALYs <- cases * DW * L
DALYs <- oapply(DALYs, cols = c("Disease"), FUN = sum)
ggplot(DALYs@output[DALYs@output$Causes == "PM2.5+" , ], aes(x = Year, weight = Result/10, fill = Response.metric)) + geom_bar() + facet_grid(FuelPolicy ~ RenovationPolicy, scale = "free_y")
############################# Plot graphs
######## Buildings in Kuopio on map
if(exists("server")) # The following code only works from Opasnet server.
{
if(language == "EN") lege2 <- "Age of building" else lege2 <- "Rakennuksen ikä"
shp <- readOGR('PG:host=localhost user=postgres dbname=spatial_db','kuopio_house')
proj4string(shp) <- ("+init=epsg:3067") # The map projection of NLS Finland.
epsg4326String <- CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
shp2 <- spTransform(shp,epsg4326String) # Convert to longitude-latitude projection.
MyRmap(shp2, plotvar = "ika", legend_title = lege2, numbins = 8, pch = 19, cex = 0.3) # Draw the map.
}
if(exists("server")) BS <- 24 else BS <- 12
### Impact graphs in ENGLISH
if(language == "EN") {
plo <- ggplot(buildings@output) + geom_bar() + facet_grid(. ~ RenovationPolicy) + theme_gray(base_size = BS) +
aes(x = Year, weight = buildingsResult/1000000, fill = Renovation) + labs(y = "Floor area (M m2)",
title = "Building impacts of renovation policy")
print(plo)
plo <- ggplot(heatingEnergy@output) + geom_bar() + facet_grid(. ~ RenovationPolicy) + theme_gray(base_size = BS) +
labs(title = "Energy impacts of renovation policy", y = "Heating energy need (GWh /a)") +
aes(x = Year, weight = heatingEnergyResult/1E+6, fill = Heating)
print(plo)
emis@output <- emis@output[emis@output$Renovation == "BAU" , ]
plo <- ggplot(emis@output) + geom_bar() + facet_grid(Pollutant ~ FuelPolicy, scales = "free_y") + theme_gray(base_size = BS) +
labs(title = "Emission impacts of biofuel policy", y = "Emissions to air (ton /a)") +
aes(x = Year, weight = Result, fill = Fuel)
print(plo)
plo <- ggplot(health@output) + geom_bar() + facet_grid(FuelPolicy ~ RenovationPolicy) + theme_gray(base_size = BS) +
labs(title = "Health impacts of fuel and renovation policy", y = "Premature deaths (# /a)") +
aes(x = Year, weight = Result, fill = Heating)
print(plo)
}
### Impact graphs in FINNISH
if(language == "FI") {
plo <- ggplot(buildings@output) + geom_bar() + facet_grid(. ~ Remonttipolitiikka) + theme_gray(base_size = BS) +
aes(x = Year, weight = buildingsResult/1000000, fill = Renovation) + labs(y = "Rakennuspinta-ala (M m2)",
title = "Remonttipolitiikan vaikutus rakennuskantaan")
print(plo)
plo <- ggplot(heatingEnergy@output) + geom_bar() + facet_grid(. ~ Remonttipolitiikka) + theme_gray(base_size = BS) +
labs(title = "Remonttipolitiikan vaikutus energiankulutukseen", y = "Lämmitysenergian tarve (GWh /a)") +
aes(x = Year, weight = heatingEnergyResult/1E+6, fill = Heating)
print(plo)
emis@output <- emis@output[emis@output$Remonttipolitiikka == "BAU" , ]
plo <- ggplot(emis@output) + geom_bar() + facet_grid(Pollutant ~ Bioenergiapolitiikka, scales = "free_y") +
theme_gray(base_size = BS) + labs(title = "Bioenergiapolitiikan päästövaikutukset", y = "Päästöt ilmaan (ton /a)") +
aes(x = Year, weight = Result, fill = Fuel)
print(plo)
plo <- ggplot(health@output) + geom_bar() + facet_grid(Bioenergiapolitiikka ~ Remonttipolitiikka) + theme_gray(base_size = BS) +
labs(title = "Remontti- ja bioenergiapolitiikan terveysvaikutukset", y = "Ennenaikaisia kuolemia (# /a)") +
aes(x = Year, weight = Result, fill = Heating)
print(plo)
}
ggplot(dose@output, aes(x = Year, y = doseResult, colour = RenovationPolicy)) + geom_point() + facet_grid(Exposcen ~ FuelPolicy)
ggplot(totcases@output, aes(x = Year, y = totcasesResult, colour = Exposcen)) + geom_point() + facet_grid(Trait ~ FuelPolicy)
ggplot(RR@output, aes(x = RRResult, colour = Trait)) + geom_density()
| |