#### Paper conservation manuscript figures and analysis
# Haley Carter

#### libraries
library(tidyverse)
library(rstatix)
library(ggpubr)
library(colorspace)
library(gtsummary)
library(kableExtra)

#### data...
setwd("/Users/haley/Documents/DataScience/R_LibraryConservation/Analyses/UpdatedAnalyses/")
colorData <- read.csv("full_colordata.csv") 
# figure 5, allStatsTable
experimental <- read.csv("full_experimental.csv") 
# plot 3
colorimetry <- read.csv("full_colorimetry.csv")  
# plot 1
cumulativeChange <- read.csv("full_cumulativechange.csv")
# plot 2
conductivity <- read.csv("full_conductivity.csv") 
# plot 4 (transformed to outliersRemovedConduct)

#### adding factor levels to data
#####
colorimetry$paper_type<-factor(colorimetry$paper_type, levels=c("Whatman Chromatography","18c","Canson Ingres", "Stonehenge Watercolor", "Alum Rosin"))
colorimetry$icondition <- factor(colorimetry$icondition, levels = c("BT-AT", "AT-AD", "BT-AD"))

conductivity$timing <- factor(conductivity$timing, levels = c("BT", "AT", "AD"))

colorData$gel<-factor(colorData$gel, levels=c("none","water","tica","kelco", "agarose"))
colorData$paper_type<-factor(colorData$paper_type, levels=c("Whatman Chromatography","18c","Canson Ingres", "Stonehenge Watercolor", "Alum Rosin"))
colorData$condition <- factor(colorData$condition, levels = c("BT", "AT", "AD"))
colorData$interleaving <- factor(colorData$interleaving, levels = c("none", "interleaving"))
colorData$gel_treatment <- factor(colorData$gel_treatment, levels = c("none_none", "water_none", "tica_none", "tica_interleaving", "kelco_none", "kelco_interleaving", "agarose_none", "agarose_interleaving"))
colorData$interleaving<-factor(colorData$interleaving, levels=c("none", "interleaving"))

cumulativeChange$gel<-factor(cumulativeChange$gel, levels=c("none","water","tica","kelco", "agarose"))
cumulativeChange$gel_treatment <- factor(cumulativeChange$gel_treatment, levels = c("none_none", "water_none", "tica_none", "tica_interleaving", "kelco_none", "kelco_interleaving", "agarose_none", "agarose_interleaving"))
cumulativeChange$interleaving<-factor(cumulativeChange$interleaving, levels=c("none", "interleaving"))
cumulativeChange$paper_labels<-factor(cumulativeChange$paper_labels, levels=c("Whatman","Antique Gelatin","Modern Gelatin", "Starch-based AKD", "Alum Rosin"))

experimental$gel<-factor(experimental$gel, levels=c("none","water","tica","kelco", "agarose"))
experimental$gel_treatment <- factor(experimental$gel_treatment, levels = c("none_none", "water_none", "tica_none", "tica_interleaving", "kelco_none", "kelco_interleaving", "agarose_none", "agarose_interleaving"))
experimental$paper_labels<-factor(experimental$paper_labels, levels=c("Whatman","Antique Gelatin","Modern Gelatin", "Starch-based AKD", "Alum Rosin"))
experimental$interleaving<-factor(experimental$interleaving, levels=c("none", "interleaving"))
experimental$condition <- factor(experimental$condition, levels = c("BT", "AT", "AD"))
#####

#### changing interleaving codes
#####
colorData$INT <- NA
for (row in 1:nrow(colorData)) {
  if (colorData[row, "interleaving"] == "interleaving") {
    colorData[row, "INT"] <-"INT"
  } else colorData[row, "INT"] <- "DA"
}
cumulativeChange$INT <- NA
for (row in 1:nrow(cumulativeChange)) {
  if (cumulativeChange[row, "interleaving"] == "interleaving") {
    cumulativeChange[row, "INT"] <-"INT"
  } else cumulativeChange[row, "INT"] <- "DA"
}
experimental$INT <- NA
for (row in 1:nrow(experimental)) {
  if (experimental[row, "interleaving"] == "interleaving") {
    experimental[row, "INT"] <-"INT"
  } else experimental[row, "INT"] <- "DA"
}
#####

