4. Output Data Processing

Process results

Data

The pseudomised data is available in the public_data folder as an example. These are NOT the original values from KBS database.

These datasets are relevant:

contrib_per_unit.xlsx
constraints.xlsx
bounds.xlsx
nutref.xlsx - this is the nutrition reference file for comparison.

The result data, result.RData is pre-computed loaded. Please refer to 3. Run Optimisation for how to reach this result.

Load data

library(readxl)
library(data.table)

# files needed during cleaning up: 
intake <- read_excel("data/bounds.xlsx") 

# contstraints
constraints <- read_excel("data/constraints.xlsx") 

# contrib per unit
cpu <- read_excel("data/contrib_per_unit.xlsx")

# nutrition reference
nutref <- read_excel("data/nutref.xlsx")

Load the results from previous step,

results <- readRDS('data/result.RData')
# examine the strucure of this R object
results |> str()

List of 4
 $ xopt  : num [1:53] 8.12 138.62 4.2 13.72 20.26 ...
 $ obj   : num -15.4
 $ status: int 1
 $ extra :List of 2
  ..$ lambda: num [1:53] 4.52 5.32e-12 -1.27e-11 1.45e-11 -3.02e-12 ...
  ..$ slack : num [1:53] -5.33e-15 -4.22e+01 1.77e+01 -1.55e+01 7.43e+01 ...

The first object in this list, results$xopt is the one we need.

Compare diets

New diet chart

Create chart to compare food composition of diets (g of each food group in new vs baseline diet, percent change between diets, lower and upper intake limits aka realism constraints)

new_intake <- round(results$xopt, digits=1) 
# New diet name (in this case 'new_intake' aka Nutrition, 
# since this scenario includes only nutritional and realism constraints)

# Baseline <- round(constraintsdata[,2],digits=1) #Baseline diet
intake_baseline <- round(intake$mean, digits=1) #Baseline diet
intake_upr <- intake$upper_bound
intake_lwr <- intake$lower_bound
foodgroup <- intake$Foodgroup

#Percent change between new diet and baseline diet
percent_change <- round((new_intake - intake_baseline)/intake_baseline, digits=3) 
percent_change

 [1] -0.905  0.194  0.105  0.384  0.829  1.529  0.212  1.095  0.000 -0.045
[11]  0.380  0.691  0.368  0.298  3.878  0.000 -0.208  0.295  0.317  1.361
[21]  0.034 -0.744  0.453 -0.903  0.437 -0.903 -0.905  0.061 -0.690 -0.737
[31] -0.231 -0.398  0.321  0.683  0.910  0.419 -0.169  0.516 -0.355 -0.254
[41]  0.189 -0.543  0.488 -0.902 -0.048 -0.048  0.369 -0.172 -1.000 -0.080
[51] -0.294  0.000 -0.903

# food chart for compact display
diet_result <- data.frame(
  foodgroup = foodgroup,
  baseline = intake_baseline, 
  new_intake = new_intake, 
  percent_change = percent_change,
  lower_limit = intake_lwr, 
  upper_limit = intake_upr
)

head(diet_result)

          foodgroup baseline new_intake percent_change lower_limit upper_limit
1        Bread fine     85.3        8.1         -0.905       8.526     244.545
2      Bread coarse    116.1      138.6          0.194      11.613     289.275
3     Flours grains      3.8        4.2          0.105       0.378      20.580
4              Rice      9.9       13.7          0.384       0.987      53.130
5             Pasta     11.1       20.3          0.829       1.113      61.635
6 Breakfast cereals     17.2       43.5          1.529       1.722      99.120

Constraint on nutrients

Code for selecting names and renaming them

