---
title: "Part 2 - linear models - Milk"
author: "Anna Edvardsson Rasmussen"
date: "2023-02-21"
output: word_document
---
# Packages and setup
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

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("writexl")
```

# Dataset base
```{r}
Part.2<-read_excel("data/Part.2.R.farm.rev.xlsx")
Part.2$Farm<-as.factor(Part.2$Farm)
Part.2$Breed<-as.factor(Part.2$Breed)
Part.2$VWP.ok<-as.factor(Part.2$VWP.ok) 
Part.2$VWP.ok.10<-as.factor(Part.2$VWP.ok.10) 
Part.2$Complete.1<-as.factor(Part.2$Complete.1)
Part.2$Trt.group<-as.factor(Part.2$Trt.group)
Part.2$Not.ins<-as.factor(Part.2$Not.ins)
Part.2$Treat<-as.factor(Part.2$Treat)
Part.2<-Part.2[which(Part.2$Breed!="Cross"),] # Cross exkluderade (n = 75)

# Calculate LL and DPL
Part.2$LL<-as.numeric(Part.2$DryDate-Part.2$Calvingdate.1)
Part.2$DPL<-as.numeric(Part.2$CInt.1-Part.2$LL)
Part.2$Pers.DO<-(Part.2$Dry.yield-Part.2$KgMilk.3)/(Part.2$Dry.DIM-Part.2$TM.DIM.3)
Part.2$Pers.3.8<-(Part.2$KgMilk.8-Part.2$KgMilk.3)/(Part.2$TM.DIM.8-Part.2$TM.DIM.3)
Part.2$DPL.kat<-as.factor(Part.2$DPL.kat)
Part.2$DPL.40<-as.factor(Part.2$DPL.40)
Part.2$DPL.70<-as.factor(Part.2$DPL.70)
Part.2$DPL.40.70<-as.factor(Part.2$DPL.40.70)
```
DPL<-Part.2[c(1,118,119)]
write_xlsx(DPL, "DPL.xlsx")



# 305d Dataset
summary(MY.305d.1) # ITT
summary(MY.305d.2) # VWP ok
```{r}
MY.305d.1<-Part.2[c(1:5,6,7,18,19,36:37,53,55,57:61,85,86,90,91,93:97,121:125, 129:130)]
MY.305d.2<-MY.305d.1[which(MY.305d.1$VWP.ok.10=="Yes"),]
MY.305d<-MY.305d.2[which(MY.305d.2$Complete.1=="Yes"),]
summary(MY.305d) # VWP ok + complete lact.1
```

## view
```{r}
view(Part.2)
```

#Subset data for each investigated variable and apply correct inclusion criteria
Subset with only complete rows  (due to different number of observations 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: random? - 3 levels  
  Farm: random - 16 levels  
  CI.Group: fix - 2 levels  
  Trt.group: fix - 3 levels
  Interaction: CI.group*trt.group - 7 levels
   
  Residual and normal QQ-plots checked for each final model

## MY dataset:
  1. CowID (n = 605)
  2. Farm (18 levels)
  3. CI.group (2 levels)
  6. Breed.group (3 levels)
  20. Inclusion criteria: VWP = ok (n =432) And complete lactation 1 (n = 342)
  36. Group: 12/12, 16/12 or 16/16
  40. 305d yield (kg) 


# 305d Model
```{r}
model.305.kg<-lmer(MY.305d.1~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.305.kg)
anova(model.305.kg)
emmeans(model.305.kg, pairwise~Treat)
plot(model.305.kg)
qqnorm(residuals(model.305.kg))
```
## output as table
```{r}
table1<-as.data.frame(emmeans(model.305.kg, pairwise~Treat))
table2<-head(table1[c(1,3,4)],3)
table3<-pivot_wider(table2, names_from = Treat,values_from = c(emmean, SE))
table4<-table3[c(1,4,2,5,3,6)]