#### conductivity outliers
#####
outliersRemovedConduct <- conductivity[conductivity$measure == "condition" & conductivity$value < 0.5, ]
outliersRemovedExp <- experimental[experimental$conductivity < 0.5, ]
outliersRemovedConduct$paper_labels<-factor(outliersRemovedConduct$paper_labels, levels=c("Whatman","Antique Gelatin","Modern Gelatin", "Starch-based AKD", "Alum Rosin"))
outliersRemovedExp$paper_labels<-factor(outliersRemovedExp$paper_labels, levels=c("Whatman","Antique Gelatin","Modern Gelatin", "Starch-based AKD", "Alum Rosin"))
#####

#### figure 1 - delta e all samples bt-at, at-ad, bt-ad
#####
plot1 <- ggplot(colorimetry, aes(y=deltae, x=icondition)) +
  geom_jitter(size = 0.7) +
  geom_boxplot(alpha = 0.7) +
  theme_classic() + 
  ylab("Delta E") + xlab("Treatment Stage") + 
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =1, ymax = 2, alpha = .08) + 
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =2, ymax = 6, alpha = .13) + 
  geom_hline(yintercept=2, linetype='dotted') + 
  geom_hline(yintercept=1, linetype='dotted') +
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 45, hjust = 1, siz = 16), axis.title = element_text(size = 16), axis.text.y = element_text(size = 16))
plot1
#####

#### figure 2 - delta e bt-ad all samples with paper type on x and gel color
#####
gelColors <- c("#7030A0", "#0000FF", "#29D573", "#FF6F31", "#FFCF37")
plot2 <- ggplot(cumulativeChange, aes(y=deltae, x=paper_labels)) + 
  geom_jitter(aes(color = gel)) +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  theme_classic() + 
  ylab("Delta E: BT-AD") + xlab("Paper type") + 
  scale_color_manual(values = gelColors) +
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =1, ymax = 2, alpha = .08) + 
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =2, ymax = 6, alpha = .13) + 
  geom_hline(yintercept=2, linetype='dotted') + 
  geom_hline(yintercept=1, linetype='dotted') +
  theme(axis.text.x = element_text(angle = 45, hjust = 1, size = 16), axis.title = element_text(size = 16), text = element_text(size = 16)) +
  labs(color = "Gel") +
  guides(colour = guide_legend(override.aes = list(size=3)))
plot2
#####

#### figure 3 (or 3 and 4 separate) - aes(time point, pH, color = paper) facet by paper, both ph and conduct
#####
plot3 <- ggplot(experimental, aes(condition, ph)) +
  geom_jitter(aes(color = paper_labels)) +
  geom_boxplot() +
  facet_grid(cols = vars(paper_labels)) +
  theme_minimal() +
  labs(x = "Time point", y = "pH", color = "Paper") +
  theme(legend.position = "bottom", text = element_text(size = 16))
plot3
plotExplore <- ggplot(experimental, aes(condition, conductivity)) +
  geom_jitter(aes(color = paper_labels)) +
  geom_boxplot() +
  facet_grid(cols = vars(paper_labels)) +
  theme_minimal() +
  labs(x = "Time point", y = "Conductivity", color = "Paper") +
  theme(legend.position = "bottom")
plotExplore
plot4 <- ggplot(outliersRemovedConduct, aes(timing, value)) +
  geom_jitter(aes(color = paper_labels)) +
  geom_boxplot() +
  facet_grid(cols = vars(paper_labels)) +
  theme_minimal() +
  labs(x = "Time point", y = "Conductivity", color = "Paper") +
  theme(legend.position = "bottom", text = element_text(size = 16))
plot4
#####

#### figure 5 (order?) - l, a, b scatterplots
#####
paper_type = "Whatman Chromatography"
whatmanA <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(a, l, color = condition)) +
  geom_point(alpha = 0.8) + theme_classic() + 
  scale_fill_gradient(low = hex(LAB(90, 0, 0)), high = hex(LAB(90, 0, 0.4)), guide = "colourbar") +
  labs(x = "a*", y = "L*") + theme(legend.position = "none", text = element_text(size = 16))

whatmanB <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(b, l, color = condition)) +
  geom_point() + theme_classic() + 
  labs(x = "b*", y = "L*") + theme(legend.position = "none", text = element_text(size = 16))

whatmanAx <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(a, 0, color = a)) +
  geom_point() +
  scale_color_gradient(low = hex(LAB(90, 0, 0.25)), high = hex(LAB(90, 0.4, 0.25))) + theme_void() + theme(legend.position = "none")

