knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 5 ) library(iomlifetR)
This vignette demonstrates impact assessment with various lag choices
For this example we will use the HRAPIE mean relative risk estiamte RR = 1.06 per 10 $\mu g/m^3$ of PM~2.5~.
Functions in this package use demographic data in a standardised format stored in demog_data
. For this example we will use the included single_year_data
to create demog_data
.
head(single_year_data) population <- subset(single_year_data, time == 2011 & sex == "Persons" & measure == "Population") population <- population[, c("age", "value")] population$age <- as.numeric(gsub(" .+", "", population$age)) colnames(population)[2] <- "population" deaths <- subset(single_year_data, time == 2011 & sex == "Persons"& measure == "Deaths") deaths <- deaths[, "value"] demog_data <- data.frame(population, deaths = deaths)
no_lag <- impact(demog_data) no_lag <- data.frame(lapply(no_lag, colSums))
Compare deaths avoided
plot(rownames(no_lag), no_lag$deaths_current, xlab = "Year", ylab = "Number", main = "Deaths avoided", type = "b") points(rownames(no_lag), no_lag$deaths_ext, col = "red") abline(h = 0) legend("topright", c("Current cohort", "Extended population"), col=c("black", "red"), pch = 21, lty=2)
Comapre life-years gained
plot(rownames(no_lag), no_lag$ly_current, xlab = "Year", ylab = "Number", main = "Life-years gained", type = "b") points(rownames(no_lag), no_lag$ly_ext, col = "red") abline(h = 0) legend("bottom", c("Current cohort", "Extended population"), col=c("black", "red"), pch = 21, lty=2)
lag <- cumsum(c(0.3, rep(0.5/4, 4), rep(0.2/15, 15))) epa_lag <- impact(demog_data, lag_structure = lag) epa_lag <- data.frame(lapply(epa_lag, colSums))
Plot a comparison of no lag and US EPA lag (Extended cohort)
plot(rownames(no_lag), no_lag$deaths_ext, xlab = "Year", ylab = "Number", main = "Deaths avoided", type = "b") points(rownames(epa_lag), epa_lag$deaths_ext, col = "red", type = "b") abline(h = 0) legend("topright", c("No lag", "US EPA lag"), col=c("black", "red"), pch = 21, lty=2)
Compare total life-years gained
plot(rownames(no_lag), no_lag$ly_ext, xlab = "Year", ylab = "Number", main = "Life-years gained" , type = "b") points(rownames(epa_lag), epa_lag$ly_ext, col = "red", type = "b") abline(h = 0) legend("bottomright", c("No lag", "US EPA lag"), col=c("black", "red"), pch = 21, lty=2)
Calculate results
pm <- seq(0.1, 1, 0.1) slow_pm <- impact(demog_data, delta_pm = pm, lag_structure = lag) slow_pm <- data.frame(lapply(slow_pm, colSums))
Compare deaths avoided
plot(rownames(no_lag), no_lag$deaths_ext, xlab = "Year", ylab = "Number", main = "Deaths avoided", type = "b") points(rownames(epa_lag), epa_lag$deaths_ext, col = "red", type = "b") points(rownames(slow_pm), slow_pm$deaths_ext, col = "blue", type = "b") abline(h = 0) legend("topright", c("No lag", "US EPA lag", "Lag and gradual fall in PM"), col=c("black", "red", "blue"), pch = 21, lty=2)
Compare life years gained
plot(rownames(no_lag), no_lag$ly_ext, xlab = "Year", ylab = "Number", main = "Life-years gained", type = "b") points(rownames(epa_lag), epa_lag$ly_ext, col = "red", type = "b") points(rownames(slow_pm), slow_pm$ly_ext, col = "blue", type = "b") abline(h = 0) legend("bottomright", c("No lag", "US EPA lag", "Lag and gradual fall in PM"), col=c("black", "red", "blue"), pch = 21, lty=2)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.