Library calls

library(correlation)
library(tidyverse)
library(broom)
library(GGally)
library(readxl)
library(ggplot2)
library(kableExtra)

theme_set(theme_classic())

Hypothesis 1 (H1)

  1. MPY/MUNY ratio is related to NUE, urine N excretion, or other indexes of N efficiency.
  1. Within cow, across diets (cow*period means for a trt)
  2. For groups of cows, across diets (group*period means)

** To be useful, the ratio must outperform MUN **

CowNFlow database

Ferreira, M., Delagarde, R., & Edouard, N. (2021). CowNflow: A dataset on nitrogen flows and balances in dairy cows fed maize forage or herbage-based diets. Data in Brief, 38, 107393. https://www.sciencedirect.com/science/article/pii/S2352340921006752

cdat <- read_excel("~/Library/CloudStorage/OneDrive-UW-Madison/Documents/2022-Fall/F22 MPY_MUNY/CowNflow_6_Cow_measurements1.xlsx") 

cdat %>% 
  group_by(exp) %>% 
  summarise(across(where(is.numeric), mean))
## # A tibble: 28 × 88
##    exp   1.-Co…¹ 2.-Co…² 2.-Bo…³ 2.-La…⁴ 2.-Ge…⁵ 3.-Fe…⁶ 3.-Fe…⁷ 3.-Fe…⁸ 3.-Fe…⁹
##    <chr>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
##  1 E01    77246.    69.1    504.    34.8   14.1    11.2    NA     NA      NA    
##  2 E02    79688.    50.4    522.    46.5   10.2    11.4    NA     NA      NA    
##  3 E03    81748.    42.5    551.    50     20.5    11.8    NA     NA      NA    
##  4 E04    79411     67.2    551     41.3    6.33   12.4    NA     NA      NA    
##  5 E05    81936.    49.2    563.    30.8   12.2     4.92    5.12   2.44    2.28 
##  6 E06    82553     62.4    642.    36.5   18.4     6.24    6.10   1.7     0.27 
##  7 E07    84081     67.3    634.    18.8    6.17   13.8     6.01   1.83    0.242
##  8 E08    84310.    75.2    628.    24.5    3.38    7.78    5.24   0.431   0.63 
##  9 E09    85297.    66.9    634.    36.8   10.2    15.1    NA     NA      NA    
## 10 E10    83836.   112.     635.    27      0       5.02    5.22   4.99    1.16 
## # … with 18 more rows, 78 more variables: `3.-Feed-5-intake-(kg-DM/day)` <dbl>,
## #   `3.-Feed-6-intake-(kg-DM/day)` <dbl>, `3.-Feed-7-intake-(kg-DM/day)` <dbl>,
## #   `3.-Feed-8-intake-(kg-DM/day)` <dbl>, `3.-Feed-9-intake-(kg-DM/day)` <dbl>,
## #   `3.-Feed-10-intake-(kg-DM/day)` <dbl>,
## #   `3.-Feed-11-intake-(kg-DM/day)` <dbl>,
## #   `3.-Forage-proportion-in-the-diet-(kg/kg,-DM-basis)` <dbl>,
## #   `3.-Maize-forage-proportion-in-the-diet-(kg/kg,-DM-basis)` <dbl>, …
cdat = cdat %>% 
  mutate(
    prot_prod_kgperd = Tru.Prot*milk_kgperd/1000,
    nue_true = prot_prod_kgperd/6.34*1000 / n_intake_gperd,
    MUN = milk_urea * 0.46,
    MUNY_gperd = MUN*milk_kgperd/100,
    MUNY_MPY = MUNY_gperd / prot_prod_kgperd*100 ,
    MPY_MUNY = prot_prod_kgperd*100 / MUNY_gperd,
    MUN_MPY = MUN / prot_prod_kgperd,
    manure_n_gperd = urine_n_gperd + `5.-Faecal-N-excretion-(g/day)`,
    diet_cp = `4.-Diet-CP-concentration-(g/kg-DM)`,
    diet_cp_c = scale(diet_cp, scale = F),
    gest_wk = `2.-Gestation-week`, 
    lact_day = `2.-Lactation-week` * 7,
    lact_day_c = scale(lact_day, scale = F),
    bw = `2.-Body-weight-(kg)`,
    lact_stage = ifelse(lact_day > 175, "late", "early"),
    trt = `1.-Treatment-ID`)