emmeans_output <- emmeans(model.305.kg, pairwise ~ Treat)
pairwise_contrasts <- summary(emmeans_output)$contrasts
table_contrasts1<-as.data.frame(pairwise_contrasts)
table_contrasts2<-table_contrasts1[c(1,6)]
table_contrasts3<-pivot_wider(table_contrasts2, names_from = c(contrast),values_from = c(p.value))
table_contrasts3

t.305d<-cbind(table4,table_contrasts3)
t.305d$model<-"305d MY"
t.305d<-t.305d[c(10,1:9)]
t.305d
```

# WL kg Model
```{r}
model.LY.kg<-lmer(LY.1~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.LY.kg)
anova(model.LY.kg)
emmeans(model.LY.kg, pairwise~Treat)
plot(model.LY.kg)
qqnorm(residuals(model.LY.kg))
```
## output as table
```{r}
table1.LY<-as.data.frame(emmeans(model.LY.kg, pairwise~Treat))
table2.LY<-head(table1.LY[c(1,3,4)],3)
table3.LY<-pivot_wider(table2.LY, names_from = Treat,values_from = c(emmean, SE))
table4.LY<-table3.LY[c(1,4,2,5,3,6)]

emmeans_output.LY <- emmeans(model.LY.kg, pairwise ~ Treat)
pairwise_contrasts.LY <- summary(emmeans_output.LY)$contrasts
table_contrasts1.LY<-as.data.frame(pairwise_contrasts.LY)
table_contrasts2.LY<-table_contrasts1.LY[c(1,6)]
table_contrasts3.LY<-pivot_wider(table_contrasts2.LY, names_from = c(contrast),values_from = c(p.value))

t.LY<-cbind(table4.LY,table_contrasts3.LY)
t.LY$model<-"LY kg"
t.LY<-t.LY[c(10,1:9)]
t.LY
```

# WL ECM Model
```{r}
model.LY.ECM<-lmer(ECM.LY.1~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.LY.ECM)
anova(model.LY.ECM)
emmeans(model.LY.ECM, pairwise~Treat)
plot(model.LY.ECM)
qqnorm(residuals(model.LY.ECM))
```

## output as table
```{r}
table1.LY.ECM<-as.data.frame(emmeans(model.LY.ECM, pairwise~Treat))
table2.LY.ECM<-head(table1.LY.ECM[c(1,3,4)],3)
table3.LY.ECM<-pivot_wider(table2.LY.ECM, names_from = Treat,values_from = c(emmean, SE))
table4.LY.ECM<-table3.LY.ECM[c(1,4,2,5,3,6)]

emmeans_output.LY.ECM<- emmeans(model.LY.ECM, pairwise ~ Treat)
pairwise_contrasts.LY.ECM<- summary(emmeans_output.LY.ECM)$contrasts
table_contrasts1.LY.ECM<-as.data.frame(pairwise_contrasts.LY.ECM)
table_contrasts2.LY.ECM<-table_contrasts1.LY.ECM[c(1,6)]
table_contrasts3.LY.ECM<-pivot_wider(table_contrasts2.LY.ECM, names_from = c(contrast),values_from = c(p.value))

t.LY.ECM<-cbind(table4.LY.ECM,table_contrasts3.LY.ECM)
t.LY.ECM$model<-"LY ECM"
t.LY.ECM<-t.LY.ECM[c(10,1:9)]
t.LY.ECM
```

# MY/day kg Model
```{r}
model.MY.day<-lmer(MY.day.1~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.MY.day)
anova(model.MY.day)
emmeans(model.MY.day, pairwise~Treat)
plot(model.MY.day)
qqnorm(residuals(model.MY.day))
```

## output as table
```{r}
table1.MY.day<-as.data.frame(emmeans(model.MY.day, pairwise~Treat))
table2.MY.day<-head(table1.MY.day[c(1,3,4)],3)
table3.MY.day<-pivot_wider(table2.MY.day, names_from = Treat,values_from = c(emmean, SE))
table4.MY.day<-table3.MY.day[c(1,4,2,5,3,6)]

