---
title: "Analysis of serial passaging experiment"
author: "Lievens et al."
date: "2019"
output: html_document
---

# PREP

Libraries:

```{r}
library(Hmisc)
library(reshape)
library(lme4)
library(car)
library(survival)
library(MuMIn)
library(RVAideMemoire)
library(MASS)
library(flexsurv)
library(multcomp)
library(lsmeans)
```

Make a dataframe that has the information on each line:

-	`Line`: evolved line ID (*1*-*24*: 2 microsporidians x 3 treatments x 4 replicates)
-	`MS`: microsporidian species (*MS1* for A. rigaudi, *MS2* for E. artemiae)
-	`EvolvedHost`: evolutionary treatment (*AF* for 'A. f. host', *AP* for 'A. p. host', *alt* for 'Alternating hosts'). Also includes a *Ctrl* level to facilitate merging with the assay datasets, but this does not represent an actual evolved line.
-	`Replicate`: replicate number of this microsporidian x treatment combination (*1*-*4*)
-	`Fridge`: *1* if the line was revived from backup spores collected after P6

*[Code not shown]*
```{r echo=FALSE}
treatments<-data.frame(Line=c("Line 1","Line 2","Line 3","Line 4","Line 5","Line 6","Line 7","Line 8","Line 9","Line 10","Line 11","Line 12","Line 13","Line 14","Line 16","Line 17","Line 18 fridge","Line 19","Line 20","Line 22","Line 23","Line 1 fridge","Line 2 fridge","Line 3 fridge","Line 4 fridge","Line 12 fridge","Ctrls 1", "Ctrls 2"),
                       MS=c("MS2","MS2","MS2","MS2","MS1","MS1","MS1","MS1","MS1","MS1","MS2","MS2","MS2","MS2","MS2","MS1","MS1","MS2","MS2","MS1","MS1","MS2","MS2","MS2","MS2","MS2","Ctrl","Ctrl"),
                       EvolvedHost=c("AF","AF","AF","AF","alt","alt","AP","AP","AP","AP","alt","alt","AP","AP","AP","AF","AF","alt","alt","AF","alt","AF","AF","AF","AF","alt","Ctrl","Ctrl"),
                       Replicate=as.factor(c(1,2,3,4,1,2,1,2,3,4,1,2,1,2,4,1,2,3,4,4,3,1,2,3,4,2,1,1)))
treatments$Fridge<-0
treatments$Fridge[grepl("fridge",treatments$Line)==T]<-1
treatments
```



# SERIAL PASSAGES

## Host survival

### Data

Variables:

-	`Line`: Line ID
-	`MS`: microsporidian species (*MS1* for A. rigaudi, *MS2* for E. artemiae)
-	`EvolvedHost`: evolutionary treatment (*AF* for 'A. f. host', *AP* for 'A. p. host', *alt* for 'Alternating hosts')
-	`Replicate`: replicate number of this microsporidian x treatment combination (*1*-*4*)
-	`Passage`: passage number after which the survival data was collected
-	`PassageHost`: host on which that passage was carried out
-	`alive`: number of surviving hosts in that passage
-	`total`: number of hosts that were originally seeded in that passage (normally 40, occasionally 45 due to counting errors)

New variables:

-	`dead`: `total` - `alive`
-	`lostline`: *1* if the parasite line was lost by this passage time (see 'Parasite population size')

*[Code not shown]*
```{r echo=F}
# Prep

surviv <- read.csv("./DATA survival during passages.csv",sep=";")

surviv$dead<-surviv$total-surviv$alive
surviv<-surviv[complete.cases(surviv),]

surviv$lostline<-0
surviv[surviv$MS=="MS2"&surviv$EvolvedHost=="AP"&surviv$Replicate==3&surviv$Passage>6,]$lostline<-1
surviv[surviv$MS=="MS1"&surviv$EvolvedHost=="AF"&surviv$Replicate==2&surviv$Passage>6,]$lostline<-1
surviv[surviv$MS=="MS1"&surviv$EvolvedHost=="AF"&surviv$Replicate==3&surviv$Passage>4,]$lostline<-1
surviv[surviv$MS=="MS1"&surviv$EvolvedHost=="alt"&surviv$Replicate==4&surviv$Passage>1,]$lostline<-1

head(surviv)
```

### Analysis

__Test: did host survival change over time?__

Model: Binomial GLMM with survival as a response variable. Because we did not maintain “control” host populations during the serial passage experiment, host survival is relative (e.g. survival in ‘Alternating hosts’ vs. ‘A. f. host’ treatments), and can only be compared within host species; we therefore analyzed it separately for A. franciscana and A. parthenogenetica.

-	Fixed effects: `Passage` (continuous), `MS`, `EvolvedHost`, interactions
-	Random effects: `Line` (to control for pseudoreplication), `Passage` (to deal with the varying background host mortality)
-	Exclude lost lines

Statistics: likelihood ratio tests on the fixed effects

*[Code not shown]*
```{r include=F,eval=F}
surviv$Passagef<-as.factor(surviv$Passage)


### AF

#data
temp<-surviv[surviv$PassageHost=="franciscana"&surviv$lostline==0,]

#models
a<-glmer(cbind(alive,dead)~EvolvedHost+MS+Passage +Passage:MS+EvolvedHost:Passage+EvolvedHost:MS+ EvolvedHost:MS:Passage+ (1|Line)+(1|Passagef),family=binomial,data=temp)
b<-glmer(cbind(alive,dead)~EvolvedHost+MS+Passage+ Passage:MS+EvolvedHost:Passage+EvolvedHost:MS+ (1|Line)+(1|Passagef),family=binomial,data=temp)
anova(a,b) #a is significantly better
summary(a)

#test sign of the rest with type III (Wald chisq)
options(contrasts=c(unordered="contr.sum",ordered="contr.poly"))
Anova(a,type="III")
options(contrasts=c(unordered="contr.treatment",ordered="contr.poly")) #reset to default

#posthoc testing
temp$dummy<-paste(temp$MS, temp$EvolvedHost)
ph<-glmer(cbind(alive,dead)~Passage +dummy+dummy:Passage+(1|Line)+(1|Passagef),family=binomial,data=temp) #doesn't converge
lsmeans(ph,pairwise~dummy:Passage,adjust="tukey")$contrasts



### AP

#data
temp<-surviv[surviv$PassageHost=="parthenogenetica"&surviv$lostline==0,]

#models
a<-glmer(cbind(alive,dead)~EvolvedHost+MS+Passage+Passage:MS+EvolvedHost:Passage+EvolvedHost:MS+EvolvedHost:MS:Passage+ (1|Line)+(1|Passagef),family=binomial,data=temp)
b<-glmer(cbind(alive,dead)~EvolvedHost+MS+Passage+Passage:MS+EvolvedHost:Passage+EvolvedHost:MS+ (1|Line)+(1|Passagef),family=binomial,data=temp)
anova(a,b) #p of +EvolvedHost:MS:Passage =0.5353, chisq(2)=0.3843
c<-glmer(cbind(alive,dead)~EvolvedHost+MS+Passage+Passage:MS+EvolvedHost:Passage+ (1|Line)+(1|Passagef),family=binomial,data=temp)
d<-glmer(cbind(alive,dead)~EvolvedHost+MS+Passage+Passage:MS+EvolvedHost:MS+ (1|Line)+(1|Passagef),family=binomial,data=temp)
e<-glmer(cbind(alive,dead)~EvolvedHost+MS+Passage+EvolvedHost:Passage+EvolvedHost:MS+ (1|Line)+(1|Passagef),family=binomial,data=temp)
anova(b,c) #p of +EvolvedHost:MS=0.2326, chisq(1) =1.4252
anova(b,d) #p of +EvolvedHost:Passage=0.003302 **, chisq(1) =8.6326
anova(b,e) #p of +Passage:MS=0.0007406 ***, chisq(1) =11.384
summary(c)

#test sign of the rest with type III (Wald chisq)
options(contrasts=c(unordered="contr.sum",ordered="contr.poly"))
Anova(c,type="III")
options(contrasts=c(unordered="contr.treatment",ordered="contr.poly")) #reset to default
```

=> AF: significant triple interaction. As the passages progressed, the survival of A. franciscana in A. rigaudi × ’A. f. host’ lines decreased as compared to the others

=> AP: significant effects of `MS` and `EvolvedHost` in interaction with `Passage`. Survival rates became progressively higher in A. rigaudi relative to E. artemiae lines, as well as in ‘Alternating hosts’ relative to ‘A. p. host’ lines.

=> Remember that all of these changes are relative to each other, not to unexposed hosts!

Model checking:

*[Code not shown]*
```{r include=F,eval=F}
### AF

temp<-surviv[surviv$PassageHost=="franciscana"&surviv$lostline==0,]
a<-glmer(cbind(alive,dead)~EvolvedHost+MS+Passage +Passage:MS+EvolvedHost:Passage+EvolvedHost:MS+ EvolvedHost:MS:Passage+ (1|Line)+(1|Passagef),family=binomial,data=temp)


## model fit

#overdispersion (for GLMs and GLMMs): overdispersed (3) not much to be done about it since there are already lots of random effects

overdisp.glmer(a)

#fitted vs. observed: fine, a bit skewed at the extremes

f <- fitted(a)
plot(temp$alive/temp$total, f)
abline(0,1)


## residuals

r <- residuals(a)

#normal: fine

hist(r)
qqnorm(r)
qqline(r)

#constant variance: fine

plot(r)
abline(0,0)
mean(r, na.rm=TRUE)

#independent of size of fitted values: fine

plot(f,r)
sm <- loess(r~f)
v <- seq(min(f,na.rm=T),max(f,na.rm=T),length=101)
lines(v,predict(sm,data.frame(f=v)),col=2)

#independent of independents: fine

plot(temp$MS,r)
plot(temp$Passage,r)
plot(temp$EvolvedHost,r)

#lag plot: fine

plot(r,r,type="n")
for (i in 1:length(r)){
  if (i==1) points(r[i]~r[length(r)]) else points(r[i]~r[i-1])
}
abline(0,0)



### AP

temp<-surviv[surviv$PassageHost=="parthenogenetica"&surviv$lostline==0,]
c<-glmer(cbind(alive,dead)~EvolvedHost+MS+Passage+Passage:MS+EvolvedHost:Passage+ (1|Line)+(1|Passagef),family=binomial,data=temp)


## model fit

#overdispersion (for GLMs and GLMMs): 1.6, ok

overdisp.glmer(c)

#fitted vs. observed: fine, a bit skewed at the extremes

f <- fitted(c)
plot(temp$alive/temp$total, f)
abline(0,1)


## residuals

r <- residuals(c)

#normal: fine

hist(r)
qqnorm(r)
qqline(r)

#constant variance: fine

plot(r)
abline(0,0)
mean(r, na.rm=TRUE)

#independent of size of fitted values: fine

plot(f,r)
sm <- loess(r~f)
v <- seq(min(f,na.rm=T),max(f,na.rm=T),length=101)
lines(v,predict(sm,data.frame(f=v)),col=2)

#independent of independents: fine

plot(temp$MS,r)
plot(temp$Passage,r)
plot(temp$EvolvedHost,r)

#lag plot: fine

plot(r,r,type="n")
for (i in 1:length(r)){
  if (i==1) points(r[i]~r[length(r)]) else points(r[i]~r[i-1])
}
abline(0,0)
```

=> AF model is quite overdispersed, but all other checks are fine.

=> AP model is fine.


## Parasite population size

### Data

Variables:

-	`Line`: Line ID
-	`MS`: microsporidian species (*MS1* for A. rigaudi, *MS2* for E. artemiae)
-	`EvolvedHost`: evolutionary treatment (*AF* for 'A. f. host', *AP* for 'A. p. host', *alt* for 'Alternating hosts')
-	`Replicate`: replicate number of this microsporidian x treatment combination (*1*-*4*)
-	`Passage`: passage number after which the survival data was collected
-	`PassageHost`: host on which that passage was carried out
-	`donors`: the number of hosts that were ground up (as a group) to count their spore content. Roughly half of all the hosts.
-	`total.spores`: the number of spores in the ground up donors

New variables:

-	`alive`: number of surviving hosts in that passage, extracted from the 'Host survival' dataset (-> PREPARE THAT DATASET FIRST!)
-	`all.spores`: (`total.spores` / `donors`) * number of hosts alive at the end of the passage
-	`all.spores.NA`: `all.spores` with 0 values replaced by NA

*[Code not shown]*
```{r echo=F}
spores <- read.csv("./DATA msp population size checks.csv",sep=";",dec=".")

temp<-surviv[surviv$Passage %in% c(1,4,7),]
temp$dummy<-as.factor(paste(temp$Line,temp$Passage))
spores$dummy<-as.factor(paste(spores$Line,spores$Passage))
spores<-merge(spores,temp[,c("dummy","alive")],all.x=TRUE)
spores$all.spores<-(spores$total.spores/spores$donors)*spores$alive

spores$all.spores.NA<-NA
spores$all.spores.NA[spores$all.spores>0]<-spores[spores$all.spores>0,]$all.spores
spores$all.spores.NA<-as.numeric(spores$all.spores.NA)

head(spores)
```

### Analysis

__Test: did parasite population size change over time?__

Model: LMM with ln+1-transformed spore counts as a response variable. A. rigaudi and E. artemiae lines were analyzed separately.

-	Fixed effects: `Passage` (continuous), `EvolvedHost`, interaction
-	Random effects: `Line` (to control for pseudoreplication)
-	Exclude lost lines

Statistics: likelihood ratio tests on the fixed effects

*[Code not shown]*
```{r include=F,eval=F}
### MS1

#data
temp<-spores[spores$MS=="MS1",]

#models
full<-lmer(log(all.spores.NA+1)~EvolvedHost*Passage+(1|Line),data=temp)
a<-lmer(log(all.spores.NA+1)~EvolvedHost+Passage+(1|Line),data=temp)
b<-lmer(log(all.spores.NA+1)~EvolvedHost+(1|Line),data=temp)
c<-lmer(log(all.spores.NA+1)~Passage+(1|Line),data=temp)
anova(full,a) #a better, p=0.326 = p of interaction, chisq(2)=2.242
anova(a,b) #a better, p=0.001605 = p of passage, chisq(1)=9.9535
anova(a,c) #a better, p=0.0272 = p of EvolvedHost, chisq(2)=7.2093

#posthoc testing
summary(glht(a, mcp(EvolvedHost="Tukey"))) #Ap sign different from Af, p = 0.02



### MS2

#data
temp<-spores[spores$MS=="MS2",]

#models
full<-lmer(log(all.spores.NA+1)~EvolvedHost*Passage+(1|Line),data=temp)
a<-lmer(log(all.spores.NA+1)~EvolvedHost+Passage+(1|Line),data=temp)
b<-lmer(log(all.spores.NA+1)~EvolvedHost+(1|Line),data=temp)
c<-lmer(log(all.spores.NA+1)~Passage+(1|Line),data=temp)
anova(full,a) #a better, p=0.1413 = p of interaction, chisq(2)=3.9132
anova(a,b) #b better, p=0.9481 = p of passage, chisq(1)=0.0042
anova(a,c) #a better, p=0.005451 = p of EvolvedHost, chisq(2)=10.424

#posthoc testing
summary(glht(b, mcp(EvolvedHost="Tukey"))) #Ap sign different from Af, p = 0.003
```

=> MS1: significant effects of `Passage` and `EvolvedHost`. Population grew over the course of the experiment, and was larger for lines evolving on AP. 

=> MS2: significant effect of `EvolvedHost`. Populations were higher in lines evolving on AF. 

Model checking:

*[Code not shown]*
```{r include=F,eval=F}
### MS1

temp<-spores[spores$MS=="MS1",]
a<-lmer(log(all.spores.NA+1)~EvolvedHost+Passage+(1|Line),data=temp)


## model fit

#overdispersion (for GLMs and GLMMs): fine

overdisp.glmer(a)

#fitted vs. observed: fine, a bit skewed at the extremes

f <- fitted(a)
plot(log(temp[!is.na(temp$all.spores.NA),]$all.spores.NA+1), f)
abline(0,1)


## residuals

r <- residuals(a)

#normal: skewed at the low extreme

hist(r)
qqnorm(r)
qqline(r)

#constant variance: fine

plot(r)
abline(0,0)
mean(r, na.rm=TRUE)

#independent of size of fitted values: fine

plot(f,r)
sm <- loess(r~f)
v <- seq(min(f,na.rm=T),max(f,na.rm=T),length=101)
lines(v,predict(sm,data.frame(f=v)),col=2)

#independent of independents: fine

plot(temp[!is.na(temp$all.spores.NA),]$Passage,r)
plot(temp[!is.na(temp$all.spores.NA),]$EvolvedHost,r)



### MS2

temp<-spores[spores$MS=="MS2",]
b<-lmer(log(all.spores.NA+1)~EvolvedHost+(1|Line),data=temp)


## model fit

#overdispersion (for GLMs and GLMMs): overdispersed

overdisp.glmer(b)

#fitted vs. observed: fine, a bit skewed at the extremes

f <- fitted(b)
plot(log(temp[!is.na(temp$all.spores.NA),]$all.spores.NA+1), f)
abline(0,1)


## residuals

r <- residuals(b)

#normal: skewed at the low extreme

hist(r)
qqnorm(r)
qqline(r)

#constant variance: fine

plot(r)
abline(0,0)
mean(r, na.rm=TRUE)

#independent of size of fitted values: fine

plot(f,r)
sm <- loess(r~f)
v <- seq(min(f,na.rm=T),max(f,na.rm=T),length=101)
lines(v,predict(sm,data.frame(f=v)),col=2)

#independent of independents: fine

plot(temp[!is.na(temp$all.spores.NA),]$EvolvedHost,r)
```

=> MS1 model is fine.

=> MS2 model is somewhat overdispersed, all other checks are fine.



# FINAL ASSAYS

## 1st (infectivity)

### Data

Variables:

-	`Line`: Line ID
-	`MS`: microsporidian species (*MS1* for A. rigaudi, *MS2* for E. artemiae)
-	`EvolvedHost`: evolutionary treatment (*AF* for 'A. f. host', *AP* for 'A. p. host', *alt* for 'Alternating hosts')
-	`Fridge`: *1* if the line was revived from backup spores collected after P6
-	`Replicate`: replicate number of this microsporidian x treatment combination (*1*-*4*)
-	`AP.infected`: number of infected hosts at the end of the assay in A. parthenogenetica 
-	`AP.total`: number of surviving hosts at the end of the assay in A. parthenogenetica 
-	`AF.infected`: number of infected hosts at the end of the assay in A. franciscana 
-	`AF.total`: number of surviving hosts at the end of the assay in A. franciscana 

*[Code not shown]*
```{r echo=FALSE}
inf <- read.csv("./DATA final assays - infectivity.csv",sep=";",dec=",")

head(inf)
```

Reshape the dataframe so there is one datapoint per row. New variables:

-	`inf`: number of infected hosts at the end of the assay
-	`TestedHost`: assay host species (*AF* for A. franciscana, *AP* for A. parthenogenetica)
-	`total`: number of surviving hosts at the end of the assay 
-	`noninf`: number of uninfected hosts at the end of the assay
-	`LineTestedHost`: concatenation of `Line` and `TestedHost`, for easier merging later on