cpu <- dplyr::select(cpu , 
                     "Foodgroup", 
                     "Energy (MJ)" = "Energy (MJ)", 
                     "Protein, g, lower" = "Protein1", 
                     "Protein, g, upper" = "Protein2", 
                     "Carbohydrates, g, lower"="Available carbohydrates1", 
                     "Carbohydrates, g, upper"="Available carbohydrates2",
                     "Added sugar, g"="Added Sugar",
                     "Dietary fiber, g"= "Dietary fiber", 
                     "Fat, g, lower"= "Fat1",  
                     "Fat, g, upper"= "Fat2",
                     "Saturated fatty acids, g"="Sum saturated fatty acids",
                     "Trans fatty acids, g"="Sum trans fatty acids", #Added
                     "n-3 fatty acids, g"="Sum n-3 fatty acids",
                     "ALA, g" = "Sum ALA",
                     "MUFA, g, lower" = "Sum monounsaturated fatty acids1",
                     "MUFA, g, upper" = "Sum monounsaturated fatty acids2",
                     "PUFA, g, lower" = "Sum polyunsaturated fatty acids1",
                     "PUFA, g, upper" = "Sum polyunsaturated fatty acids2",
                     "Vitamin A, RE µg"= "Vitamin A",
                     "Vitamin E, alfa-TE"= "Vitamin E",
                     "Thiamin (Vitamin B1), mg"= "Thiamin (Vitamin B1)",
                     "Riboflavin (Vitamin B2), mg"= "Riboflavin (Vitamin B2)",
                     "Niacin, NE"="Niacin equivalent",
                     "Vitamin B6, mg"="Vitamin B6",
                     "Folate, µg"="Folate",
                     "Vitamin B12, µg"="Vitamin B12",
                     "Vitamin C, mg"="Vitamin C",
                     "Vitamin D, µg" = "Vitamin D",
                     "Sodium, mg"="Sodium",
                     "Potassium, mg"="Potassium",
                     "Calcium, mg"="Calcium",
                     "Magnesium, mg"= "Magnesium",
                     "Phosphorus, mg"="Phosphorus",
                     "Iron, mg"="Iron",
                     "Zinc, mg"="Zinc",
                     "Iodine, µg"="Iodine",
                     "Selenium, µg"="Selenium",
                     "Copper, µg"="Copper",
                     "Alcohol, g" = "Alcohol",
                     "GHGE, kg CO2-eq" = "GHGE",
                     "FE, g P-eq" = "FE",
                     "ME, g N-eq" = "ME",
                     "TA, g SO2-eq" = "ACID",
                     "WU, m2" = "WU",
                     "LU, m3a" = "LU",
                     "Whole grains, g" = "Whole grains",
                     "Total fruit, g" = "Fruit",
                     "Total vegetables, g" = "Vegetables",
                     "Total dairy, minimum, g" = "Dairy1",
                     "Total dairy, maximum, g" = "Dairy2",
                     "Total fish, minimum, g" = "Fish1",
                     "Total fish, maximum, g" = "Fish2",
                     "Total red meat, maximum, g" = "Red meat",
                     "Total white meat, maximum, g" = "White meat" 
)

cstr <- nutref |> 
  dplyr::select("Energy (MJ)",
         "Protein, g, lower" = "Protein1", 
         "Protein, g, upper" = "Protein2", 
         "Carbohydrates, g, lower"="Available carbohydrates1", 
         "Carbohydrates, g, upper"="Available carbohydrates2",
         "Added sugar, g"="Added Sugar",
         "Dietary fiber, g"= "Dietary fiber", 
         "Fat, g, lower"= "Fat1",  
         "Fat, g, upper"= "Fat2",
         "Saturated fatty acids, g"="Sum saturated fatty acids",
         "Trans fatty acids, g"="Sum trans fatty acids", #Added
         "n-3 fatty acids, g"="Sum n-3 fatty acids",
         "ALA, g" = "Sum ALA",
         "MUFA, g, lower" = "Sum monounsaturated fatty acids1",
         "MUFA, g, upper" = "Sum monounsaturated fatty acids2",
         "PUFA, g, lower" = "Sum polyunsaturated fatty acids1",
         "PUFA, g, upper" = "Sum polyunsaturated fatty acids2",
         "Vitamin A, RE µg"= "Vitamin A",
         "Vitamin E, alfa-TE"= "Vitamin E",
         "Thiamin (Vitamin B1), mg"= "Thiamin (Vitamin B1)",
         "Riboflavin (Vitamin B2), mg"= "Riboflavin (Vitamin B2)",
         "Niacin, NE"="Niacin equivalent",
         "Vitamin B6, mg"="Vitamin B6",
         "Folate, µg"="Folate",
         "Vitamin B12, µg"="Vitamin B12",
         "Vitamin C, mg"="Vitamin C",
         "Vitamin D, µg" = "Vitamin D",
         "Sodium, mg"="Sodium",
         "Potassium, mg"="Potassium",
         "Calcium, mg"="Calcium",
         "Magnesium, mg"= "Magnesium",
         "Phosphorus, mg"="Phosphorus",
         "Iron, mg"="Iron",
         "Zinc, mg"="Zinc",
         "Iodine, µg"="Iodine",
         "Selenium, µg"="Selenium",
         "Copper, µg"="Copper",
         "Alcohol, g" = "Alcohol",
         "GHGE, kg CO2-eq" = "GHGE",
         "FE, g P-eq" = "FE",
         "ME, g N-eq" = "ME",
         "TA, g SO2-eq" = "ACID",
         "WU, m2" = "WU",
         "LU, m3a" = "LU",
         "Whole grains, g" = "Whole grains",
         "Total fruit, g" = "Fruit",
         "Total vegetables, g" = "Vegetables",
         "Total dairy, minimum, g" = "Dairy1",
         "Total dairy, maximum, g" = "Dairy2",
         "Total fish, minimum, g" = "Fish1",
         "Total fish, maximum, g" = "Fish2",
         "Total red meat, maximum, g" = "Red meat",
         "Total white meat, maximum, g" = "White meat"  
  )

Compute the nutrient for the new diet, and compare with the upper and lower constraints. Make sure the names match!

