---
title: "Manus 1 - Milk final models 2022-06-01"
author: "Anna Edvardsson Rasmussen"
date: '2022-06-01'
output: word_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

```{r echo=TRUE, include=FALSE}
library("reshape2")
library("plyr")
library("readxl")
library("magrittr")
library("data.table")
library("lubridate")
library("broom")
library("tidyverse")
library("tidylog")
library("lme4")
library("lmerTest")
library("emmeans")
library("car")
library("xtable")
library("sjPlot")
```

  These are the models used for the results section for manus 1 milk.
  Cross are excluded and ECM yield corrected and MY/2CI day corrected.

# Lactation 1 dataset and preparation
```{r echo=TRUE, include=FALSE}
setwd(".")
Lact.1.full.ny<- read_excel("Lact.1.R.ny.xlsx")
Lact.1.full.ny$CI.plan<-as.factor(Lact.1.full.ny$CI.plan)
Lact.1.full.ny$Farm<-as.factor(Lact.1.full.ny$Farm)
Lact.1.full.ny$Breed<-as.factor(Lact.1.full.ny$Breed)
Lact.1.full.ny$DP.cat<-as.factor(Lact.1.full.ny$DP.cat.1)
Lact.1.full.ny$DP.40<-as.factor(Lact.1.full.ny$DP.40.1)
Lact.1.full.ny$DP.40.70<-as.factor(Lact.1.full.ny$DP.40.70.1)
Lact.1.full.ny$DP.70<-as.factor(Lact.1.full.ny$DP.70.1)
Lact.1.full.ny$Excl.O<-as.factor(Lact.1.full.ny$Excl.O)
Lact.1.full.ny$VWP.complete.1.ok<-as.factor(Lact.1.full.ny$VWP.complete.1.ok)
Lact.1.full.ny$VWP.compl.MY.OK.1<-as.factor(Lact.1.full.ny$VWP.compl.MY.OK.1)
Lact.1.full.ny$VWP.notOK.MY.OK.1<-as.factor(Lact.1.full.ny$VWP.notOK.MY.OK.1)
Lact.1.full.ny$VWP.compl.MY.DO.OK.1<-as.factor(Lact.1.full.ny$VWP.compl.MY.DO.OK.1)
Lact.1.full.ny$VWP.complete.2lact<-as.factor(Lact.1.full.ny$VWP.complete.2lact)
```

```{r echo=FALSE, include=FALSE}
View(Lact.1.full.ny)
```

# Subset data for each investigated variable and apply correct inclusion criteria
Subset with only complete rows  (due to different number of observation for each parameter). 

## Subsets with only VWPok+only complete rows
  1. MY.305d.cI => 305d yield (kg + ECM), 
  2. LY.cI => Whole lactation yield (kg + ECM)
  3. LY.CI.cI => MY/CI day (kg + ECM)
  4. MYCId.2lact.cI => MY/2CI day (kg + ECM)
  5. CI.cI => CI
  
## Subsets with only VWPok+ daily MYok +only complete rows
  6. LL.cI => LL
  7. TM.50.cI => MY TM near 50 d before dry off
  
## Subset with only VWPok + complete lact + daily MYok + DOY = ok + only complete rows
  8.DOY.cI => DOY (kg)
  
## Model description:
  Breed: fix - 2 levels  
  Farm: random - 16 levels  
  Group: fix - 2 levels  
  Breed:group: fix  

  If the interaction was not significant a second model without interaction was used
  Residual and normal QQ-plots checked for each final model

# 305 d yield dataset:
  1. CowID (n = 347)
  2. Farm (16 levels)
  3. VWP group (2 levels)
  6. Breed group (2 levels)
  20. Inclusion criteria: VWP = ok and complete lactation 1
  40. 305d yield (kg) 
  41. 305d yield (ECM) 
```{r}
MY.305d.v<-c(1:3,6,20,40,41)
MY.305d<-Lact.1.full.ny[MY.305d.v]
MY.305d.c<-MY.305d[complete.cases(MY.305d),]
MY.305d.cI<-MY.305d.c[which(MY.305d.c$VWP.complete.1.ok=="Yes"),]
summary(MY.305d.cI)
```

