# Code to accompany my blog post entitled
#  "Some problems in a high-profile study of processed vs unprocessed diets".
# By Nick Brown, January 2021.
# License: CC-0.

library(sas7bdat)             # if necessary: install.packages("sas7bdat")

ucfirst <- function (s)
{
  paste0(substr(toupper(s), 1, 1), substr(s, 2, nchar(s)))
}

adir <- "."                   # put a relative or absolute path here if your data files are not in the same directory as the code
#adir <- "./ADLDataSAScode1"  # (example for a subdirectory)

# Change 0 to 1 in the next line to export the SAS data files in CSV format.
if (0) {
  files <- list.files(path=adir, pattern="*.sas7bdat", full.names=FALSE)
  lapply (files, function(x) {
    sasfile <- (paste0(adir, "/", x))
    csvfile <- paste0(adir, "/", gsub("\\.sas7bdat", ".csv", x))
    sas <- read.sas7bdat(sasfile)
    write.csv(sas, csvfile, row.names=FALSE)
  })
}

df.baseline <- read.sas7bdat(paste0(adir, "/", "baseline.sas7bdat"))
names(df.baseline)[1] <- "StudyID"     # for consistency with other tables; file contains "subjectID"
names(df.baseline)[6] <- "Day1BW"      # file contains "Bw", which will be confusing later
df.baseline <- df.baseline[order(df.baseline$StudyID),]

df.dailyintake <- read.sas7bdat(paste0(adir, "/", "dailyintake.sas7bdat"))
df.dailyintake <- df.dailyintake[order(df.dailyintake$StudyID, df.dailyintake$Period, df.dailyintake$DaysOnDiet),]
ndaily <- nrow(df.dailyintake)

df.deltabc <- read.sas7bdat(paste0(adir, "/", "deltabc.sas7bdat"))
df.deltabc <- df.deltabc[order(df.deltabc$StudyID),]

df.deltabcadj14 <- read.sas7bdat(paste0(adir, "/", "deltabcadj14.sas7bdat"))
df.deltabcadj14 <- df.deltabcadj14[order(df.deltabcadj14$StudyID),]

df.intakebymeal <- read.sas7bdat(paste0(adir, "/", "intakebymeal.sas7bdat"))
df.intakebymeal <- df.intakebymeal[order(df.intakebymeal$StudyID, df.intakebymeal$Period, df.intakebymeal$DaysOnDiet, df.intakebymeal$Meal),]

for (i in unique(df.intakebymeal$Day)) {      # correct erroneous diet for ADL002
  adl002 <- (
      (df.intakebymeal$StudyID == "ADL002")
    & (df.intakebymeal$Day == i)
  )
  df.intakebymeal[((adl002) & (df.intakebymeal$Meal != "Snack")),]$Period <-
    df.intakebymeal[((adl002) & (df.intakebymeal$Meal == "Snack")),]$Period
}

df.chamber <- read.sas7bdat(paste0(adir, "/", "chamber.sas7bdat"))
df.chamber <- df.chamber[order(df.chamber$SubID, df.chamber$Diet, df.chamber$Visit),]

# Compare daily total EI with macronutrient calories.
macro.fields <- c("Protein_Kcal", "fat_Kcal", "CHO_Kcal")
kcal.fields <- c("EI", macro.fields)
kcal.field.titles <- paste(ucfirst(gsub("_Kcal", "", kcal.fields)), collapse="\t")

df.dailyintake$calc.EI <- 0
df.dailyintake$diff.EI <- 0

for (i in 1:ndaily) {
  daily <- df.dailyintake[i,]

  calc.EI <- 0
  for (field in macro.fields) {
    calc.EI <- calc.EI + daily[[field]]
  }
  df.dailyintake[i,]$calc.EI <- calc.EI
  df.dailyintake[i,]$diff.EI <- daily$EI - calc.EI
}

ei.fields <- c("EI", "calc.EI", "diff.EI")
cat("Differences between daily EI field and calculated EI from macronutrient kcal data\n")
cat("StudyID\tPeriod\tDietDay\t", paste(ei.fields, collapse="\t"), "\n", sep="")

