---
title: "Manus 1 - Milk - LM lact 2"
author: "Anna Edvardsson Rasmussen"
date: '2022-08-23'
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")
```


# Lactation 2 dataset and preparation, Cross + farm X = excluded.
```{r}
setwd(".")
Lact.2.O<- read_excel("Lact.2.R.no.O.xlsx")
Lact.2.O$CI.plan<-as.factor(Lact.2.O$CI.plan)
Lact.2.O$Farm<-as.factor(Lact.2.O$Farm)
Lact.2.O$Breed<-as.factor(Lact.2.O$Breed)
Lact.2.O$DP.cat.2<-as.factor(Lact.2.O$DP.cat.2)
Lact.2.O$DP.40.2<-as.factor(Lact.2.O$DP.40.2)
Lact.2.O$DP.40.70.2<-as.factor(Lact.2.O$DP.40.70.2)
Lact.2.O$DP.70.2<-as.factor(Lact.2.O$DP.70.2)
Lact.2.O$VWP.compl.2.OK<-as.factor(Lact.2.O$VWP.compl.2.OK)
Lact.2.O$VWP.compl.MY.OK.2<-as.factor(Lact.2.O$VWP.compl.MY.OK.2)
Lact.2.O$VWP.notOK.MY.OK.2<-as.factor(Lact.2.O$VWP.notOK.MY.OK.2)
Lact.2.O$VWP.compl.MY.DO.OK.2<-as.factor(Lact.2.O$VWP.compl.MY.DO.OK.2)
Lact.2.O$VWP.acc.plan<-as.factor(Lact.2.O$VWP.acc.plan)
Lact.2.O$Complete.lact.2<-as.factor(Lact.2.O$Complete.lact.2)
Lact.2.O$MY.acc.plan.2<-as.factor(Lact.2.O$MY.acc.plan.2)
Lact.2.O$Min.DIM.OK<-as.factor(Lact.2.O$Min.DIM.OK)
Lact.2.O$Missing.OK<-as.factor(Lact.2.O$Missing.OK)
Lact.2.O$Max.DIM.OK<-as.factor(Lact.2.O$Max.DIM.OK)
```

```{r echo=FALSE, include=FALSE}
View(Lact.2.O)
```

# 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. CI.cI => CI
  
## Subsets with only VWPok+ daily MYok +only complete rows
  5. LL.cI => LL
  6. TM.50.cI => MY TM near 50 d before dry off
  
## Subset with only VWPok + complete lact + daily MYok + DOY = ok + only complete rows
  7.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 = 221)
  2. Farm (16 levels)
  3. VWP group (2 levels)
  6. Breed group (2 levels)
  20. Inclusion criteria: VWP = ok and complete lactation 1
  43. 305d yield (kg) 
  44. 305d yield (ECM) 
```{r}
Lact.2.cows<-c(1:3,6,14,16,18,20:22,25, 81:83)
Lact.2.cows<-Lact.2.O[Lact.2.cows]
summary(Lact.2.cows)
```
```{r}
table(Lact.2.cows$CI.plan, Lact.2.cows$VWP.acc.plan)
```

```{r}
table(Lact.2.cows$CI.plan, Lact.2.cows$Complete.lact.2)
```

```{r}
table(Lact.2.cows$CI.plan, Lact.2.cows$VWP.compl.2.OK)
```

```{r}
table(Lact.2.cows$CI.plan, Lact.2.cows$MY.acc.plan.2)
```

```{r}
table(Lact.2.cows$CI.plan, Lact.2.cows$ VWP.compl.MY.OK.2)
```

```{r}
table(Lact.2.cows$CI.plan, Lact.2.cows$VWP.acc.plan)
```

```{r}
table(Lact.2.cows$CI.plan, Lact.2.cows$Min.DIM.OK)
```

```{r}
table(Lact.2.cows$CI.plan, Lact.2.cows$Missing.OK)
```

```{r}
table(Lact.2.cows$CI.plan, Lact.2.cows$Max.DIM.OK)
```

```{r}
table(Lact.2.cows$CI.plan, Lact.2.cows$Breed)
```

```{r}
MY.305d.2.v<-c(1:3,6,20,43,44)
MY.305d.2<-Lact.2.O[MY.305d.2.v]
MY.305d.2.c<-MY.305d.2[complete.cases(MY.305d.2),]
MY.305d.2.cI<-MY.305d.2.c[which(MY.305d.2.c$VWP.compl.2.OK=="Yes"),]
summary(MY.305d.2.cI)
```

```{r}
table(MY.305d.2.cI$CI.plan, MY.305d.2.cI$VWP.compl.2.OK)
```

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

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

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

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

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

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

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

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

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

```{r}
model.LY.2.ECM<-lm(LactY.ECM.TM.2~CI.plan, data = LY.2.cI)
summary(model.LY.2.ECM)
anova(model.LY.2.ECM)
emmeans(model.LY.2.ECM, pairwise~CI.plan)
plot(model.LY.2.ECM)
qqnorm(residuals(model.LY.2.ECM))
```


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

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

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

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

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

```{r}
model.LY.CI.2.ECM<-lm(ECM.CI.day.2~CI.plan, data = LY.CI.2.cI)
summary(model.LY.CI.2.ECM)
anova(model.LY.CI.2.ECM)
emmeans(model.LY.CI.2.ECM, pairwise~CI.plan)
plot(model.LY.CI.2.ECM)
qqnorm(residuals(model.LY.CI.2.ECM))
```

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

# Models from DOY dataset
```{r}
model.DOY.2<-lmer(DO.MY2~CI.plan+Breed+(1|Farm)+CI.plan*Breed, data = DOY.2.cI)
summary(model.DOY.2)
anova(model.DOY.2)
emmeans(model.DOY.2, pairwise~CI.plan+Breed)
```  
Tendency to interaction..

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