# Models from 305.d dataset
```{r}
model.305.kg<-lmer(Lact.305.kg.VXA1~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = MY.305d.cI)
summary(model.305.kg)
anova(model.305.kg)
emmeans(model.305.kg, pairwise~CI.plan)
emmeans(model.305.kg, pairwise~Breed)
plot(model.305.kg)
qqnorm(residuals(model.305.kg))
```
Interaction not significant => model 2

```{r}
model.305.kg2<-lmer(Lact.305.kg.VXA1~CI.plan+Breed+(1|Farm), data = MY.305d.cI)
summary(model.305.kg2)
anova(model.305.kg2)
emmeans(model.305.kg2, pairwise~CI.plan)
emmeans(model.305.kg2, pairwise~Breed)
plot(model.305.kg2)
qqnorm(residuals(model.305.kg2))
```

```{r}
model.305.ECM<-lmer(Lact.305.ECM.VXA1~CI.plan+Breed+(1|Farm)+ CI.plan*Breed, data = MY.305d.cI)
summary(model.305.ECM)
anova(model.305.ECM)
emmeans(model.305.ECM, pairwise~CI.plan)
emmeans(model.305.ECM, pairwise~Breed)
plot(model.305.ECM)
qqnorm(residuals(model.305.ECM))
```
Interaction not significant => model 2

```{r}
model.305.ECM2<-lmer(Lact.305.ECM.VXA1~CI.plan+Breed+(1|Farm), data = MY.305d.cI)
summary(model.305.ECM2)
anova(model.305.ECM2)
emmeans(model.305.ECM2, pairwise~CI.plan)
emmeans(model.305.ECM2, pairwise~Breed)
plot(model.305.ECM2)
qqnorm(residuals(model.305.ECM2))
```

# Whole lactation yield dataset:
  1. CowID (n = 349)
  2. Farm (16 levels)
  3. VWP group (2 levels)
  6. Breed group (2 levels)
  20. Inclusion criteria: VWP = ok and complete lactation 1
  42. Whole lactation yield (kg) from TM
  43. Whole lactation yield (ECM) from TM
```{r}
LY.v<-c(1:3,6,20,42,43)
LY<-Lact.1.full.ny[LY.v]
LY.c<-LY[complete.cases(LY),]
LY.cI<-LY.c[which(LY.c$VWP.complete.1.ok=="Yes"),]
summary(LY.cI)
```

# Models from LY dataset
```{r}
model.LY.kg<-lmer(LactY.kg.TM.1~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = LY.cI)
summary(model.LY.kg)
anova(model.LY.kg)
emmeans(model.LY.kg, pairwise~CI.plan)
emmeans(model.LY.kg, pairwise~Breed)
plot(model.LY.kg)
qqnorm(residuals(model.LY.kg))
```
Interaction not significant => model 2

```{r}
model.LY.kg2<-lmer(LactY.kg.TM.1~CI.plan+Breed+(1|Farm), data = LY.cI)
summary(model.LY.kg2)
anova(model.LY.kg2)
emmeans(model.LY.kg2, pairwise~CI.plan)
emmeans(model.LY.kg2, pairwise~Breed)
plot(model.LY.kg2)
qqnorm(residuals(model.LY.kg2))
```

```{r}
model.LY.ECM2<-lmer(LactY.ECM.TM1~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = LY.cI)
summary(model.LY.ECM2)
anova(model.LY.ECM2)
emmeans(model.LY.ECM2, pairwise~CI.plan+Breed)
emmeans(model.LY.ECM2, pairwise~CI.plan)
emmeans(model.LY.ECM2, pairwise~Breed)
plot(model.LY.ECM2)
qqnorm(residuals(model.LY.ECM2))
```
Interaction not significant => model 2

```{r}
model.LY.ECM<-lmer(LactY.ECM.TM1~CI.plan+Breed+(1|Farm), data = LY.cI)
summary(model.LY.ECM)
anova(model.LY.ECM)
emmeans(model.LY.ECM, pairwise~CI.plan)
emmeans(model.LY.ECM, pairwise~Breed)
plot(model.LY.ECM)
qqnorm(residuals(model.LY.ECM))
```

# Average yeld per CI day dataset:
  1. CowID (n = 349)
  2. Farm (16 levels)
  3. VWP group (2 levels)
  6. Breed group (2 levels)
  20. Inclusion criteria: VWP = ok and complete lactation 1
  44. Whole lactation yield (ECM) from TM/days in the first CI
  45. Whole lactation yield (kg) from TM/days in the first CI