only_cross_3trt = cdat %>% 
  group_by(exp, cow) %>% 
  summarise(n = n()) %>% 
  filter(n >2)
## `summarise()` has grouped output by 'exp'. You can override using the `.groups`
## argument.
# filter out fresh herbage diets
cdat1 = cdat %>% 
  right_join(only_cross_3trt) %>% 
  filter(`3.-Diet-type` != "Fresh_herbage") %>% 
  filter(`3.-Diet-type` != "Maize_Fresh_herbage")
## Joining, by = c("exp", "cow")

Show structure of this dataframe and SIGNIFICANT missingness for certain variables.

cr = cdat1 %>% 
  dplyr::select(cow, exp, nue_true, milk_kgperd, MUN, MUNY_gperd, MPY_MUNY, MUN_MPY, manure_n_gperd, diet_cp, lact_stage) 

str(cr)
## tibble [220 × 11] (S3: tbl_df/tbl/data.frame)
##  $ cow           : chr [1:220] "E07_C1" "E07_C2" "E07_C3" "E07_C4" ...
##  $ exp           : chr [1:220] "E07" "E07" "E07" "E07" ...
##  $ nue_true      : num [1:220] 0.317 0.289 0.268 0.258 0.27 ...
##  $ milk_kgperd   : num [1:220] 39.2 37.3 36.4 30.5 38.4 39.1 33.3 30.9 33.6 29.8 ...
##  $ MUN           : num [1:220] NA NA NA NA NA NA NA NA NA NA ...
##  $ MUNY_gperd    : num [1:220] NA NA NA NA NA NA NA NA NA NA ...
##  $ MPY_MUNY      : num [1:220] NA NA NA NA NA NA NA NA NA NA ...
##  $ MUN_MPY       : num [1:220] NA NA NA NA NA NA NA NA NA NA ...
##  $ manure_n_gperd: num [1:220] 360 379 395 351 414 345 361 387 348 377 ...
##  $ diet_cp       : num [1:220] 155 165 173 167 158 159 164 166 169 161 ...
##  $ lact_stage    : chr [1:220] "early" "early" "early" "early" ...
summary(cr)
##      cow                exp               nue_true        milk_kgperd   
##  Length:220         Length:220         Min.   :0.08225   Min.   : 5.50  
##  Class :character   Class :character   1st Qu.:0.22965   1st Qu.:21.93  
##  Mode  :character   Mode  :character   Median :0.26707   Median :27.45  
##                                        Mean   :0.26721   Mean   :26.68  
##                                        3rd Qu.:0.30255   3rd Qu.:32.30  
##                                        Max.   :0.42649   Max.   :47.00  
##                                                                         
##       MUN           MUNY_gperd        MPY_MUNY         MUN_MPY      
##  Min.   : 2.760   Min.   :0.7121   Min.   : 15.70   Min.   : 3.065  
##  1st Qu.: 7.130   1st Qu.:1.6916   1st Qu.: 23.65   1st Qu.: 7.883  
##  Median : 8.901   Median :2.6037   Median : 34.85   Median :10.419  
##  Mean   : 9.841   Mean   :2.8461   Mean   : 37.90   Mean   :11.382  
##  3rd Qu.:12.420   3rd Qu.:3.9500   3rd Qu.: 45.63   3rd Qu.:14.179  
##  Max.   :21.528   Max.   :5.9848   Max.   :126.45   Max.   :28.919  
##  NA's   :128      NA's   :128      NA's   :128      NA's   :128     
##  manure_n_gperd     diet_cp       lact_stage       
##  Min.   :168.0   Min.   :107.0   Length:220        
##  1st Qu.:265.8   1st Qu.:140.0   Class :character  
##  Median :317.5   Median :151.0   Mode  :character  
##  Mean   :323.1   Mean   :153.4                     
##  3rd Qu.:376.2   3rd Qu.:167.2                     
##  Max.   :661.0   Max.   :241.0                     
##