whatmany <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(0, l, color = l)) +
  geom_point() + scale_color_gradient(low = hex(LAB(90, 0, 0.25)), high = hex(LAB(93, 0, 0.25))) + theme_void() + theme(legend.position = "none")

whatmanBx <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(b, 0, color = b)) +
  geom_point() +
  scale_color_gradient(low = hex(LAB(90, 0, 0.25)), high = hex(LAB(90, 0, 2.5))) + theme_void() + theme(legend.position = "none")

whatmanAlab <- ggarrange(whatmany, whatmanA, NA, whatmanAx, ncol = 2, nrow = 2, widths = c(0.1, 1), heights = c(1, 0.05), align = "hv")
whatmanBlab <- ggarrange(whatmany, whatmanB, NA, whatmanBx, ncol = 2, nrow = 2, widths = c(0.1, 1), heights = c(1, 0.05), align = "hv")

antiqueA <- ggplot(colorData[colorData$paper_type =="18c" & colorData$condition == "BT" | colorData$paper_type =="18c" & colorData$condition == "AD", ], aes(a, l, color = condition)) +
  geom_point(alpha = 0.8) + theme_classic() + 
  labs(x = "a*", y = "L*") + theme(legend.position = "none", text = element_text(size = 16))
antiqueB <- ggplot(colorData[colorData$paper_type =="18c" & colorData$condition == "BT" | colorData$paper_type =="18c" & colorData$condition == "AD", ], aes(b, l, color = condition)) +
  geom_point(alpha = 0.8) + theme_classic() + 
  labs(x = "b*", y = "L*") + theme(legend.position = "none", text = element_text(size = 16))

antiqueAx <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(a, 0, color = a)) +
  geom_point() +
  scale_color_gradient(low = hex(LAB(75, 1, 10)), high = hex(LAB(75, 6.5, 10))) + theme_void() + theme(legend.position = "none")

antiquey <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(0, l, color = l)) +
  geom_point() + scale_color_gradient(low = hex(LAB(75, 1, 10)), high = hex(LAB(90, 1, 10))) + theme_void() + theme(legend.position = "none")

antiqueBx <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(b, 0, color = b)) +
  geom_point() +
  scale_color_gradient(low = hex(LAB(75, 1, 10)), high = hex(LAB(75, 1, 20))) + theme_void() + theme(legend.position = "none")

antiqueAlab <- ggarrange(antiquey, antiqueA, NA, antiqueAx, ncol = 2, nrow = 2, widths = c(0.1, 1), heights = c(1, 0.05), align = "hv")
antiqueBlab <- ggarrange(antiquey, antiqueB, NA, antiqueBx, ncol = 2, nrow = 2, widths = c(0.1, 1), heights = c(1, 0.05), align = "hv")

paper_type = "Canson Ingres"
modernA <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(a, l, color = condition)) +
  geom_point() + theme_classic() + 
  labs(x = "a*", y = "L*") + theme(legend.position = "none", text = element_text(size = 16))
modernB <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(b, l, color = condition)) +
  geom_point() + theme_classic() + 
  labs(x = "b*", y = "L*") + theme(legend.position = "none", text = element_text(size = 16))

modernAx <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(a, 0, color = a)) +
  geom_point() +
  scale_color_gradient(low = hex(LAB(91, 0.23, 0.5)), high = hex(LAB(91, 0.62, 0.5))) + theme_void() + theme(legend.position = "none")

moderny <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(0, l, color = l)) +
  geom_point() + scale_color_gradient(low = hex(LAB(91, 0.23, 0.5)), high = hex(LAB(93.5, 0.23, 0.5))) + theme_void() + theme(legend.position = "none")

modernBx <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(b, 0, color = b)) +
  geom_point() +
  scale_color_gradient(low = hex(LAB(91, 0.23, 0.5)), high = hex(LAB(91, 0.23, 4.25))) + theme_void() + theme(legend.position = "none")

modernAlab <- ggarrange(moderny, modernA, NA, modernAx, ncol = 2, nrow = 2, widths = c(0.1, 1), heights = c(1, 0.05), align = "hv")
modernBlab <- ggarrange(moderny, modernB, NA, modernBx, ncol = 2, nrow = 2, widths = c(0.1, 1), heights = c(1, 0.05), align = "hv")

paper_type = "Stonehenge Watercolor"
starchA <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(a, l, color = condition)) +
  geom_point() + theme_classic() +
  labs(x = "a*", y = "L*") + theme(legend.position = "none", text = element_text(size = 16))