emmeans_output.MY.day<- emmeans(model.MY.day, pairwise ~ Treat)
pairwise_contrasts.MY.day<- summary(emmeans_output.MY.day)$contrasts
table_contrasts1.MY.day<-as.data.frame(pairwise_contrasts.MY.day)
table_contrasts2.MY.day<-table_contrasts1.MY.day[c(1,6)]
table_contrasts3.MY.day<-pivot_wider(table_contrasts2.MY.day, names_from = c(contrast),values_from = c(p.value))

t.MY.day<-cbind(table4.MY.day,table_contrasts3.MY.day)
t.MY.day$model<-"MY/day kg"
t.MY.day<-t.MY.day[c(10,1:9)]
t.MY.day
```
# MY/day ECM Model
```{r}
model.MY.d.ECM<-lmer(ECMY.day.1~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.MY.d.ECM)
anova(model.MY.d.ECM)
emmeans(model.MY.d.ECM, pairwise~Treat)
plot(model.MY.d.ECM)
qqnorm(residuals(model.MY.d.ECM))
```

## output as table
```{r}
table1.MY.d.ECM<-as.data.frame(emmeans(model.MY.d.ECM, pairwise~Treat))
table2.MY.d.ECM<-head(table1.MY.d.ECM[c(1,3,4)],3)
table3.MY.d.ECM<-pivot_wider(table2.MY.d.ECM, names_from = Treat,values_from = c(emmean, SE))
table4.MY.d.ECM<-table3.MY.d.ECM[c(1,4,2,5,3,6)]

emmeans_output.MY.d.ECM<- emmeans(model.MY.d.ECM, pairwise ~ Treat)
pairwise_contrasts.MY.d.ECM<- summary(emmeans_output.MY.d.ECM)$contrasts
table_contrasts1.MY.d.ECM<-as.data.frame(pairwise_contrasts.MY.d.ECM)
table_contrasts2.MY.d.ECM<-table_contrasts1.MY.d.ECM[c(1,6)]
table_contrasts3.MY.d.ECM<-pivot_wider(table_contrasts2.MY.d.ECM, names_from = c(contrast),values_from = c(p.value))

t.MY.d.ECM<-cbind(table4.MY.d.ECM,table_contrasts3.MY.d.ECM)
t.MY.d.ECM$model<-"MY/Day ECM"
t.MY.d.ECM<-t.MY.d.ECM[c(10,1:9)]
t.MY.d.ECM
```
# CInt Model
```{r}
model.CInt<-lmer(CInt.1~(Treat)+(1|Farm)+Breed, data = MY.305d)
summary(model.CInt)
anova(model.CInt)
emmeans(model.CInt, pairwise~Treat)
plot(model.CInt)
qqnorm(residuals(model.CInt))
```

## output as table
```{r}
table1.CInt<-as.data.frame(emmeans(model.CInt, pairwise~Treat))
table2.CInt<-head(table1.CInt[c(1,3,4)],3)
table3.CInt<-pivot_wider(table2.CInt, names_from = Treat,values_from = c(emmean, SE))
table4.CInt<-table3.CInt[c(1,4,2,5,3,6)]

emmeans_output.CInt<- emmeans(model.CInt, pairwise ~ Treat)
pairwise_contrasts.CInt<- summary(emmeans_output.CInt)$contrasts
table_contrasts1.CInt<-as.data.frame(pairwise_contrasts.CInt)
table_contrasts2.CInt<-table_contrasts1.CInt[c(1,6)]
table_contrasts3.CInt<-pivot_wider(table_contrasts2.CInt, names_from = c(contrast),values_from = c(p.value))

