MAW636 ADSA Analysis
M.G. Erickson
1/24/2022
1 Set-Up Workspace
#.Set Working Directory----------------------------------------------------------------------------------------------------
setwd("~/Cow/S21-MAW636 Analysis/Production Data Analysis/")
# Load packages----------------------------------------------------------------------------------------------------.
library(plotly)
library(readxl)
library(tidyverse)
library(lubridate)
library(lmerTest)
library(lme4)
library(broom.mixed)
library(emmeans)
library(ggplot2)
library(gtsummary)
library(afex)
library(broom.mixed)
# Upload data, by cow*period, filter production square, n = 32--------------------------------------------------------------------------------
# treatment schedule
trtsched_upload <- read_excel("MAW636 Treatment Schedule.xlsx") %>%
pivot_longer(cols = -Cow, names_to = "period", values_to = "trt") %>%
rename(cow = "Cow") %>%
mutate(cow = as.character(cow)) %>%
na_if("NA") %>%
drop_na() %>%
separate(trt, c("osc_lvl", "prot_lvl"), sep = cumsum(c(1, 2))) %>%
mutate(trt = paste0(osc_lvl, prot_lvl)) # separate into factorial trt
# MILK -----------------------------------------------------------------------------
milk_comp <- read.csv("milk_comp_wt_summary.csv") %>%
mutate(cow = as.character(cow)) %>%
select(-X) %>%
filter(period != "P5") %>%
filter(!(cow == "8993" & period == "P2")) %>%
filter(!(cow == "8993" & period == "P3")) %>%
filter(!(cow == "8993" & period == "P4")) %>%
filter(!(cow == "9067" & period == "P4"))
milk_comp_long <- milk_comp %>%
left_join(trtsched_upload) %>%
mutate(
osc_lvl = if_else(str_detect(trt, "^S"), "S", "O"),
prot_lvl = if_else(str_detect(trt, "HP"), "HP", "LP"),
prot_MUNY_gperg_ratio = prot_prod_kgperd * 1000 / MUNY_gperd
) %>%
pivot_longer(
cols = -c(cow, period, osc_lvl, prot_lvl, trt),
names_to = "var", values_to = "val"
) %>%
split(., .$var) # 10 variables = 720
# DMI -----------------------------------------------------------------------------
DMI_dat <- read.csv("DMI_dat_period.csv") %>%
select(-X) %>%
filter(period != "P5") %>%
mutate(cow = as.character(cow)) %>% # 62 obs, checks out.
left_join(trtsched_upload) %>%
mutate(
osc_lvl = if_else(str_detect(trt, "^S"), "S", "O"),
prot_lvl = if_else(str_detect(trt, "HP"), "HP", "LP")
)
DMI_dat_osc <- read.csv("DMI_dat_osc.csv") %>%
select(-X) %>%
filter(period != "P5") %>%
mutate(cow = as.character(cow)) %>%
left_join(trtsched_upload) %>%
mutate(
osc_lvl = if_else(str_detect(trt, "^S"), "S", "O"),
prot_lvl = if_else(str_detect(trt, "HP"), "HP", "LP")
)
# Feed Efficiency - combine dmi and milk datasets
fe_dat <- milk_comp %>%
left_join(DMI_dat) %>%
mutate(fe = fpcm / DMI_kg) # kg/kg
# Checks
unique(milk_comp_long$cow) # check that we have 16 cows## NULLunique(milk_comp_long$period)## NULL# DAIRYLAND FEED ----------------------------------------------------------------------
feed_dairyland <- read.csv("MAW636 Feed Dairyland.csv") %>%
rename(Index = Description1) %>%
mutate(Index = as.character(Index))
feed_index <- read.csv("MAW636 Feed Index Dairyland.csv") %>%
mutate(Index = as.character(Index))
feed_dy <- feed_dairyland %>%
left_join(feed_index)
rm(feed_dairyland)
# pivot long
feed_dy_long <- feed_dy[, 13:42] %>%
pivot_longer(cols = -c(Diet, Period, OscDay), names_to = "var", values_to = "val")
feed_summary <- feed_dy_long %>%
filter(Period != "P5") %>%
group_by(Diet, var) %>%
summarise(
meanval = mean(val),
sdval = sd(val)
) %>%
pivot_wider(names_from = Diet, values_from = c(meanval, sdval)) %>%
mutate(across(where(is.numeric), round, 2))
write.csv(feed_summary, "adsa_dairyland_feed_summary.csv")
# MARY FEED ---------------------------------------------------------------------------
mary_CP <- read_xlsx("MAW636 Feed USDA MGE.xlsx", sheet = "CP") %>%
mutate(rep = rep(1:3, each = 48),
Index = as.character(Index))
mary_DM <- read_xlsx("MAW636 Feed USDA MGE.xlsx", sheet = "DM") %>%
mutate(rep = rep(1:3, each = 40),
Index = as.character(Index))
mary_NDF_ADF <- read_xlsx("MAW636 Feed USDA MGE.xlsx", sheet = "NDF ADF") %>%
mutate(rep = rep(1:3, each = 48),
Index = as.character(Index))
# includes ADF reruns, corrected for DM
mary_feed <- mary_CP %>%
left_join(mary_DM, by = c("Index", "rep")) %>%
left_join(mary_NDF_ADF, by = c("Index", "rep")) %>%
filter(Index != "Blank" | Index != "EDTA") %>%
left_join(feed_index, by = "Index") %>%
mutate(Index = as.numeric(Index)) %>% # force numeric
filter(!is.na(Index)) %>%
mutate(
CPcorr = CP * (100 / DM),
NDFcorr = NDF * (100 / DM),
ADFcorr = ADF * (100 / DM)
) %>%
select(-c(CP, ADF, NDF, DM))
# pivot long
mary_feed_long <- mary_feed %>%
pivot_longer(cols = -c(Index, rep, Period, OscDay, Diet),
names_to = "var", values_to = "val")
rm(mary_CP)
rm(mary_DM)
rm(mary_NDF_ADF)2 Visualize diet composition period by period
The haylage changed between P1 and P2. ## Mary feed analysis
pmary <- mary_feed_long %>%
mutate(
Diet = factor(Diet,
levels = c("A", "B", "C", "D", "Corn Silage", "Haylage")
),
rep = as.character(rep)
) %>%
group_by(Diet, Period, var) %>%
ggplot(aes(y = val, x = Period, col = var, group = var)) +
facet_wrap(~Diet, nrow = 1) +
labs(title = "Mary feed samples, n = 3 analytical replicates per n = 1 period*diet sample") +
theme_minimal() +
geom_smooth(alpha = 0.2, lwd = 0.2) +
stat_smooth(geom='line', alpha=0.5, se=FALSE) +
scale_shape_manual(values=c(0, 2, 3))+
geom_point(aes(shape = rep))
pmary
ggsave("pmary.tiff", pmary, dpi = 300, height = 4, width = 9, device = "tiff")2.1 Dairyland feed analysis by period
p <- feed_dy_long %>%
mutate(Diet = factor(Diet,
levels = c("A", "B", "C", "D", "Corn Silage", "Haylage")
)) %>%
group_by(var) %>%
ggplot(aes(y = val, x = Period, col = var, group = var)) + geom_point(size = 0.5)+
facet_wrap(~Diet, nrow = 1) +
labs(title = "Dairyland feed samples, 1 sample per period*diet", x = NULL, y = "% of DM") +
theme_minimal() +
geom_line()
ggplotly(p, width = 800, height = 400) %>% layout(yaxis = list(autorange = TRUE)) # 3 Data cleaning for production variables
Includes n = 62 observations, 16 cows x 4 periods. See Milk Data Organization.Rmd for details.
4 Set-Up for Modeling
Set orthogonal (sum) contrasts and define any convenience functions to make output of mixed() more readable/workable. The mixed() function is a wrapper for lmer(). Whereas lmer() by default outputs dummy (treatment) contrasts, mixed() outputs orthogonal (sum) contrasts with F tests based on Type III sums of squares.
# Set contrast options ----------------------------------------------------
afex::set_sum_contrasts()
anova_proc <- function(model) {
result = tidy(model$anova_table)
}5 Model expression
In the ANOVA-style notation, a general model can be expressed for the values of each dependent variable: \[ value_{ijkl} = \mu + period_{j} + protlevel_k + feedpattern_l + (protlevel_k \times feedpattern_l) + cow_i + \epsilon_{ijkl}\] where
\(\mu\) is the grand mean of the observations
\(period\) is the fixed effect of period (\(j\) = P1, P3, P4, P5)
\(protlevel\) is the fixed main effect of protein level (\(k\) = low, high)
\(feedpattern\) is the fixed main effect of feeding pattern (\(l\) = O, S)
\(protlevel_k \times feedpattern_l\) is the interaction of feeding pattern and protein level.
\(cow\) is the random effect of cow (\(i\) = 1,…,8)
6 Production Models
For production models, data frames contain 32 observations, representing 8 cows in 4 periods. When no significant treatment effects, get overall mean of treatment means. This is different from the grand mean.
6.1 DMI
dmi_mod <- mixed(DMI_kg ~ period + prot_lvl * osc_lvl + (1 | cow), data = DMI_dat)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]dmi_mod## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: DMI_kg ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: DMI_dat
## Effect df F p.value
## 1 period 3, 40.38 24.94 *** <.001
## 2 prot_lvl 1, 40.17 0.37 .549
## 3 osc_lvl 1, 40.87 0.31 .584
## 4 prot_lvl:osc_lvl 1, 40.14 0.48 .495
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1plot(emmeans(dmi_mod,~ period ))
emmeans(dmi_mod, "1") # emmean overall, on intercept. ## 1 emmean SE df lower.CL upper.CL
## overall 25.9 0.575 15.8 24.7 27.1
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956.1.1 Oscillation Phase differences?
Out of curiosity, check for differences in DMI due to oscillation phase.
# check for any DMI differences in oscillating cows based on osc_phase
DMI_dat_lmer <- DMI_dat_osc %>%
filter(osc_lvl == "O") %>%
lmer(DMI_kg ~ OscDay + period + (1|cow), data = .)
summary(DMI_dat_lmer) ## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: DMI_kg ~ OscDay + period + (1 | cow)
## Data: .
##
## REML criterion at convergence: 268.5
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.85789 -0.61934 -0.01746 0.57398 2.09109
##
## Random effects:
## Groups Name Variance Std.Dev.
## cow (Intercept) 5.601 2.367
## Residual 2.738 1.655
## Number of obs: 62, groups: cow, 16
##
## Fixed effects:
## Estimate Std. Error df t value Pr(>|t|)
## (Intercept) 26.06950 0.62905 15.12336 41.442 < 2e-16 ***
## OscDay1 -0.07809 0.21013 42.56008 -0.372 0.712
## period1 2.83833 0.46552 50.96122 6.097 1.45e-07 ***
## period2 -2.38993 0.47231 49.22582 -5.060 6.23e-06 ***
## period3 -0.59088 0.44718 48.89147 -1.321 0.193
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) OscDy1 perid1 perid2
## OscDay1 0.000
## period1 -0.026 0.000
## period2 0.029 0.000 -0.297
## period3 -0.010 0.000 -0.169 -0.5696.2 Feed efficiency
fe_mod <- mixed(fe ~ period + prot_lvl * osc_lvl + (1 | cow), data = fe_dat)## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]fe_mod## Mixed Model Anova Table (Type 3 tests, S-method)
##
## Model: fe ~ period + prot_lvl * osc_lvl + (1 | cow)
## Data: fe_dat
## Effect df F p.value
## 1 period 3, 39.34 11.61 *** <.001
## 2 prot_lvl 1, 39.14 0.13 .720
## 3 osc_lvl 1, 39.80 0.72 .401
## 4 prot_lvl:osc_lvl 1, 39.11 0.03 .870
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1emmeans(fe_mod, "1")## 1 emmean SE df lower.CL upper.CL
## overall 1.46 0.0334 15.9 1.39 1.53
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.956.3 Milk
Run models in bulk with purrr::map and list of data.frames in milk_comp_long.
milk_mods <- milk_comp_long %>%
purrr::map(., ~ mixed(val ~ 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]
## 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]
## 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]
## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]
## Fitting one lmer() model. [DONE]
## Calculating p-values. [DONE]names <- names(milk_mods)
# ANOVA ---------------------------------------------------------------------------
anova_milk <- milk_mods %>%
purrr::map(., ~ anova_proc(.)) %>%
bind_rows() %>%
mutate(var = rep(names, each = 4)) %>%
relocate(var, 1) %>%
mutate(across(where(is.numeric), round, 3)) %>%
mutate(
sign = ifelse(p.value < 0.05, "*", ""),
tend = ifelse(.05 < p.value & p.value < 0.10, "^", "")
) %>%
filter(term != "period")
# Reformat ---------------------------------------------------------------------------
anova_milk_1 <- anova_milk %>%
mutate(p.paste = paste0(p.value, sign, tend)) %>%
select(var, term, p.paste) %>%
pivot_wider(names_from = "term", values_from = "p.paste") %>%
mutate(var = factor(var, levels = c("milk_kgperd", "fpcm", "fat_prod_kgperd", "prot_prod_kgperd", "lactose_prod_kgperd", "MUNY_gperd", "prot_MUNY_gperg_ratio", "Fat", "Tru.Prot", "Lactose", "Urea"))) %>%
arrange(., var)
write.csv(anova_milk_1, "anova_milk_1.csv")7 Model-implied trt means for anova with >1 sig. F-test
Regardless of whether the interaction is significant, presenting all four treatment means for consistency/transparency.
emmeans(milk_mods$MUNY_gperd, ~ prot_lvl * osc_lvl)## prot_lvl osc_lvl emmean SE df lower.CL upper.CL
## HP O 4.54 0.161 35.2 4.21 4.86
## LP O 3.31 0.157 33.5 3.00 3.63
## HP S 4.68 0.160 35.2 4.36 5.01
## LP S 3.53 0.157 33.4 3.21 3.85
##
## Results are averaged over the levels of: period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95emmeans(milk_mods$Tru.Prot, ~ prot_lvl * osc_lvl)## prot_lvl osc_lvl emmean SE df lower.CL upper.CL
## HP O 2.93 0.0382 23.7 2.85 3.01
## LP O 2.88 0.0377 22.7 2.80 2.96
## HP S 2.89 0.0382 23.7 2.81 2.97
## LP S 2.90 0.0378 22.8 2.83 2.98
##
## Results are averaged over the levels of: period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95emmeans(milk_mods$Urea, ~ prot_lvl * osc_lvl)## prot_lvl osc_lvl emmean SE df lower.CL upper.CL
## HP O 11.93 0.321 43.3 11.28 12.58
## LP O 8.70 0.311 42.0 8.07 9.32
## HP S 12.46 0.321 43.4 11.81 13.11
## LP S 9.19 0.312 41.6 8.56 9.82
##
## Results are averaged over the levels of: period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95emmeans(milk_mods$prot_MUNY_gperg_ratio, ~ prot_lvl * osc_lvl)## prot_lvl osc_lvl emmean SE df lower.CL upper.CL
## HP O 258 9.79 50.5 238 278
## LP O 338 9.47 49.9 319 357
## HP S 239 9.78 50.6 220 259
## LP S 321 9.47 49.6 302 340
##
## Results are averaged over the levels of: period
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.958 Model-implied overall (intercept) means for anova with non-significant F-tests
emmeans(milk_mods$milk_kgperd, "1")## 1 emmean SE df lower.CL upper.CL
## overall 38.5 1.15 15.9 36.1 40.9
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95emmeans(milk_mods$fpcm, "1")## 1 emmean SE df lower.CL upper.CL
## overall 37.6 0.916 15.8 35.7 39.5
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95emmeans(milk_mods$fat_prod_kgperd, "1")## 1 emmean SE df lower.CL upper.CL
## overall 1.57 0.0386 15.7 1.48 1.65
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95emmeans(milk_mods$prot_prod_kgperd, "1")## 1 emmean SE df lower.CL upper.CL
## overall 1.12 0.0311 15.9 1.05 1.18
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95emmeans(milk_mods$lactose_prod_kgperd, "1")## 1 emmean SE df lower.CL upper.CL
## overall 1.8 0.0604 16 1.67 1.93
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95emmeans(milk_mods$Fat, "1")## 1 emmean SE df lower.CL upper.CL
## overall 4.11 0.0903 15.7 3.92 4.3
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95emmeans(milk_mods$Lactose, "1")## 1 emmean SE df lower.CL upper.CL
## overall 4.66 0.0331 16 4.59 4.73
##
## Results are averaged over the levels of: period, prot_lvl, osc_lvl
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95##%######################################################%## # # #### Miscellaneous and visualizations #### # # ##%######################################################%##
9 Miscellaneous and visualizations
anova_proc <- function(model) {
result = tidy(model$anova_table)
#depvar = as.character(terms(model$full_model)[[2]])
#depvarc = rep(depvar, each = 4)
#cbind(depvarc, result)
}
em_proc <- function(model) {
result <- data.frame(emmeans(model, ~ prot_lvl * osc_lvl)) %>%
mutate(across(where(is.numeric), round, 3))
return(result)
#depvar = as.character(terms(model$full_model)[[2]])
#depvarc = rep(depvar, each = 4)
#cbind(depvarc, result) # depvar is "val" in list of df
}
em_proc(dmi_mod)
#-----------------------------------------------------------------------------
a_mix <- function(mod) {
data.frame(mod$anova_table) %>%
rownames_to_column() %>%
select(rowname, Pr..F.) %>%
mutate(P = round(Pr..F., 3)) %>%
select(-Pr..F.)
}
e_mix <- function(mod) {
data.frame(emmeans(mod, ~ prot_lvl * osc_lvl)) %>%
mutate(across(where(is.numeric), round, 3))
} # Not used in 2022.01.17 due to computing overall mu when F test not significant
#-----------------------------------------------------------------------------
em_milk <- milk_mods %>%
purrr::map(., ~em_proc(.)) %>%
bind_rows() %>%
mutate(var = rep(names, each = 4)) %>%
relocate(var, 1)
em_milk # maybe wait to do emmeans til after the anova()
ADSA_milk_em_01 <- milk_mods %>%
purrr::map(., ~emmeans(., ~osc_lvl*prot_lvl))
ADSA_milk_em_01
ADSA_milk_em_01_dat <- ADSA_milk_em_01 %>%
purrr::map(., ~as.data.frame(.)) %>%
bind_rows(.id = "column_label") %>%
select(-df, -lower.CL, -upper.CL) %>%
pivot_wider(names_from = c(osc_lvl, prot_lvl), values_from = c(emmean, SE)) %>%
mutate(SEM = pmax(SE_O_HP, SE_S_HP, SE_O_LP, SE_S_LP)) %>% # keep maximum of SE values
select(-contains("SE_")) %>%
mutate(across(where(is.numeric), round, 2)) %>%
filter(column_label != "milkwt") %>%
mutate(column_label = factor(column_label, levels = c("milk_kgperd", "fpcm", "fat_prod_kgperd", "prot_prod_kgperd", "lactose_prod_kgperd", "MUNY_gperd", "prot_MUNY_gperg_ratio" ,"Fat", "Tru.Prot", "Lactose", "Urea"))) %>%
arrange(., column_label)
write.csv(ADSA_milk_em_01_dat, "ADSA_milk_em_01_dat.csv")
ADSA_milk_em_02 <- milk_mods %>%
purrr::map(., ~emmeans(., ~1))
ADSA_milk_em_02
ADSA_milk_em_02_dat <- ADSA_milk_em_01 %>%
purrr::map(., ~as.data.frame(.)) %>%
bind_rows(.id = "column_label")%>%
select(-df, -lower.CL, -upper.CL) %>%
mutate(column_label = factor(column_label, levels = c("milk_kgperd", "fpcm", "fat_prod_kgperd", "prot_prod_kgperd", "lactose_prod_kgperd", "MUNY_gperd", "prot_MUNY_gperg_ratio" ,"Fat", "Tru.Prot", "Lactose", "Urea"))) %>%
arrange(., column_label)
write.csv(ADSA_milk_em_02_dat, "ADSA_milk_em_02_dat.csv")
#-----------------------------------------------------------------------------
write.csv(anova_milk, "")
library(gridExtra)
milk_mods %>%
purrr::map(., ~em_proc(.))
p1 <- em_milk %>%
filter(var == "prot_prod_kgperd") %>%
ggplot(aes(
x = interaction(osc_lvl, prot_lvl),
y = emmean, group = var
)) +
geom_point() +
geom_line() +
geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2) +
labs(title = "prot_prod_kgperd")
p2 <- em_milk %>%
filter(var == "MUNY_gperd") %>%
ggplot(aes(
x = interaction(osc_lvl, prot_lvl),
y = emmean, group = var
)) +
geom_point() +
geom_line() +
geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2) +
labs(title = "MUNY")
p3 <- em_milk %>%
filter(var == "Tru.Prot") %>%
ggplot(aes(
x = interaction(osc_lvl, prot_lvl),
y = emmean, group = var
)) +
geom_point() +
geom_line() +
geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2) +
labs(title = "TruProt")
p4 <- em_milk %>%
filter(var == "Urea") %>%
ggplot(aes(
x = interaction(osc_lvl, prot_lvl),
y = emmean, group = var
)) +
geom_point() +
geom_line() +
geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2) +
labs(title = "Urea")
p4 <- em_milk %>%
filter(var == "prot_MUNY_gperg_ratio") %>%
ggplot(aes(
x = interaction(osc_lvl, prot_lvl),
y = emmean, group = var
)) +
geom_point() +
geom_line() +
geom_errorbar(aes(ymin = emmean - SE, ymax = emmean + SE), width = 0.2) +
labs(title = "prot_MUNY_gperg_ratio")
grid.arrange(p1, p2, p3, p4)
# Less urea with oscillating...more incorporated into MCP? Intake differences day by day?