H1 procedures

Test correlations between variables

220 treatment means from 56 cows in crossover experiments. Methods do not account for pseudoreplication caused by within-cow repeated measurements.

98 treatment means for most varaibles due to missingness.

crs = cdat1 %>% 
  dplyr::select(nue_true, milk_kgperd, MUN, MUNY_gperd, MPY_MUNY, MUN_MPY, manure_n_gperd, diet_cp, lact_stage, cow) 

crs$cow %>% unique %>%  length
## [1] 56
crs %>% 
  dplyr::select(-cow, -lact_stage) %>% 
  ggpairs(lower = list(continuous = wrap(ggally_points, size = .2, color = "black"))) + theme_bw() + ggtitle("Fig. 1. Pairwise relations between variables in trt means (n = 98) \n in crossover experiments in CowNFlow. ")

Distribution gets skewed for MPY_MUNY. This is because the ratio involves the product of MPY * (MUNY^-1).

Regressions

Still using CowNFlow treatment means. In regressing NUE on MUN, MPY, and the MUN:MPY and MPY:MUN ratios, MUN was the best predictor. Conceptually, using a ratio as a predictor is the same as including main effects and an interaction between MPY and the inverse of MUN (1/MUN) and fixing the coefficients to equality. Therefore, I tested models including different combinations of MUN, MP%, MUNY, MPY, milk yield, to predict urinary N and NUE.

not shown, for conciseness

H1 conclusion

MUN is more predictive of NUE and urine N excretion than any of the MPY:MUNY ratios. Additionally, some of the ratios impose non-linear relationships and non-normal distributions, reducing their convenience compared to MUN.

Hypothesis 2 (H2)

  1. A certain level of MPY/MUNY can indicate the breakpoint where the milk protein yield response (Y) to additional dietary CP as a % of DM (X) is zero, i.e., where the slope is zero.

H2 procedures

Find data for a group of cows fed different dietary CP for at least 3-4 levels, where it would be possible to detect a quadratic trend in milk protein yield vs. dietary CP %.

I tried subsetting the crossover trials in the CowNFlow database, but in most cases the observations did not span broad enough CP to show a clear linear and quadratic response for the peak dietary CP. I have not tested this with covariates. Lactation stage or intial production may be important.

H2 conclusion

Inconclusive - unsure where to get enough data preferably within cow across several dietary CP levels to test this. Neither GIZ13 nor MAW636 works (?).

Hypothesis 3 (H3)

For GIZ13, test if rates of change in MUN or other related variables correlate with rates of change in NUE across experimental days within lactation stage. The implication would be that dynamic changes in certain variables could indicate declining NUE even without measuring DMI.

H3 procedures

For GIZ13, using the entire dataset (day by day).

dat <- read_excel("SAS.xlsx", sheet = "SAS")

trt = dat %>% 
  select(AnID, Lact, Trt) %>% 
mutate(AnID = as.character(AnID)) %>% 
  distinct(AnID, Lact, Trt)

# covariate
cvt = dat %>% 
  filter(Expwk == 1 | Expwk == 2) %>% 
  group_by(AnID) %>% 
  summarise(across(where(is.numeric), mean, .names = "cvt_{col}"))

dat = dat %>% 
  filter(!(Expwk == 1) ) %>% 
  filter(!(Expwk == 2)) %>% 
  mutate(MPY_MUNY_ratio = PrtYield*1000/MUNY,
         DIM_c = scale(DIM, center = T, scale = F)) %>% 
  group_by(AnID) %>% 
  mutate( DIM_1 = min(DIM),
          expday_c = scale(expday, scale = F),
          BW_c = scale(BW, scale = F)) %>% 
  filter(NUE < 50) %>% 
  left_join(cvt) %>% 
  mutate(cvt_BW_c = scale(cvt_BW, scale = F))