t.CInt<-cbind(table4.CInt,table_contrasts3.CInt)
t.CInt$model<-"CInt"
t.CInt<-t.CInt[c(10,1:9)]
t.CInt
```

# CInt Model (months)
```{r}
model.CInt<-lmer((CInt.1/30.4368499)~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.CInt)
anova(model.CInt)
emmeans(model.CInt, pairwise~Treat)
plot(model.CInt)
qqnorm(residuals(model.CInt))
```

## output as table
```{r}
table1.CInt<-as.data.frame(emmeans(model.CInt, pairwise~Treat))
table2.CInt<-head(table1.CInt[c(1,3,4)],3)
table3.CInt<-pivot_wider(table2.CInt, names_from = Treat,values_from = c(emmean, SE))
table4.CInt<-table3.CInt[c(1,4,2,5,3,6)]

emmeans_output.CInt<- emmeans(model.CInt, pairwise ~ Treat)
pairwise_contrasts.CInt<- summary(emmeans_output.CInt)$contrasts
table_contrasts1.CInt<-as.data.frame(pairwise_contrasts.CInt)
table_contrasts2.CInt<-table_contrasts1.CInt[c(1,6)]
table_contrasts3.CInt<-pivot_wider(table_contrasts2.CInt, names_from = c(contrast),values_from = c(p.value))

t.CInt<-cbind(table4.CInt,table_contrasts3.CInt)
t.CInt$model<-"CInt"
t.CInt<-t.CInt[c(10,1:9)]
t.CInt
```

# LL Model
```{r}
model.LL<-lmer(LL~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.LL)
anova(model.LL)
emmeans(model.LL, pairwise~Treat)
plot(model.LL)
qqnorm(residuals(model.LL))
```

## output as table
```{r}
table1.LL<-as.data.frame(emmeans(model.LL, pairwise~Treat))
table2.LL<-head(table1.LL[c(1,3,4)],3)
table3.LL<-pivot_wider(table2.LL, names_from = Treat,values_from = c(emmean, SE))
table4.LL<-table3.LL[c(1,4,2,5,3,6)]

emmeans_output.LL<- emmeans(model.LL, pairwise ~ Treat)
pairwise_contrasts.LL<- summary(emmeans_output.LL)$contrasts
table_contrasts1.LL<-as.data.frame(pairwise_contrasts.LL)
table_contrasts2.LL<-table_contrasts1.LL[c(1,6)]
table_contrasts3.LL<-pivot_wider(table_contrasts2.LL, names_from = c(contrast),values_from = c(p.value))

t.LL<-cbind(table4.LL,table_contrasts3.LL)
t.LL$model<-"LL"
t.LL<-t.LL[c(10,1:9)]
t.LL
```

# DPL Model
```{r}
model.DPL<-lmer(DPL~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.DPL)
anova(model.DPL)
emmeans(model.DPL, pairwise~Treat)
plot(model.DPL)
qqnorm(residuals(model.DPL))
```

## output as table
```{r}
table1.DPL<-as.data.frame(emmeans(model.DPL, pairwise~Treat))
table2.DPL<-head(table1.DPL[c(1,3,4)],3)
table3.DPL<-pivot_wider(table2.DPL, names_from = Treat,values_from = c(emmean, SE))
table4.DPL<-table3.DPL[c(1,4,2,5,3,6)]

emmeans_output.DPL<- emmeans(model.DPL, pairwise ~ Treat)
pairwise_contrasts.DPL<- summary(emmeans_output.DPL)$contrasts
table_contrasts1.DPL<-as.data.frame(pairwise_contrasts.DPL)
table_contrasts2.DPL<-table_contrasts1.DPL[c(1,6)]
table_contrasts3.DPL<-pivot_wider(table_contrasts2.DPL, names_from = c(contrast),values_from = c(p.value))