#Create chart and assign column names and row names
# colnames(output_newdiet)
# colnames(cstr)
# data.frame(colnames(output_newdiet), colnames(cstr))
# make sure the match the names!
output_newdiet <- t(as.matrix(new_intake)) %*%  as.matrix(cpu[,2:54])

Combine the results and add tags

const_result <- t(rbind(output_newdiet, cstr))
colnames(const_result) <- c('new_diet','const_lwr', 'const_upr')
head(const_result)

                         new_diet const_lwr const_upr
Energy (MJ)              10.49982      10.0      10.0
Protein, g, lower       107.21292      59.9     119.8
Protein, g, upper       107.21292      59.9     119.8
Carbohydrates, g, lower 299.17897     269.5     359.3
Carbohydrates, g, upper 299.17897     269.5     359.3
Added sugar, g           51.97298       0.0      59.9

It can be convenient to add some indicators for us to read easily.

# extract nutrition names
nutnames <- constraints$tag_outcome

const_result <- data.table(const_result)
const_result[, nutrient := nutnames]

#Check if results meet constraints
const_result[, is_ok := 'Yes']
const_result[new_diet < const_lwr, is_ok := 'beyond lower']
const_result[new_diet > const_upr, is_ok := 'beyond upper']
const_result[, relative_dev := 0]
const_result[is_ok == 'beyond lower', relative_dev := round((new_diet - const_lwr)/const_lwr, 3)]
const_result[is_ok == 'beyond upper', relative_dev := round((new_diet - const_upr)/const_upr, 3)]

head(const_result)

    new_diet const_lwr const_upr                nutrient        is_ok
       <num>     <num>     <num>                  <char>       <char>
1:  10.49982      10.0      10.0             Energy (MJ) beyond upper
2: 107.21292      59.9     119.8       Protein, g, lower          Yes
3: 107.21292      59.9     119.8       Protein, g, upper          Yes
4: 299.17897     269.5     359.3 Carbohydrates, g, lower          Yes
5: 299.17897     269.5     359.3 Carbohydrates, g, upper          Yes
6:  51.97298       0.0      59.9          Added sugar, g          Yes
   relative_dev
          <num>
1:         0.05
2:         0.00
3:         0.00
4:         0.00
5:         0.00
6:         0.00

Departure

Diet departure is customized indicator for comparing new and current diet. They are based on customized functions, so you need to be careful with the formula and implementation.

Technical note

Diet departure was computed from a formula where each diet element is manually selected and compared, then summed together. This can be compactly written as below, utilizing the vectorized computation. Please compare if this is the object you need!

currentintake <- intake_baseline
head(diet_result)

          foodgroup baseline new_intake percent_change lower_limit upper_limit
1        Bread fine     85.3        8.1         -0.905       8.526     244.545
2      Bread coarse    116.1      138.6          0.194      11.613     289.275
3     Flours grains      3.8        4.2          0.105       0.378      20.580
4              Rice      9.9       13.7          0.384       0.987      53.130
5             Pasta     11.1       20.3          0.829       1.113      61.635
6 Breakfast cereals     17.2       43.5          1.529       1.722      99.120

# sum(diet_result$percent_change^2) # slightly different numbers, due to rounding

# departure
diet_departure <- ((diet_result$new_intake - diet_result$baseline)/diet_result$baseline)^2 |> sum()
diet_departure

[1] 33.4085

Red meat departure, need to select a subset that are relevant for red meat first:

# red meat departure
# select necessary data rows 
redmeat <- diet_result[c(23, 24, 25, 26, 27),]
redmeat

                    foodgroup baseline new_intake percent_change lower_limit
23 Ruminant meat, unprocessed     23.2       33.7          0.453      2.3205
24   Ruminant meat, processed     24.7        2.4         -0.903      2.4675
25          Pork, unprocessed     21.5       30.9          0.437      2.1525
26            Pork, processed     23.6        2.3         -0.903      2.3625
27  Red meat mixed, processed     34.6        3.3         -0.905      3.4650
   upper_limit
23     153.405
24     138.600
25     139.650
26      86.100
27     151.725

redmeat_departure <- ((sum(redmeat$new_intake) - sum(redmeat$baseline))/sum(redmeat$baseline))*100
redmeat_departure

[1] -43.10345

It is convenient to use a function instead of computing from formula

# define function 
f_meat_departure <- function(meat_data){
  
  meat_departure <- ((sum(meat_data$new_intake) - sum(meat_data$baseline))/sum(meat_data$baseline))*100
  return(meat_departure)
}

# use the function 
f_meat_departure(meat_data = redmeat)

[1] -43.10345

This is giving the same value as above. Now use the formula for the other departures. First select the relevant foods, then pass them into the function.

# total meat departure
totalmeat <- diet_result[c(23, 24, 25, 26, 27, 28, 29, 30, 31),]
f_meat_departure(meat_data = totalmeat)

[1] -34.68208

# ruminant meat departure
ruminantmeat <- diet_result[c(23, 24),]
f_meat_departure(meat_data = ruminantmeat)

[1] -24.63466