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

library("readxl")
library("lubridate")
library("emmeans")
library("car")
library("glmmTMB")
library("tidyverse")
library("DHARMa")
library("writexl")
library("dplyr")
```

# prep.
```{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$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)
```

# Dataset grouped analysis
  1. CowID (n = 513)
  2. Farm (18 levels)
  3. CI.group (2 levels)
  5. VWP = ok +/- 10 d
  6. Breed.group (3 levels)
  72. n.AI (lakt 1)
  77. n.dis (lakt)
  85. time lact 1
  87. out lact 1, Yes = culled
```{r}
Grouped.1<-Part.2[c(1:3,5,6,7,23,19,36,37,72,77,85,87,89,124)]
Grouped.1b<-Grouped.1[which(Grouped.1$Farm!="B"),]
Grouped.2<-Grouped.1[which(Grouped.1$VWP.ok.10=="Yes"),]
Grouped.2b<-Grouped.2[which(Grouped.2$Farm!="B"),] # För - farm B har inga kor i ConvConv gruppen
Grouped.3<-Grouped.2b[which(Grouped.2b$Complete.1=="Yes"),]
summary(Grouped.2b)
```

```{r}
summary(Grouped.2)
```

```{r}
summary(Grouped.1b)
```

# ITT - trt. group / randomized to 12 or 16
```{r}
table(Part.2$Trt.group,Part.2$CI.plan)
```

# Grouping of the cows per farm and VWP group
```{r}
#aggregate(cull.z.VWP$Out.lact.1, list(cull.z.VWP$CI.plan), FUN=sum)
#aggregate(cull.z.VWP$Time.lact.1, list(cull.z.VWP$CI.plan), FUN=sum)
```


# Culling rate models - glmmTMB - random farm utan breed - Grouped by farm and CI group - negative binomial model
  1. ITT analysis: dataset Grouped.1
  2. VWP ok analysis: dataset Grouped.2
  
## Culling model  - ITT - "intention to treat" analysis
```{r}
cull.p2.ITT<-Grouped.1b%>%dplyr::group_by(Farm, Treat)%>%dplyr::summarize(sumout=sum(Out.b1),sumtime=sum(Time.l1))
mc.ITT<-glmmTMB(sumout~Treat+(1|Farm)+offset(log(sumtime)), family=nbinom2, data=cull.p2.ITT)
summary(mc.ITT)
car::Anova(mc.ITT)
emmeans(mc.ITT, specs=~Treat,type="response",  offset=log(365*100))
qqnorm(residuals(mc.ITT))
simulationOutput <- simulateResiduals(fittedModel =mc.ITT)
plot(simulationOutput)
emmeans(mc.ITT, pairwise~Treat, type="response",  offset=log(365*100))
```

## Culling model  - VWP ok - "per protocol" analysis
```{r}
cull.p2<-Grouped.2b%>%dplyr::group_by(Farm, Treat)%>%dplyr::summarize(sumout=sum(Out.b1),sumtime=sum(Time.l1))
mc.vwp<-glmmTMB(sumout~Treat+(1|Farm)+offset(log(sumtime)), family=nbinom2, data=cull.p2)
summary(mc.vwp)
car::Anova(mc.vwp)
emmeans(mc.vwp, specs=~Treat,type="response",  offset=log(365*100))
qqnorm(residuals(mc.vwp))
simulationOutput <- simulateResiduals(fittedModel =mc.vwp)
plot(simulationOutput)
emmeans(mc.vwp, pairwise~Treat, type="response",  offset=log(365*100))
```

## Disease model  - VWP ok - "per protocol" analysis - diseases from 33 DIM
```{r}
cull.p2<-Grouped.2b%>%dplyr::group_by(Farm, Treat)%>%dplyr::summarize(sumdis=sum(ndis.33),sumtime=sum(Time.l1))
md<-glmmTMB(sumdis~Treat+(1|Farm)+offset(log(sumtime)), family=nbinom2, data=cull.p2)
summary(md)
car::Anova(md)
emmeans(md, specs=~Treat,type="response",  offset=log(365*100))
qqnorm(residuals(md))
simulationOutput <- simulateResiduals(fittedModel =md)
plot(simulationOutput)
emmeans(md, pairwise~Treat, type="response",  offset=log(365*100))
```

# Disease model  - ITT - "per protocol" analysis - diseases from 33 DIM
```{r}
cull.p2<-Grouped.1b%>%dplyr::group_by(Farm, Treat)%>%dplyr::summarize(sumdis=sum(ndis.33),sumtime=sum(Time.l1))
md<-glmmTMB(sumdis~Treat+(1|Farm)+offset(log(sumtime)), family=nbinom2, data=cull.p2)
summary(md)
car::Anova(md)
emmeans(md, specs=~Treat,type="response",  offset=log(365*100))
qqnorm(residuals(md))
simulationOutput <- simulateResiduals(fittedModel =md)
plot(simulationOutput)
emmeans(md, pairwise~Treat, type="response",  offset=log(365*100))
```

# Number of disease events  from 33 DIM - VWP ok
```{r}
table(Grouped.2$ndis.33,Grouped.2$Treat)
```

# Number of disease events  from 33 DIM - ITT
```{r}
table(Part.2$ndis.33,Part.2$Treat)
```

## Disease model  - VWP ok - "per protocol" analysis - total number of diseases
```{r}
cull.p2<-Grouped.2b%>%dplyr::group_by(Farm, Treat)%>%dplyr::summarize(sumout=sum(n.dis),sumtime=sum(Time.l1))
md<-glmmTMB(sumout~Treat+(1|Farm)+offset(log(sumtime)), family=nbinom2, data=cull.p2)
summary(md)
car::Anova(md)
emmeans(md, specs=~Treat,type="response",  offset=log(365*100))
qqnorm(residuals(md))
simulationOutput <- simulateResiduals(fittedModel =md)
plot(simulationOutput)
emmeans(md, pairwise~Treat, type="response",  offset=log(365*100))
```

# NINS model  - VWP ok - "per protocol" analysis
```{r}
cull.preg <- Grouped.2b %>%
  dplyr::group_by(Farm, Treat) %>%
  dplyr::summarize(dplyr::across(c(n.AI, Preg.calv), sum, .names = "sum{.col}"))
mn<-glmmTMB(sumn.AI~Treat+(1|Farm)+offset(log(sumPreg.calv)), family=nbinom2, data=cull.preg)
summary(mn)
car::Anova(mn)
emmeans(mn, specs=~Treat,type="response",  offset=log(1))
qqnorm(residuals(mn))
simulationOutput <- simulateResiduals(fittedModel =mn)
plot(simulationOutput)
emmeans(mn, pairwise~Treat, type="response",  offset=log(1))
```