t.DPL<-cbind(table4.DPL,table_contrasts3.DPL)
t.DPL$model<-"DPL"
t.DPL<-t.DPL[c(10,1:9)]
t.DPL
```


# DOY Model
```{r}
MY.305d$Dry.yield<-ifelse(MY.305d$Dry.yield==0,NA,MY.305d$Dry.yield)
model.DOY<-lmer(Dry.yield~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.DOY)
anova(model.DOY)
emmeans(model.DOY, pairwise~Treat)
plot(model.DOY)
qqnorm(residuals(model.DOY))
```

## output as table
```{r}
table1.DOY<-as.data.frame(emmeans(model.DOY, pairwise~Treat))
table2.DOY<-head(table1.DOY[c(1,3,4)],3)
table3.DOY<-pivot_wider(table2.DOY, names_from = Treat,values_from = c(emmean, SE))
table4.DOY<-table3.DOY[c(1,4,2,5,3,6)]

emmeans_output.DOY<- emmeans(model.DOY, pairwise ~ Treat)
pairwise_contrasts.DOY<- summary(emmeans_output.DOY)$contrasts
table_contrasts1.DOY<-as.data.frame(pairwise_contrasts.DOY)
table_contrasts2.DOY<-table_contrasts1.DOY[c(1,6)]
table_contrasts3.DOY<-pivot_wider(table_contrasts2.DOY, names_from = c(contrast),values_from = c(p.value))

t.DOY<-cbind(table4.DOY,table_contrasts3.DOY)
t.DOY$model<-"DOY"
t.DOY<-t.DOY[c(10,1:9)]
t.DOY
```

# 4.33 Model
```{r}
MY.305d$MY.4.33<-ifelse(MY.305d$MY.4.33==0,NA,MY.305d$MY.4.33)
model.4.33<-lmer(MY.4.33~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.4.33)
anova(model.4.33)
emmeans(model.4.33, pairwise~Treat)
plot(model.4.33)
qqnorm(residuals(model.4.33))
```

## output as table
```{r}
table1.4.33<-as.data.frame(emmeans(model.4.33, pairwise~Treat))
table2.4.33<-head(table1.4.33[c(1,3,4)],3)
table3.4.33<-pivot_wider(table2.4.33, names_from = Treat,values_from = c(emmean, SE))
table4.4.33<-table3.4.33[c(1,4,2,5,3,6)]

emmeans_output.4.33<- emmeans(model.4.33, pairwise ~ Treat)
pairwise_contrasts.4.33<- summary(emmeans_output.4.33)$contrasts
table_contrasts1.4.33<-as.data.frame(pairwise_contrasts.4.33)
table_contrasts2.4.33<-table_contrasts1.4.33[c(1,6)]
table_contrasts3.4.33<-pivot_wider(table_contrasts2.4.33, names_from = c(contrast),values_from = c(p.value))

t.4.33<-cbind(table4.4.33,table_contrasts3.4.33)
t.4.33$model<-"4.33"
t.4.33<-t.4.33[c(10,1:9)]
t.4.33
```

# PI Model
```{r}
MY.305d$PI.nov<-ifelse(MY.305d$PI.nov==0,NA,MY.305d$PI.nov)
model.PI<-lmer(PI.nov~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.PI)
anova(model.PI)
emmeans(model.PI, pairwise~Treat)
plot(model.PI)
qqnorm(residuals(model.PI))
```

## output as table
```{r}
table1.PI<-as.data.frame(emmeans(model.PI, pairwise~Treat))
table2.PI<-head(table1.PI[c(1,3,4)],3)
table3.PI<-pivot_wider(table2.PI, names_from = Treat,values_from = c(emmean, SE))
table4.PI<-table3.PI[c(1,4,2,5,3,6)]

emmeans_output.PI<- emmeans(model.PI, pairwise ~ Treat)
pairwise_contrasts.PI<- summary(emmeans_output.PI)$contrasts
table_contrasts1.PI<-as.data.frame(pairwise_contrasts.PI)
table_contrasts2.PI<-table_contrasts1.PI[c(1,6)]
table_contrasts3.PI<-pivot_wider(table_contrasts2.PI, names_from = c(contrast),values_from = c(p.value))

