1 Load packages

library(ggplot2)
library(GGally)
library(tidyverse)
library(gridExtra)
library(knitr)
library(ggpubr)
library(ggthemes)
theme_set(theme_classic())# ggpubr::theme_pubr())

2 Upload data

We will use two datasets, at two levels of aggregation.

  • The urine_dat dataset is the highest resolution at the timepoint level. It includes experimental design factors and urine variables only. This will enable us to look at diurnal variation, and variation related to different sampling timepoints.

  • The urine_dat_d dataset is at the day level of aggregation and includes experimental design factors AND both urine and milk related data. This will ease interpretation of day-to-day values which are commonly reported in the literature.

urine_dat_d <- read.csv("00_input/urine.dat.d.export.csv") %>%
  select(-X) %>%
  rename(MUN = Urea) %>%
  mutate(cow = as.character(cow))

urine_dat <- read.csv("00_input/urine.dat.export.csv") %>%
  select(-X) %>%
  mutate(hr = factor(hr, levels = c("4AM", "12PM", "8PM"))) %>%
  mutate(day_hr = paste0(day, "_", hr)) %>%
  mutate(day_hr = factor(day_hr, levels = c(
    "Day1_4AM", "Day1_12PM", "Day1_8PM",
    "Day2_4AM", "Day2_12PM", "Day2_8PM",
    "Day3_4AM", "Day3_12PM", "Day3_8PM",
    "Day4_4AM", "Day4_12PM", "Day4_8PM"
  ))) %>%
  mutate(cow = as.character(cow))

tc_n_imp <- read.csv("01_organize/tc_n_imp.csv") %>% 
  mutate(hr = factor(hr, levels = c("4AM", "12PM", "8PM"))) %>%
  mutate(day_hr = paste0(day, "_", hr)) %>%
  mutate(day_hr = factor(day_hr, levels = c(
    "Day1_4AM", "Day1_12PM", "Day1_8PM",
    "Day2_4AM", "Day2_12PM", "Day2_8PM",
    "Day3_4AM", "Day3_12PM", "Day3_8PM",
    "Day4_4AM", "Day4_12PM", "Day4_8PM"
  ))) #%>% 
 # dplyr::select(-X, -X.1)

2.1 Timepoint-level (urine_dat)

kable(summary(urine_dat)[1:3, 1:7])
period day hr cow osc_lvl prot_lvl trt
Length:360 Length:360 4AM :120 Length:360 Length:360 Length:360 Length:360
Class :character Class :character 12PM:120 Class :character Class :character Class :character Class :character
Mode :character Mode :character 8PM :120 Mode :character Mode :character Mode :character Mode :character 
kable(summary(urine_dat)[1:6, 8:12])
urine_kg urine_ph urine_n_gpertp urine_urea_n_gpertp urine_creatinine_mgpertp
Min. : 1.251 Min. :-0.01139 Min. : 13.81 Min. : 1.831 Min. :1136
1st Qu.: 7.151 1st Qu.: 1.89450 1st Qu.: 45.57 1st Qu.:28.299 1st Qu.:5024
Median : 8.501 Median : 2.54830 Median : 57.96 Median :38.256 Median :5844
Mean : 8.954 Mean : 2.84197 Mean : 59.78 Mean :39.599 Mean :5740
3rd Qu.:10.213 3rd Qu.: 3.61050 3rd Qu.: 73.10 3rd Qu.:50.267 3rd Qu.:6544
Max. :19.451 Max. : 6.57900 Max. :122.07 Max. :93.238 Max. :9375

2.2 Day-level (urine_dat_d)

