Difference between revisions of "Compound intake estimator"

library(OpasnetUtils)
library(ggplot2)

################## SETUP MODEL

openv.setN(100) # Sets the number of model iterations to 100.

# Fetches the ovariables needed in the model from [[Compound intake estimator]].
# Ovariables downloaded: guidance.values, compound.intake, food.risk, extraction, concplot.

objects.latest("Op_en6193", code_name="initiate") 

# compound.conc is the concentration of compounds in the critical constituent.
# Data comes from [[Compound intake estimator]].

compound.conc <- Ovariable(name = "compound.conc", ddata = "Op_en6193.example_data")

################## USER INPUTS

# name is the name of product. It is used in the output graphs.

# food.amount is the amount that the consumer consumes the defined product (g /day)

food.amount <- Ovariable("food.amount", data = data.frame(Result = food.amount))

# fr.const is the fraction of a critical constituent in the food or food supplement (given as percentage)

fr.const <- Ovariable("fr.const", data = data.frame(Result = fr.const))

# oil.amount is the amount that the consumer consumes essential oil in addition to the product (g /day)

oil.amount <- Ovariable("oil.amount", data = data.frame(Result = oil.amount))

# extraction is a list of extraction techniques to consider. (by default, all extraction techniques)

if(exists('extracted')) {
	extraction@data <- extraction@data[extracted , ] # Order of extraction techniques must not change!
} else {
	extraction@data <- data.frame(extractionResult = 1) # Remove column Extraction.techniques
}

# addcomp are user-defined additional compounds for a supplement, if they are not in the table.
## addcompound is the name and addconc is the concentration of the additional compound.
# addcomp and addcompound can be vectors.

addcompound <- ifelse(! exists('addcompound') || addcompound == ' ', NA, addcompound)
addconc <- ifelse(! exists('addconc') || addconc == ' ', 0, addconc)

addcomp <- Ovariable("addcomp", data = data.frame(Compound = addcompound, Result = addconc))


############## ACTUAL MODEL

food.risk <- EvalOutput(food.risk)

################ PRINTING MODEL OUTPUTS

if(intermediates) {

	cat("Concentrations of compounds in ", name, " (mg /kg DW)\n")

	oprint(summary(compound.conc, marginals = c("Compound")))

	cat("Intakes of different compounds (mg /kg BW /d)\n")

	oprint(summary(compound.intake, marginals = c("Compound")), digits = 3)

}

cat("Risks relative to guidance values (ADI, TDI, BMDL10, or TTC). Values < 1 are of minimal concern and are not shown in the graph. Note! BMDL10 is used to calculate margin of exposure (MOE), but instead of MOE, (10000 / MOE) is shown on this graph for consistency.\n")

oprint(summary(food.risk, marginals = c("Compound")))

test <- unique(food.risk@output[food.risk@output$food.riskResult >= 1 , "Compound"])

temp <- food.risk@output[food.risk@output$Compound %in% test , ]

if(nrow(temp) == 0) {
	cat("There are no compounds that exceed a guidance value.\n")
	stop("There are no compounds that exceed a guidance value.\n")
}

ggplot(temp, aes(x = Compound, y = food.riskResult, fill = Guidance)) + 
	geom_boxplot() + 
	theme_grey(base_size = 24) +
	labs(title = paste("Exposure to compounds in", name), y = "Ratio (Exposure /Guidance)") +
	theme(axis.text.x = element_text(angle = 90, hjust = 1))

### Load some R packages needed to implement the code.

library(OpasnetUtils)
library(ggplot2)

################ SETUP THE MODEL

openv.setN(100) # Set model iterations

# Fetches the ovariables needed in the model from [[Compound intake estimator]].
# Ovariables downloaded: guidance.values, compound.intake, food.risk, extraction, concplot

objects.latest("Op_en6193", code_name="initiate") 

### Download the data from the database. You can also see the raw data on this page below,
### or access the database via this link: http://heande.opasnet.org/wiki/Special:Opasnet_Base?id=Heande3857

# compound.conc is the concentration of compounds in the critical constituent.
# Data comes from [[Compound intake estimator]].

compound.conc <- Ovariable(name = "compound.conc", ddata = "Op_en6193.example_data")

################ USER INPUTS

### Select only those rows that fulfill the criteria of compound and extraction method set by the user (if any).

dat <- compound.conc@data
if(!allcompounds) {dat <- dat[dat$Compound %in% compound , ]}
if(!allextractions) {dat <- dat[dat$Extraction.technique %in% extraction , ]}
compound.conc@data <- dat

### Create an information object (called ovariable) out of the data for easy handling.

compound.conc <- EvalOutput(compound.conc)

### Print a summary table of the output.

oprint(summary(compound.conc), digits = 3)

### Use the graph formatting for several graphs. The graphs only differ by the range of concentration that is shown on Y axis.

concplot(compound.conc@output$compound.concResult >= limit)
concplot(compound.conc@output$compound.concResult <= 10 & compound.conc@output$compound.concResult > 1)
concplot(compound.conc@output$compound.concResult <= 1 )

library(OpasnetUtils)

################## PRE-DEFINED INPUTS