t.PI<-cbind(table4.PI,table_contrasts3.PI)
t.PI$model<-"PI"
t.PI<-t.PI[c(10,1:9)]
t.PI
```

# TM8.TM3 Model
```{r}
MY.305d$TM8.TM3<-ifelse(MY.305d$TM8.TM3==0,NA,MY.305d$TM8.TM3)
model.TM8.TM3<-lmer(TM8.TM3~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.TM8.TM3)
anova(model.TM8.TM3)
emmeans(model.TM8.TM3, pairwise~Treat)
plot(model.TM8.TM3)
qqnorm(residuals(model.TM8.TM3))
```

## output as table
```{r}
table1.TM8.TM3<-as.data.frame(emmeans(model.TM8.TM3, pairwise~Treat))
table2.TM8.TM3<-head(table1.TM8.TM3[c(1,3,4)],3)
table3.TM8.TM3<-pivot_wider(table2.TM8.TM3, names_from = Treat,values_from = c(emmean, SE))
table4.TM8.TM3<-table3.TM8.TM3[c(1,4,2,5,3,6)]

emmeans_output.TM8.TM3<- emmeans(model.TM8.TM3, pairwise ~ Treat)
pairwise_contrasts.TM8.TM3<- summary(emmeans_output.TM8.TM3)$contrasts
table_contrasts1.TM8.TM3<-as.data.frame(pairwise_contrasts.TM8.TM3)
table_contrasts2.TM8.TM3<-table_contrasts1.TM8.TM3[c(1,6)]
table_contrasts3.TM8.TM3<-pivot_wider(table_contrasts2.TM8.TM3, names_from = c(contrast),values_from = c(p.value))

t.TM8.TM3<-cbind(table4.TM8.TM3,table_contrasts3.TM8.TM3)
t.TM8.TM3$model<-"TM8.TM3"
t.TM8.TM3<-t.TM8.TM3[c(10,1:9)]
t.TM8.TM3
```

# DrY.TM3 Model
```{r}
model.DrY.TM3<-lmer(Pers.DO~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.DrY.TM3)
anova(model.DrY.TM3)
emmeans(model.DrY.TM3, pairwise~Treat)
plot(model.DrY.TM3)
qqnorm(residuals(model.DrY.TM3))
```

## output as table
```{r}
table1.DrY.TM3<-as.data.frame(emmeans(model.DrY.TM3, pairwise~Treat))
table2.DrY.TM3<-head(table1.DrY.TM3[c(1,3,4)],3)
table3.DrY.TM3<-pivot_wider(table2.DrY.TM3, names_from = Treat,values_from = c(emmean, SE))
table4.DrY.TM3<-table3.DrY.TM3[c(1,4,2,5,3,6)]

emmeans_output.DrY.TM3<- emmeans(model.DrY.TM3, pairwise ~ Treat)
pairwise_contrasts.DrY.TM3<- summary(emmeans_output.DrY.TM3)$contrasts
table_contrasts1.DrY.TM3<-as.data.frame(pairwise_contrasts.DrY.TM3)
table_contrasts2.DrY.TM3<-table_contrasts1.DrY.TM3[c(1,6)]
table_contrasts3.DrY.TM3<-pivot_wider(table_contrasts2.DrY.TM3, names_from = c(contrast),values_from = c(p.value))

t.DrY.TM3<-cbind(table4.DrY.TM3,table_contrasts3.DrY.TM3)
t.DrY.TM3$model<-"DrY.TM3"
t.DrY.TM3<-t.DrY.TM3[c(10,1:9)]
t.DrY.TM3
```

# Persistency TM8.TM3 Model
```{r}
model.TM8.3<-lmer(Pers.3.8~Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.TM8.3)
anova(model.TM8.3)
emmeans(model.TM8.3, pairwise~Treat)
plot(model.TM8.3)
qqnorm(residuals(model.TM8.3))
```

## output as table
```{r}
table1.TM8.3<-as.data.frame(emmeans(model.TM8.3, pairwise~Treat))
table2.TM8.3<-head(table1.TM8.3[c(1,3,4)],3)
table3.TM8.3<-pivot_wider(table2.TM8.3, names_from = Treat,values_from = c(emmean, SE))
table4.TM8.3<-table3.TM8.3[c(1,4,2,5,3,6)]