starchB <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(b, l, color = condition)) +
  geom_point() + theme_classic() + 
  labs(x = "b*", y = "L*") + theme(legend.position = "none", text = element_text(size = 16))

starchAx <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(a, 0, color = a)) +
  geom_point() +
  scale_color_gradient(low = hex(LAB(91.5, 0, 0)), high = hex(LAB(91.5, 0.6, 0))) + theme_void() + theme(legend.position = "none")

starchy <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(0, l, color = l)) +
  geom_point() + scale_color_gradient(low = hex(LAB(91.5, 0, 0)), high = hex(LAB(97, 0, 0))) + theme_void() + theme(legend.position = "none")

starchBx <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(b, 0, color = b)) +
  geom_point() +
  scale_color_gradient(low = hex(LAB(91.5, 0, 0)), high = hex(LAB(91.5, 0, 8.5))) + theme_void() + theme(legend.position = "none")

starchAlab <- ggarrange(starchy, starchA, NA, starchAx, ncol = 2, nrow = 2, widths = c(0.1, 1), heights = c(1, 0.05), align = "v")
starchBlab <- ggarrange(starchy, starchB, NA, starchBx, ncol = 2, nrow = 2, widths = c(0.1, 1), heights = c(1, 0.05), align = "v")

paper_type = "Alum Rosin"
alumA <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(a, l, color = condition)) +
  geom_point() + theme_classic() + 
  labs(x = "a*", y = "L*") + theme(legend.position = "none", text = element_text(size = 16))
alumB <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(b, l, color = condition)) +
  geom_point() + theme_classic() +
  labs(x = "b*", y = "L*") + theme(legend.position = "none", text = element_text(size = 16))

alumAx <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(a, 0, color = a)) +
  geom_point() +
  scale_color_gradient(low = hex(LAB(87.5, 0.4, 0.5)), high = hex(LAB(87.5, 0.4, 1.1))) + theme_void() + theme(legend.position = "none")

alumy <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(0, l, color = l)) +
  geom_point() + scale_color_gradient(low = hex(LAB(87.5, 0.4, 0.5)), high = hex(LAB(91.5, 0.4, 0.5))) + theme_void() + theme(legend.position = "none")

alumBx <- ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(b, 0, color = b)) +
  geom_point() +
  scale_color_gradient(low = hex(LAB(87.5, 0.4, 0.5)), high = hex(LAB(87.5, 0.4, 7))) + theme_void() + theme(legend.position = "none")

alumAlab <- ggarrange(alumy, alumA, NA, alumAx, ncol = 2, nrow = 2, widths = c(0.1, 1), heights = c(1, 0.05), align = "hv")
alumBlab <- ggarrange(alumy, alumB, NA, alumBx, ncol = 2, nrow = 2, widths = c(0.1, 1), heights = c(1, 0.05), align = "hv")

colorLegend <- get_legend(ggplot(colorData[colorData$paper_type ==paper_type & colorData$condition == "BT" | colorData$paper_type ==paper_type & colorData$condition == "AD", ], aes(a, l, color = condition)) +
                            geom_point() + theme_classic() + theme(legend.position = "bottom", text = element_text(size = 16)))
plot5 <- ggarrange(whatmanAlab, whatmanBlab, antiqueAlab, antiqueBlab, modernAlab, modernBlab, starchAlab, starchBlab, alumAlab, alumBlab, colorLegend, nrow = 6, ncol = 2, labels = c("Whatman", "Whatman", "Antique Gelatin", "Antique Gelatin", "Modern Gelatin", "Modern Gelatin", "Starch-based AKD", "Starch-based AKD", "Alum Rosin", "Alum Rosin"), label.y = 0, vjust = -1, hjust = c(-0.5, -0.5, -0.4, -0.4, -0.4, -0.4, -0.3, -0.3, -0.5, -0.5), font.label = list(size = 20), heights = c(1, 1, 1, 1, 1, 0.3))

#####

#### calculate stats l, a, b BT to AD
##### 
whatman <- colorData[colorData$paper_type == "Whatman Chromatography", ]
whatman <- whatman[whatman$condition != "AT", ]

parameters <- c("l", "a", "b")
BTmeans <- list()
AAmeans <- list()
pvalues <- list()
for (p in parameters) {
  BTmeans[p] <- mean(whatman[whatman$condition == "BT", p])
  AAmeans[p] <- mean(whatman[whatman$condition == "AD", p])
  ktests <- kruskal.test(whatman[, p] ~ whatman$condition)
  pvalues[p] <- ktests$p.value
}

