---
title: "Manus 1 - No Cross Binomial milk lact 1"
author: "Anna Edvardsson Rasmussen"
date: '2022-05-09'
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")
```

# Cross = excluded.
```{r echo=TRUE, include=FALSE}
setwd(".")
Lact.1.full<- read_excel("Lact.1.R.ny.xlsx")
Lact.1.full$CI.plan<-as.factor(Lact.1.full$CI.plan)
Lact.1.full$Farm<-as.factor(Lact.1.full$Farm)
Lact.1.full$Breed<-as.factor(Lact.1.full$Breed)
Lact.1.full$DP.cat<-as.factor(Lact.1.full$DP.cat.1)
Lact.1.full$DP.40.1<-as.factor(Lact.1.full$DP.40.1)
Lact.1.full$DP.40.70.1<-as.factor(Lact.1.full$DP.40.70.1)
Lact.1.full$DP.70.1<-as.factor(Lact.1.full$DP.70.1)
Lact.1.full$DO.MY.less.12.5.kg.1<-as.factor(Lact.1.full$DO.MY.less.12.5.kg.1)
Lact.1.full$VWP.complete.1.ok<-as.factor(Lact.1.full$VWP.complete.1.ok)
Lact.1.full$VWP.compl.MY.OK.1<-as.factor(Lact.1.full$VWP.compl.MY.OK.1)
Lact.1.full$VWP.notOK.MY.OK.1<-as.factor(Lact.1.full$VWP.notOK.MY.OK.1)
Lact.1.full$VWP.compl.MY.DO.OK.1<-as.factor(Lact.1.full$VWP.compl.MY.DO.OK.1)
```

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

# Dry period length dataset
  1. CowID (n = 320)
  2. Farm (16 levels)
  3. VWP group (2 levels)
  6. Breed group (3 levels)
  21. Inclusion criteria: VWP.compl.MY.OK.1 = VWP acc to plan, complete lactation + daily MY = ok
  30. Dry period of less than 40 days (2 levels: Yes/No)
  31. Dry period between 40-70 days (2 levels: Yes/No)
  32. Dry period of more than 70 days (2 levels: Yes/No)

# 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). 
```{r echo=TRUE, include=FALSE}
DPL.class.v<-c(1:3,6,21,30:32)
DPL.class<-Lact.1.full[DPL.class.v]
DPL.class.c<-DPL.class[complete.cases(DPL.class),]
DPL.class.cI<-DPL.class.c[which(DPL.class.c$VWP.compl.MY.OK.1=="Yes"),]
```

# Summary of dataset
```{r}
summary(DPL.class.cI)
```


```{r}
DPL.40.cI<-DPL.class.cI[which((DPL.class.cI$DP.40.1=="Yes") | (DPL.class.cI$DP.40.70.1=="Yes")),]
table(DPL.40.cI$Farm, DPL.40.cI$DP.40.1)
table(DPL.class.cI$Farm, DPL.class.cI$DP.40.70.1)
DPL.70.cI<-DPL.class.cI[which((DPL.class.cI$DP.70.1=="Yes") | (DPL.class.cI$DP.40.70.1=="Yes")),]
table(DPL.70.cI$Farm, DPL.70.cI$DP.70.1)
```
# DP-length: 
Binary variables:
  Dry period < 40 d
  Dry period 40-70 d
  Dry period > 70 d
  
# Model description:
  Breed: fix - 3 levels  
  Farm: random - 16 levels  
  Group: fix - 2 levels  
  Breed:group interaction: not included due to perfect separation problem (too small subgroups)  
  
# Short dry period dataset
Dataset for only short and "normal" dry period. (The tables shows the problem with perfect separation)
```{r}
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)
```

# DPL 40 logistic mixed model: Test without interaction is
```{r}
model.DPL.40.1<-glm(formula = DP.40.1~CI.plan+Breed+Farm+CI.plan*Breed,
                    family = binomial, data = DPL.40.cI)
summary(model.DPL.40.1)
anova(model.DPL.40.1)
emmeans(model.DPL.40.1, pairwise~CI.plan)
emmeans(model.DPL.40.1, pairwise~Breed)
plot(model.DPL.40.1)
qqnorm(residuals(model.DPL.40.1))
car::Anova(model.DPL.40.1)
```
Still singularity problem with equalizer.. Not enough variation between farms (to few observations?)
```{r}
model.DPL.40.1<-glmer(formula = DP.40.1~CI.plan+Breed+(1|Farm)+ CI.plan*Breed,
                      family = binomial, control = glmerControl(optimizer="bobyqa"),  data = DPL.40.cI)
