---
title: "Part 2 - Binomial models - Fertility"
author: "Anna Edvardsson Rasmussen"
date: "2023-03-07"
output: word_document
---

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

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

# Dataset
```{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$FSCR<-as.factor(Part.2$FSCR)
Part.2$Pregloss.5<-as.factor(Part.2$Pregloss.5)
Part.2$EI.1<-as.factor(Part.2$EI.1)
Part.2$Still.lact.1<-as.factor(Part.2$Still.lact.1)
Part.2<-Part.2[which(Part.2$Still.lact.1=="No"),] # Excluded cows still lact (n = 17)
Part.2<-Part.2[which(Part.2$Breed!="Cross"),] # Cross exkluderade (n = 75)
```

# table A1 - manus 4:
```{r}
table(Part.2$Breed, Part.2$Treat)
table(Part.2$Trt.group, Part.2$Treat)
table(Part.2$VWP.ok.10, Part.2$Treat)
table(Part.2$Not.ins, Part.2$Treat)
table(Part.2$Complete.1, Part.2$Treat)
table(Part.2$Complete.1, Part.2$VWP.ok.10, Part.2$Treat)
```


```{r}
summary(Part.2[c(3,5,6,7,36,37)])
```

# Fertility - Binomial models:
  1. Complience - ITT vs VWP ok +/- 10
  2. FSCR
  3. EI
  4. Pregloss
  5. SCC - last TM <  100,000 cells/mL

# Fertility - linear model dataset
  1. CowID (n = 513)
  2. Farm (18 levels)
  3. CI.group (2 levels)
  6. Breed.group (3 levels)
  19. Inclusion criteria: VWP = ok (n =432) 
  36. Group: 12/12, 16/12 or 16/16
  63. FSCR
  65. Pregloss
  71. EI
```{r}
Fert.LM.1<-Part.2[c(1:5,6,7,19,36,63,65,71,73,74)] # ITT
Fert.LM<-Fert.LM.1[which(Fert.LM.1$VWP.ok.10=="Yes"),]
summary(Fert.LM) # VWP ok +/- 10 d
```

# Complience - Binomial model

## Glmer binomial mixed model; Complience
Breed is significant => should we run an interaction? Trt*Breed?
```{r}
Model.compl<-glmer(formula = VWP.ok.10~Treat+(1|Farm)+Breed,family = binomial, data = Fert.LM.1)
summary(Model.compl)
anova(Model.compl)
emmeans(Model.compl, pairwise~Treat)
plot(Model.compl)
qqnorm(residuals(Model.compl))
car::Anova(Model.compl)
```
# FSCR - Binomial model


```{r}
table(Fert.LM$Treat,Fert.LM$FSCR)
```

## Glmer binomial mixed model: FSCR
  - inclusion criteria: VWP ok = only inseminated animals.
```{r}
Model.FSCR<-glmer(formula = FSCR~Treat+(1|Farm)+Breed,family = binomial, data = Fert.LM)
summary(Model.FSCR)
anova(Model.FSCR)
emmeans(Model.FSCR, pairwise~Treat)
plot(Model.FSCR)
qqnorm(residuals(Model.FSCR))
car::Anova(Model.FSCR)
```

# Pregloss - Binomial model

## Glmer binomial mixed model: Pregloss
  - inclusion criteria: VWP ok = only inseminated animals.
  - is a zero-inflation model more appropriate?
```{r}
Model.pregloss<-glmer(formula = Pregloss.5~Treat+(1|Farm)+Breed,family = binomial, data = Fert.LM)
summary(Model.pregloss)
anova(Model.pregloss)
emmeans(Model.pregloss, pairwise~Treat)
plot(Model.pregloss)
qqnorm(residuals(Model.pregloss))
car::Anova(Model.pregloss)
```

# EI dataset
```{r}
Fert.LM$EI.0.2<-ifelse(Fert.LM$EI.1== "0"|Fert.LM$EI.1== "1"|Fert.LM$EI.1== "2","1","0")
Fert.LM$EI.0.2<-as.factor(Fert.LM$EI.0.2)
Fert.LM$EI.3<-ifelse(Fert.LM$EI.1== "3","1","0")
Fert.LM$EI.3<-as.factor(Fert.LM$EI.3)
Fert.LM$EI.4.5<-ifelse(Fert.LM$EI.1== "4"|Fert.LM$EI.1== "5","1","0")
Fert.LM$EI.4.5<-as.factor(Fert.LM$EI.4.5)
Fert.LM<-Fert.LM[which(Fert.LM$Farm!="B"),] # Excluded 6 farms
Fert.LM<-Fert.LM[which(Fert.LM$Farm!="D"),] # Excluded 6 farms
Fert.LM<-Fert.LM[which(Fert.LM$Farm!="J"),] # Excluded 6 farms
Fert.LM<-Fert.LM[which(Fert.LM$Farm!="K"),] # Excluded 6 farms
Fert.LM<-Fert.LM[which(Fert.LM$Farm!="U"),] # Excluded 6 farms
Fert.LM<-Fert.LM[which(Fert.LM$Farm!="Y"),] # Excluded 6 farms
summary(Fert.LM)
```