whatmanStats <- cbind(parameters, as.numeric(BTmeans), as.numeric(AAmeans), as.numeric(pvalues))
antique <- colorData[colorData$paper_type == "18c", ]
antique <- antique[antique$condition != "AT", ]

parameters <- c("l", "a", "b")
BTmeans <- list()
AAmeans <- list()
pvalues <- list()
for (p in parameters) {
  BTmeans[p] <- mean(antique[antique$condition == "BT", p])
  AAmeans[p] <- mean(antique[antique$condition == "AD", p])
  ktests <- kruskal.test(antique[, p] ~ antique$condition)
  pvalues[p] <- ktests$p.value
}

antiqueStats <- cbind(parameters, as.numeric(BTmeans), as.numeric(AAmeans), as.numeric(pvalues))

modern <- colorData[colorData$paper_type == "Canson Ingres", ]
modern <- modern[modern$condition != "AT", ]

parameters <- c("l", "a", "b")
BTmeans <- list()
AAmeans <- list()
pvalues <- list()
for (p in parameters) {
  BTmeans[p] <- mean(modern[modern$condition == "BT", p])
  AAmeans[p] <- mean(modern[modern$condition == "AD", p])
  ktests <- kruskal.test(modern[, p] ~ modern$condition)
  pvalues[p] <- ktests$p.value
}

modernStats <- cbind(parameters, as.numeric(BTmeans), as.numeric(AAmeans), as.numeric(pvalues))

starch <- colorData[colorData$paper_type == "Stonehenge Watercolor", ]
starch <- starch[starch$condition != "AT", ]

parameters <- c("l", "a", "b")
BTmeans <- list()
AAmeans <- list()
pvalues <- list()
for (p in parameters) {
  BTmeans[p] <- mean(starch[starch$condition == "BT", p])
  AAmeans[p] <- mean(starch[starch$condition == "AD", p])
  ktests <- kruskal.test(starch[, p] ~ starch$condition)
  pvalues[p] <- ktests$p.value
}

starchStats <- cbind(parameters, as.numeric(BTmeans), as.numeric(AAmeans), as.numeric(pvalues))

alum <- colorData[colorData$paper_type == "Alum Rosin", ]
alum <- alum[alum$condition != "AT", ]

parameters <- c("l", "a", "b")
BTmeans <- list()
AAmeans <- list()
pvalues <- list()
for (p in parameters) {
  BTmeans[p] <- mean(alum[alum$condition == "BT", p])
  AAmeans[p] <- mean(alum[alum$condition == "AD", p])
  ktests <- kruskal.test(alum[, p] ~ alum$condition)
  pvalues[p] <- ktests$p.value
}

alumStats <- cbind(parameters, as.numeric(BTmeans), as.numeric(AAmeans), as.numeric(pvalues))


#### stats table
#####
allStats <- rbind(whatmanStats, antiqueStats, modernStats, starchStats, alumStats)
papertypelabels <- c("Whatman", "Whatman", "Whatman", "Antique Gelatin", "Antique Gelatin", "Antique Gelatin", "Modern Gelatin", "Modern Gelatin", "Modern Gelatin", "Starch-based AKD", "Starch-based AKD", "Starch-based AKD", "Alum Rosin", "Alum Rosin", "Alum Rosin")
allStats <- cbind(papertypelabels, allStats)
allStats <- as.data.frame(allStats)
names(allStats) <- c("Paper sizing", "Color parameter", "BT mean", "AD mean", "p value (Kruskal-Wallis test)")
allStats$`Color parameter` <- c("L*", "a*", "b*", "L*", "a*", "b*", "L*", "a*", "b*", "L*", "a*", "b*", "L*", "a*", "b*")
allStats$`BT mean` <- round(as.numeric(allStats$`BT mean`), 2)
allStats$`AD mean` <- round(as.numeric(allStats$`AD mean`), 2)
allStats$`p value (Kruskal-Wallis test)` <- format(as.numeric(allStats$`p value (Kruskal-Wallis test)`), digits = 2, justify = "right")

allStatsTable <- allStats %>% 
  kbl() %>% 
  kable_classic()
#####