```{r}
LY.CI.v<-c(1:3,6,20,44,45)
LY.CI<-Lact.1.full.ny[LY.CI.v]
LY.CI.c<-LY.CI[complete.cases(LY.CI),]
LY.CI.cI<-LY.CI.c[which(LY.CI.c$VWP.complete.1.ok=="Yes"),]
summary(LY.CI.cI)
```

# Models from LY/CI dataset
```{r}
model.LY.CI.kg<-lmer(Kg.CI.day.TM.1~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = LY.CI.cI)
summary(model.LY.CI.kg)
anova(model.LY.CI.kg)
emmeans(model.LY.CI.kg, pairwise~CI.plan)
emmeans(model.LY.CI.kg, pairwise~Breed)
plot(model.LY.CI.kg)
qqnorm(residuals(model.LY.CI.kg))
```
Interaction not significant => model 2

```{r}
model.LY.CI.kg2<-lmer(Kg.CI.day.TM.1~CI.plan+Breed+(1|Farm), data = LY.CI.cI)
summary(model.LY.CI.kg2)
anova(model.LY.CI.kg2)
emmeans(model.LY.CI.kg2, pairwise~CI.plan)
emmeans(model.LY.CI.kg2, pairwise~Breed)
plot(model.LY.CI.kg2)
qqnorm(residuals(model.LY.CI.kg2))
```

```{r}
model.LY.CI.ECM2<-lmer(ECM.CI.day.TM.1~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = LY.CI.cI)
summary(model.LY.CI.ECM2)
anova(model.LY.CI.ECM2)
emmeans(model.LY.CI.ECM2, pairwise~CI.plan+Breed)
emmeans(model.LY.CI.ECM2, pairwise~CI.plan)
emmeans(model.LY.CI.ECM2, pairwise~Breed)
emmeans(model.LY.CI.ECM2, pairwise~Breed+CI.plan)
```
Interaction not significant => model 2

```{r}
model.LY.CI.ECM<-lmer(ECM.CI.day.TM.1~CI.plan+Breed+(1|Farm), data = LY.CI.cI)
summary(model.LY.CI.ECM)
anova(model.LY.CI.ECM)
emmeans(model.LY.CI.ECM, pairwise~CI.plan)
emmeans(model.LY.CI.ECM, pairwise~Breed)
plot(model.LY.CI.ECM)
qqnorm(residuals(model.LY.CI.ECM))
```

# Average yeld per day in 2 CI dataset:
  1. CowID (n = 245)
  2. Farm (16 levels)
  3. VWP group (2 levels)
  6. Breed group (2 levels)
  65. Inclusion criteria: VWP = ok and complete lactation 2
  46. Sum of 2 lactations yield (kg) from TM/days in both CI
  47. Sum of 2 lactations yield (ECM) from TM/days in both CI
```{r}
MYCId.2lact.v2.ny<-c(1:3,6,65,46,47)
MYCId.2lact2.ny<-Lact.1.full.ny[MYCId.2lact.v2.ny]
MYCId.2lact.c2.ny<-MYCId.2lact2.ny[complete.cases(MYCId.2lact2.ny),]
MYCId.2lact.cI2.ny<-MYCId.2lact.c2.ny[which(MYCId.2lact.c2.ny$VWP.complete.2lact=="Yes"),]
summary(MYCId.2lact.cI2.ny)
```

# Models from 2LY/2CI dataset
```{r}
model.2LY.2CI.ECM2<-lmer(LY2.ECM.2CI.d~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = MYCId.2lact.cI2.ny)
summary(model.2LY.2CI.ECM2)
anova(model.2LY.2CI.ECM2)
emmeans(model.2LY.2CI.ECM2, pairwise~CI.plan+ Breed)
```

```{r}
model.2LY.2CI.ECM<-lmer(LY2.ECM.2CI.d~CI.plan+Breed+(1|Farm), data = MYCId.2lact.cI2.ny)
summary(model.2LY.2CI.ECM)
anova(model.2LY.2CI.ECM)
emmeans(model.2LY.2CI.ECM, pairwise~CI.plan)
emmeans(model.2LY.2CI.ECM, pairwise~Breed)
plot(model.2LY.2CI.ECM)
qqnorm(residuals(model.2LY.2CI.ECM))
```