car::Anova(model.DPL.40.1)
```

# Final DP 40 model
  This one is the final model despite Singular warning - the results are not much different - not worth it to change the model
```{r}
model.DPL.40.1.2<-glmer(formula = DP.40.1~CI.plan+Breed+(1|Farm),
                        family = binomial, data = DPL.40.cI)
summary(model.DPL.40.1.2)
anova(model.DPL.40.1.2)
car::Anova(model.DPL.40.1.2)
```


# "Normal" dry period data set
No perfect separation - but want to do all models the same - excluding interaction
```{r}
summary(DPL.class.cI)
table(DPL.class.cI$Breed, DPL.class.cI$DP.40.70.1)
table(DPL.class.cI$CI.plan, DPL.class.cI$DP.40.70.1)
```

# DPL 40 - 70 logistic mixed model
```{r}
model.DPL.40.70.1<-glmer(formula = DP.40.70.1~(1|Farm)+CI.plan+Breed+CI.plan*Breed,
                         family = binomial, data = DPL.class.cI)
summary(model.DPL.40.70.1)
anova(model.DPL.40.70.1)
car::Anova(model.DPL.40.70.1)
```
No significant interaction => model 2

```{r}
model.DPL.40.70.1.2<-glmer(formula = DP.40.70.1~CI.plan+Breed+(1|Farm),
                           family = binomial, data = DPL.class.cI)
summary(model.DPL.40.70.1.2)
anova(model.DPL.40.70.1.2)
car::Anova(model.DPL.40.70.1.2)
```

# Long dry period
```{r}
summary(DPL.70.cI)
table(DPL.70.cI$Breed, DPL.70.cI$DP.70.1)
table(DPL.70.cI$CI.plan, DPL.70.cI$DP.70.1)
```

# DPL 70 logistic mixed model - med felmeddelande pga overdispersion
```{r}
model.DPL.70.1<-glmer(formula = DP.70.1~(1|Farm)+CI.plan+Breed+CI.plan*Breed,
                      family = binomial, data = DPL.70.cI)
summary(model.DPL.70.1)
anova(model.DPL.70.1)
car::Anova(model.DPL.70.1)
```

No significant interaction => try model 2
```{r}
model.DPL.70.1.2<-glmer(formula = DP.70.1~(1|Farm)+CI.plan+Breed,
                        family = binomial, data = DPL.70.cI)
summary(model.DPL.70.1.2)
anova(model.DPL.70.1.2)
car::Anova(model.DPL.70.1.2)
```

# MY < 12,5 kg dataset
  1. CowID (n = 435)
  2. Farm (16 levels)
  3. VWP group (2 levels)
  6. Breed group (3 levels)
  22. Inclusion criteria: with only VWPok+only complete lact + Daily MY OK + DO date OK included
  34. DO.MY.less.12.5.kg.1: Dry off MY less than 12,5 kg: (2 levels (Yes/No)) 
  
## Subset with only complete cases: Binomial parameter dry off yield less than 12,5 kg (Yes/No)
```{r echo=TRUE, include=FALSE}
DOY.less.12.5.v<-c(1:3,6,22,34)
DOY.less.12.5<-Lact.1.full[DOY.less.12.5.v]
DOY.less.12.5.c<-DOY.less.12.5[complete.cases(DOY.less.12.5),]
DOY.less.12.5.cI<-DOY.less.12.5.c[which(DOY.less.12.5.c$VWP.compl.MY.DO.OK.1=="Yes"),]
```

```{r}
summary(DOY.less.12.5.cI)
table(DOY.less.12.5.cI$Breed, DOY.less.12.5.cI$DO.MY.less.12.5.kg.1)
table(DOY.less.12.5.cI$CI.plan, DOY.less.12.5.cI$DO.MY.less.12.5.kg.1)
```


# Problem with perfect separation again => interaction is excluded from the model
```{r}
summary(DOY.less.12.5.cI)
table(DOY.less.12.5.cI$Breed, DOY.less.12.5.cI$DO.MY.less.12.5.kg.1)
table(DOY.less.12.5.cI$CI.plan, DOY.less.12.5.cI$DO.MY.less.12.5.kg.1)
table(DOY.less.12.5.cI$CI.plan, DOY.less.12.5.cI$Breed, DOY.less.12.5.cI$DO.MY.less.12.5.kg.1)
table(DOY.less.12.5.cI$Farm, DOY.less.12.5.cI$DO.MY.less.12.5.kg.1)
```

# 12,5 kg generalized logistic mixed effect binomial model
```{r}
model.DOy.12.5<-glmer(formula = DO.MY.less.12.5.kg.1~(1|Farm)+CI.plan+Breed,
                      family = binomial, data = DOY.less.12.5.cI)
summary(model.DOy.12.5)
car::Anova(model.DOy.12.5)
```