emmeans_output.TM8.3<- emmeans(model.TM8.3, pairwise ~ Treat)
pairwise_contrasts.TM8.3<- summary(emmeans_output.TM8.3)$contrasts
table_contrasts1.TM8.3<-as.data.frame(pairwise_contrasts.TM8.3)
table_contrasts2.TM8.3<-table_contrasts1.TM8.3[c(1,6)]
table_contrasts3.TM8.3<-pivot_wider(table_contrasts2.TM8.3, names_from = c(contrast),values_from = c(p.value))

t.TM8.3<-cbind(table4.TM8.3,table_contrasts3.TM8.3)
t.TM8.3$model<-"TM8.3"
t.TM8.3<-t.TM8.3[c(10,1:9)]
t.TM8.3
```

```{r}
z<-rbind(t.TM8.3,t.DrY.TM3)
write_xlsx(z, "Persistency.xlsx")
```

# PI Subgroup dataset
```{r}
MY.305d$PI.nov<-ifelse(MY.305d$PI.nov==0,NA,MY.305d$PI.nov)
PIs<-MY.305d[which(MY.305d$Trt.group=="PI"),] # subdataset with only the PI group
summary(PIs)
```

# PI Subgroup Models
```{r}
model.PIs<-lmer(Pers.3.8~Treat+(1|Farm)+Breed, data = PIs)
summary(model.PIs)
anova(model.PIs)
emmeans(model.PIs, pairwise~Treat)
plot(model.PIs)
qqnorm(residuals(model.PIs))
```

```{r}
model.PIs.DO<-lmer(Pers.DO~Treat+(1|Farm)+Breed, data = PIs)
summary(model.PIs.DO)
anova(model.PIs.DO)
emmeans(model.PIs.DO, pairwise~Treat)
plot(model.PIs.DO)
qqnorm(residuals(model.PIs.DO))
```
```{r}
model.PIs.MY.d<-lmer(MY.day.1~Treat+(1|Farm)+Breed, data = PIs)
summary(model.PIs.MY.d)
anova(model.PIs.MY.d)
emmeans(model.PIs.MY.d, pairwise~Treat)
plot(model.PIs.MY.d)
qqnorm(residuals(model.PIs.MY.d))
```
```{r}
model.PIs.ECMY.d<-lmer(ECMY.day.1~Treat+(1|Farm)+Breed, data = PIs)
summary(model.PIs.ECMY.d)
anova(model.PIs.ECMY.d)
emmeans(model.PIs.ECMY.d, pairwise~Treat)
plot(model.PIs.ECMY.d)
qqnorm(residuals(model.PIs.ECMY.d))
```
```{r}
model.LM<-lmer(Pers.3.8~PI.nov+Treat+(1|Farm)+Breed, data = MY.305d)
summary(model.LM)
```
```{r}
model.LM<-lm(Pers.3.8~PI.nov, data = MY.305d)
summary(model.LM)
```


## output as table
```{r}
table1.PIs<-as.data.frame(emmeans(model.PIs, pairwise~Treat))
table2.PIs<-head(table1.PIs[c(1,3,4)],3)
table3.PIs<-pivot_wider(table2.PIs, names_from = Treat,values_from = c(emmean, SE))
table4.PIs<-table3.PIs[c(1,4,2,5,3,6)]

emmeans_output.PIs<- emmeans(model.PIs, pairwise ~ Treat)
pairwise_contrasts.PIs<- summary(emmeans_output.PIs)$contrasts
table_contrasts1.PIs<-as.data.frame(pairwise_contrasts.PIs)
table_contrasts2.PIs<-table_contrasts1.PIs[c(1,6)]
table_contrasts3.PIs<-pivot_wider(table_contrasts2.PIs, names_from = c(contrast),values_from = c(p.value))