## Joining, by = "AnID"
split = dat %>% 
  split(~.$AnID)

Show one element of “split” (measurements for cow 5062)

str(split$`5062`)
## gropd_df [20 × 63] (S3: grouped_df/tbl_df/tbl/data.frame)
##  $ AnID          : num [1:20] 5062 5062 5062 5062 5062 ...
##  $ DIM           : num [1:20] 244 245 251 252 258 259 265 266 267 268 ...
##  $ Preg          : num [1:20] 19 20 26 27 33 34 40 41 42 43 ...
##  $ Lact          : chr [1:20] "4_Late" "4_Late" "4_Late" "4_Late" ...
##  $ Trt           : chr [1:20] "C" "C" "C" "C" ...
##  $ Expwk         : num [1:20] 3 3 4 4 5 5 6 6 6 6 ...
##  $ expday        : num [1:20] 16 17 23 24 30 31 37 38 39 40 ...
##  $ DailyMY       : num [1:20] 41.7 42 43.2 43 39.9 ...
##  $ FatYield      : num [1:20] 1.82 1.63 1.65 1.79 1.61 ...
##  $ PrtYield      : num [1:20] 1.31 1.35 1.4 1.4 1.33 ...
##  $ LactY         : num [1:20] 1.77 1.79 1.8 1.78 1.71 ...
##  $ SNFY          : num [1:20] 3.53 3.6 3.68 3.66 3.49 ...
##  $ SCC           : num [1:20] 17.1 27.4 58 36.4 27.2 ...
##  $ MUN           : num [1:20] 12.4 12.6 11.9 12.1 12.2 ...
##  $ FatPct        : num [1:20] 4.36 3.87 3.81 4.16 4.03 ...
##  $ PrtPct        : num [1:20] 3.13 3.2 3.24 3.25 3.34 ...
##  $ LactPct       : num [1:20] 4.25 4.25 4.16 4.13 4.28 ...
##  $ SNFPct        : num [1:20] 8.46 8.56 8.52 8.5 8.74 ...
##  $ MUNY          : num [1:20] 5.17 5.29 5.13 5.19 4.88 ...
##  $ DMI           : num [1:20] 30.7 28.6 28.4 32.1 33.1 ...
##  $ FPCM          : num [1:20] 42.9 40.9 41.9 43.5 40 ...
##  $ ECM           : num [1:20] 47.2 45.1 46.2 47.9 44.1 ...
##  $ NEL           : num [1:20] 31 29.5 30.1 31.3 29 ...
##  $ NI            : num [1:20] 0.823 0.765 0.761 0.86 0.886 ...
##  $ NUE           : num [1:20] 24.9 27.6 28.8 25.5 23.6 ...
##  $ FE            : num [1:20] 1.36 1.47 1.52 1.34 1.21 ...
##  $ FPCME         : num [1:20] 1.4 1.43 1.47 1.36 1.21 ...
##  $ ECME          : num [1:20] 1.54 1.58 1.63 1.49 1.33 ...
##  $ NELE          : num [1:20] 1.01 1.034 1.058 0.974 0.876 ...
##  $ BW            : num [1:20] 753 745 775 763 785 ...
##  $ MPY_MUNY_ratio: num [1:20] 252 255 273 270 273 ...
##  $ DIM_c         : num [1:20, 1] 47.7 48.7 54.7 55.7 61.7 ...
##   ..- attr(*, "scaled:center")= num 196
##  $ DIM_1         : num [1:20] 244 244 244 244 244 244 244 244 244 244 ...
##  $ expday_c      : num [1:20, 1] -27.5 -26.5 -20.5 -19.5 -13.5 -12.5 -6.5 -5.5 -4.5 -3.5 ...
##  $ BW_c          : num [1:20, 1] -30.19 -37.44 -7.51 -19.3 2.47 ...
##  $ cvt_DIM       : num [1:20] 236 236 236 236 236 ...
##  $ cvt_Preg      : num [1:20] 10.8 10.8 10.8 10.8 10.8 ...
##  $ cvt_Expwk     : num [1:20] 1.67 1.67 1.67 1.67 1.67 ...
##  $ cvt_expday    : num [1:20] 7.83 7.83 7.83 7.83 7.83 ...
##  $ cvt_DailyMY   : num [1:20] 40 40 40 40 40 ...
##  $ cvt_FatYield  : num [1:20] 1.72 1.72 1.72 1.72 1.72 ...
##  $ cvt_PrtYield  : num [1:20] 1.21 1.21 1.21 1.21 1.21 ...
##  $ cvt_LactY     : num [1:20] 1.74 1.74 1.74 1.74 1.74 ...
##  $ cvt_SNFY      : num [1:20] 3.39 3.39 3.39 3.39 3.39 ...
##  $ cvt_SCC       : num [1:20] 31.6 31.6 31.6 31.6 31.6 ...
##  $ cvt_MUN       : num [1:20] 11 11 11 11 11 ...
##  $ cvt_FatPct    : num [1:20] 4.3 4.3 4.3 4.3 4.3 ...
##  $ cvt_PrtPct    : num [1:20] 3.02 3.02 3.02 3.02 3.02 ...
##  $ cvt_LactPct   : num [1:20] 4.34 4.34 4.34 4.34 4.34 ...
##  $ cvt_SNFPct    : num [1:20] 8.49 8.49 8.49 8.49 8.49 ...
##  $ cvt_MUNY      : num [1:20] 4.41 4.41 4.41 4.41 4.41 ...
##  $ cvt_DMI       : num [1:20] 32.6 32.6 32.6 32.6 32.6 ...
##  $ cvt_FPCM      : num [1:20] 40.5 40.5 40.5 40.5 40.5 ...
##  $ cvt_ECM       : num [1:20] 44.6 44.6 44.6 44.6 44.6 ...
##  $ cvt_NEL       : num [1:20] 29.4 29.4 29.4 29.4 29.4 ...
##  $ cvt_NI        : num [1:20] 0.845 0.845 0.845 0.845 0.845 ...
##  $ cvt_NUE       : num [1:20] 22.5 22.5 22.5 22.5 22.5 ...
##  $ cvt_FE        : num [1:20] 1.23 1.23 1.23 1.23 1.23 ...
##  $ cvt_FPCME     : num [1:20] 1.24 1.24 1.24 1.24 1.24 ...
##  $ cvt_ECME      : num [1:20] 1.37 1.37 1.37 1.37 1.37 ...
##  $ cvt_NELE      : num [1:20] 0.903 0.903 0.903 0.903 0.903 ...
##  $ cvt_BW        : num [1:20] 732 732 732 732 732 ...
##  $ cvt_BW_c      : num [1:20, 1] 0 0 0 0 0 0 0 0 0 0 ...
##  - attr(*, "groups")= tibble [1 × 2] (S3: tbl_df/tbl/data.frame)
##   ..$ AnID : num 5062
##   ..$ .rows: list<int> [1:1] 
##   .. ..$ : int [1:20] 1 2 3 4 5 6 7 8 9 10 ...
##   .. ..@ ptype: int(0) 
##   ..- attr(*, ".drop")= logi TRUE
names = names(split)
mod_res = split %>% 
  purrr::map(., ~lm(cbind(NUE, MUN, MUNY, MPY_MUNY_ratio, DailyMY, FPCM, PrtPct)~expday_c, .)) %>% 
  purrr::map(., ~tidy(.)) %>% 
  bind_rows() %>% 
  mutate(AnID = rep(names, each = length(.$estimate)/length(names))) %>% 
 # filter(term!="(Intercept)") %>% 
  left_join(trt)