```{r}
model.2LY.2CI.kg2<-lmer(LY2.kg.2CI.d~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = MYCId.2lact.cI2.ny)
summary(model.2LY.2CI.kg2)
anova(model.2LY.2CI.kg2)
emmeans(model.2LY.2CI.kg2, pairwise~CI.plan+Breed)
```

```{r}
model.2LY.2CI.kg<-lmer(LY2.kg.2CI.d~CI.plan+Breed+(1|Farm), data = MYCId.2lact.cI2.ny)
summary(model.2LY.2CI.kg)
anova(model.2LY.2CI.kg)
emmeans(model.2LY.2CI.kg, pairwise~CI.plan)
emmeans(model.2LY.2CI.kg, pairwise~Breed)
plot(model.2LY.2CI.kg)
qqnorm(residuals(model.2LY.2CI.kg))
```


# Calving interval length dataset:
  1. CowID (n = 349)
  2. Farm (16 levels)
  3. VWP group (2 levels)
  6. Breed group (2 levels)
  20. Inclusion criteria: VWP = ok and complete lactation 1
  49. First calving interval
```{r}
CI.v<-c(1:3,6,20,49)
CI<-Lact.1.full.ny[CI.v]
CI.c<-CI[complete.cases(CI),]
CI.cI<-CI.c[which(CI.c$VWP.complete.1.ok=="Yes"),]
summary(CI.cI)
```

# Models from CI dataset
```{r}
model.CI<-lmer(CI.1~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = CI.cI)
summary(model.CI)
anova(model.CI)
emmeans(model.CI, pairwise~CI.plan+ Breed)
```
Interaction not significant => model 2

```{r}
model.CI2<-lmer(CI.1~CI.plan+Breed+(1|Farm), data = CI.cI)
summary(model.CI2)
anova(model.CI2)
emmeans(model.CI2, pairwise~CI.plan)
emmeans(model.CI2, pairwise~Breed)
plot(model.CI2)
qqnorm(residuals(model.CI2))
```

# Lactation length dataset:
  1. CowID (n = 320)
  2. Farm (16 levels)
  3. VWP group (2 levels)
  6. Breed group (2 levels)
  21. Inclusion criteria: VWP = ok, complete lactation 1 and daily MY = ok
  27. Lactation length first lactation
  28. Dry period length first lactation
```{r}
LL.v<-c(1:3,6,21,27,28)
LL<-Lact.1.full.ny[LL.v]
LL.c<-LL[complete.cases(LL),]
LL.cI<-LL.c[which(LL.c$VWP.compl.MY.OK.1=="Yes"),] 
summary(LL.cI)
```

# Models from LL dataset
```{r}
model.LL<-lmer(Lact.lenght.1~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = LL.cI)
summary(model.LL)
anova(model.LL)
emmeans(model.LL, pairwise~CI.plan)
emmeans(model.LL, pairwise~Breed)
```
Interaction not significant => model 2

```{r}
model.LL2<-lmer(Lact.lenght.1~CI.plan+Breed+(1|Farm), data = LL.cI)
summary(model.LL2)
anova(model.LL2)
emmeans(model.LL2, pairwise~CI.plan)
emmeans(model.LL2, pairwise~Breed)
plot(model.LL2)
qqnorm(residuals(model.LL2))
```

```{r}
model.DPL<-lmer(Dry.period1~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = LL.cI)
summary(model.DPL)
anova(model.DPL)
emmeans(model.DPL, pairwise~CI.plan)
emmeans(model.DPL, pairwise~Breed)
plot(model.DPL)
qqnorm(residuals(model.DPL))
```
Interaction not significant => model 2

```{r}
model.DPL2<-lmer(Dry.period1~CI.plan+Breed+(1|Farm), data = LL.cI)
summary(model.DPL2)
anova(model.DPL2)
emmeans(model.DPL2, pairwise~CI.plan)
emmeans(model.DPL2, pairwise~Breed)
plot(model.DPL2)
qqnorm(residuals(model.DPL2))
```

# Test milking 50 d before dry off dataset:
  1. CowID (n = 285)
  2. Farm (16 levels)
  3. VWP group (2 levels)
  6. Breed group (2 levels)
  21. Inclusion criteria: VWP = ok, complete lactation 1 and daily MY = ok
  36. TM between 50 and 20 before dry off