# compound.conc is the concentration of a compound in the critical constituent (or food, if fr.const == 1) 
# fr.const is the fraction of the critical constituent in the whole food
# extraction is a criteria to select specific extraction techniques
# addcomp is an additional compound not in the product. Data about its concentration in food is given by the user.

# compound.intake is the intakes of the compounds via food or food supplement, each given as mg /kg /d.

compound.intake <- Ovariable("compound.intake", 
	dependencies = data.frame(Names = c(
		"food.amount", # amount of food supplement consumed
		"compound.conc", # concentration of compounds in essential oil
		"oil.amount", # amount of essential oil consumed in addition to the supplement product
		"fr.const", # fraction of essential oil in the supplement product
		"extraction", # Which extraction methods the user wants to consider
		"addcomp" # concentrations of compounds in the essential oil (missing from default data)
		# conc is the temporary ovariable for compound.conc used in calculations
	)), 
	formula = function(...) {

		conc <- compound.conc * 10 * extraction 
		# remove extractions that are not considered.
		# compound.conc from w/w% to mg/g

#		oapply would be an alternative to remove extra indices, but there is a bizarre error.
#		compound.conc <- oapply(compound.conc, cols = "Extraction.technique", fun = mean)

		test <- colnames(conc@output) != "Extraction.technique"
		conc@marginal	<- conc@marginal[test]
		conc@output 	<- conc@output[test] # unit: mg /g
		# removes index Extraction.technique, which is no longer needed.

		conc@output <- orbind(conc, addcomp * 1)
		# addcomp must not contain columns that are not in conc. It must contain Compound.
		# It is multiplied by 1 to make the result column name "Result".

		out <- (food.amount * fr.const / 100 + oil.amount) * conc / 60 
		# calculates the oil intake from food supplement and then adds direct oil intake.
		# scale to 60 kg person
		# fr.const from % to fraction

		test2 <- unique(out@output$Compound[result(out) > 0]) # remove all compounds where intake = 0
		out@output <- out@output[out@output$Compound %in% test2 , ]

		return(out)
	}
)

# guidance.values is the guidance value (ADI, BMDL10, TTC or similar) for compounds. unit: mg /kg /d
# cramer is cramer classes of compounds (temporary)
# cramer.values are the thresholds of toxicological concern for each Cramer class (temporary)

cramer <- tidy(opbase.data("Op_en6193", subset = "cramer_classes_of_compounds"))
cramer$Result <- NULL
cramer$Guidance <- "TTC" # Threshold for toxicological concern

cramer.values <- tidy(opbase.data("Op_en6193", subset = "cramer_values"))
cramer.values <- merge(cramer, cramer.values)
cramer.values[["Cramer class"]] <- NULL
cramer.values$Result <- cramer.values$Result / 1000 # from ug /kg /d to mg /kg /d

guidance.values <- tidy(opbase.data("Op_en6193", subset = "guidance_values"))
guidance.values <- orbind(cramer.values, guidance.values)
guidance.values <- Ovariable("guidance.values", data = guidance.values)

safety.factor <- Ovariable("safety.factor", ddata = "Op_en6193.safety_factors") #  [[Compound intake estimator]]

# food.risk calculates "hazard quotients" or risks relative to guidance values. Value < 1 are of minimal concern

food.risk <- Ovariable(name = "food.risk",
	dependencies = data.frame(Name = c("compound.intake", "guidance.values", "safety.factor")),
	formula = function(...) {

		# guidance.risk is intake compared with guidance values.

		out <- compound.intake / (guidance.values / safety.factor)

		removepoorguidance <- function(dat, goodg) { # Remove inferior guidance values

			# Compounds that have good guidance values (goodg).
			hasgoodg <- unique(dat[dat$Guidance == goodg , "Compound"])

			# The row contains data about a chemical with a good guidance value
			# but the data is about an inferior guidance value.

			out <- dat[!(dat$Compound %in% hasgoodg & dat$Guidance != goodg) , ]

			return(out)
		}

		out <- removepoorguidance(out@output, "ADI")
		out <- removepoorguidance(out, "TDI")
		out <- removepoorguidance(out, "BMDL10")

		test2 <- unique(out$Compound[out$Result > 0]) # remove all compounds where risk = 0
		out <- out[out$Compound %in% test2 , ]

		return(out)
	}
)

extraction <- Ovariable("extraction", ddata = "Op_en6193.extraction_techniques")
extraction@data$extractionResult <- 1

### Define formatting for the concentration graphs that will be shown.

concplot <- function(condition) {
	dat <- compound.conc@output[condition, ]
	if(nrow(dat) != 0) {
		out <- ggplot(dat, 
			aes(
				x = Compound, 
				y = compound.concResult, 
				colour = Extraction.technique
			)) + 
			geom_point(size=5) +
			theme_grey(base_size = 24) +
			theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
			labs(y = "Concentration (%)")
		return(out)
	} else {
		 warning("No data to show")
	}
}

objects.store(compound.intake, guidance.values, food.risk, extraction, concplot, safety.factor)

cat("The following objects were successfully stored:
compound.intake, guidance.values, food.risk, extraction, concplot, safety.factor\n")