1 Load packages

library(ggplot2)
library(GGally)
library(tidyverse)
library(gridExtra)
library(knitr)
library(ggpubr)
library(ggthemes)
library(afex)
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.

  • The milk_comp_wt_summary_tp dataset has milk composition at each (AM and PM) timepoint.

urine_dat_d <- read.csv("urine.dat.d.export.csv") %>%
  select(-X) %>%
  rename(MUN = Urea) %>%
  mutate(cow = as.character(cow))
urine_dat <- read.csv("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))

    

milk_comp_wt_summary_tp <- read.csv("milk_comp_wt_summary_tp.csv") %>% 
  mutate(day_hr = paste0(day, "_", AMPM))  %>%
  mutate(day_hr_n = case_when(
    day_hr == "Day1_AM" ~ 1,
    day_hr == "Day1_PM" ~ 2,
    day_hr == "Day2_AM" ~ 3,
    day_hr == "Day2_PM" ~ 4,
    day_hr == "Day3_AM" ~ 5,
    day_hr == "Day3_PM" ~ 6,
    day_hr == "Day4_AM" ~ 7,
    day_hr == "Day4_PM" ~ 8,
  )) %>% 
  mutate(day_hr = factor(day_hr, levels = c(
    "Day1_AM", "Day1_PM", 
    "Day2_AM", "Day2_PM", 
    "Day3_AM", "Day3_PM", 
    "Day4_AM", "Day4_PM"
  ))) %>% 
  filter(period != "P5") #%>%
 # filter(!(cow == "8993" & period == "P2")) %>%
#  #filter(!(cow == "8993" & period == "P3")) %>%
#  filter(!(cow == "8993" & period == "P4")) %>%
#  filter(!(cow == "9067" & period == "P4")) ##### %>%
#  filter(!(cow == 145)) %>%
#  filter(!(cow == 148)) %>%
#  filter(!(cow == 164)) %>%
#  filter(!(cow == 7961)) %>%
#  filter(!(cow == 8156)) %>%
#  filter(!(cow == 8254)) %>%
#  filter(!(cow == 8359)) %>%
#  filter(!(cow == 9140))

3 Timepoint Level

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

3.1 Milk (n = 16 cows)

3.1.1 Milk production (kg)

There are minimal differences visible related to the oscillation pattern.

ymin = min(milk_comp_wt_summary_tp$milk_kgpertp)

milk_comp_wt_summary_tp %>% 
  ggplot(aes(x = day_hr, y = milk_kgpertp)) +
  labs(title = "Milk 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 = 1.25, xmax = 5.25, ymin = -Inf, ymax = Inf,
           alpha = .1, fill = "grey50") + 
  annotate("text", x = 3.5, y = 1.07*ymin, label = "High Diet Fed", size = 3) + 
  annotate("text", x = 7, y = 1.07*ymin, label = "Low Diet Fed", size = 3)

3.1.2 Milk protein (kg)

There are minimal differences visible related to the oscillation pattern.

ymin = min(milk_comp_wt_summary_tp$prot_prod_kgpertp)

milk_comp_wt_summary_tp %>% 
  ggplot(aes(x = day_hr, y = prot_prod_kgpertp)) +
  labs(title = "Milk protein (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 = 1.25, xmax = 5.25, ymin = -Inf, ymax = Inf,
           alpha = .1, fill = "grey50") + 
  annotate("text", x = 3.5, y = 1.07*ymin, label = "High Diet Fed", size = 3) + 
  annotate("text", x = 7, y = 1.07*ymin, label = "Low Diet Fed", size = 3)

3.1.3 Milk fat production (kg)

ymin = min(milk_comp_wt_summary_tp$fat_prod_kgpertp)

milk_comp_wt_summary_tp %>% 
  ggplot(aes(x = day_hr, y = fat_prod_kgpertp)) +
  labs(title = "Milk fat (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 = 1.25, xmax = 5.25, ymin = -Inf, ymax = Inf,
           alpha = .1, fill = "grey50") + 
  annotate("text", x = 3.5, y = 1.07*ymin, label = "High Diet Fed", size = 3) + 
  annotate("text", x = 7, y = 1.07*ymin, label = "Low Diet Fed", size = 3)

3.1.4 MUNY (g)

There is a curve in the oscillating treatment here.

ymin = min(milk_comp_wt_summary_tp$MUNY_gpertp)

milk_comp_wt_summary_tp %>% 
  ggplot(aes(x = day_hr, y = MUNY_gpertp)) +
  labs(title = "MUNY (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("rect", xmin = 1.25, xmax = 5.25, ymin = -Inf, ymax = Inf,
           alpha = .1, fill = "grey50") + 
  annotate("text", x = 3.5, y = 1.07*ymin, label = "High Diet Fed", size = 3) + 
  annotate("text", x = 7, y = 1.07*ymin, label = "Low Diet Fed", size = 3)

3.1.5 Milk nitrogen (g)

There is a curve in the oscillating treatment here.

ymin = min(milk_comp_wt_summary_tp$milk_n_gpertp)

milk_comp_wt_summary_tp %>% 
  ggplot(aes(x = day_hr, y = milk_n_gpertp)) +
  labs(title = "Milk 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("rect", xmin = 1.25, xmax = 5.25, ymin = -Inf, ymax = Inf,
           alpha = .1, fill = "grey50") + 
  annotate("text", x = 3.5, y = 1.07*ymin, label = "High Diet Fed", size = 3) + 
  annotate("text", x = 7, y = 1.07*ymin, label = "Low Diet Fed", size = 3)

3.2 Urine (n = 8 cows)

3.2.1 Urine pH

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

3.2.2 Urine production (kg)

ymin = min(urine_dat$urine_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 = 1.5, xmax = 7.5, ymin = -Inf, ymax = Inf,
           alpha = .1, fill = "grey50") + 
  annotate("text", x = 4.5, y = 1.2*ymin, label = "High Diet Fed", size = 3) + 
  annotate("text", x = 10, y = 1.2*ymin, label = "Low Diet Fed", size = 3)

3.2.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.2.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.

ymin = min(urine_dat$urine_urea_n_gpertp)

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 = 1.5, xmax = 7.5, ymin = -Inf, ymax = Inf,
           alpha = .1, fill = "grey50") + 
  annotate("text", x = 4.5, y = 1.2*ymin, label = "High Diet Fed", size = 3) + 
  annotate("text", x = 10, y = 1.2*ymin, label = "Low Diet Fed", size = 3)

3.2.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.

ymin = min(urine_dat$urine_n_gpertp)

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("rect", xmin = 1.5, xmax = 7.5, ymin = -Inf, ymax = Inf,
           alpha = .1, fill = "grey50") + 
  annotate("text", x = 4.5, y = 1.2*ymin, label = "High Diet Fed", size = 3) + 
  annotate("text", x = 10, y = 1.2*ymin, label = "Low Diet Fed", size = 3)

3.2.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.

ymin = min(urine_dat$urine_creatinine_mgpertp)

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) + 
  annotate("rect", xmin = 1.5, xmax = 7.5, ymin = -Inf, ymax = Inf,
           alpha = .1, fill = "grey50") + 
  annotate("text", x = 4.5, y = 1.2*ymin, label = "High Diet Fed", size = 3) + 
  annotate("text", x = 10, y = 1.2*ymin, label = "Low Diet Fed", size = 3)

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))

4.2 This version does not include polynomial models. See D02 to continue.