```{r}
table(Fert.LM$Treat,Fert.LM$EI.0.2)
```

```{r}
table(Fert.LM$Treat,Fert.LM$EI.3)
```


```{r}
table(Fert.LM$Treat,Fert.LM$EI.4.5)
```

# EI - Binomial models

## Glmer binomial mixed model: EI 0-2
  - inclusion criteria: VWP ok = only inseminated animals.
  - is a zero-inflation model more appropriate?
```{r}
Model.EI.0.2<-glmer(formula = EI.0.2~Treat+(1|Farm)+Breed,family = binomial, data = Fert.LM)
summary(Model.EI.0.2)
anova(Model.EI.0.2)
emmeans(Model.EI.0.2, pairwise~Treat)
plot(Model.EI.0.2)
qqnorm(residuals(Model.EI.0.2))
car::Anova(Model.EI.0.2)
```

## Glmer binomial mixed model: EI 3
  - inclusion criteria: VWP ok = only inseminated animals.
  - is a zero-inflation model more appropriate?
```{r}
Model.EI.3<-glmer(formula = EI.3~Treat+(1|Farm)+Breed,family = binomial, data = Fert.LM)
summary(Model.EI.3)
anova(Model.EI.3)
emmeans(Model.EI.3, pairwise~Treat)
plot(Model.EI.3)
qqnorm(residuals(Model.EI.3))
car::Anova(Model.EI.3)
```

## Glmer binomial mixed model: EI 4.5
  - inclusion criteria: VWP ok = only inseminated animals.
  - is a zero-inflation model more appropriate?
```{r}
Model.EI.4.5<-glmer(formula = EI.4.5~Treat+(1|Farm)+Breed,family = binomial, data = Fert.LM)
summary(Model.EI.4.5)
anova(Model.EI.4.5)
emmeans(Model.EI.4.5, pairwise~Treat)
plot(Model.EI.4.5)
qqnorm(residuals(Model.EI.4.5))
car::Anova(Model.EI.4.5)
```

# SCC - Dataset
  1. CowID (n = 530)
  2. Farm (18 levels)
  3. Breed.group (2 levels)
  5. VWP.ok.10 Inclusion criteria: VWP = ok (n = 432)
  6.  and complete lactation (n = 338)
  36. Treatment group (3 levels)
  73. SCC at first TM
  74. SCC at last TM
  Ny1. Cows with SCC < 100,000 cells/mL at first TM ="No"
  Ny2. Cows with SCC < 100,000 cells/mL at last TM ="No"
```{r}
Fert.LM.2<-Fert.LM[which(Fert.LM$Complete.1=="Yes"),]
Fert.LM.2$Ny1<-ifelse(Fert.LM.2$SCC.first<100,"1","0")
Fert.LM.2$Ny2<-ifelse(Fert.LM.2$SCC.last<100,"1","0")
Fert.LM.2$Ny1<-as.factor(Fert.LM.2$Ny1)
Fert.LM.2$Ny2<-as.factor(Fert.LM.2$Ny2)
summary(Fert.LM.2)
```
## Glmer binomial mixed model: SCC first TM
  - inclusion criteria: VWP ok + complete lactation
```{r}
Model.Ny1<-glmer(formula = Ny1~Treat+(1|Farm)+Breed,family = binomial, data = Fert.LM.2)
summary(Model.Ny1)
anova(Model.Ny1)
emmeans(Model.Ny1, pairwise~Treat)
plot(Model.Ny1)
qqnorm(residuals(Model.Ny1))
car::Anova(Model.Ny1)
```


## Glmer binomial mixed model: SCC last TM
  - inclusion criteria: VWP ok + complete lactation
```{r}
Model.Ny2<-glmer(formula = Ny2~Treat+(1|Farm)+Breed,family = binomial, data = Fert.LM.2)
summary(Model.Ny2)
anova(Model.Ny2)
emmeans(Model.Ny2, pairwise~Treat)
plot(Model.Ny2)
qqnorm(residuals(Model.Ny2))
car::Anova(Model.Ny2)
```

```{r}
table(Fert.LM.2$Treat,Fert.LM.2$Ny1)
```

```{r}
table(Fert.LM.2$Treat,Fert.LM.2$Ny2)
```