#' @author Gerald C. Nelson, \email{nelson.gerald.c@@gmail.com}
#' @keywords utilities, nutrient data, IMPACT food commodities nutrient lookup
# Intro -------------------------------------------------------------------
#Copyright (C) 2015 Gerald C. Nelson, except where noted
# This program is free software: you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by the Free
# Software Foundation, either version 3 of the License, or (at your option)
# any later version.
#
# This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
# for more details at http://www.gnu.org/licenses/.
#' @description
# Intro -------------------------------------------------------------------
#' This script reads in IMPACT nutrient lookup table, does some manipulations of the data,
#and writes out results to an excel spreadsheet
#Copyright (C) 2015 Gerald C. Nelson, except where noted
#Important code contributions from Brendan Power.
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details at http://www.gnu.org/licenses/.
source(file = "R/setup.global.R")
# Load functions for the rest of the script
# source(file="nutrientFunctions.R") # may not be needed now (Nov 13)
# Read in data ----------
# source(file ="IMPACTdataLoading.R") replaced with the readRDS line
# this file is generated by dataPrep.IMPACT.R.) Only needs to be redone when new IMPACT data arrives
dt.df1 <- getNewestVersionIMPACT("IMPACTfood")
dt.fish <- getNewestVersion(("fishScenarios"))
dt.alc <- getNewestVersion(("alcScenarios"))
# read in nutrients data ----
# rds file generated in dataPrep.nutrients.R
nutrients <- getNewestVersion("nutrients")
dt.nutrients <- as.data.table(nutrients)
#list of columns with cooking retention values
#' @param cookretn.cols
cookretn.cols <-
colnames(nutrients.clean[, grep('_cr', names(nutrients.clean))])
# a switch to turn on or off the use of the cooking retention values
#' @param applyCookingRetention
applyCookingRetention <- "yes"
if (applyCookingRetention %in% "yes") {
for (i in 1:length(cookretn.cols)) {
nutrientName <-
substr(x = cookretn.cols[i], 1, nchar(cookretn.cols[i]) - 3)
nutColName <- paste ("nutrients$", nutrientName, sep = "")
# print(nutColName)
nutRetentColName <-
paste ("nutrients$", cookretn.cols[i], sep = "")
# print(nutRetentColName)
temp <-
as.data.frame(eval(parse(text = nutRetentColName)) * eval(parse(text = nutColName)) /
100)
colnames(temp) <- nutrientName
nutrients[, nutrientName] <- temp
}
}
#generated in SSPPopExtract.R because it is adjusted to the SSP age and gender groups
req.EAR <- getNewestVersion("req.EAR")
req.RDA.vits <- getNewestVersion("req.RDA.vits")
req.RDA.minrls <- getNewestVersion("req.RDA.minrls")
req.RDA.macro <- getNewestVersion("req.RDA.macro")
req.UL.vitsFile <- getNewestVersion("req.UL.vitsFile")
req.UL.minrls <- getNewestVersion("req.UL.minrls")
# calculate the share of per capita income spent on IMPACT commodities
setkey(dt.df1, "scenario", "region", "year")
dt.df1[, budget.Pw := sum(food * Pw / 365 / 1000), by = key(dt.df1)]
dt.df1[, budget.Pcon := sum(food * Pc * (1 - CSE) / 365 / 1000), by = key(dt.df1)]
dt.df1[, budget.Pc := sum(food * Pc / 365 / 1000), by = key(dt.df1)]
dt.budget <-
unique(dt.df1[, c(key(dt.df1), "budget.Pw", "budget.Pc", "budget.Pcon", "pcGDP"), with =
F])
#incomeShare <-join(t1.pcGDP,budget)
dt.budget[, incSharePw := budget.Pw / ((pcGDP * 1000) / 365)]
dt.budget[, incSharePc := budget.Pc / ((pcGDP * 1000) / 365)]
dt.budget[, incSharePcon := budget.Pcon / ((pcGDP * 1000) / 365)]
#nutrient stuff
#include only columns that are needed; IMPACT code and nutrient names
#choices are
# - all - all relevant nutrients, 23
# - macro - "energy", "protein", "fat", "carbohydrate", "fiber", "sugar"
# - minerals - "calcium", "iron", "potassium", "sodium", "zinc"
# - vitamins - "vitamin_c", "thiamin", "riboflavin", "niacin", "vitamin_b6", "folate", "vitamin_b12",
# "vitamin_a_RAE", "vitamin_e", "vitamin_d2_3"
# - fattyAcids - "ft_acds_tot_sat", "ft_acds_plyunst"
# next line is where the choice is made
# short.name <- c("vitamins")
# nut.list <- eval(parse(text = short.name))
# includes <- c("IMPACT_code", nut.list)
# nuts.reduced <- nutrients[, (names(nutrients) %in% includes)]
# setkey(dt.df1,"IMPACT_code")
# temp <- as.data.table(foodGroupsInfo)
# setkey(temp,"IMPACT_code")
# dt.d2 <- dt.df1[temp]
keepList <- c("scenario", "IMPACT_code", "region", "year", "food")
dt.food <- dt.df1[, names(dt.df1) %in% keepList, with = FALSE]
#convert from annual availability to daily availability. 'food' variable is in kgs/person/year
dt.food[, food := food / 365]
setkey(dt.food, "IMPACT_code")
#a list of the requirements types. Each has a different set of nutrients. These are a subset
# of what are in the nutrients requirements tables from IOM. They are the nutrients common to
# both the IOM and nutrient content lookup spreadsheet
reqList <-
c(
"req.EAR",
"req.RDA.vits" ,
"req.RDA.minrls",
"req.RDA.macro",
"req.UL.vits",
"req.UL.minrls"
)
for (i in 1:length(reqList)) {
#get list of nutrients from the requirement df
nut.list <-
unique(eval(parse(text = paste(
reqList[i], "$nutrient", sep = ""
))))
#nuts.reduced <- nutrients[, (names(nutrients) %in% nut.list)]
#keep only nutrients in the requirements being worked on
dt.nuts.reduced <-
dt.nutrients[, c("IMPACT_code", nut.list, "staple.code", "food.group.code"), with = FALSE]
setkey(dt.nuts.reduced, "IMPACT_code")
dt.d3 <- dt.food[dt.nuts.reduced]
# get sum of nutrient intake by food group
setkey(dt.d3, "scenario", "region", "food.group.code", "year")
nut.list.Q <- paste(nut.list, "Q", sep = ".")
nut.list.Q.sum <- paste(nut.list, "Q.sum", sep = ".")
#splitting up these operations seems to make them faster than doing it in one go
#note - * 10 is need because nutrients are per 100 g and food is in kgs
temp2 <-
dt.d3[, (nut.list.Q) := lapply(.SD, function(x)
(x * dt.d3[['food']] * 10)), .SDcols = nut.list]
temp3 <-
temp2[, (nut.list.Q.sum) := lapply(.SD, sum), .SDcols = nut.list.Q, by = key(dt.d3)]
dt.food.group.sum <-
unique(temp3[, c("scenario",
"region",
"food.group.code",
"year",
nut.list.Q.sum), with = F])
# get sum of nutrient intake by staple
setkey(dt.d3, "scenario", "region", "staple.code", "year")
temp3 <-
temp2[, (nut.list.Q.sum) := lapply(.SD, sum), .SDcols = nut.list.Q, by = key(dt.d3)]
dt.staples.sum <-
unique(temp3[, c("scenario", "region", "staple.code", "year", nut.list.Q.sum), with =
F])
# get sum of nutrient intake by country
setkey(dt.d3, "scenario", "region", "year")
temp3 <-
temp2[, (nut.list.Q.sum) := lapply(.SD, sum), .SDcols = nut.list.Q, by = key(dt.d3)]
dt.food.sum <-
unique(temp3[, c("scenario", "region", "year", nut.list.Q.sum), with =
F])
# convert the requirements dt from wide (years as columns) to long (add year and value columns)
dt.reqs.long <- melt(
eval(parse(text = reqList[i])),
id.vars = c("region", "nutrient"),
measure.vars = keepYearList,
variable.name = "year",
value.name = "nut_req"
)
dt.reqs.long[, year := as.character(year)]
# .. and then back to wide with nutrients as columns, to match the structure of the sum DTs.
dt.reqs <-
dcast.data.table(dt.reqs.long, region + year ~ nutrient, value.var = "nut_req")
setkey(dt.reqs, "region", "year")
setkey(dt.food.sum, "region", "year", "scenario")
dt.food.tot <- dt.reqs[dt.food.sum]
nut.list.ratio <- paste(nut.list, "ratio", sep = ".")
#loop through all the nutrients in this requirement set and calculate the ratio of consumption to requirement
for (j in 1:length(nut.list)) {
#calculate the share of the requirement from IMPACT commodities
dt.food.tot[, (nut.list.ratio[j]) := eval(parse(text = nut.list[j])) /
eval(parse(text = nut.list.Q.sum[j]))]
}
dt.food.ratio <- dt.food.tot[, nut.list.ratio, with = FALSE]
colMax <- function(data)
lapply(data, max, na.rm = TRUE)
colMin <- function(data)
lapply(data, min, na.rm = TRUE)
cMax <- as.data.table(colMax(dt.food.ratio))
cMin <- as.data.table(colMin(dt.food.ratio))
# dt.maxRows = dt.food.tot[0]
# rbind[dt.maxRows, dt.food.tot[maxRowNum,]]
#
saveRDS(dt.food.tot,
file = paste("results/", reqList[i], ".results.", Sys.Date(), ".rds", sep =
""))
#reshape the results to get years in columns
dt.food.tot.long <- melt(
dt.food.tot,
id.vars = c("scenario", "region", "year"),
measure.vars = c(nut.list, nut.list.Q.sum, nut.list.ratio),
variable.name = "nutrient",
value.name = "nut_req"
)
# dt.food.tot.long[,nutrient:= as.character(nutrient)]
dt.food.tot.wide <-
dcast.data.table(dt.food.tot.long,
scenario + region + nutrient ~ year,
value.var = "nut_req")
dt.food.tot.wide[, nutrient := as.character(nutrient)]
#set up initial worksheets in the spreadsheet for the ith set of requirements
tmp <- f.createGeneralWorksheet()
# this structure is needed to get two data frames out of the function
wbGeneral <- tmp[[1]]
wbInfoGeneral <- tmp[[2]]
#add a worksheet describing the nutrients in this requirement
addWorksheet(wbGeneral, sheetName = reqList[i])
writeData(
wbGeneral,
eval(parse(text = reqList[i])),
sheet = reqList[i],
startRow = 1,
startCol = 1,
rowNames = FALSE
)
addStyle(
wbGeneral,
sheet = reqList[i],
style = numStyle,
rows = 1:nrow(eval(parse(text = reqList[i]))) + 1,
cols = 2:ncol(eval(parse(text = reqList[i]))),
gridExpand = TRUE
)
setColWidths(
wbGeneral,
sheet = reqList[i],
cols = 1:ncol(eval(parse(text = reqList[i]))),
widths = "auto"
)
wbInfoGeneral[(nrow(wbInfoGeneral) + 1), ] <-
c(reqList[i], paste("metadata on nutrients included in ", reqList[i]))
# add data to the spreadsheet -----
for (j in 1:length(nut.list)) {
# add an info sheet for the jth nutrient in the ith set of requirements
temp <- dt.food.tot.wide[nutrient %in% nut.list.ratio[j]]
# add the results to the spreadsheet
addWorksheet(wbGeneral, sheetName = nut.list[j])
writeData(
wbGeneral,
temp,
sheet = nut.list[j],
startRow = 1,
startCol = 1,
rowNames = FALSE,
withFilter = TRUE
)
setColWidths(
wbGeneral,
sheet = nut.list[j],
cols = 1,
widths = "auto",
ignoreMergedCells = FALSE
)
addStyle(
wbGeneral,
sheet = nut.list[j],
style = numStyle,
rows = 1:nrow(temp),
cols = 2:ncol(temp),
gridExpand = TRUE
)
# code to create a graph of the nutrient ratios for each of the nutrients in the current requirements list for a single scenario
scenario.name <- "SSP2-GFDL"
nutrient.name <- sub("^(.*)_.*$", "\\1", nut.list[j])
temp <-
dt.food.tot.long[nutrient %like% "ratio" &
nutrient %like% nutrient.name & scenario %in% scenario.name, ]
temp$year <- gsub("X", "", temp$year)
temp$year <- as.numeric(temp$year)
gg <- ggplot(data = temp,
aes(x = year, y = nut_req, color = region)) + geom_line() +
theme(
axis.text.x = element_text(
color = "black",
angle = 0,
size = 9,
vjust = 0.5
),
axis.text.y = element_text(
color = "black",
size = 9,
vjust = 0.5
),
axis.title.y = element_text(
color = "black",
size = 9,
vjust = 0.5
),
legend.text = element_text(
color = "black",
size = 7,
vjust = 0.5
),
plot.title = element_text(
color = "black",
face = "bold",
size = 11,
hjust = 0.5,
vjust = 1
),
panel.background = element_blank(),
panel.border = element_rect(
linetype = "solid",
colour = "black",
fill = NA
),
legend.position = "bottom",
legend.title = element_blank(),
legend.key = element_rect(fill = "white"),
legend.background = element_rect(fill = NA)
) +
xlab("year") + ylab(nutrient.name) +
ggtitle(paste(
"Cons. share of require., ",
nutrient.name,
", ",
unique(temp$scenario),
sep = ""
)) +
guides(color = guide_legend(nrow = 5, byrow = TRUE))
# add plot to the spreadsheet
plotSheet <- paste("plot_", nutrient.name, sep = "")
print(plotSheet)
addWorksheet(wbGeneral, sheetName = plotSheet)
print(gg) #needs to be showing for the next step to work
insertPlot(
wbGeneral,
sheet = plotSheet,
width = 9,
height = 5,
fileType = "png",
units = "in",
startRow = 5
)
wbInfoGeneral[(nrow(wbInfoGeneral) + 1), ] <-
c(nut.list[j], paste("Ratio results for ", nutrient.name))
#write the rows that have the max and min values for the current nutrient
# eval parse needed here to get the which to run. uggh. This displays the row num of the row with the max/min value
maxRowNum <-
which(dt.food.ratio[, eval(parse(text = colnames(cMax)[j]))] == cMax[1, eval(parse(text = colnames(cMax)[j]))])
minRowNum <-
which(dt.food.ratio[, eval(parse(text = colnames(cMin)[j]))] == cMin[1, eval(parse(text = colnames(cMin)[j]))])
writeData(
wbGeneral,
dt.food.tot[maxRowNum, ],
sheet = plotSheet,
startRow = 1,
startCol = 2,
rowNames = FALSE,
withFilter = FALSE
)
writeData(
wbGeneral,
paste(
"Country and year with max value for the ",
nutrient.name,
" ratio",
sep = ""
),
sheet = plotSheet,
startRow = 1,
startCol = 1,
rowNames = FALSE,
withFilter = FALSE
)
writeData(
wbGeneral,
dt.food.tot[minRowNum, ],
sheet = plotSheet,
startRow = 3,
startCol = 2,
rowNames = FALSE,
withFilter = FALSE
)
writeData(
wbGeneral,
paste("row with min value for the ", nutrient.name, " ratio", sep = ""),
sheet = plotSheet,
startRow = 3,
startCol = 1,
rowNames = FALSE,
withFilter = FALSE
)
addStyle(
wbGeneral,
sheet = plotSheet,
style = numStyle,
rows = 1:4,
cols = 2:ncol(dt.food.tot),
gridExpand = TRUE
)
setColWidths(
wbGeneral,
sheet = plotSheet,
cols = 1:ncol(dt.food.tot[minRowNum, ]),
widths = "auto",
ignoreMergedCells = FALSE
)
wbInfoGeneral[(nrow(wbInfoGeneral) + 1), ] <-
c(plotSheet,
paste("Max/min rows and graph of results for ", nutrient.name))
}
#----------------------------
short.name <- reqList[i]
f.finalizeWB(wbGeneral, wbInfoGeneral, short.name)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.