MAW636-Analyze-Intensive-Rumen-V01.Rmd
M.G. Erickson
08/29/2022
1 Load packages
# Load packages-----------------------------------------------------------------
library(readxl)
library(tidyverse)
library(lmerTest)
library(lme4)
library(emmeans)
library(ggplot2)
library(afex)
library(kableExtra)
library(nlme)
library(MAW636)
afex::set_sum_contrasts()
theme_set(theme_bw())
# Expect 640 obs
# # Some issues with numeric day and hr with both data.frames, no effect on current
dat <- read.csv("01_organize/rumen_imp.csv")
nh3dat <- read.csv("01_organize/rumen_nh3_imp.csv")
#convert to factors## day is already char
dat$cow <- as.factor(dat$cow)
dat$time_n_cat <- as.factor(dat$time_n)
dat$hr_cat <- as.factor(dat$hr)
nh3dat$hr_cat = as.factor(nh3dat$hr_n)
# Summarize to period level ### TO BE MOVED TO ORG FILE #####
# Expect 32 obs
dat_period <- dat %>%
group_by(cow, period, prot_lvl, osc_lvl ) %>%
summarise(rumen_ph_mn = mean(rumen_ph))%>%
mutate(cow = as.character(cow))
nh3dat_period = nh3dat %>%
group_by(cow, period, prot_lvl, osc_lvl ) %>%
summarise(ammonia_mM_mn = mean(ammonia_mM),
taa_mM_mn = mean(taa_mM)) %>%
mutate(cow = as.character(cow))
rum_per_list = dat_period %>%
left_join(nh3dat_period) %>%
pivot_longer(., cols = -c(cow, period, prot_lvl, osc_lvl)) %>%
split(.$name)2 Aggregated Analysis
rum_mods <- rum_per_list %>%
purrr::map(., ~ mixed(value ~ period + prot_lvl * osc_lvl + (1 | cow), data = .))## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]
## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]
## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]rum_anova <- anova_out(rum_mods)
rum_em = em_df(rum_mods)
rum_em_table = em_table(rum_em)
rum_pval = p_vals(rum_mods)
rum_em_table %>%
cbind(., rum_pval[,2:4]) %>%
kable(caption = "Rumen pH Model") %>%
kable_styling() | name | emmean_OLP | emmean_SLP | emmean_OHP | emmean_SHP | max_err | prot_lvl | osc_lvl | prot_lvl:osc_lvl |
|---|---|---|---|---|---|---|---|---|
| ammonia_mM_mn | 1.25 | 1.38 | 1.64 | 1.68 | 0.10 | 0.000 | 0.238 | 0.524 |
| rumen_ph_mn | 6.35 | 6.38 | 6.37 | 6.38 | 0.04 | 0.728 | 0.414 | 0.624 |
| taa_mM_mn | 2.70 | 2.82 | 2.56 | 2.50 | 0.21 | 0.258 | 0.883 | 0.654 |
3 Over Time Analysis
3.1 pH
Significant effects and tendencies
- Effect of period is significant but not of interest.
- Sampling hour (hr) is significant.
- There is a tendency for an interaction between CP level (prot_lvl) and sampling hour (hr)
3.1.1 pH Model
rumen_ph_mod1 <- lme(rumen_ph ~ period + osc_lvl * prot_lvl * day * hr_cat,
random = ~ 1 | cow / period / day, correlation = corAR1(form = ~ 1 | cow / period / day), data = dat
)
#summary(rumen_ph_mod1)
df = anova(rumen_ph_mod1) %>%
data.frame() %>%
mutate(across(where(is.numeric), round, 3))
df %>%
kable(caption = "Rumen pH Model") %>%
kable_styling() %>%
row_spec(which(df$p.value < 0.05), bold = T, background = "lightyellow") %>%
row_spec(which(df$p.value > 0.05 & df$p.value < 0.10), bold = T, background = "grey25")| numDF | denDF | F.value | p.value | |
|---|---|---|---|---|
| (Intercept) | 1 | 448 | 37724.621 | 0.000 |
| period | 3 | 17 | 12.871 | 0.000 |
| osc_lvl | 1 | 17 | 0.761 | 0.395 |
| prot_lvl | 1 | 17 | 0.086 | 0.772 |
| day | 3 | 84 | 0.965 | 0.413 |
| hr_cat | 4 | 448 | 310.654 | 0.000 |
| osc_lvl:prot_lvl | 1 | 17 | 0.259 | 0.617 |
| osc_lvl:day | 3 | 84 | 0.380 | 0.768 |
| prot_lvl:day | 3 | 84 | 0.354 | 0.786 |
| osc_lvl:hr_cat | 4 | 448 | 0.401 | 0.808 |
| prot_lvl:hr_cat | 4 | 448 | 2.295 | 0.058 |
| day:hr_cat | 12 | 448 | 1.040 | 0.410 |
| osc_lvl:prot_lvl:day | 3 | 84 | 1.188 | 0.319 |
| osc_lvl:prot_lvl:hr_cat | 4 | 448 | 0.784 | 0.536 |
| osc_lvl:day:hr_cat | 12 | 448 | 0.884 | 0.564 |
| prot_lvl:day:hr_cat | 12 | 448 | 1.316 | 0.205 |
| osc_lvl:prot_lvl:day:hr_cat | 12 | 448 | 1.031 | 0.419 |
3.1.2 pH Model-implied means
e_rumen_ph = emmeans(rumen_ph_mod1, ~ osc_lvl*prot_lvl*day*hr_cat)
em_rumen_ph = e_rumen_ph %>% data.frame() %>% mutate(hr = as.numeric(as.character(hr_cat)))
# get estimates for marginal effects but don't save them
e1 = emmeans(rumen_ph_mod1, ~ hr_cat) ; e1## hr_cat emmean SE df lower.CL upper.CL
## 6 6.69 0.0375 8 6.61 6.78
## 8 6.75 0.0375 8 6.66 6.83
## 10 6.34 0.0375 8 6.26 6.43
## 12 6.14 0.0375 8 6.05 6.23
## 19 5.92 0.0375 8 5.83 6.00
##
## Results are averaged over the levels of: period, osc_lvl, prot_lvl, day
## Degrees-of-freedom method: containment
## Confidence level used: 0.95e2 = emmeans(rumen_ph_mod1, ~ prot_lvl*hr_cat) ; e2## prot_lvl hr_cat emmean SE df lower.CL upper.CL
## HP 6 6.65 0.0434 8 6.55 6.75
## LP 6 6.74 0.0434 8 6.64 6.84
## HP 8 6.78 0.0434 8 6.68 6.88
## LP 8 6.72 0.0434 8 6.62 6.82
## HP 10 6.34 0.0434 8 6.24 6.44
## LP 10 6.34 0.0434 8 6.24 6.44
## HP 12 6.17 0.0434 8 6.07 6.27
## LP 12 6.11 0.0434 8 6.01 6.21
## HP 19 5.93 0.0434 8 5.83 6.03
## LP 19 5.91 0.0434 8 5.81 6.01
##
## Results are averaged over the levels of: period, osc_lvl, day
## Degrees-of-freedom method: containment
## Confidence level used: 0.95# show the sampling hr and prot_lvl main effects and interaction
e2 %>%
data.frame() %>%
mutate(hr = as.numeric(as.character(hr_cat))) %>%
ggplot(aes(x = hr, y = emmean, col = prot_lvl, group = prot_lvl)) +
geom_point() +
geom_line() +
scale_x_continuous(breaks = c(6, 8, 10, 12, 19)) +
ggtitle(paste0("Rumen pH by CP level (LP = 13.8%, HP = 15.5% CP) at 5x/d sampling \n during 4-d sampling periods marginalized across 4 exp. periods for n = 8 cows")) +
labs(x = "Sampling hour", y = "Model-implied mean")
3.2 Ammonia (NH3)
Significant effects and tendencies
- Effect of period is significant but not of interest
- CP level (prot_lvl) significant
- Sampling hour (hr) significant
- CP level x sampling hour significant
- Day x sampling hour (hr) effect is significant
- Feeding pattern (osc_lvl) x day x sampling hour 3-way interaction is a tendency**
3.2.1 NH3 Model
rumen_nh3_mod1 <- lme(ammonia_mM ~ period + osc_lvl * prot_lvl * day * hr_cat,
random = ~ 1 | cow / period / day, correlation = corAR1(form = ~ 1 | cow / period / day), data = nh3dat
)
#summary(rumen_nh3_mod1)
df = anova(rumen_nh3_mod1) %>%
data.frame() %>%
mutate(across(where(is.numeric), round, 3))
df %>%
kable(caption = "Rumen NH3 Model") %>%
kable_styling() %>%
row_spec(which(df$p.value < 0.05), bold = T, background = "lightyellow") %>%
row_spec(which(df$p.value > 0.05 & df$p.value < 0.10), bold = T, background = "grey25")| numDF | denDF | F.value | p.value | |
|---|---|---|---|---|
| (Intercept) | 1 | 448 | 386.421 | 0.000 |
| period | 3 | 17 | 5.874 | 0.006 |
| osc_lvl | 1 | 17 | 1.251 | 0.279 |
| prot_lvl | 1 | 17 | 25.544 | 0.000 |
| day | 3 | 84 | 0.361 | 0.782 |
| hr_cat | 4 | 448 | 21.530 | 0.000 |
| osc_lvl:prot_lvl | 1 | 17 | 0.356 | 0.559 |
| osc_lvl:day | 3 | 84 | 0.964 | 0.413 |
| prot_lvl:day | 3 | 84 | 1.835 | 0.147 |
| osc_lvl:hr_cat | 4 | 448 | 0.876 | 0.478 |
| prot_lvl:hr_cat | 4 | 448 | 4.316 | 0.002 |
| day:hr_cat | 12 | 448 | 4.106 | 0.000 |
| osc_lvl:prot_lvl:day | 3 | 84 | 0.805 | 0.495 |
| osc_lvl:prot_lvl:hr_cat | 4 | 448 | 1.398 | 0.234 |
| osc_lvl:day:hr_cat | 12 | 448 | 1.667 | 0.071 |
| prot_lvl:day:hr_cat | 12 | 448 | 1.313 | 0.208 |
| osc_lvl:prot_lvl:day:hr_cat | 12 | 448 | 1.084 | 0.372 |
3.2.2 NH3 Model-implied means
e3 = e_rumen_nh3 = emmeans(rumen_nh3_mod1, ~ hr_cat) ; e3## hr_cat emmean SE df lower.CL upper.CL
## 6 1.29 0.0899 8 1.08 1.49
## 8 1.43 0.0899 8 1.23 1.64
## 10 1.93 0.0899 8 1.72 2.13
## 12 1.36 0.0899 8 1.16 1.57
## 19 1.42 0.0899 8 1.21 1.63
##
## Results are averaged over the levels of: period, osc_lvl, prot_lvl, day
## Degrees-of-freedom method: containment
## Confidence level used: 0.95e4 = e_rumen_nh3 = emmeans(rumen_nh3_mod1, ~ prot_lvl) ; e4## prot_lvl emmean SE df lower.CL upper.CL
## HP 1.66 0.0828 8 1.47 1.85
## LP 1.31 0.0828 8 1.12 1.50
##
## Results are averaged over the levels of: period, osc_lvl, day, hr_cat
## Degrees-of-freedom method: containment
## Confidence level used: 0.95e5 = e_rumen_nh3 = emmeans(rumen_nh3_mod1, ~ prot_lvl*hr_cat) ; e5## prot_lvl hr_cat emmean SE df lower.CL upper.CL
## HP 6 1.52 0.108 8 1.273 1.77
## LP 6 1.05 0.108 8 0.801 1.30
## HP 8 1.44 0.108 8 1.189 1.69
## LP 8 1.43 0.108 8 1.178 1.68
## HP 10 2.03 0.108 8 1.786 2.28
## LP 10 1.82 0.108 8 1.570 2.07
## HP 12 1.58 0.108 8 1.328 1.83
## LP 12 1.15 0.108 8 0.902 1.40
## HP 19 1.72 0.108 8 1.473 1.97
## LP 19 1.12 0.108 8 0.870 1.37
##
## Results are averaged over the levels of: period, osc_lvl, day
## Degrees-of-freedom method: containment
## Confidence level used: 0.95e5 %>%
data.frame() %>%
mutate(hr = as.numeric(as.character(hr_cat))) %>%
ggplot(aes(x = hr, y = emmean, col = prot_lvl, group = prot_lvl)) +
geom_point() +
geom_line() +
scale_x_continuous(breaks = c(6, 8, 10, 12, 19)) +
ggtitle(paste0("Rumen ammonia by CP level (LP = 13.8%, HP = 15.5% CP) at 5x/d sampling \n during 4-d sampling periods marginalized across 4 exp. periods for n = 8 cows")) +
labs(x = "Sampling hour", y = "Model-implied mean")
e6 <- e_rumen_nh3 <- emmeans(rumen_nh3_mod1, ~ osc_lvl * day * hr_cat)
e6## osc_lvl day hr_cat emmean SE df lower.CL upper.CL
## O Day1 6 0.718 0.182 8 0.299 1.14
## S Day1 6 1.105 0.182 8 0.686 1.52
## O Day2 6 1.471 0.182 8 1.051 1.89
## S Day2 6 1.741 0.182 8 1.322 2.16
## O Day3 6 1.520 0.182 8 1.101 1.94
## S Day3 6 1.309 0.182 8 0.890 1.73
## O Day4 6 0.961 0.182 8 0.542 1.38
## S Day4 6 1.461 0.182 8 1.042 1.88
## O Day1 8 0.977 0.182 8 0.558 1.40
## S Day1 8 1.530 0.182 8 1.111 1.95
## O Day2 8 1.368 0.182 8 0.948 1.79
## S Day2 8 1.293 0.182 8 0.874 1.71
## O Day3 8 1.659 0.182 8 1.240 2.08
## S Day3 8 1.385 0.182 8 0.965 1.80
## O Day4 8 1.651 0.182 8 1.232 2.07
## S Day4 8 1.598 0.182 8 1.178 2.02
## O Day1 10 2.174 0.182 8 1.755 2.59
## S Day1 10 2.069 0.182 8 1.650 2.49
## O Day2 10 1.937 0.182 8 1.518 2.36
## S Day2 10 1.960 0.182 8 1.540 2.38
## O Day3 10 1.730 0.182 8 1.311 2.15
## S Day3 10 1.878 0.182 8 1.459 2.30
## O Day4 10 1.622 0.182 8 1.203 2.04
## S Day4 10 2.045 0.182 8 1.626 2.46
## O Day1 12 1.605 0.182 8 1.186 2.02
## S Day1 12 1.328 0.182 8 0.908 1.75
## O Day2 12 1.597 0.182 8 1.178 2.02
## S Day2 12 1.461 0.182 8 1.042 1.88
## O Day3 12 1.253 0.182 8 0.834 1.67
## S Day3 12 1.543 0.182 8 1.124 1.96
## O Day4 12 0.991 0.182 8 0.572 1.41
## S Day4 12 1.131 0.182 8 0.711 1.55
## O Day1 19 1.672 0.182 8 1.253 2.09
## S Day1 19 1.406 0.182 8 0.987 1.83
## O Day2 19 1.337 0.182 8 0.918 1.76
## S Day2 19 1.244 0.182 8 0.824 1.66
## O Day3 19 1.316 0.182 8 0.897 1.74
## S Day3 19 1.290 0.182 8 0.870 1.71
## O Day4 19 1.370 0.182 8 0.951 1.79
## S Day4 19 1.727 0.182 8 1.308 2.15
##
## Results are averaged over the levels of: period, prot_lvl
## Degrees-of-freedom method: containment
## Confidence level used: 0.95e6 %>%
data.frame() %>%
mutate(hr = as.numeric(as.character(hr_cat))) %>%
ggplot(aes(x = hr, y = emmean, col = osc_lvl, group = osc_lvl)) +
geom_point() +
geom_line() +
scale_color_manual(values = c("purple", "grey"))+
scale_x_continuous(breaks = c(6, 8, 10, 12, 19)) +
facet_wrap(~day, nrow = 1) +
ggtitle(paste0("Rumen ammonia by feeding pattern (O = oscillating 48-hr, S = static) \n at 5x/d sampling during 4-d sampling periods \n marginalized across 4 exp. periods for n = 8 cows")) +
labs(x = "Sampling hour", y = "Model-implied mean")
3.3 Total Amino Acid (TAA)
Significant effects and tendencies
- Effect of period is significant but not of interest
- Sampling hour (hr) is significant
- CP level (prot_lvl) x sampling hour (hr) tendency
- day x sampling hour (hr) tendency
- Feeding pattern (osc_lvl) x day x sampling hour 3-way interaction is a tendency
3.3.1 TAA Model
rumen_taa_mod1 <- lme(taa_mM ~ period + osc_lvl * prot_lvl * day * hr_cat,
random = ~ 1 | cow / period / day, correlation = corAR1(form = ~ 1 | cow / period / day), data = nh3dat
)
#summary(rumen_taa_mod1)
df = anova(rumen_taa_mod1) %>%
data.frame() %>%
mutate(across(where(is.numeric), round, 3))
df %>%
kable(caption = "Rumen TAA Model") %>%
kable_styling() %>%
row_spec(which(df$p.value < 0.05), bold = T, background = "lightyellow") %>%
row_spec(which(df$p.value > 0.05 & df$p.value < 0.10), bold = T, background = "grey25")| numDF | denDF | F.value | p.value | |
|---|---|---|---|---|
| (Intercept) | 1 | 448 | 577.333 | 0.000 |
| period | 3 | 17 | 6.076 | 0.005 |
| osc_lvl | 1 | 17 | 0.022 | 0.883 |
| prot_lvl | 1 | 17 | 1.493 | 0.238 |
| day | 3 | 84 | 2.814 | 0.044 |
| hr_cat | 4 | 448 | 42.566 | 0.000 |
| osc_lvl:prot_lvl | 1 | 17 | 0.222 | 0.643 |
| osc_lvl:day | 3 | 84 | 0.631 | 0.597 |
| prot_lvl:day | 3 | 84 | 1.558 | 0.206 |
| osc_lvl:hr_cat | 4 | 448 | 0.571 | 0.684 |
| prot_lvl:hr_cat | 4 | 448 | 2.226 | 0.065 |
| day:hr_cat | 12 | 448 | 2.100 | 0.016 |
| osc_lvl:prot_lvl:day | 3 | 84 | 0.941 | 0.425 |
| osc_lvl:prot_lvl:hr_cat | 4 | 448 | 0.134 | 0.970 |
| osc_lvl:day:hr_cat | 12 | 448 | 1.618 | 0.083 |
| prot_lvl:day:hr_cat | 12 | 448 | 0.336 | 0.982 |
| osc_lvl:prot_lvl:day:hr_cat | 12 | 448 | 1.032 | 0.418 |
3.3.2 TAA Model-implied means
e7 = emmeans(rumen_taa_mod1, ~ hr_cat) ; e7## hr_cat emmean SE df lower.CL upper.CL
## 6 1.19 0.197 8 0.731 1.64
## 8 2.40 0.197 8 1.941 2.85
## 10 4.50 0.197 8 4.047 4.96
## 12 2.41 0.197 8 1.960 2.87
## 19 2.74 0.197 8 2.282 3.19
##
## Results are averaged over the levels of: period, osc_lvl, prot_lvl, day
## Degrees-of-freedom method: containment
## Confidence level used: 0.95e8 = emmeans(rumen_taa_mod1, ~ prot_lvl*hr_cat) ; e8## prot_lvl hr_cat emmean SE df lower.CL upper.CL
## HP 6 1.16 0.274 8 0.531 1.80
## LP 6 1.21 0.274 8 0.574 1.84
## HP 8 1.79 0.274 8 1.162 2.43
## LP 8 3.00 0.274 8 2.363 3.63
## HP 10 4.43 0.274 8 3.796 5.06
## LP 10 4.57 0.274 8 3.940 5.21
## HP 12 2.43 0.274 8 1.802 3.07
## LP 12 2.39 0.274 8 1.760 3.03
## HP 19 2.83 0.274 8 2.195 3.46
## LP 19 2.65 0.274 8 2.013 3.28
##
## Results are averaged over the levels of: period, osc_lvl, day
## Degrees-of-freedom method: containment
## Confidence level used: 0.95e8 %>%
data.frame() %>%
mutate(hr = as.numeric(as.character(hr_cat))) %>%
ggplot(aes(x = hr, y = emmean, col = prot_lvl, group = prot_lvl)) +
geom_point() +
geom_line() +
scale_x_continuous(breaks = c(6, 8, 10, 12, 19)) +
ggtitle(paste0("Rumen total amino acid by CP level (LP = 13.8%, HP = 15.5% CP) at 5x/d sampling \n during 4-d sampling periods marginalized across 4 exp. periods for n = 8 cows")) +
labs(x = "Sampling hour", y = "Model-implied mean")
e9 <- emmeans(rumen_taa_mod1, ~ osc_lvl * day * hr_cat); e9## osc_lvl day hr_cat emmean SE df lower.CL upper.CL
## O Day1 6 0.727 0.53 8 -0.49547 1.95
## S Day1 6 0.564 0.53 8 -0.65766 1.79
## O Day2 6 1.047 0.53 8 -0.17487 2.27
## S Day2 6 1.001 0.53 8 -0.22088 2.22
## O Day3 6 1.562 0.53 8 0.34004 2.78
## S Day3 6 1.899 0.53 8 0.67681 3.12
## O Day4 6 1.220 0.53 8 -0.00165 2.44
## S Day4 6 1.459 0.53 8 0.23678 2.68
## O Day1 8 1.770 0.53 8 0.54751 2.99
## S Day1 8 2.634 0.53 8 1.41150 3.86
## O Day2 8 0.894 0.53 8 -0.32768 2.12
## S Day2 8 1.547 0.53 8 0.32459 2.77
## O Day3 8 3.199 0.53 8 1.97727 4.42
## S Day3 8 2.470 0.53 8 1.24801 3.69
## O Day4 8 3.842 0.53 8 2.61945 5.06
## S Day4 8 2.806 0.53 8 1.58417 4.03
## O Day1 10 5.164 0.53 8 3.94153 6.39
## S Day1 10 3.408 0.53 8 2.18622 4.63
## O Day2 10 4.552 0.53 8 3.33010 5.77
## S Day2 10 5.107 0.53 8 3.88462 6.33
## O Day3 10 3.801 0.53 8 2.57911 5.02
## S Day3 10 5.496 0.53 8 4.27355 6.72
## O Day4 10 3.533 0.53 8 2.31087 4.76
## S Day4 10 4.949 0.53 8 3.72724 6.17
## O Day1 12 2.641 0.53 8 1.41856 3.86
## S Day1 12 2.481 0.53 8 1.25901 3.70
## O Day2 12 2.592 0.53 8 1.37023 3.81
## S Day2 12 2.418 0.53 8 1.19594 3.64
## O Day3 12 2.537 0.53 8 1.31506 3.76
## S Day3 12 3.003 0.53 8 1.78099 4.23
## O Day4 12 2.171 0.53 8 0.94920 3.39
## S Day4 12 1.467 0.53 8 0.24487 2.69
## O Day1 19 2.364 0.53 8 1.14160 3.59
## S Day1 19 1.783 0.53 8 0.56072 3.01
## O Day2 19 2.835 0.53 8 1.61314 4.06
## S Day2 19 2.722 0.53 8 1.50027 3.94
## O Day3 19 3.014 0.53 8 1.79221 4.24
## S Day3 19 2.696 0.53 8 1.47413 3.92
## O Day4 19 3.161 0.53 8 1.93885 4.38
## S Day4 19 3.317 0.53 8 2.09506 4.54
##
## Results are averaged over the levels of: period, prot_lvl
## Degrees-of-freedom method: containment
## Confidence level used: 0.95e9 %>%
data.frame() %>%
mutate(hr = as.numeric(as.character(hr_cat))) %>%
ggplot(aes(x = hr, y = emmean, col = osc_lvl, group = osc_lvl)) +
geom_point() +
geom_line() +
scale_color_manual(values = c("purple", "grey"))+
scale_x_continuous(breaks = c(6, 8, 10, 12, 19)) +
facet_wrap(~day, nrow = 1) +
ggtitle(paste0("Rumen total amino acid by feeding pattern (O = oscillating 48-hr, S = static) at 5x/d sampling \n during 4-d sampling periods marginalized across 4 exp. periods for n = 8 cows")) +
labs(x = "Sampling hour", y = "Model-implied mean")
##%######################################################%##
# #
#### End ####
# #
##%######################################################%##