for (i in 1:ndaily) {
  daily <- df.dailyintake[i,]

  cat(daily$StudyID, daily$Period, daily$DaysOnDiet, sep="\t")
  for (field in ei.fields) {
    cat("\t", sprintf("%.2f", daily[[field]]), sep="")
  }
  cat("\n")
}
cat("------------------------------------------------------------\n")
cat("Maximum difference")
cat("\t", sprintf("%.2f", max(df.dailyintake$diff.EI)), sep="")
cat("\n")

cat("Minimum difference")
cat("\t", sprintf("%.2f", min(df.dailyintake$diff.EI)), sep="")
cat("\n")

cat("Minimum abs. diff.")
cat("\t", sprintf("%.2f", min(abs(df.dailyintake$diff.EI))), sep="")
cat("\n")

cat("------------------------------------------------------------\n")
cat("\n")

# Compare daily intake with sum of meal intakes.
for (field in kcal.fields) {
  diff.field <- paste0("diff.", field)
  df.dailyintake[diff.field] <- 0
}

for (i in 1:ndaily) {
  daily <- df.dailyintake[i,]

  meal.period <- daily$Period
  meals <- df.intakebymeal[
        (df.intakebymeal$StudyID == daily$StudyID)
      & (df.intakebymeal$Period == daily$Period)
      & (df.intakebymeal$DaysOnDiet == daily$DaysOnDiet)
    ,]
  for (field in kcal.fields) {
    diff.field <- paste0("diff.", field)
    df.dailyintake[i, diff.field] <- daily[[field]] - sum(meals[[field]])
  }
}

cat("Differences between dailyintake and intakebymeal data\n")
cat("StudyID\tPeriod\tDietDay\t", kcal.field.titles, "\n", sep="")
for (i in 1:ndaily) {
  daily <- df.dailyintake[i,]

  cat(daily$StudyID, daily$Period, daily$DaysOnDiet, sep="\t")
  for (field in kcal.fields) {
    diff.field <- paste0("diff.", field)
    cat("\t", sprintf("%.2f", daily[[diff.field]]), sep="")
  }
  cat("\n")
}

cat("------------------------------------------------------------\n")
cat("Maximum difference")
for (field in kcal.fields) {
  diff.field <- paste0("diff.", field)
  cat("\t", sprintf("%.2f", max(df.dailyintake[[diff.field]])), sep="")
}
cat("\n")

cat("Minimum difference")
for (field in kcal.fields) {
  diff.field <- paste0("diff.", field)
  cat("\t", sprintf("%.2f", min(df.dailyintake[[diff.field]])), sep="")
}
cat("\n")

cat("Minimum abs. diff.")
for (field in kcal.fields) {
  diff.field <- paste0("diff.", field)
  cat("\t", sprintf("%.2f", min(abs(df.dailyintake[[diff.field]]))), sep="")
}
cat("\n")

cat("------------------------------------------------------------\n")
cat("\n")

# Compare meal intake macronutrient calories with total energy input.
# Also examine relation between implied mass difference and macronutrient calorie/EI discrepancy.

nmeals <- nrow(df.intakebymeal)
cals.gram <- c(4, 9, 4)
names(cals.gram) <- macro.fields

df.intakebymeal$calc.EI <- NA
df.intakebymeal$calc.mass <- NA

for (i in 1:nmeals) {
  meal <- df.intakebymeal[i,]

  calc.EI <- 0
  for (field in macro.fields) {
    calc.EI <- calc.EI + meal[[field]]
  }
  df.intakebymeal[i,]$calc.EI <- calc.EI

  calc.mass <- meal$Water
  for (field in macro.fields) {
    calc.mass <- calc.mass + (meal[[field]] / cals.gram[field])
  }
  df.intakebymeal[i,]$calc.mass <- calc.mass
}

df.intakebymeal$diff.EI <- df.intakebymeal$calc.EI - df.intakebymeal$EI
df.intakebymeal$diff.mass <- df.intakebymeal$calc.mass - df.intakebymeal$Mass

min.diff.EI <- 1.0
cat("Differences >", min.diff.EI, " kcal between total EI and sum of macronutrient calories, per meal\n", sep="")
cat("StudyID\tPeriod\tDietDay\tMeal\t\t", paste(ei.fields, collapse="\t"), "\n", sep="")