kable(summary(urine_dat_d)[1:3, 1:7])
period day cow osc_lvl prot_lvl trt n_intake_gperd
Length:120 Length:120 Length:120 Length:120 Length:120 Length:120 Min. :372.6
Class :character Class :character Class :character Class :character Class :character Class :character 1st Qu.:524.1
Mode :character Mode :character Mode :character Mode :character Mode :character Mode :character Median :562.1 
kable(summary(urine_dat_d)[, 7:14])
n_intake_gperd urine_kgperd urine_n_gperd urine_urea_n_gperd urine_creatinine_mgperd milk_n_gperd MUNY_gperd MUN
Min. :372.6 Min. :15.10 Min. : 78.74 Min. : 44.54 Min. : 7514 Min. :122.6 Min. :1.472 Min. : 4.520
1st Qu.:524.1 1st Qu.:22.38 1st Qu.:145.23 1st Qu.: 89.33 1st Qu.:15063 1st Qu.:155.3 1st Qu.:3.285 1st Qu.: 8.489
Median :562.1 Median :24.98 Median :173.85 Median :112.86 Median :17639 Median :175.6 Median :3.845 Median : 9.977
Mean :588.8 Mean :26.86 Mean :179.35 Mean :118.80 Mean :17220 Mean :177.1 Mean :3.953 Mean :10.272
3rd Qu.:640.0 3rd Qu.:29.95 3rd Qu.:216.75 3rd Qu.:147.52 3rd Qu.:19339 3rd Qu.:195.0 3rd Qu.:4.791 3rd Qu.:11.819
Max. :963.4 Max. :51.00 Max. :293.03 Max. :212.29 Max. :22755 Max. :250.1 Max. :7.437 Max. :17.028

3 Timepoint Level

Univariate plots of timepoint-level data showing raw values over time.

3.0.1 Urine pH

Nothing remarkable in urine pH plots by day and timepoint. Plots omitted.

3.0.2 Urine production (kg)