#### full page figures (5 of them) - all AD responses for each papertype
# whatman
#####
paper_type = "Whatman Chromatography"
paper_label = "Whatman"
whatmanPLOT <- ggplot(data=cumulativeChange[cumulativeChange$paper_type ==  "Whatman Chromatography", ], aes(y=deltae, x=gel_treatment)) + 
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =1, ymax = 2, alpha = .08) + 
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =2, ymax = 6.2, alpha = .13) + 
  geom_hline(yintercept=2, linetype='dotted') + 
  geom_hline(yintercept=1, linetype='dotted') +
  ylab("Delta E") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
whatmanLplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, l)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("L*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
whatmanAplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, a)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("a*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
whatmanBplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, b)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("b*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
whatmanPHplot <- ggplot(experimental[experimental$paper_labels == paper_label & experimental$condition == "AD", ], aes(INT, ph)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("pH") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
whatmanCplot <- ggplot(experimental[experimental$paper_labels == paper_label & experimental$condition == "AD", ], aes(INT, conductivity)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("Conductivity") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.title.x = element_blank(), text = element_text(size = 16))

whatmanPAGE <- ggarrange(whatmanPLOT, whatmanLplot, whatmanAplot, whatmanBplot, whatmanPHplot, whatmanCplot, ncol = 1, common.legend = T, align = "hv")
#####
# antique
#####
paper_label = "Antique Gelatin"
deltaePLOT <- ggplot(data=cumulativeChange[cumulativeChange$paper_labels ==  paper_label, ], aes(y=deltae, x=gel_treatment)) + 
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =1, ymax = 2, alpha = .08) + 
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =2, ymax = 6.2, alpha = .13) + 
  geom_hline(yintercept=2, linetype='dotted') + 
  geom_hline(yintercept=1, linetype='dotted') +
  ylab("Delta E") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Lplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, l)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("L*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Aplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, a)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("a*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Bplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, b)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("b*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
PHplot <- ggplot(experimental[experimental$paper_labels == paper_label & experimental$condition == "AD", ], aes(INT, ph)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("pH") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Cplot <- ggplot(experimental[experimental$paper_labels == paper_label & experimental$condition == "AD", ], aes(INT, conductivity)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("Conductivity") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.title.x = element_blank(), text = element_text(size = 16))

antiquePAGE <- ggarrange(deltaePLOT, Lplot, Aplot, Bplot, PHplot, Cplot, ncol = 1, common.legend = T, align = "hv")
antiquePAGE
#####
# modern
#####
paper_label = "Modern Gelatin"
deltaePLOT <- ggplot(data=cumulativeChange[cumulativeChange$paper_labels ==  paper_label, ], aes(y=deltae, x=gel_treatment)) + 
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =1, ymax = 2, alpha = .08) + 
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =2, ymax = 6.2, alpha = .13) + 
  geom_hline(yintercept=2, linetype='dotted') + 
  geom_hline(yintercept=1, linetype='dotted') +
  ylab("Delta E") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Lplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, l)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("L*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Aplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, a)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("a*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Bplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, b)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("b*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
PHplot <- ggplot(experimental[experimental$paper_labels == paper_label & experimental$condition == "AD", ], aes(INT, ph)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("pH") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Cplot <- ggplot(experimental[experimental$paper_labels == paper_label & experimental$condition == "AD", ], aes(INT, conductivity)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("Conductivity") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.title.x = element_blank(), text = element_text(size = 16))

modernPAGE <- ggarrange(deltaePLOT, Lplot, Aplot, Bplot, PHplot, Cplot, ncol = 1, common.legend = T, align = "hv")
#####
# starch
#####
paper_label = "Starch-based AKD"
deltaePLOT <- ggplot(data=cumulativeChange[cumulativeChange$paper_labels ==  paper_label, ], aes(y=deltae, x=gel_treatment)) + 
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =1, ymax = 2, alpha = .08) + 
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =2, ymax = 6.2, alpha = .13) + 
  geom_hline(yintercept=2, linetype='dotted') + 
  geom_hline(yintercept=1, linetype='dotted') +
  ylab("Delta E") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Lplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, l)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("L*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Aplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, a)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("a*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Bplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, b)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("b*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
PHplot <- ggplot(experimental[experimental$paper_labels == paper_label & experimental$condition == "AD", ], aes(INT, ph)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("pH") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Cplot <- ggplot(experimental[experimental$paper_labels == paper_label & experimental$condition == "AD", ], aes(INT, conductivity)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("Conductivity") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.title.x = element_blank(), text = element_text(size = 16))

starchPAGE <- ggarrange(deltaePLOT, Lplot, Aplot, Bplot, PHplot, Cplot, ncol = 1, common.legend = T, align = "hv")
#####
# alum
#####
paper_label = "Alum Rosin"
deltaePLOT <- ggplot(data=cumulativeChange[cumulativeChange$paper_labels ==  paper_label, ], aes(y=deltae, x=gel_treatment)) + 
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =1, ymax = 2, alpha = .08) + 
  ggplot2::annotate("rect",xmin = -Inf, xmax = Inf, ymin =2, ymax = 6.2, alpha = .13) + 
  geom_hline(yintercept=2, linetype='dotted') + 
  geom_hline(yintercept=1, linetype='dotted') +
  ylab("Delta E") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Lplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, l)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("L*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Aplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, a)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("a*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Bplot <- ggplot(colorData[colorData$paper_labels == paper_label & colorData$condition == "AD", ], aes(INT, b)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("b*") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
PHplot <- ggplot(experimental[experimental$paper_labels == paper_label & experimental$condition == "AD", ], aes(INT, ph)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("pH") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.text.x = element_blank(), axis.title.x = element_blank(), text = element_text(size = 16))
Cplot <- ggplot(experimental[experimental$paper_labels == paper_label & experimental$condition == "AD", ], aes(INT, conductivity)) +
  theme_classic() + 
  geom_jitter(width = 0.3, aes(color = INT))  +
  geom_boxplot(na.rm=TRUE, alpha = 0.7) + 
  ylab("Conductivity") +
  facet_grid(cols = vars(gel), scales = "free") +
  theme(legend.position = "bottom", axis.title.x = element_blank(), text = element_text(size = 16))

alumPAGE <- ggarrange(deltaePLOT, Lplot, Aplot, Bplot, PHplot, Cplot, ncol = 1, common.legend = T, align = "hv")
#####

#### linear models for AD each paper type
# delta e
#####
whatmanModel <- lm(deltae~gel_treatment, cumulativeChange[cumulativeChange$paper_labels == "Whatman", ])
antiqueModel <- lm(deltae~gel_treatment, cumulativeChange[cumulativeChange$paper_labels == "Antique Gelatin", ])
modernModel <- lm(deltae~gel_treatment, cumulativeChange[cumulativeChange$paper_labels == "Modern Gelatin", ])
starchModel <- lm(deltae~gel_treatment, cumulativeChange[cumulativeChange$paper_labels == "Starch-based AKD", ])
alumModel <- lm(deltae~gel_treatment, cumulativeChange[cumulativeChange$paper_labels == "Alum Rosin", ])

#####
# L*
#####
whatmanLModel <- lm(l~gel_treatment, colorData[colorData$paper_labels == "Whatman" & colorData$condition == "AD", ])
antiqueLModel <- lm(l~gel_treatment, colorData[colorData$paper_labels == "Antique Gelatin" & colorData$condition == "AD", ])
modernLModel <- lm(l~gel_treatment, colorData[colorData$paper_labels == "Modern Gelatin" & colorData$condition == "AD", ])
starchLModel <- lm(l~gel_treatment, colorData[colorData$paper_labels == "Starch-based AKD" & colorData$condition == "AD", ])
alumLModel <- lm(l~gel_treatment, colorData[colorData$paper_labels == "Alum Rosin" & colorData$condition == "AD", ])
#####
# b*
#####
whatmanBModel <- lm(b~gel_treatment, colorData[colorData$paper_labels == "Whatman" & colorData$condition == "AD", ])
antiqueBModel <- lm(b~gel_treatment, colorData[colorData$paper_labels == "Antique Gelatin" & colorData$condition == "AD", ])
modernBModel <- lm(b~gel_treatment, colorData[colorData$paper_labels == "Modern Gelatin" & colorData$condition == "AD", ])
starchBModel <- lm(b~gel_treatment, colorData[colorData$paper_labels == "Starch-based AKD" & colorData$condition == "AD", ])
alumBModel <- lm(b~gel_treatment, colorData[colorData$paper_labels == "Alum Rosin" & colorData$condition == "AD", ])
#####
# a*
#####
whatmanAModel <- lm(a~gel_treatment, colorData[colorData$paper_labels == "Whatman" & colorData$condition == "AD", ])
antiqueAModel <- lm(a~gel_treatment, colorData[colorData$paper_labels == "Antique Gelatin" & colorData$condition == "AD", ])
modernAModel <- lm(a~gel_treatment, colorData[colorData$paper_labels == "Modern Gelatin" & colorData$condition == "AD", ])
starchAModel <- lm(a~gel_treatment, colorData[colorData$paper_labels == "Starch-based AKD" & colorData$condition == "AD", ])
alumAModel <- lm(a~gel_treatment, colorData[colorData$paper_labels == "Alum Rosin" & colorData$condition == "AD", ])
#####
# pH and conductivity
#####
whatmanPHModel <- lm(ph~gel_treatment, experimental[experimental$paper_labels == "Whatman" & experimental$condition == "AD", ])
antiquePHModel <- lm(ph~gel_treatment, experimental[experimental$paper_labels == "Antique Gelatin" & experimental$condition == "AD", ])
modernPHModel <- lm(ph~gel_treatment, experimental[experimental$paper_labels == "Modern Gelatin" & experimental$condition == "AD", ])
starchPHModel <- lm(ph~gel_treatment, experimental[experimental$paper_labels == "Starch-based AKD" & experimental$condition == "AD", ])
alumPHModel <- lm(ph~gel_treatment, experimental[experimental$paper_labels == "Alum Rosin" & experimental$condition == "AD", ])

whatmanCModel <- lm(conductivity~gel_treatment, experimental[experimental$paper_labels == "Whatman" & experimental$condition == "AD", ])
antiqueCModel <- lm(conductivity~gel_treatment, experimental[experimental$paper_labels == "Antique Gelatin" & experimental$condition == "AD", ])
modernCModel <- lm(conductivity~gel_treatment, experimental[experimental$paper_labels == "Modern Gelatin" & experimental$condition == "AD", ])
starchCModel <- lm(conductivity~gel_treatment, experimental[experimental$paper_labels == "Starch-based AKD" & experimental$condition == "AD", ])
alumCModel <- lm(conductivity~gel_treatment, experimental[experimental$paper_labels == "Alum Rosin" & experimental$condition == "AD", ])
#####

#linear model tables
#####
whatmanDeltaE <- whatmanModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
whatmanL <- whatmanLModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
whatmanA <- whatmanAModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
whatmanB <- whatmanBModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
whatmanpH <- whatmanPHModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
whatmanC <- whatmanCModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()

antiqueDeltaE <- antiqueModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
antiqueL <- antiqueLModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
antiqueA <- antiqueAModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
antiqueB <- antiqueBModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
antiquepH <- antiquePHModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
antiqueC <- antiqueCModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()

modernDeltaE <- modernModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
modernL <- modernLModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
modernA <- modernAModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
modernB <- modernBModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
modernpH <- modernPHModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
modernC <- modernCModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()

starchDeltaE <- starchModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
starchL <- starchLModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
starchA <- starchAModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
starchB <- starchBModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
starchpH <- starchPHModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
starchC <- starchCModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()

alumDeltaE <- alumModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
alumL <- alumLModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
alumA <- alumAModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
alumB <- alumBModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
alumpH <- alumPHModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
alumC <- alumCModel %>% tbl_regression() %>% 
  bold_p(t = 0.10) %>%
  bold_labels() %>%
  italicize_levels()
#####

# save the plots
#####
tiff(filename = "figures/plot1.tiff", height = 800, width = 800)
plot1
dev.off()

tiff(filename = "figures/plot2.tiff", height = 800, width = 800)
plot2
dev.off()

tiff(filename = "figures/plot3.tiff", height =800, width = 800)
plot3
dev.off()

tiff(filename = "figures/plot4.tiff", height =800, width = 800)
plot4
dev.off()

tiff(filename = "figures/plot5.tiff", height = 1000, width = 1000)
plot5
dev.off()

tiff(filename = "figures/plot6Whatman.tiff", height = 800, width = 480)
whatmanPAGE
dev.off()

tiff(filename = "figures/plot7antique.tiff", height = 800, width = 480)
antiquePAGE
dev.off()

tiff(filename = "figures/plot8modern.tiff", height = 800, width = 480)
modernPAGE
dev.off()

tiff(filename = "figures/plot9starch.tiff", height = 800, width = 480)
starchPAGE
dev.off()

tiff(filename = "figures/plot10alum.tiff", height = 800, width = 480)
alumPAGE
dev.off()
#####