df.intakebymeal.diff.EI.min <- df.intakebymeal[(abs(df.intakebymeal$diff.EI) >= min.diff.EI),]
for (i in 1:nrow(df.intakebymeal.diff.EI.min)) {
  meal <- df.intakebymeal.diff.EI.min[i,]
  cat(meal$StudyID, meal$Period, meal$DaysOnDiet, sprintf("%-12s", meal$Meal), sep="\t")
  for (ei.field in ei.fields) {
    cat("\t", sprintf("%.2f", meal[ei.field]), sep="")
  }
  cat("\n")
}

# Correlate missing mass with missing calories.
mass.cal.cor <- cor(df.intakebymeal.diff.EI.min$diff.mass, df.intakebymeal.diff.EI.min$diff.EI)
cat("\n")
cat("Correlation between missing calories and missing mass=", sprintf("%.3f", mass.cal.cor), "\n", sep="")

# Identify chamber days.
cat("\n")
cat("Day numbers within each 7-day meal cycle where energy intake corresponds to respiratory chamber data\n")
cat("StudyID\tPeriod\tMealCycleDay\n")

chamberDays <- rep(0, 7)
notSameDay <- 0
first <- 0
for (i in 1:nrow(df.chamber)) {
  cham <- df.chamber[i,]

  di <- df.dailyintake[
      (df.dailyintake$StudyID == cham$SubID)
    & (df.dailyintake$EI == cham$Intake)
  ,]

  chamberDay <- "?"
  dietDay <- 99
  if (nrow(di) == 1) {
    chamberDay <- di$DaysOnDiet
    dietDay <- ((chamberDay - 1) %% 7) + 1
    chamberDays[dietDay] <- chamberDays[dietDay] + 1
  }
  else if (nrow(di) > 2) {
    chamberDay <- "multiple"
  }

  if (first == 0) {
    first <- dietDay      # two records per participant per diet
  }
  else {
    if ((first != 99) && (dietDay != 99) && (first != dietDay)) {
      notSameDay <- notSameDay + 1
    }

    first <- 0
  }

  cat(di$StudyID, di$Period, chamberDay, "\n", sep="\t")
}

# Calculate mean sodium intake per kcal for each diet x day.
cat("\n")
cat("Mean sodium intake per kcal\n")
cat("Period\tDietDay\tSodium (mg)\n")

df.dailyintake$DietDay <- (((df.dailyintake$DaysOnDiet - 1) %% 7) + 1)
for (period in c("PROC", "UNPROC")) {
  p.sum <- 0
  for (dday in 1:7) {
    daily <- df.dailyintake[((df.dailyintake$Period == period) & (df.dailyintake$DietDay == dday)),]
    sodium.EI <- mean(daily$Sodium / daily$EI)
    p.sum <- p.sum + sodium.EI
    cat(period, "\t", dday, "\t", sprintf("%.2f", sodium.EI), "\n", sep="")
  }
  cat(period, "\t", "Mean", "\t", sprintf("%.2f", (p.sum / 7)), "\n", sep="")
}

# Compare unadjusted and adjusted weight and fat/fat-free mass records.
# First, copy records from original into adjusted data frame,
# and give the variables new names based on the adjusted names.
adj.names <- names(df.deltabcadj14)[-1]
old.names <- names(df.deltabc)[c(8, 7, 9, 10, 13, 14)]

for (i in 1:length(adj.names)) {
  adj.name <- adj.names[i]
  old.name <- old.names[i]

  orig.name <- paste0("orig.", adj.name)
  df.deltabcadj14[orig.name] <- df.deltabc[old.name]

  ratio.name <- paste0("ratio.", adj.name)
  df.deltabcadj14[ratio.name] <- df.deltabcadj14[adj.name] / df.deltabcadj14[orig.name]

  diff.name <- paste0("diff.", adj.name)
  df.deltabcadj14[diff.name] <- df.deltabcadj14[adj.name] - df.deltabcadj14[orig.name]
}

cat("\n")
cat("Largest (in magnitude) differences between adjusted and original data\n")
cat("Explore df.deltabcadj14 for more\n")
for (weight in c("BW", "FM", "FFM")) {
  for (period in c("Proc", "UnProc")) {
    var <- paste0("diff.delta", weight, period)
    sortv <- order(-abs(df.deltabcadj14[var]))
    show <- c("StudyID", var)
    cat("\n")
    print(head(df.deltabcadj14[sortv,][show], 5))
  }
}