```{r echo=FALSE}
inf_melt<-melt(inf,c("MS","Line","EvolvedHost","Fridge","Replicate","AP.total","AF.total"),c("AP.infected","AF.infected"))
names(inf_melt)[9]<-"inf"

inf_melt$TestedHost<-NA
  inf_melt$TestedHost[inf_melt$variable=="AP.infected"]<-"AP"
  inf_melt$TestedHost[inf_melt$variable=="AF.infected"]<-"AF"

inf_melt$total<-NA
  inf_melt$total[inf_melt$TestedHost=="AP"]<-inf_melt[inf_melt$TestedHost=="AP",]$AP.total
  inf_melt$total[inf_melt$TestedHost=="AF"]<-inf_melt[inf_melt$TestedHost=="AF",]$AF.total

inf_melt$noninf<-inf_melt$total-inf_melt$inf
inf_melt<-droplevels(inf_melt)  

inf_melt$LineTestedHost<-paste(inf_melt$Line,inf_melt$TestedHost)

head(inf_melt)
```

### Analysis

__Pretest: are surviving and revived lines different?__

Model: Binomial GLMM with infection as a response variable. A. rigaudi and E. artemiae lines were analyzed separately.

-	Fixed effects: `Fridge`, `EvolvedHost`, `TestedHost`, interactions
-	Random effects: `Line` (to control for pseudoreplication)
-	Exclude treatments for which there were no revived lines.

Statistics: likelihood ratio tests on the fixed effects

*[Code not shown]*
```{r include=FALSE,eval=FALSE}
### MS1

#data
with(inf_melt[inf_melt$MS=="MS1",],table(EvolvedHost,TestedHost,Fridge))
temp<-inf_melt[inf_melt$MS=="MS1"&inf_melt$EvolvedHost=="AF",]

#models
fridge_full<-glmer(cbind(inf,noninf)~Fridge*TestedHost+(1|Line),data=temp,family=binomial)
fridge_a<-glmer(cbind(inf,noninf)~Fridge+TestedHost+(1|Line),data=temp,family=binomial)
fridge_b<-glmer(cbind(inf,noninf)~Fridge+(1|Line),data=temp,family=binomial)
fridge_c<-glmer(cbind(inf,noninf)~TestedHost+(1|Line),data=temp,family=binomial)
anova(fridge_full,fridge_a) #a better with p=0.898 = p of interaction
anova(fridge_a,fridge_c) #c better with p=0.207 = p of fridge, chisq(1)=1.5922



### MS2

#data
with(inf_melt[inf_melt$MS=="MS2",],table(EvolvedHost,TestedHost,Fridge))
temp<-inf_melt[inf_melt$MS=="MS2"&inf_melt$EvolvedHost!="AP",]

#models
fridge_full<-glmer(cbind(inf,noninf)~Fridge*TestedHost*EvolvedHost+(1|Line),data=temp,family=binomial)
fridge_a<-glmer(cbind(inf,noninf)~Fridge+TestedHost+EvolvedHost+Fridge:TestedHost+Fridge:EvolvedHost+TestedHost:EvolvedHost+(1|Line),data=temp,family=binomial)
fridge_b<-glmer(cbind(inf,noninf)~Fridge+TestedHost+EvolvedHost+Fridge:TestedHost+Fridge:EvolvedHost+(1|Line),data=temp,family=binomial)
fridge_c<-glmer(cbind(inf,noninf)~Fridge+TestedHost+EvolvedHost+Fridge:TestedHost+TestedHost:EvolvedHost+(1|Line),data=temp,family=binomial)
fridge_d<-glmer(cbind(inf,noninf)~Fridge+TestedHost+EvolvedHost+Fridge:EvolvedHost+TestedHost:EvolvedHost+(1|Line),data=temp,family=binomial)
anova(fridge_full,fridge_a) #a better with p=0.1354 = p of interaction
anova(fridge_a,fridge_b) #a better with p=0.01402 = p of TestedHost:EvolvedHost
anova(fridge_a,fridge_c) #c better with p=0.6337 = p of fridge:EvolvedHost
anova(fridge_a,fridge_d) #d better with p=0.8641 = p of TestedHost:fridge
fridge_e<-glmer(cbind(inf,noninf)~Fridge+TestedHost+EvolvedHost+TestedHost:EvolvedHost+(1|Line),data=temp,family=binomial)
fridge_f<-glmer(cbind(inf,noninf)~TestedHost+EvolvedHost+TestedHost:EvolvedHost+(1|Line),data=temp,family=binomial)
anova(fridge_e,fridge_f) #e better with p=0.03288 = p of fridge,chisq(1)=4.5518
summary(fridge_e) #fridge lines poorer
```

=> MS1: no effect of revival

=> MS2: significant effect of revival. Revived lines do worse.

__Test: effect of passaging treatment__

Model: Binomial GLMM with infection as a response variable. A. rigaudi and E. artemiae lines were analyzed separately.

-	Fixed effects: `TestedHost`, `EvolvedHost`, interaction
-	Random effects: `Line` (to control for pseudoreplication)
-	Include revived lines if they were not different in the pretest.

Statistics: likelihood ratio tests on the fixed effects

*[Code not shown]*
```{r include=FALSE,eval=FALSE}
### MS1

#data
temp<-inf_melt[inf_melt$MS=="MS1",]

#models
full<-glmer(cbind(inf,noninf)~EvolvedHost*TestedHost+(1|Line),data=temp,family=binomial)
a<-glmer(cbind(inf,noninf)~EvolvedHost+TestedHost+(1|Line),data=temp,family=binomial)
b<-glmer(cbind(inf,noninf)~EvolvedHost+(1|Line),data=temp,family=binomial)
c<-glmer(cbind(inf,noninf)~TestedHost+(1|Line),data=temp,family=binomial)
anova(full,a) #a better with p=0.7305 = p of interaction, chisq(2)=0.628
anova(a,b) #b better with p=0.05453 = p of TestedHost, chisq(1)=3.6964
anova(a,c) #c better with p=0.2465 = p of EvolvedHost, chisq(2)=2.8008

#posthoc
inf_melt$dummy<-paste(inf_melt$TestedHost,inf_melt$EvolvedHost)
ph<-glmer(cbind(inf,noninf)~dummy+(1|Line),data=inf_melt[inf_melt$MS=="MS1",],family=binomial)
summary(glht(ph, mcp(dummy="Tukey"))) # no sign interactions



### MS2

#data
temp<-inf_melt[inf_melt$MS=="MS2"&inf_melt$Fridge==0,]

#models
full<-glmer(cbind(inf,noninf)~EvolvedHost*TestedHost+(1|Line),data=temp,family=binomial)
a<-glmer(cbind(inf,noninf)~EvolvedHost+TestedHost+(1|Line),data=temp,family=binomial)
anova(full,a) #full model better with p=0.01783 = p of interaction, chisq(2)=8.0532
summary(full)

#type III tests for the other factors
options(contrasts=c(unordered="contr.sum",ordered="contr.poly"))
Anova(full,type="III")
options(contrasts=c(unordered="contr.treatment",ordered="contr.poly")) #reset to default

#posthoc
inf_melt$dummy<-paste(inf_melt$TestedHost,inf_melt$EvolvedHost)
ph<-glmer(cbind(inf,noninf)~dummy+(1|Line),data=inf_melt[inf_melt$MS=="MS2"&inf_melt$Fridge==0,],family=binomial)
summary(glht(ph, mcp(dummy="Tukey")))
```

=> MS1: no effect of `EvolvedHost`, a marginally non-sign. effect of `TestedHost`.

=> MS2: significant interaction effect. MS2 lines always infect AF better, and lines that encountered AP in their evolution have higher infectivity in AP than those that only encountered AF.

Model checking:

*[Code not shown]*
```{r include=FALSE,eval=FALSE} 
### MS1

temp<-inf_melt[inf_melt$MS=="MS1",]
d<-glmer(cbind(inf,noninf)~1+(1|Line),data=temp,family=binomial)


## model fit

#overdispersion (for GLMs and GLMMs): underdispersed

overdisp.glmer(d)

#fitted vs. observed: fine

f <- fitted(d)
with(temp,plot(I(inf/total), f))
abline(0,1)


## residuals

r <- residuals(d)

#normal: not normal but that's ok because the regression is binomial

hist(r)
qqnorm(r)
qqline(r)

#constant variance: decreases towards higher values

plot(r)
abline(0,0)
mean(r, na.rm=TRUE)

#independent of size of fitted values: fine

plot(f,r)
sm <- loess(r~f)
v <- seq(min(f,na.rm=T),max(f,na.rm=T),length=101)
lines(v,predict(sm,data.frame(f=v)),col=2)

#independent of independents: NA



### MS2

temp<-inf_melt[inf_melt$MS=="MS2"&inf_melt$Fridge==0,]
full<-glmer(cbind(inf,noninf)~EvolvedHost*TestedHost+(1|Line),data=temp,family=binomial)


## model fit

#overdispersion (for GLMs and GLMMs): underdispersed

overdisp.glmer(full)

#fitted vs. observed: fine

f <- fitted(full)
with(temp,plot(I(inf/total), f))
abline(0,1)


## residuals

r <- residuals(full)

#normal: not normal but that's ok because the regression is binomial

hist(r)
qqnorm(r)
qqline(r)

#constant variance: fine

plot(r)
abline(0,0)
mean(r, na.rm=TRUE)

#independent of size of fitted values: fine

plot(f,r)
sm <- loess(r~f)
v <- seq(min(f,na.rm=T),max(f,na.rm=T),length=101)
lines(v,predict(sm,data.frame(f=v)),col=2)

#independent of independents: fine

plot(r~as.factor(temp$TestedHost))
plot(r~temp$EvolvedHost)
```

=> MS1 model is fine.

=> MS2 model is fine.


## 2nd (virulence & transmission)

### Data

The final assay was run in 3 blocks:

- Block 1 was exposed from day 0-1 and incubated from day 2-91.
- Block 2 was exposed from day 0-1 and incubated from day 2-91.
- Block 3 was exposed from day 1-2 and incubated from day 3-92.

Virulence dataset, containing one datapoint for each assay host. This dataset is expanded to include information on the infection status, resulting in the following final dataset:

-	`LineTestedHost`: individual ID, a concatenation of `Line` and `TestedHost`
-	`Line`: Line ID
-	`TestedHost`: assay host species (*AF* for A. franciscana, *AP* for A. parthenogenetica)
-	`ID`: assay host ID within each line x assay host combination
-	`Sex`: sex of assay host (*m* for male, *f* for female)
-	`Died`: *1* if the host died during the final assay
-	`Death.day`: day on which the host died (*NA* if it didn't die during the assay)
-	`Last.day.seen`: last day on which the host was seen alive (*91*/*92* if the host didn't die during the assay)
-	`Start.day`: day on which the host was first exposed to the parasite (*0* for hosts in block 1 or 2, *1* for hosts in block 3)
-	`MS`: microsporidian species (*MS1* for A. rigaudi, *MS2* for E. artemiae)
-	`EvolvedHost`: evolutionary treatment (*AF* for 'A. f. host', *AP* for 'A. p. host', *alt* for 'Alternating hosts')
-	`Replicate`: replicate number of this microsporidian x treatment combination (*1*-*4*)
-	`Fridge`: *1* if the line was revived from backup spores collected after P6
-	`Death.day.H`: `Death.day` standardized across blocks (`Death.day` - `Start.day`)
-	`Last.day.seen.H`: `Last.day.seen` standardized across blocks (`Last.day.seen` - `Start.day`)
-	`PCR_MS1`: did the individual test positive for A. rigaudi DNA?
-	`PCR_MS2`: did the individual test positive for E. artemiae DNA?
-	`poop_MS1`: did the individual have A. rigaudi spores in its faecal sample?
-	`poop_MS2`: did the individual have E. artemiae spores in its faecal sample?
-	`infected`: taking the PCR and poop information into account, do we consider the individual to be infected?
    + exposed individuals that were found dead before or on day 2 are not PCR-tested & have no poop sample: PCR=*NA*, poop=*NA* -> infection=*NA*
    + exposed individuals that were found dead before or on day 10 are PCR-tested, but didn't have a poop sample yet:
        o PCR=*0*, poop=*NA* -> infection=*NA* (PCR could be a false negative due to degraded DNA)
        o PCR=*1*, poop=*NA* -> infection=*1*
        o PCR=*NA*, poop=*NA* -> infection=*NA* (if the body was lost)
    + exposed individuals that died after day 10 are PCR-tested and had a poop sample if they tested negative:
        o PCR=*0*, poop=*0* -> infection=*0*
        o PCR=*0*, poop=*1* -> infection=*1* (PCR was a false negative due to degraded DNA or reduced parasite load towards the end of the experiment)
        o PCR=*1*, poop=*NA* -> infection=*1* (no need to check the poop if the PCR was positive)
        o PCR=*NA*, poop=*0* -> infection=*0* (if the body was lost but there were stored poop samples)
        o PCR=*NA*, poop=*1* -> infection=*1* (if the body was lost but there were stored poop samples)
    + control individuals are presumed uninfected until proven otherwise; there were no poop samples taken
        o PCR=*NA*, poop=*NA* -> infection=*0*
        o PCR=*0*, poop=*NA* -> infection=*0*
        o PCR=*1*, poop=*NA* -> infection=*1*

*[Code not shown]*
```{r echo=FALSE}
dat<- read.csv("./DATA final assays - virtrans.csv", sep=";",dec=",")

dat<-dat[!is.na(dat$ID),] #remove missing individuals

dat<-merge(dat,treatments,by="Line",all.x=T) #add line data

dat$LineTestedHostID<-paste(dat$Line,dat$TestedHost,dat$ID,sep=" ") #add an individual identifier

dat$Death.day.H<-dat$Death.day-dat$Start.day #homogenize with respect to Start.day
dat$Last.day.seen.H<-dat$Last.day.seen-dat$Start.day #homogenize with respect to Start.day


### determine who was infected or not

## add data on who was infected or not, based on PCRs (susceptible to false 0s)

PCRs<- read.csv("./DATA final assays - PCRs plates 1-29.csv", sep=",",dec=".")

MS1<-PCRs[which(PCRs$msp==1),c("infected","DNA")]
names(MS1)[1]<-"PCR_MS1"
MS2<-PCRs[which(PCRs$msp==2),c("infected","DNA")]
names(MS2)[1]<-"PCR_MS2"

dat<-merge(dat,MS1,by="DNA",all.x=T,sort=F)
dat<-merge(dat,MS2,by="DNA",all.x=T,sort=F)


## add data on who was infected or not, based on checks of poop samples

#link location on the poop collection plate (spore_pl_col_row) to individual ID (LineTestedHostID)
sporecoll <- read.csv("./DATA final assays - virtrans - spore collection.csv", sep=";",dec=",")
sporecoll$LineTestedHostID<-paste(sporecoll$LineTestedHost,sporecoll$ID,sep=" ")
sporecoll$spore_pl_col_row<-paste(sporecoll$spore_pl,sporecoll$spore_col,sporecoll$spore_row,sep="_")
sporecoll<-sporecoll[,c("LineTestedHostID","spore_pl_col_row")]

#link location on the poop collection plate (spore_pl_col_row) to presence/absence of msp in the checks
check <- read.csv("./DATA final assays - virtrans - checks of infection in spore samples.csv", sep=";")
check<-check[,c("spore_pl_col_row","msp1","msp2")]
names(check)[names(check)=="msp1"]<-"poop_MS1"
names(check)[names(check)=="msp2"]<-"poop_MS2"

#merge so that presence/absence of msp in the checks is linked to individual ID (LineTestedHostID)
check<-merge(check,sporecoll,by="spore_pl_col_row",all.x=T,all.y=F)

#merge with overall dataset
dat<-merge(dat,check[,c("poop_MS1","poop_MS2","LineTestedHostID")],by="LineTestedHostID",all.x=T)


## conclude which of the exposed individuals were infected or not based on PCR and poop data

dat$infected<-NA

dat[dat$MS=="MS1",]$infected<-with(dat[dat$MS=="MS1",],ifelse(PCR_MS1==1|poop_MS1==1,1,NA)) #if either one has a positive result, infection is confirmed
dat[which(dat$MS=="MS1"&dat$PCR_MS1==0&dat$poop_MS1==0),]$infected<-0 #confirmed un-infected if there was both a PCR test and a poop sample
dat[which(dat$MS=="MS1"&is.na(dat$PCR_MS1)&dat$poop_MS1==0),]$infected<-0 #confirmed un-infected if the body was lost but there was a poop sample
# with(dat[dat$MS=="MS1",],table(PCR_MS1,poop_MS1,infected,useNA="ifany")) #check that all is well

dat[dat$MS=="MS2",]$infected<-with(dat[dat$MS=="MS2",],ifelse(PCR_MS2==1|poop_MS2==1,1,NA)) #if either one has a positive result, infection is confirmed
dat[which(dat$MS=="MS2"&dat$PCR_MS2==0&dat$poop_MS2==0),]$infected<-0 #confirmed un-infected if there was both a PCR test and a poop sample
dat[which(dat$MS=="MS2"&is.na(dat$PCR_MS2)&dat$poop_MS2==0),]$infected<-0 #confirmed un-infected if the body was lost but there was a poop sample
# with(dat[dat$MS=="MS2",],table(PCR_MS2,poop_MS2,infected,useNA="ifany")) #check that all is well



### check controls for contamination & mark them as (un-)infected accordingly

# with(dat[dat$MS=="Ctrl",],table(PCR_MS1,PCR_MS2,useNA="ifany")) #one Ctrl is infected with MS2 (don't need to check the poop collection columns because I didn't check the poop samples of any Ctrls); 9 have NAs (these died too early for PCR)
dat[which(dat$MS=="Ctrl"&dat$PCR_MS1==0&dat$PCR_MS2==0),]$infected<-0 #any controls that didn't test positive are uninfected
dat[which(dat$MS=="Ctrl"&is.na(dat$PCR_MS1)&is.na(dat$PCR_MS2)),]$infected<-0 #any controls that didn't test positive are uninfected
dat[which(dat$MS=="Ctrl"&dat$PCR_MS2==1),]$infected<-1 #the one that tested positive for MS2 is infected

head(dat)
```

Spore production dataset, containing datapoints at the line x assay host level. This dataset is reworked and reduced to include only the pertinent information:

-	`LineTestedHost`: group ID, a concatenation of `Line` and `TestedHost`
-	`StandardizedCountAve`: a standardized spore count for each `LineTestedHost`, as counted in 9 circles of a Quick-Read Urinalysis Slide (Globe Scientific). Multiplying this count by 350 gives the number of spores in the spore sample of 1 infected host.
-	`Line`: Line ID
-	`Fridge`: *1* if the line was revived from backup spores collected after P6
-	`TestedHost`: assay host species (*AF* for A. franciscana, *AP* for A. parthenogenetica)
-	`nb_inf`: number of infected hosts, i.e. number of hosts that contributed a fecal sample to the pooled sample for that line x assay host
-	`MS`: microsporidian species (*MS1* for A. rigaudi, *MS2* for E. artemiae)
-	`EvolvedHost`: evolutionary treatment (*AF* for 'A. f. host', *AP* for 'A. p. host', *alt* for 'Alternating hosts')

*[Code not shown]*
```{r echo=FALSE}
counts <- read.csv("./DATA final assays - virtrans - transmission week 3.csv",sep=";",dec=",")

counts$LineTestedHost<-paste(counts$Line,counts$TestedHost) #create an identifier


## calculate spores per host and spores in all hosts

# counts$SporesPerHost<-NA
# counts[counts$CountType=="preliminary",]$SporesPerHost<-with(counts[counts$CountType=="preliminary",],(((SporeCount/CirclesCounted)*9)*(Volume*10))/5) #preliminary spore counts only included individuals who tested positive for infection by PCR. For these counts, the total spore mix contained 200uL fecal sample from each host. I took a 1mL subsample of this (~ samples from 5 hosts), concentrated it, and counted the number of spores inside.
# counts[counts$CountType=="final",]$SporesPerHost<-with(counts[counts$CountType=="final",],(((SporeCount/CirclesCounted)*9)*(Volume*10))/2) #final spore counts included spore samples from all infected hosts.For these counts, the total spore mix contained 500uL fecal sample from each host. I took a 1mL subsample of this (~ samples from 2 hosts), concentrated it, and counted the number of spores inside. 

# counts$SporesAllHosts<-counts$SporesPerHost*counts$nb_inf


## checks

# #check whether preliminary and final spore counts are similar
# finalcountsavg<-with(counts[counts$CountType=="final",],tapply(log(SporesAllHosts+1),list(Line,TestedHost),mean,na.rm=T))
# prelimcountsavg<-with(counts[counts$CountType=="preliminary",],tapply(log(SporesAllHosts+1),list(Line,TestedHost),mean,na.rm=T))
# plot(prelimcountsavg[,1],finalcountsavg[,1]) #AF
# cor.test(prelimcountsavg[,1],finalcountsavg[,1]) #AF: cor=0.60, p=0.001382 -> low value due to the one outlier?
# cor.test(prelimcountsavg[1:24,1],finalcountsavg[1:24,1]) #AF without outlier: cor=0.88, p=1.888e-08 -> yes, low value due to outlier
# points(prelimcountsavg[,2],finalcountsavg[,2],col="red") #AP
# cor.test(prelimcountsavg[,2],finalcountsavg[,2]) #AF: cor=0.87, p=1.787e-08
# abline(0,1)

# #check repeatability of the final spore counts
# plot(with(counts[counts$CountType=="final",],table(LineTestedHost))) #all have 2 measures
# temp<-counts[counts$CountType=="final",c("LineTestedHost","SporesAllHosts")]
# temp<-temp[order(temp$LineTestedHost),]
# plot(log(temp$SporesAllHosts[seq(1,99,2)]+1),log(temp$SporesAllHosts[seq(2,100,2)]+1))
# abline(0,1)
# cor.test(log(temp$SporesAllHosts[seq(1,99,2)]+1),log(temp$SporesAllHosts[seq(2,100,2)]+1)) #cor=0.90; p<2.2e-16


## subset for final counts only

counts<-counts[counts$CountType=="final",]


## manually add line 7 (it didn't have any infected individuals)

levels(counts$Line)<-c(levels(counts$Line),"Line 7")
counts<-rbind(counts,
              c("Line 7",0,"AP","final",NA,NA,0,70,"Line 7 AP",NA,NA), #NA for spore production: we don't know what it would have produced in infected individuals
              c("Line 7",0,"AF","final",NA,NA,0,70,"Line 7 AF",NA,NA)) #NA for spore production: we don't know what it would have produced in infected individuals


## classes

# counts$SporesPerHost<-as.numeric(counts$SporesPerHost)
# counts$SporesAllHosts<-as.numeric(counts$SporesAllHosts)
counts$nb_inf<-as.numeric(counts$nb_inf)
counts$TestedHost<-as.factor(counts$TestedHost)
counts$Fridge<-as.numeric(counts$Fridge)
counts$SporeCount<-as.numeric(counts$SporeCount)
counts$CirclesCounted<-as.numeric(counts$CirclesCounted)
counts$Volume<-as.numeric(counts$Volume)


## get line descriptions

counts<-merge(counts,treatments)


## add a standardized spore count to keep the error structure as correct as possible

counts$StandardizedCount<-(counts$SporeCount/counts$CirclesCounted)*9*(counts$Volume/70) #this is the actual count, standardized to 9 circles counted and an original volume of 70uL (this probably does distort the errors a bit, but i can't figure out a way to do it better). dividing this by 2 gives the spore count per host (since this was a count for a 1ml subsample, and each host contributed 500ul, see above)


## average the standardized spore counts (I counted each LineTestedHost sample twice), then round this to make a proper count measure

StandardizedCountAve<-round(with(counts,tapply(StandardizedCount,LineTestedHost,mean)))
temp<-data.frame(StandardizedCountAve)
temp$LineTestedHost<-row.names(temp)

temp2<-counts[,c("LineTestedHost","Line","Fridge","TestedHost","nb_inf","MS","EvolvedHost")]
temp2<-temp2[!duplicated(temp2),]

counts_simplified<-merge(temp,temp2)
head(counts_simplified)
```

### Infection rates

What are the overall rates of detected infection? Only includes individuals that died after at least 11 days (i.e. had a detectable infection, see Virulence dataset).

```{r echo=FALSE}
# with(dat[dat$Last.day.seen.H>10,],tapply(infected,list(EvolvedHost,TestedHost,MS),mean,na.rm=T))
with(dat[dat$Last.day.seen.H>10,],tapply(infected,list(TestedHost,MS),mean,na.rm=T)) 
```

### Analysis: virulence

__Pretest: which parametric survival model should we use to analyse the survival curves?__

*[Code not shown]*
```{r include=FALSE,eval=FALSE}
## AP

temp<-dat[dat$TestedHost=="AP"&dat$EvolvedHost!="Ctrl"&dat$infected==1&dat$Last.day.seen.H>10,]
temp<-rbind(temp,dat[dat$TestedHost=="AP"&dat$EvolvedHost=="Ctrl"&dat$infected==0&dat$Last.day.seen.H>10,])
temp<-droplevels(temp[!is.na(temp$LineTestedHostID),])
temp$Last.day.seen.H10<-temp$Last.day.seen.H-10
temp<-droplevels(temp)

weib<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="weibull")
exp<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="exponential")
extreme<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="extreme")
gaussian<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="gaussian")
log<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="logistic")
lognorm<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="lognormal")
loglog<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="loglogistic")
ray<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="rayleigh")

comparison<-AICc(weib,exp,extreme,gaussian,log,lognorm,loglog,ray)
comparison[order(comparison$AICc),]
  #                df     AICc
  # loglog   8.001204 3565.106
  # lognorm  8.001126 3567.370
  # weib     8.000155 3582.785
  # exp      7.000136 3582.895
  # gaussian 8.000117 3797.960
  # ray      8.492613 3807.131
  # log      8.058034 3821.182
  # extreme  8.060147 3898.493


## AF

temp<-dat[dat$TestedHost=="AF"&dat$EvolvedHost!="Ctrl"&dat$infected==1&dat$Last.day.seen.H>10,]
temp<-rbind(temp,dat[dat$TestedHost=="AF"&dat$EvolvedHost=="Ctrl"&dat$infected==0&dat$Last.day.seen.H>10,])
temp<-droplevels(temp[!is.na(temp$LineTestedHostID),])
temp$Last.day.seen.H10<-temp$Last.day.seen.H-10
temp<-droplevels(temp)

weib<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="weibull")
exp<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="exponential")
extreme<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="extreme")
gaussian<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="gaussian")
log<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="logistic")
lognorm<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="lognormal")
loglog<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="loglogistic")
ray<-survreg(Surv(Last.day.seen.H10,Died)~MS*EvolvedHost+frailty(Line,method="aic"),data=temp,dist="rayleigh")

comparison<-AICc(weib,exp,extreme,gaussian,log,lognorm,loglog,ray)
comparison[order(comparison$AICc),]
  #                df     AICc
  # weib     8.003078 3602.680
  # loglog   8.002319 3608.657
  # ray      7.005512 3633.265
  # lognorm  8.001689 3637.923
  # exp      7.001376 3650.899
  # gaussian 8.068430 3660.626
  # log      8.066501 3684.916
  # extreme  8.084279 3704.847
```

=> Try loglogistic: it's the best for AP and second-best (though pretty far from the best) for AF.

__Pretest: are surviving and revived lines different?__

Model: log-logistic survival models, stratified over `TestedHost` (this allows the host species to have a different baseline survival shape, https://stat.ethz.ch/R-manual/R-devel/library/survival/html/survreg.html). A. rigaudi and E. artemiae lines were analyzed separately.

-	Fixed effects: `Fridge`, `EvolvedHost`, `TestedHost`, interactions
-	Random effects: `Line` (to control for pseudoreplication)
-	Exclude treatments for which there were no revived lines.
- Exclude exposed-but-uninfected hosts, as well as hosts that died before day 11 (infection couldn't be reliably detected before day 11, see Virulence dataset; cf. Lievens et al. 2018)

Statistics: likelihood ratio tests on the fixed effects

*[Code not shown]*
```{r include=FALSE,eval=FALSE}
### MS1

#data
temp<-dat[dat$MS=="MS1"&dat$infected==1&dat$Last.day.seen.H>10,] #only those that were successfully infected with MS1 during the final assays
temp<-rbind(temp,dat[dat$MS=="Ctrl"&dat$infected==0&dat$Last.day.seen.H>10,]) #adds the controls that weren't contaminated
temp<-droplevels(temp[!is.na(temp$LineTestedHostID),])
temp$Last.day.seen.H10<-temp$Last.day.seen.H-10
temp<-droplevels(temp)
temp$EvolvedHost<-factor(temp$EvolvedHost,levels(temp$EvolvedHost)[c(4,1:3)]) #reorders so 'Ctrl' becomes the first level of EvolvedHost, ie so that it is the intercept

#models
with(temp,table(Fridge,EvolvedHost))
now<-temp[temp$EvolvedHost=="AF",]
fridge_full<-survreg(Surv(Last.day.seen.H10,Died)~TestedHost*Fridge+strata(TestedHost),data=now,dist="loglogistic") #removed +frailty(Line,method="aic") because model doesn't work with it in
fridge_a<-survreg(Surv(Last.day.seen.H10,Died)~TestedHost+Fridge+strata(TestedHost),data=now,dist="loglogistic")
anova(fridge_full,fridge_a) #interaction: p = 0.4364855
fridge_a<-survreg(Surv(Last.day.seen.H10,Died)~TestedHost+Fridge+strata(TestedHost)+frailty(Line,method="aic"),data=now,dist="loglogistic") #re-included frailty(Line,method="aic") to test fridge effect
fridge_b<-survreg(Surv(Last.day.seen.H10,Died)~TestedHost+strata(TestedHost)+frailty(Line,method="aic"),data=now,dist="loglogistic")
anova(fridge_a,fridge_b) # fridge effect: X2(-0.9999335)=-0.2149281, p = 0.6429212



### MS2

#data
temp<-dat[dat$MS=="MS2"&dat$infected==1&dat$Last.day.seen.H>10,] #only those that were successfully infected with MS1 during the final assays
temp<-rbind(temp,dat[dat$MS=="Ctrl"&dat$infected==0&dat$Last.day.seen.H>10,]) #adds the controls that weren't contaminated
temp<-droplevels(temp[!is.na(temp$LineTestedHostID),])
temp$Last.day.seen.H10<-temp$Last.day.seen.H-10
temp<-droplevels(temp)
temp$EvolvedHost<-factor(temp$EvolvedHost,levels(temp$EvolvedHost)[c(4,1:3)]) #reorders so 'Ctrl' becomes the first level of EvolvedHost, ie so that it is the intercept

#models
with(temp,table(Fridge,EvolvedHost))
now<-temp[temp$EvolvedHost!="AP"&temp$EvolvedHost!="Ctrl",]
fridge_full<-survreg(Surv(Last.day.seen.H10,Died)~TestedHost*Fridge+strata(TestedHost)+frailty(Line,method="aic"),data=now,dist="loglogistic")
fridge_a<-survreg(Surv(Last.day.seen.H10,Died)~TestedHost+Fridge+strata(TestedHost)+frailty(Line,method="aic"),data=now,dist="loglogistic")
fridge_b<-survreg(Surv(Last.day.seen.H10,Died)~TestedHost+strata(TestedHost)+frailty(Line,method="aic"),data=now,dist="loglogistic")
anova(fridge_full,fridge_a) #interaction: p = 0.6597018
anova(fridge_a,fridge_b) # fridge effect: X2(1.1531)=0.9288465,p = 0.3856883
```

=> MS1: no effect of revival

=> MS2: no effect of revival

__Test: effect of passaging treatment on host survival__

Model: log-logistic survival models, stratified over `TestedHost` (this allows the host species to have a different baseline survival shape, https://stat.ethz.ch/R-manual/R-devel/library/survival/html/survreg.html). A. rigaudi and E. artemiae lines were analyzed separately.

-	Fixed effects: `TestedHost`, `EvolvedHost`, interaction
-	Frailty effects: `Line` (to control for pseudoreplication)
-	Include revived lines if they were not different in the pretest.
- Exclude exposed-but-uninfected hosts, as well as hosts that died before day 11 (infection couldn't be reliably detected before day 11, see Virulence dataset; cf. Lievens et al. 2018)

Statistics: likelihood ratio tests on the fixed effects

*[Code not shown]*
```{r include=FALSE,eval=FALSE}
### MS1

#data
temp<-dat[dat$MS=="MS1"&dat$infected==1&dat$Last.day.seen.H>10,] #only those that were successfully infected with MS1 during the final assays
temp<-rbind(temp,dat[dat$MS=="Ctrl"&dat$infected==0&dat$Last.day.seen.H>10,]) #adds the controls that weren't contaminated
temp<-droplevels(temp[!is.na(temp$LineTestedHostID),])
temp$Last.day.seen.H10<-temp$Last.day.seen.H-10
temp<-droplevels(temp)
temp$EvolvedHost<-factor(temp$EvolvedHost,levels(temp$EvolvedHost)[c(4,1:3)]) #reorders so 'Ctrl' becomes the first level of EvolvedHost, ie so that it is the intercept

#LRT
surv1_full<-survreg(Surv(Last.day.seen.H10,Died)~TestedHost*EvolvedHost+strata(TestedHost)+frailty(Line,method="aic"),data=temp,dist="loglogistic")
surv1_a<-survreg(Surv(Last.day.seen.H10,Died)~TestedHost+EvolvedHost+strata(TestedHost)+frailty(Line,method="aic"),data=temp,dist="loglogistic")
anova(surv1_full,surv1_a) #interaction: X2(-3.000052)=-10.43558, p=0.01520504
summary(surv1_full)

#type III tests for the others:
options(contrasts=c(unordered="contr.sum",ordered="contr.poly"))
Anova(surv1_full,type="III")
options(contrasts=c(unordered="contr.treatment",ordered="contr.poly")) #reset to default

#posthoc
lsmeans(surv1_full,pairwise~TestedHost*EvolvedHost,adjust="tukey")$contrasts



### MS2

#data
temp<-dat[dat$MS=="MS2"&dat$infected==1&dat$Last.day.seen.H>10,] #only those that were successfully infected with MS1 during the final assays
temp<-rbind(temp,dat[dat$MS=="Ctrl"&dat$infected==0&dat$Last.day.seen.H>10,]) #adds the controls that weren't contaminated
temp<-droplevels(temp[!is.na(temp$LineTestedHostID),])
temp$Last.day.seen.H10<-temp$Last.day.seen.H-10
temp<-droplevels(temp)
temp$EvolvedHost<-factor(temp$EvolvedHost,levels(temp$EvolvedHost)[c(4,1:3)]) #reorders so 'Ctrl' becomes the first level of EvolvedHost, ie so that it is the intercept

#LRT
surv2_full<-survreg(Surv(Last.day.seen.H10,Died)~TestedHost*EvolvedHost+strata(TestedHost)+frailty(Line,method="aic"),data=temp,dist="loglogistic")
surv2_a<-survreg(Surv(Last.day.seen.H10,Died)~TestedHost+EvolvedHost+strata(TestedHost)+frailty(Line,method="aic"),data=temp,dist="loglogistic")
surv2_b<-survreg(Surv(Last.day.seen.H10,Died)~TestedHost+strata(TestedHost)+frailty(Line,method="aic"),data=temp,dist="loglogistic")
surv2_c<-survreg(Surv(Last.day.seen.H10,Died)~EvolvedHost+strata(TestedHost)+frailty(Line,method="aic"),data=temp,dist="loglogistic")
anova(surv2_full,surv2_a) #intreaction: X2(-4.487883)=-5.785175,p=0.2686614
anova(surv2_a,surv2_b) #evolvedhost: X2(3.551824)=3.719976,p=0.3768402
anova(surv2_a,surv2_c) #testedhost: X2(-4.00538)=-58.11918,p=7.261798e-12
summary(surv2_b) #APs die more
```

=> MS1: significant interaction effect. Lines evolved on AP are less virulent when tested on AP; APs also die more quickly in general.

=> MS2: significant `TestedHost` effect. APs die more quickly.

Model checking:

*[Code not shown]*
```{r include=FALSE,eval=FALSE}
### MS1

par(mfrow=c(1,2))

# AF
plot(survfit(Surv(Last.day.seen.H-10,Died)~1,data=dat[dat$TestedHost=="AF"&dat$MS=="Ctrl"&dat$infected==0&dat$Last.day.seen.H>10,]),conf.int=F,col="gray",lwd=2,main="AF infected with MS1")
# # all lines
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="alt",]),col=adjustcolor("purple",.5),conf.int=F,lwd=1.8)
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="AF"&temp$Fridge==0,]),col=adjustcolor("red",.5),conf.int=F,lwd=1.8)
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="AF"&temp$Fridge==1,]),col=adjustcolor("red",.5),conf.int=F,lwd=1.8,lty="dashed")
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="AP",]),col=adjustcolor("blue",.5),conf.int=F,lwd=1.8)
# means
lines(survfit(Surv(Last.day.seen.H10,Died)~1,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="alt",]),col="purple",conf.int=F,lwd=2)
lines(survfit(Surv(Last.day.seen.H10,Died)~1,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="AF",]),col="red",conf.int=F,lwd=2)
lines(survfit(Surv(Last.day.seen.H10,Died)~1,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="AP",]),col="blue",conf.int=F,lwd=2)
# predictions
lines(flexsurvreg(Surv(Last.day.seen.H10,Died)~EvolvedHost,data=temp[temp$TestedHost=="AF",],dist="llogis"),col="black")

# AP
plot(survfit(Surv(Last.day.seen.H-10,Died)~1,data=dat[dat$TestedHost=="AP"&dat$MS=="Ctrl"&dat$infected==0&dat$Last.day.seen.H>10,]),conf.int=F,col="gray",lwd=2,main="AP infected with MS1")
# # all lines
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="alt",]),col=adjustcolor("purple",.5),conf.int=F,lwd=1.8)
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="AF"&temp$Fridge==0,]),col=adjustcolor("red",.5),conf.int=F,lwd=1.8)
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="AF"&temp$Fridge==1,]),col=adjustcolor("red",.5),conf.int=F,lwd=1.8,lty="dashed")
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="AP",]),col=adjustcolor("blue",.5),conf.int=F,lwd=1.8)
# means
lines(survfit(Surv(Last.day.seen.H10,Died)~1,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="alt",]),col="purple",conf.int=F,lwd=2)
lines(survfit(Surv(Last.day.seen.H10,Died)~1,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="AF",]),col="red",conf.int=F,lwd=2)
lines(survfit(Surv(Last.day.seen.H10,Died)~1,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="AP",]),col="blue",conf.int=F,lwd=2)
# predictions
lines(flexsurvreg(Surv(Last.day.seen.H10,Died)~EvolvedHost,data=temp[temp$TestedHost=="AP",],dist="llogis"),col="black")



### MS2


par(mfrow=c(1,2))

# AF
plot(survfit(Surv(Last.day.seen.H-10,Died)~1,data=dat[dat$TestedHost=="AF"&dat$MS=="Ctrl"&dat$infected==0&dat$Last.day.seen.H>10,]),conf.int=F,col="gray",lwd=2,main="AF infected with MS2")
# # all lines
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="alt",]),col=adjustcolor("purple",.5),conf.int=F,lwd=1.8)
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="AF"&temp$Fridge==0,]),col=adjustcolor("red",.5),conf.int=F,lwd=1.8)
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="AF"&temp$Fridge==1,]),col=adjustcolor("red",.5),conf.int=F,lwd=1.8,lty="dashed")
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="AP",]),col=adjustcolor("blue",.5),conf.int=F,lwd=1.8)
# means
lines(survfit(Surv(Last.day.seen.H10,Died)~1,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="alt",]),col="purple",conf.int=F,lwd=2)
lines(survfit(Surv(Last.day.seen.H10,Died)~1,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="AF",]),col="red",conf.int=F,lwd=2)
lines(survfit(Surv(Last.day.seen.H10,Died)~1,data=temp[temp$TestedHost=="AF"&temp$EvolvedHost=="AP",]),col="blue",conf.int=F,lwd=2)
# predictions
lines(flexsurvreg(Surv(Last.day.seen.H10,Died)~EvolvedHost,data=temp[temp$TestedHost=="AF",],dist="llogis"),col="black")

# AP
plot(survfit(Surv(Last.day.seen.H-10,Died)~1,data=dat[dat$TestedHost=="AP"&dat$MS=="Ctrl"&dat$infected==0&dat$Last.day.seen.H>10,]),conf.int=F,col="gray",lwd=2,main="AP infected with MS2")
# # all lines
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="alt",]),col=adjustcolor("purple",.5),conf.int=F,lwd=1.8)
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="AF"&temp$Fridge==0,]),col=adjustcolor("red",.5),conf.int=F,lwd=1.8)
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="AF"&temp$Fridge==1,]),col=adjustcolor("red",.5),conf.int=F,lwd=1.8,lty="dashed")
# lines(survfit(Surv(Last.day.seen.H10,Died)~Line,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="AP",]),col=adjustcolor("blue",.5),conf.int=F,lwd=1.8)
# means
lines(survfit(Surv(Last.day.seen.H10,Died)~1,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="alt",]),col="purple",conf.int=F,lwd=2)
lines(survfit(Surv(Last.day.seen.H10,Died)~1,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="AF",]),col="red",conf.int=F,lwd=2)
lines(survfit(Surv(Last.day.seen.H10,Died)~1,data=temp[temp$TestedHost=="AP"&temp$EvolvedHost=="AP",]),col="blue",conf.int=F,lwd=2)
# predictions
lines(flexsurvreg(Surv(Last.day.seen.H10,Died)~EvolvedHost,data=temp[temp$TestedHost=="AP",],dist="llogis"),col="black")
```

=> MS1: loglog works well, except to fit the AF Ctrls' survival curve

=> MS2: loglog works well, except to fit the AF Ctrls' survival curve

### Analysis: spore production

__Pretest: are surviving and revived lines different?__

Model: Negative binomial GLMM with the standardized spore count as a response variable. A. rigaudi and E. artemiae lines were analyzed separately.

-	Fixed effects: `Fridge`, `EvolvedHost`, `TestedHost`, interactions
-	Random effects: `Line` (to control for pseudoreplication)
-	Exclude treatments for which there were no revived lines.
- Exposed-but-uninfected hosts are automatically excluded, because spore samples were only collected from infected hosts.

Statistics: likelihood ratio tests on the fixed effects

*[Code not shown]*
```{r include=FALSE,eval=FALSE}
### MS1

#data
temp<-counts_simplified[counts_simplified$MS=="MS1",]

#models
with(temp,table(EvolvedHost,TestedHost,Fridge))
fridge_full<-glmer.nb(StandardizedCountAve~Fridge*TestedHost+(1|Line),data=temp[temp$EvolvedHost=="AF",]) #warning
fridge_a<-glmer.nb(StandardizedCountAve~Fridge+TestedHost+(1|Line),data=temp[temp$EvolvedHost=="AF",]) #warning
anova(fridge_full,fridge_a) #a better with X2(1)=1.9515, p=0.1624 = p of interaction
fridge_b<-glmer.nb(StandardizedCountAve~TestedHost+(1|Line),data=temp[temp$EvolvedHost=="AF",]) #warning
anova(fridge_a,fridge_b) #effect of fridge X2(1)=0.1622, p=0.6871



### MS2

#data
temp<-counts_simplified[counts_simplified$MS=="MS2",]

#models
with(temp,table(EvolvedHost,TestedHost,Fridge))
fridge_full<-glmer.nb(StandardizedCountAve~Fridge*TestedHost*EvolvedHost+(1|Line),data=temp[temp$EvolvedHost!="AP",]) #warning
fridge_a<-glmer.nb(StandardizedCountAve~Fridge+TestedHost+EvolvedHost+Fridge:TestedHost+Fridge:EvolvedHost+TestedHost:EvolvedHost+(1|Line),data=temp[temp$EvolvedHost!="AP",])
anova(fridge_full,fridge_a) #a better with X2(1)=0.2092, p=0.6474 = p of interaction
fridge_b<-glmer.nb(StandardizedCountAve~Fridge+TestedHost+EvolvedHost+Fridge:TestedHost+Fridge:EvolvedHost+(1|Line),data=temp[temp$EvolvedHost!="AP",])
fridge_c<-glmer.nb(StandardizedCountAve~Fridge+TestedHost+EvolvedHost+Fridge:TestedHost+TestedHost:EvolvedHost+(1|Line),temp[temp$EvolvedHost!="AP",])
fridge_d<-glmer.nb(StandardizedCountAve~Fridge+TestedHost+EvolvedHost+Fridge:EvolvedHost+TestedHost:EvolvedHost+(1|Line),data=temp[temp$EvolvedHost!="AP",])
anova(fridge_a,fridge_b) #b better with X2(1)=0.0519,p=0.8199 = p of TestedHost:EvolvedHost
anova(fridge_a,fridge_c) #c better with p=0.5855 = p of fridge:EvolvedHost
anova(fridge_a,fridge_d) #d better with p=0.8743 = p of TestedHost:fridge
fridge_e<-glmer.nb(StandardizedCountAve~Fridge+TestedHost+EvolvedHost+(1|Line),data=temp[temp$EvolvedHost!="AP",]) #warning
fridge_f<-glmer.nb(StandardizedCountAve~TestedHost+EvolvedHost+(1|Line),data=temp[temp$EvolvedHost!="AP",])
fridge_g<-glmer.nb(StandardizedCountAve~Fridge+EvolvedHost+(1|Line),data=temp[temp$EvolvedHost!="AP",])
fridge_h<-glmer.nb(StandardizedCountAve~Fridge+TestedHost+(1|Line),data=temp[temp$EvolvedHost!="AP",]) #warning
anova(fridge_e,fridge_f) #e better with p=0.3818 = p of fridge, chisq(1)=0.7649
anova(fridge_e,fridge_g) #e better with p=8.52e-05 = p of TestedHost
anova(fridge_e,fridge_h) #e better with p=0.8621 = p of EvolvedHost
```

=> MS1: no effect of revival

=> MS2: no effect of revival

__Test: effect of passaging treatment on spore production__

Model: Negative binomial GLMM with the standardized spore count as a response variable. A. rigaudi and E. artemiae lines were analyzed separately. Observations are not weighted by the number of infected individuals in the trial (the more individuals there were, the more precise the spore sample is expected to be). I checked the effects of including/not weights and they were negligble; it seems like the variation introduced by the number of hosts was not that important. Therefore, I choose to run the less heavy models.

-	Fixed effects: `TestedHost`, `EvolvedHost`, interaction
-	Frailty effects: `Line` (to control for pseudoreplication)
-	Include revived lines if they were not different in the pretest.
- Exposed-but-uninfected hosts are automatically excluded, because spore samples were only collected from infected hosts.

Statistics: likelihood ratio tests on the fixed effects

*[Code not shown]*
```{r include=FALSE,eval=FALSE}
### MS1

#data
temp<-counts_simplified[counts_simplified$MS=="MS1",]

# #poisson or negbin?
# full<-glmer(StandardizedCountAve~EvolvedHost*TestedHost+(1|Line),data=temp,family=poisson)
# overdisp.glmer(full) #1.794
# full_nb<-glmer.nb(StandardizedCountAve~EvolvedHost*TestedHost+(1|Line),data=temp)
# AICc(full,full_nb) #poisson fits better, but go with nb anyway (slightt overdispersion + nb is more robust)

#LRT
full<-glmer.nb(StandardizedCountAve~EvolvedHost*TestedHost+(1|Line),data=temp)#warning, but SEs are ok
a<-glmer.nb(StandardizedCountAve~EvolvedHost+TestedHost+(1|Line),data=temp)
b<-glmer.nb(StandardizedCountAve~EvolvedHost+(1|Line),data=temp)
c<-glmer.nb(StandardizedCountAve~TestedHost+(1|Line),data=temp)
anova(full,a) #p=0.3581 = p of interaction, chisq(2)=2.0542
anova(a,b) #p=0.0001211 = p of TestedHost, chisq(1)=14.776
anova(a,c) #p=0.9868 = p of EvolvedHost, chsiq(2)=0.0266

#Tukey
temp$dummy<-paste(temp$EvolvedHost,temp$TestedHost)
ph<-glmer.nb(StandardizedCountAve~dummy+(1|Line),data=temp)
summary(glht(ph, mcp(dummy="Tukey")))



### MS2

#data
temp<-counts_simplified[counts_simplified$MS=="MS2",]

# #poisson or negbin?
# full<-glmer(StandardizedCountAve~EvolvedHost*TestedHost+(1|Line),data=temp,family=poisson)
# overdisp.glmer(full) #4.341
# full_nb<-glmer.nb(StandardizedCountAve~EvolvedHost*TestedHost+(1|Line),data=temp)
# AICc(full,full_nb) #neg bin fits better

#LRT
full<-glmer.nb(StandardizedCountAve~EvolvedHost*TestedHost+(1|Line),data=temp) #warning, but SEs are ok
a<-glmer.nb(StandardizedCountAve~EvolvedHost+TestedHost+(1|Line),data=temp)
b<-glmer.nb(StandardizedCountAve~EvolvedHost+(1|Line),data=temp) #warning, but SEs are ok
c<-glmer.nb(StandardizedCountAve~TestedHost+(1|Line),data=temp)
anova(full,a) #p=0.4314 = p of interaction, chisq(2)=1.6816
anova(a,b) #p=6.174e-05 = p of TestedHost, chisq(1)=16.048
anova(a,c) #p=0.7109 = p of EvolvedHost, chsiq(2)=0.6824
summary(c) #higher spore production in AF
```

=> MS1: significant `TestedHost` effect. Spore production is higher in AP.

=> MS2: significant `TestedHost` effect. Spore production is higher in AF.

Model checking:

*[Code not shown]*
```{r include=FALSE,eval=FALSE}
### MS1

temp<-counts_simplified[counts_simplified$MS=="MS1",]
c<-glmer.nb(StandardizedCountAve~TestedHost+(1|Line),data=temp)


## model fit

#overdispersion (for GLMs and GLMMs): fine

overdisp.glmer(c)

#fitted vs. observed: fine

f <- fitted(c)
plot(temp$StandardizedCountAve[!is.na(temp$StandardizedCountAve)], f)
abline(0,1)


## residuals

r <- residuals(c)

#normal: fine

hist(r)
qqnorm(r)
qqline(r)

#constant variance: fine

plot(r)
abline(0,0)
mean(r, na.rm=TRUE)

#independent of size of fitted values: not great

plot(f,r)
sm <- loess(r~f)
v <- seq(min(f,na.rm=T),max(f,na.rm=T),length=101)
lines(v,predict(sm,data.frame(f=v)),col=2)

#independent of independents: fine

plot(r~temp$TestedHost[!is.na(temp$StandardizedCountAve)])



### MS2

temp<-counts_simplified[counts_simplified$MS=="MS2",]
c<-glmer.nb(StandardizedCountAve~TestedHost+(1|Line),data=temp)


## model fit

#overdispersion (for GLMs and GLMMs): fine

overdisp.glmer(c)

#fitted vs. observed: fine

f <- fitted(c)
plot(temp$StandardizedCount, f)
abline(0,1)


## residuals

r <- residuals(c)

#normal: fine

hist(r)
qqnorm(r)
qqline(r)

#constant variance: fine

plot(r)
abline(0,0)
mean(r, na.rm=TRUE)

#independent of size of fitted values: fine

plot(f,r)
sm <- loess(r~f)
v <- seq(min(f,na.rm=T),max(f,na.rm=T),length=101)
lines(v,predict(sm,data.frame(f=v)),col=2)

#independent of independents: fine

plot(r~temp$TestedHost)
```

=> MS1: overall fine

=> MS2: fine


## Fitness

### Data

New variables appended to the 'Spore production' dataset, requiring the 'Infectivity' dataset (-> PREPARE THOSE DATASETS FIRST!).

-	`inf`: number of infected hosts at the end of the assay
-	`total`: number of surviving hosts at the end of the assay 
-	`infectiousness`: the probability of a single spore to start an infection, calculated based on an independent action model with birth-death processes (Schmid-Hempel 2011, pg. 225-6).
-	`Fitness`: the projected number of infections occurring if the line were passaged onto a new set of susceptible hosts, calculated as the total number of spores produced by the surviving individuals over a two-day period after three weeks of incubation (thus virulence is implicit), multiplied by the infectiousness of a single spore.

*[Code not shown]*
```{r echo=FALSE}
counts_simplified<-merge(counts_simplified,inf_melt[,c("inf","total","LineTestedHost")])

counts_simplified$infectiousness<-with(counts_simplified,(-1)*log((total-inf)/total)/750) #dose 750 spores
counts_simplified$Fitness<-with(counts_simplified,infectiousness*(StandardizedCountAve*10*70*.5)*nb_inf) #number of new infections that can be started in a two-day period by all the infected hosts together
counts_simplified$Fitness[counts_simplified$Line=="Line 7"]<-0 #fitness = 0 because no hosts were infected

head(counts_simplified)
```

### Analysis

__Test: effect of passaging treatment on spore production__

Model: LMM with ln+1-transformed fitness as a response variable. A. rigaudi and E. artemiae lines were analyzed separately. There is an issue with error propagation (each measurement has its error) that I wasn't able to solve, but the results seem to be robust.

-	Fixed effects: `TestedHost`, `EvolvedHost`, interaction
-	Frailty effects: `Line` (to control for pseudoreplication)
-	Exclude revived lines for MS2 (shown to have an effect for infectivity)

Statistics: likelihood ratio tests on the fixed effects

*[Code not shown]*
```{r include=FALSE,eval=FALSE}
### MS1

#data
temp<-counts_simplified[counts_simplified$MS=="MS1",]
temp$lnFitness<-log(temp$Fitness+1)

#LRT
full<-lmer(lnFitness~EvolvedHost*TestedHost+(1|Line),data=temp)
a<-lmer(lnFitness~EvolvedHost+TestedHost+(1|Line),data=temp)
b<-lmer(lnFitness~EvolvedHost+(1|Line),data=temp)
c<-lmer(lnFitness~TestedHost+(1|Line),data=temp)
anova(full,a) #a better with p=0.2926 = p of interaction, chisq(2)=2.458
anova(a,b) #a better with p=0.0005755 = p of TestedHost, chisq(1)=11.854
anova(a,c) #c better with p=0.3356 = p of EvolvedHost, chsiq(2)=2.1834
summary(c)



### MS2

#data
temp<-counts_simplified[counts_simplified$MS=="MS2"&counts_simplified$Fridge==0,]
temp$lnFitness<-log(temp$Fitness+1)

#LRT
full<-lmer(lnFitness~EvolvedHost*TestedHost+(1|Line),data=temp)
a<-lmer(lnFitness~EvolvedHost+TestedHost+(1|Line),data=temp)
anova(full,a) #full better with p=0.0414 = p of interaction, chisq(2)=6.369
summary(full)

#type III tests for the others:
options(contrasts=c(unordered="contr.sum",ordered="contr.poly"))
Anova(full,type="III")
options(contrasts=c(unordered="contr.treatment",ordered="contr.poly")) #reset to default
```

=> MS1: significant `TestedHost` effect. Fitness is higher in AP.
=> MS2: significant interaction effect. Fitness is higher in AF but lines evolved on AP or alternating hosts do better when tested on AP.

Model checking:

*[Code not shown]*
```{r include=FALSE,eval=FALSE}
### MS1

temp<-counts_simplified[counts_simplified$MS=="MS1",]
temp$lnFitness<-log(temp$Fitness+1)
c<-lmer(lnFitness~TestedHost+(1|Line),data=temp)


## model fit

#overdispersion (for GLMs and GLMMs): underdispersed

overdisp.glmer(c)

#fitted vs. observed: fine

f <- fitted(c)
with(temp,plot(lnFitness, f,col=as.numeric(as.factor(paste(EvolvedHost,TestedHost)))))
abline(0,1) #good


## residuals

r <- residuals(c)

#normal: skewed at low end

hist(r)
qqnorm(r)
qqline(r)

#constant variance: decreases a little towards higher values, but still ok

plot(r)
abline(0,0)
mean(r, na.rm=TRUE)

#independent of size of fitted values: fine

plot(f,r)
sm <- loess(r~f)
v <- seq(min(f,na.rm=T),max(f,na.rm=T),length=101)
lines(v,predict(sm,data.frame(f=v)),col=2)

#independent of independents: fine

plot(r~temp$TestedHost)



### MS2

temp<-counts_simplified[counts_simplified$MS=="MS2"&counts_simplified$Fridge==0,]
temp$lnFitness<-log(temp$Fitness+1)
full<-lmer(lnFitness~EvolvedHost*TestedHost+(1|Line),data=temp)


## model fit

#overdispersion (for GLMs and GLMMs): underdispersed

overdisp.glmer(full)

#fitted vs. observed: fine

f <- fitted(full)
with(temp[temp$Fridge==0,],plot(lnFitness, f,col=as.numeric(as.factor(paste(EvolvedHost,TestedHost)))))
abline(0,1) #good


## residuals

r <- residuals(full)

#normal: fine

hist(r)
qqnorm(r)
qqline(r)

#constant variance: fine

plot(r)
abline(0,0)
mean(r, na.rm=TRUE)

#independent of size of fitted values: fine

plot(f,r)
sm <- loess(r~f)
v <- seq(min(f,na.rm=T),max(f,na.rm=T),length=101)
lines(v,predict(sm,data.frame(f=v)),col=2)

#independent of independents: fine

plot(r~temp$TestedHost)
plot(r~temp$EvolvedHost)
```

=> MS1: slightly off, but still passable

=> MS2: fine