```{r}
TM.50.v<-c(1:3,6,21,36)
TM.50<-Lact.1.full.ny[TM.50.v]
TM.50.c<-TM.50[complete.cases(TM.50),]
TM.50.cI<-TM.50.c[which(TM.50.c$VWP.compl.MY.OK.1=="Yes"),] 
summary(TM.50.cI)
```

# Models from TM.50 dataset
```{r}
model.TM.50<-lmer(TM.near.50.1~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = TM.50.cI)
summary(model.TM.50)
anova(model.TM.50)
emmeans(model.TM.50, pairwise~CI.plan)
emmeans(model.TM.50, pairwise~Breed)
plot(model.TM.50)
qqnorm(residuals(model.TM.50))
```
Interaction not significant => model 2

```{r}
model.TM.502<-lmer(TM.near.50.1~CI.plan+Breed+(1|Farm), data = TM.50.cI)
summary(model.TM.502)
anova(model.TM.502)
emmeans(model.TM.502, pairwise~CI.plan)
emmeans(model.TM.502, pairwise~Breed)
plot(model.TM.502)
qqnorm(residuals(model.TM.502))
```

# MY day 4.33 for cows with VWP NOT ok
  1. CowID (n = 88)
  2. Farm (16 levels)
  3. VWP group (2 levels)
  6. Breed group (2 levels)
  25. Inclusion criteria: VWP = NOT ok
  48. Mean MY day 4-33 in lactation
```{r}
MY.4.33.v<-c(1:3,6,25,48)
MY.4.33<-Lact.1.full.ny[MY.4.33.v]
MY.4.33.c<-MY.4.33[complete.cases(MY.4.33),]
MY.4.33.cI<-MY.4.33.c[which(MY.4.33.c$VWP.notOK.MY.OK.1=="Yes"),] 
summary(MY.4.33.cI)
```

# Models from MY.4.33 dataset
```{r}
model.MY.4.33<-lmer(MY.4.33.1~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = MY.4.33.cI)
summary(model.MY.4.33)
anova(model.MY.4.33)
emmeans(model.MY.4.33, pairwise~CI.plan)
emmeans(model.MY.4.33, pairwise~Breed)
emmeans(model.MY.4.33, pairwise~Breed+CI.plan)
plot(model.MY.4.33)
qqnorm(residuals(model.MY.4.33))
```
Interaction not significant => model 2

```{r}
model.MY.4.33.2<-lmer(MY.4.33.1~CI.plan+Breed+(1|Farm), data = MY.4.33.cI)
summary(model.MY.4.33.2)
anova(model.MY.4.33.2)
emmeans(model.MY.4.33.2, pairwise~CI.plan)
emmeans(model.MY.4.33.2, pairwise~Breed)
plot(model.MY.4.33.2)
qqnorm(residuals(model.MY.4.33.2))
```

# Dry off yield:
  1. CowID (n = 166)
  2. Farm (16 levels)
  3. VWP group (2 levels)
  6. Breed group (2 levels)
  22. Inclusion criteria: VWP = ok, complete lactation 1, daily MY = ok and dry off date = ok
  33. Dry off yield
```{r}
DOY.v<-c(1:3,6,22,33)
DOY<-Lact.1.full.ny[DOY.v]
DOY.c<-DOY[complete.cases(DOY),]
DOY.cI<-DOY.c[which(DOY.c$VWP.compl.MY.DO.OK.1=="Yes"),]
summary(DOY.cI)
```

# Models from DOY dataset
```{r}
model.DOY<-lmer(DO.MY.1~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = DOY.cI)
summary(model.DOY)
anova(model.DOY)
emmeans(model.DOY, pairwise~CI.plan)
emmeans(model.DOY, pairwise~Breed)
plot(model.DOY)
qqnorm(residuals(model.DOY))
```  

```{r}
model.DOY2<-lmer(DO.MY.1~CI.plan+Breed+(1|Farm), data = DOY.cI)
summary(model.DOY2)
anova(model.DOY2)
emmeans(model.DOY2, pairwise~CI.plan)
emmeans(model.DOY2, pairwise~Breed)
plot(model.DOY2)
qqnorm(residuals(model.DOY2))
```