t.PIs<-cbind(table4.PIs,table_contrasts3.PIs)
t.PIs$model<-"PIs"
t.PIs<-t.PIs[c(8,1:4,7)]
t.PIs 
```

```{r}
write_xlsx(PIs, "PIs.xlsx")
```

## Dry period cathegory

```{r}
table(MY.305d$DPL.kat,MY.305d$Treat)
```

# DPL kategory binomial models

DPL.40.cI<-DPL.class.cI[which((DPL.class.cI$DP.40.1=="Yes") | (DPL.class.cI$DP.40.70.1=="Yes")),]
summary(DPL.40.cI)
table(DPL.40.cI$Breed, DPL.40.cI$DP.40.1)
table(DPL.40.cI$CI.plan, DPL.40.cI$DP.40.1)

```{r}
Data.40.cI<-MY.305d[which((MY.305d$DPL.40=="Yes") | (MY.305d$DPL.40.70=="Yes")),]
table(Data.40.cI$DPL.40, Data.40.cI$Treat)
```


```{r}
Model.DPL.40<-glmer(formula = DPL.40~Treat+(1|Farm)+Breed,family = binomial, data = Data.40.cI)
summary(Model.DPL.40)
anova(Model.DPL.40)
emmeans(Model.DPL.40, pairwise~Treat)
plot(Model.DPL.40)
qqnorm(residuals(Model.DPL.40))
car::Anova(Model.DPL.40)
```


```{r}
Data.70.cI<-MY.305d[which((MY.305d$DPL.70=="Yes") | (MY.305d$DPL.40.70=="Yes")),]
table(Data.70.cI$DPL.70, Data.70.cI$Treat)
```


```{r}
Model.DPL.70<-glmer(formula = DPL.70~Treat+(1|Farm)+Breed,family = binomial, data = Data.70.cI)
summary(Model.DPL.70)
anova(Model.DPL.70)
emmeans(Model.DPL.70, pairwise~Treat)
plot(Model.DPL.70)
qqnorm(residuals(Model.DPL.70))
car::Anova(Model.DPL.70)
```

```{r}
Model.DPL.40.70<-glmer(formula = DPL.40.70~Treat+(1|Farm)+Breed,family = binomial, data = MY.305d)
summary(Model.DPL.40.70)
anova(Model.DPL.40.70)
emmeans(Model.DPL.40.70, pairwise~Treat)
plot(Model.DPL.40.70)
qqnorm(residuals(Model.DPL.40.70))
car::Anova(Model.DPL.40.70)
```

# Results table
```{r}
x<-rbind(t.305d,t.LY,t.LY.ECM,t.MY.day,t.MY.d.ECM,t.CInt,t.LL,t.DPL,t.DOY,t.4.33,t.PI,t.TM8.TM3,t.DrY.TM3)
y<-round(x[c(2:7)],2)
q<-round(x[c(8:10)],3)
z<-cbind(x[c(1)],y,q)
z
```

```{r}
write_xlsx(z, "LM.part.2.xlsx")
```

# Basics - to appendix
```{r}
MY.305d.1<-Part.2[c(1:5,6,7,19,36:37,53,55,57:61,85,86)]
MY.305d.2<-MY.305d.1[which(MY.305d.1$VWP.ok.10=="Yes"),]
MY.305d<-MY.305d.2[which(MY.305d.2$Complete.1=="Yes"),]
summary(MY.305d.1) # ITT
summary(MY.305d.2) # VWP ok
summary(MY.305d) # VWP ok + complete lact.1
```

```{r}
table(MY.305d.1$Breed, MY.305d.1$Treat)
```
```{r}
table(MY.305d.1$VWP.ok.10, MY.305d.1$Treat)
```
```{r}
table(MY.305d.1$Complete.1, MY.305d.1$Treat)
```
```{r}
table(MY.305d.1$Not.ins, MY.305d.1$Treat)
```
```{r}
table(MY.305d.2$Not.ins, MY.305d.2$Treat)
```