## Joining, by = "AnID"
# dat$expday_c %>% summary

day = -30:24

AnID = dat$AnID %>% unique

days = expand.grid(day = day, AnID = AnID) %>% mutate(AnID = as.character(AnID))
mod_res_y = mod_res %>% 
  dplyr::select(term, response, estimate, AnID, Lact, Trt) %>% 
  pivot_wider(id_cols = c(AnID, response, Trt, Lact), names_from = term, values_from = estimate) %>% 
  rename(int = `(Intercept)`, 
         sl = expday_c) %>% 
  left_join(days) %>% 
  mutate(y = sl*day + int)
## Joining, by = "AnID"
mod_res_y %>% 
  ggplot(aes(x = day, y = y, color = Lact, group = Lact)) + 
  geom_point(size = .1, alpha = .2) + 
  facet_wrap(~response, scales = "free") + 
  #theme(legend.position = "none") + 
  geom_smooth(method = "lm") + 
  ggtitle("Fig. 2. By lactation stage,  individual cow and mean linear rates of change across \n experimental days for selected variables in GIZ13")
## `geom_smooth()` using formula 'y ~ x'

mod_res_y %>% 
  ggplot(aes(x = day, y = y, color = Trt, group = Trt)) + 
  geom_point(size = .1, alpha = .2) + 
  facet_wrap(~response, scales = "free") + 
  #theme(legend.position = "none") + 
  geom_smooth(method = "lm")+ 
  ggtitle("Fig. 3. By CP level, individual cow and mean linear rates of change across \n experimental days for selected variables in GIZ13")
