---
title: "Part 2 - linear models - Fertility"
author: "Anna Edvardsson Rasmussen"
date: "2023-03-07"
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")
library("survival")
library("coxme")
library("survminer")
```

# 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$DPL.kat<-as.factor(Part.2$DPL.kat)
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)
```

# Fertility - linear models:
  1. CFI
  2. IPL

# Fertility - linear model dataset
  1. CowID (n = 530)
  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
  66. CFI
  69. IPL
```{r}
Fert.LM<-Part.2[c(1:5,6,7,19,36,66,69,88,97,125)]
Fert.LM<-Fert.LM[which(Fert.LM$VWP.ok.10=="Yes"),]
summary(Fert.LM) # VWP ok +/- 10 d
```

# CFI - linear model
```{r}
model.CFI<-lmer(CFI~Treat+(1|Farm)+Breed, data = Fert.LM)
summary(model.CFI)
anova(model.CFI)
emmeans(model.CFI, pairwise~Treat)
plot(model.CFI)
qqnorm(residuals(model.CFI))
```

# IPL - linear model
```{r}
model.IPL<-lmer(IPL~Treat+(1|Farm)+Breed, data = Fert.LM)
summary(model.IPL)
anova(model.IPL)
emmeans(model.IPL, pairwise~Treat)
plot(model.IPL)
qqnorm(residuals(model.IPL))
```

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

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

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

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

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

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


```{r}
x<-rbind(t.CFI,t.IPL)
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.fert.xlsx")
```

# Kaplan Meier - IPL - Part 2
```{r}
model_fit <- survfit(Surv(IPL.cens, Preg) ~ Treat, data = Fert.LM)

ggsurvplot(
  model_fit,
  xlim = c(0, 80),
  break.x.by = 20,
  xlab = "Insemination period length (days)",
  ylab = "Non-pregnant cows (%)\n",
  title = "IPL, Part 2\n",
  legend = "right",
  data = Fert.LM,
  font.x = 20,
  font.y = 20,
  font.tickslab = 18,
  font.title = 24,
  vjust = 1,
  font.legend = 18,
  theme = theme(plot.width = 3, plot.height = 10, text = element_text(size = 18)),
  surv.scale = "percent",
  palette = c("lightgrey","darkgrey","black"),
  legend.title = "VWP group",
  legend.labs = c("ConvConv","ExtConv","ExtExt")  # Add a placeholder label
)
```
# culling causes - Part 2 - ITT
```{r}
cull.2<-Part.2[c(1:3,5,6,7,36,86,126,128)]
cull.2$cull.cause.1<-as.factor(cull.2$cull.cause.1)
summary(cull.2)
```
# cull cause ITT - part 2
```{r}
table(cull.2$Treat,cull.2$cull.cause.1)
```
# culling causes - Part 2 - "per protocol"
```{r}
cull.2<-Part.2[c(1:3,5,6,7,36,86,126,128)]
cull.2$cull.cause.1<-as.factor(cull.2$cull.cause.1)
cull.2<-cull.2[which(cull.2$VWP.ok.10 == "Yes"),]
summary(cull.2)
```

# cull cause "per protocol" - part 2
```{r}
table(cull.2$Treat,cull.2$cull.cause.1)
```