urine_dat %>%
  ggplot(aes(x = day_hr, y = urine_kg)) +
  labs(title = "Urine production (kg) per timepoint") +
  geom_point(
    alpha = 0.5, position = position_jitterdodge(jitter.width = 0.15),
    aes(col = prot_lvl, group = prot_lvl)
  ) +
  geom_boxplot(aes(color = prot_lvl), fill = NA, 
               position = position_dodge(width = 0.75), outlier.shape = NA) +
  geom_violin(aes(fill = prot_lvl), color = NA, alpha = 0.2, 
              position = position_dodge(width = 0.75)) +
  facet_grid("osc_lvl") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  theme(axis.title.x = element_blank()) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1)+
  annotate("rect", xmin = -Inf, xmax = 1.5, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") +
  annotate("rect", xmin = 7.5, xmax = Inf, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") + 
  annotate("text", x = 1.5, y = 15, label = "FED HIGH", angle = 90, size = 2)+ 
  annotate("text", x = 4.5, y = 15, label = "FED HIGH", angle = 90, size = 2)+
  annotate("text", x = 7.5, y = 15, label = "FED LOW", angle = 90, size = 2)+ 
  annotate("text", x = 10.5, y = 15, label = "FED LOW", angle = 90, size = 2)

3.0.2.1 Highlight 7182 urine production (kg)

Looking at the same data as above, with boxplots removed and a single cow’s observations highlighted. The urine production (kg) for 7182 is consistently greater than the other cows. This cow is also the largest bodyweight of the group.

urine_dat %>%
  mutate(is7182 = ifelse(cow == "7182", "7182", "Not7182")) %>%
  ggplot(aes(x = day_hr, y = urine_kg)) +
  labs(title = "Urine production (kg) per timepoint") +
  geom_point(
    alpha = 0.5, position = position_jitterdodge(jitter.width = 0.15),
    aes(col = is7182, group = prot_lvl)
  ) +
  facet_grid("osc_lvl") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  theme(axis.title.x = element_blank()) +
  scale_colour_manual(values = c("7182" = "orange")) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1)

3.0.3 Urine Urea N

THere is a slight increase in urine urea N is visible in Day 2-3, which may be a delayed effect of the HIGH-CP intake in Day 1-2. There also may be some smaller diurnal or sampling-time-related variation visible in the static data, where the 4 AM timepoint is greater than 12 PM in somem cases.

urine_dat %>%
  ggplot(aes(x = day_hr, y = urine_urea_n_gpertp)) +
  labs(title = "Urine Urea N (g) per timepoint") +
  geom_point(
    alpha = 0.5, position = position_jitterdodge(jitter.width = 0.15),
    aes(col = prot_lvl, group = prot_lvl)
  ) +
  geom_boxplot(aes(color = prot_lvl), fill = NA, 
               position = position_dodge(width = 0.75), outlier.shape = NA) +
  geom_violin(aes(fill = prot_lvl), color = NA, alpha = 0.2, 
              position = position_dodge(width = 0.75)) +
  facet_grid("osc_lvl", scales = "free") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  geom_hline(yintercept = -1) +
  scale_y_continuous(limits = c(0, 100)) +
  theme(axis.title.x = element_blank(), axis.line = element_line()) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1)+
  annotate("rect", xmin = -Inf, xmax = 1.5, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") +
  annotate("rect", xmin = 7.5, xmax = Inf, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") + 
  annotate("text", x = 1.5, y = 100, label = "FED HIGH", angle = 90, size = 2)+ 
  annotate("text", x = 4.5, y = 100, label = "FED HIGH", angle = 90, size = 2)+
  annotate("text", x = 7.5, y = 100, label = "FED LOW", angle = 90, size = 2)+ 
  annotate("text", x = 10.5, y = 100, label = "FED LOW", angle = 90, size = 2)

3.0.4 Urine N

Urine N appears to be correlated with urine urea N, as expected. Again, a slight increase in urine N is visible in Day 2-3, which may be a delayed effect of the HIGH-CP intake in Day 1-2.

urine_dat %>%
  ggplot(aes(x = day_hr, y = urine_n_gpertp)) +
  labs(title = "Urine N (g) per timepoint") +
  geom_point(
    alpha = 0.5, position = position_jitterdodge(jitter.width = 0.15),
    aes(col = prot_lvl, group = prot_lvl)
  ) +
  geom_boxplot(aes(color = prot_lvl), fill = NA, position = position_dodge(width = 0.75), outlier.shape = NA) +
  geom_violin(aes(fill = prot_lvl), color = NA, alpha = 0.2, position = position_dodge(width = 0.75)) +
  facet_grid("osc_lvl") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  theme(axis.title.x = element_blank()) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1)+ 
  annotate("rect", xmin = -Inf, xmax = 1.5, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") +
  annotate("rect", xmin = 7.5, xmax = Inf, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") + 
  annotate("text", x = 1.5, y = 100, label = "FED HIGH", angle = 90, size = 2)+ 
  annotate("text", x = 4.5, y = 100, label = "FED HIGH", angle = 90, size = 2)+
  annotate("text", x = 7.5, y = 100, label = "FED LOW", angle = 90, size = 2)+ 
  annotate("text", x = 10.5, y = 100, label = "FED LOW", angle = 90, size = 2)

3.0.5 Urine creatinine

Creatinine is expected to be 27.3 mg/kg of BW based on Tebbe (2018). With mean 667 kg BW, the average should be around 6,069 per timepoint (shown with horizontal line in figure). Creatinine appears to be mostly stable across time regardless of oscillation phase.

urine_dat %>%
  ggplot(aes(x = day_hr, y = urine_creatinine_mgpertp)) +
  labs(title = "Urine creatinine (mg) per timepoint") +
  geom_point(
    alpha = 0.5, position = position_jitterdodge(jitter.width = 0.15),
    aes(col = prot_lvl, group = prot_lvl)
  ) +
  geom_boxplot(aes(color = prot_lvl), fill = NA, position = position_dodge(width = 0.75), alpha = .2, outlier.shape = NA) +
  geom_violin(aes(fill = prot_lvl), color = NA, alpha = 0.2, position = position_dodge(width = 0.75)) +
  facet_grid("osc_lvl") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  theme(axis.title.x = element_blank()) +
  geom_hline(yintercept = 6069, color = "blue", lwd = 0.5, linetype = "solid") +
  geom_text(aes(9, 9400, label = "Line shows expected creatinine, 27.3 mg/kg BW"), size = 3, check_overlap = TRUE, color = "blue") +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1)

3.1 Feces output (kg)

tc_n_imp %>%
  ggplot(aes(x = day_hr, y = feces_kg)) +
  labs(title = "Feces output (kg)") +
  geom_point(
    alpha = 0.5, position = position_jitterdodge(jitter.width = 0.15),
    aes(col = prot_lvl, group = prot_lvl)
  ) +
  geom_boxplot(aes(color = prot_lvl), fill = NA, position = position_dodge(width = 0.75), outlier.shape = NA) +
  geom_violin(aes(fill = prot_lvl), color = NA, alpha = 0.2, position = position_dodge(width = 0.75)) +
  facet_grid("osc_lvl") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  theme(axis.title.x = element_blank()) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1)+ 
  annotate("rect", xmin = -Inf, xmax = 1.5, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") +
  annotate("rect", xmin = 7.5, xmax = Inf, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") + 
  annotate("text", x = 1.5, y = 30, label = "FED HIGH", angle = 90, size = 2)+ 
  annotate("text", x = 4.5, y = 30, label = "FED HIGH", angle = 90, size = 2)+
  annotate("text", x = 7.5, y = 30, label = "FED LOW", angle = 90, size = 2)+ 
  annotate("text", x = 10.5, y = 30, label = "FED LOW", angle = 90, size = 2)

3.2 Feces DM output (kg)

tc_n_imp %>%
  ggplot(aes(x = day_hr, y = feces_dm_kg)) +
  labs(title = "Feces DM output (kg)") +
  geom_point(
    alpha = 0.5, position = position_jitterdodge(jitter.width = 0.15),
    aes(col = prot_lvl, group = prot_lvl)
  ) +
  geom_boxplot(aes(color = prot_lvl), fill = NA, position = position_dodge(width = 0.75), outlier.shape = NA) +
  geom_violin(aes(fill = prot_lvl), color = NA, alpha = 0.2, position = position_dodge(width = 0.75)) +
  facet_grid("osc_lvl") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  theme(axis.title.x = element_blank()) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1)+
  annotate("rect", xmin = -Inf, xmax = 1.5, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") +
  annotate("rect", xmin = 7.5, xmax = Inf, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") + 
  annotate("text", x = 1.5, y = 3, label = "FED HIGH", angle = 90, size = 2)+ 
  annotate("text", x = 4.5, y = 3, label = "FED HIGH", angle = 90, size = 2)+
  annotate("text", x = 7.5, y = 3, label = "FED LOW", angle = 90, size = 2)+ 
  annotate("text", x = 10.5, y = 3, label = "FED LOW", angle = 90, size = 2)

3.3 Fecal N (g)

tc_n_imp %>%
  ggplot(aes(x = day_hr, y = feces_n_gpertp)) +
  labs(title = "Feces N output (g)") +
  geom_point(
    alpha = 0.5, position = position_jitterdodge(jitter.width = 0.15),
    aes(col = prot_lvl, group = prot_lvl)
  ) +
  geom_boxplot(aes(color = prot_lvl), fill = NA, position = position_dodge(width = 0.75), outlier.shape = NA) +
  geom_violin(aes(fill = prot_lvl), color = NA, alpha = 0.2, position = position_dodge(width = 0.75)) +
  facet_grid("osc_lvl") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  theme(axis.title.x = element_blank()) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1)+
  annotate("rect", xmin = -Inf, xmax = 1.5, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") +
  annotate("rect", xmin = 7.5, xmax = Inf, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") + 
  annotate("text", x = 1.5, y = 90, label = "FED HIGH", angle = 90, size = 2)+ 
  annotate("text", x = 4.5, y = 90, label = "FED HIGH", angle = 90, size = 2)+
  annotate("text", x = 7.5, y = 90, label = "FED LOW", angle = 90, size = 2)+ 
  annotate("text", x = 10.5, y = 90, label = "FED LOW", angle = 90, size = 2)

3.4 Manure output (kg)

tc_n_imp %>%
  ggplot(aes(x = day_hr, y = manure_kgpertp)) +
  labs(title = "Manure output (kg)") +
  geom_point(
    alpha = 0.5, position = position_jitterdodge(jitter.width = 0.15),
    aes(col = prot_lvl, group = prot_lvl)
  ) +
  geom_boxplot(aes(color = prot_lvl), fill = NA, position = position_dodge(width = 0.75), outlier.shape = NA) +
  geom_violin(aes(fill = prot_lvl), color = NA, alpha = 0.2, position = position_dodge(width = 0.75)) +
  facet_grid("osc_lvl") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  theme(axis.title.x = element_blank()) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1)+
  annotate("rect", xmin = -Inf, xmax = 1.5, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") +
  annotate("rect", xmin = 7.5, xmax = Inf, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") + 
  annotate("text", x = 1.5, y = 45, label = "FED HIGH", angle = 90, size = 2)+ 
  annotate("text", x = 4.5, y = 45, label = "FED HIGH", angle = 90, size = 2)+
  annotate("text", x = 7.5, y = 45, label = "FED LOW", angle = 90, size = 2)+ 
  annotate("text", x = 10.5, y = 45, label = "FED LOW", angle = 90, size = 2)

3.5 Manure N output

tc_n_imp %>%
  ggplot(aes(x = day_hr, y = manure_n_gpertp)) +
  labs(title = "Manure N output (g)") +
  geom_point(
    alpha = 0.5, position = position_jitterdodge(jitter.width = 0.15),
    aes(col = prot_lvl, group = prot_lvl)
  ) +
  geom_boxplot(aes(color = prot_lvl), fill = NA, position = position_dodge(width = 0.75), outlier.shape = NA) +
  geom_violin(aes(fill = prot_lvl), color = NA, alpha = 0.2, position = position_dodge(width = 0.75)) +
  facet_grid("osc_lvl") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  theme(axis.title.x = element_blank()) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1)+
  annotate("rect", xmin = -Inf, xmax = 1.5, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") +
  annotate("rect", xmin = 7.5, xmax = Inf, ymin = -Inf, ymax = Inf, alpha = .3, fill = "lightblue") + 
  annotate("text", x = 1.5, y = 180, label = "FED HIGH", angle = 90, size = 2)+ 
  annotate("text", x = 4.5, y = 180, label = "FED HIGH", angle = 90, size = 2)+
  annotate("text", x = 7.5, y = 180, label = "FED LOW", angle = 90, size = 2)+ 
  annotate("text", x = 10.5, y = 180, label = "FED LOW", angle = 90, size = 2)

4 Day Level

4.1 Univariate and bivariate

Visualize univariate with density, and bivariate with scatter plot matrices.

4.1.1 Pairs by period

Shortening variable names (shown in code) so they fit on figure. We can see from the density plots on the diagonal Period 1 was somewhat different than P2 through P4 for several of the variables here. However, most of the relationships appear to be similar across periods.

urine_dat_d <- urine_dat_d %>% 
  rename(U_creat_mg = urine_creatinine_mgperd,
        UN_gperd = urine_n_gperd,
         UUN_gperd = urine_urea_n_gperd)

ggpairs(urine_dat_d[, 7:14], aes(col = urine_dat_d$period),
  lower = list(
    continuous = wrap("points", alpha = .8, shape = 23, size = .1),
    combo = wrap("dot", alpha = 0.4, size = 0.2)
  ),
  upper = list(continuous = wrap("cor", size = 3)),
  diag = list(continuous = wrap("densityDiag", alpha = 0.5))
) + theme(strip.text = element_text(size = 8)) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1) +
  annotate("segment", x = -Inf, xend = -Inf, y = -Inf, yend = Inf)

4.1.2 Pairs by prot_lvl

Here, we can see the pairwise relationships between variables with colors indicating the protein level. The orange high-CP dots and teal low-CP dots are often in distinct groups, yet in the upper diagonal we can see that the directionality of the correlations is similar across HP and LP oscillation levels.

ggpairs(urine_dat_d[, 7:14], aes(col = urine_dat_d$prot_lvl),
  lower = list(
    continuous = wrap("points", alpha = .8, shape = 23, size = .1),
    combo = wrap("dot", alpha = 0.4, size = 0.2)
  ),
  upper = list(continuous = wrap("cor", size = 3)),
  diag = list(continuous = wrap("densityDiag", alpha = 0.5))
) +theme( strip.text = element_text(size = 8)) +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1) + 
  annotate("segment", x = -Inf, xend = -Inf, y = -Inf, yend = Inf)

4.1.3 Pairs by osc_lvl

p <- ggpairs(urine_dat_d[, 7:14], aes(col = urine_dat_d$osc_lvl),
  lower = list(
    continuous = wrap("points", alpha = .8, shape = 23, size = .1),
    combo = wrap("dot", alpha = 0.4, size = 0.2)
  ),
  upper = list(continuous = wrap("cor", size = 3)),
  diag = list(continuous = wrap("densityDiag", alpha = 0.5))
)   +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1)) +
  annotate("segment", x = -Inf, xend = Inf, y = -Inf, yend = -Inf, lwd = 1) + 
  annotate("segment", x = -Inf, xend = -Inf, y = -Inf, yend = Inf)

# Loop through each plot changing relevant scales
for(i in 1:p$nrow) {
  for(j in 1:p$ncol){
    p[i,j] <- p[i,j] + 
        scale_fill_manual(values=c("purple", "darkgrey")) +
        scale_color_manual(values=c("purple", "darkgrey"))  
  }
}
p +theme( strip.text = element_text(size = 8))