## `geom_smooth()` using formula 'y ~ x'

mod_res_y %>% 
  ggplot(aes(x = day, y = y, color = Trt, group = Trt)) + 
  geom_point(size = .1, alpha = .2) + 
  facet_grid(response~Lact, scales = "free") + 
  #theme(legend.position = "none") + 
  geom_smooth(method = "lm")+theme_classic() +
  ggtitle("Fig. 4. By CP level (color) and lactation stage (facet), \n individual cow and mean linear rate of change across experimental days for \n selected variables in GIZ13")
## `geom_smooth()` using formula 'y ~ x'

H3 conclusion

Delta NUE appears related to Delta MUNY in late lactation. At all stages of lactation Delta NUE is related to Delta milk yield. Therefore, the appropriate time to change diets (to optimize NUE) may be when milk yield starts to drop. Delta MUNY and Delta MUN may be useful indicators of dropping NUE in late lactation…more data would be useful to confirm this.

Hypothesis 4 (H4)

Using GIZ13 data (and possibly combining it with other similar RCBD studies such as Barros et al. 2017), we can demonstrate that grouping and separately feeding late lactation cows lower protein results in significant environmental benefits, i.e., reduction in urine urea N, manure volatile solids, predicted ammonia emissions, GHG emissions.

H4 procedures

Extension of RCBD designs in GIZ13 and other recent trials.

We could report some measured values (e.g., urine N), and some predicted (e.g., methane). This would describe the rate of increase in environmental impact across stages of lactation. For GIZ13 where we have BW change, we can also consider Productive N After accounting for change in BW and gestation, although these are only a small fraction of predicted N use.

H4 conclusions

The conclusions would recommend optimal late lactation feeding strategies for certain groups of cows to minimize environmental impact. For example, when to switch diets and how to do it (literature review or simulated diets in NASEM).

Other

Can now run NASEM in R on “batches” of simulated conditions with Excel input. For example, can simulate increasing stage of lactation or gestation day by day, or changing diets day-by-day.

inputs = read_excel("FuncInputsBWchange.xlsx")
inputs[4:15, 1:5] %>% 
  kable() %>% 
  row_spec(11:12, background = "yellow")
DietID “test1” “test2” “test3” “test4”
An_BW_mature 700 700 700 700
Trg_FrmGain 0 0 0 0
Trg_RsrvGain 0 0 0 0
An_Parity_rl 2 2 2 2
An_BCS 3.5 3.5 3.5 3.5
An_GestLength 283 283 283 283
An_AgeConcept1st NA NA NA NA
An_DIMConcept 100 101 102 103
Fet_BWbrth 40.909090909090907 40.909090909090907 40.909090909090907 40.909090909090907
An_StatePhys “Lactating Cow” “Lactating Cow” “Lactating Cow” “Lactating Cow”
An_LactDay 150 151 152 153
An_GestDay 50 51 52 53