---
title: "Statistical Analysis"
author: "Claire A. McLean"
date: ""
output:
  pdf_document: default
  html_document: default
---

```{r libraries, echo=FALSE, message=FALSE,warning=FALSE}
library(lme4)
library(lmerTest)
library(MuMIn)
library(knitr)
```
# Divergent male and female mate preferences do not explain incipient speciation between lizard lineages
## Claire A. McLean ^1^^,^ ^2^, Richard A. Bartle ^1^, Caroline M. Dong ^1^^,^ ^2^, Katrina J. Rankin ^1^ & Devi Stuart-Fox ^1^
###^1^ The University of Melbourne, ^2^ Museums Victoria

 

##**Female-male behavioural trials**

```{r, echo=FALSE}
female_male_data <- read.csv("female_male_data.csv", header=TRUE, sep=",", dec=".", strip.white=TRUE)
female_male_data$pairing_number<-factor(female_male_data$pairing_number)
```
####The dataset consists of `r nrow(female_male_data)` behavioural trials involving `r length(unique(female_male_data$female_ID))` female and `r length(unique(female_male_data$male_ID))` male tawny dragons (*Ctenophorus decresii*).

```{r, echo=FALSE}
options(contrasts = c("contr.sum","contr.poly"))
```

&nbsp;

###*Model 1: Number of copulation attempts*
####Generalised linear mixed model with female lineage, male lineage, female lineage * male lineage as fixed effects and female ID, male ID and pairing number (female's first or second trial) as random effects.

```{r, message=FALSE, warning=FALSE}
lmer.copulation.attempts.model= lmer(cop_attempts~female_lineage+male_lineage+female_lineage*male_lineage+(1|female_ID)+(1|male_ID)+(1|pairing_number),data=female_male_data, REML = FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```
####Check the normality of the residuals
```{r}
plot(lmer.copulation.attempts.model)
qqnorm(resid(lmer.copulation.attempts.model))
qqline(resid(lmer.copulation.attempts.model))
shapiro.test(resid(lmer.copulation.attempts.model))
```
####log transform the data to meet model assumptions of normality
```{r, message=FALSE, warning=FALSE}
log.lmer.copulation.attempts.model= lmer(log1p(cop_attempts)~female_lineage+male_lineage+female_lineage*male_lineage+(1|female_ID)+(1|male_ID)+(1|pairing_number),data=female_male_data, REML = FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(log.lmer.copulation.attempts.model)
anova(log.lmer.copulation.attempts.model, type = 3)
confint(log.lmer.copulation.attempts.model)
lsmeansLT(log.lmer.copulation.attempts.model, test.effs = NULL)
difflsmeans(log.lmer.copulation.attempts.model, test.effs = NULL)
r.squaredGLMM(log.lmer.copulation.attempts.model)
```

&nbsp;

###*Model 2: Number of male head-bobs*
####Generalised linear mixed model with female lineage, male lineage, female lineage * male lineage as fixed effects and female ID, male ID and pairing number (female's first or second trial) as random effects.

```{r, message=FALSE, warning=FALSE}
lmer.male.headbob.model= lmer(male_headbob~female_lineage+male_lineage+female_lineage*male_lineage+(1|female_ID)+(1|male_ID)+(1|pairing_number),data=female_male_data, REML = FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```
####Check the normality of the residuals
```{r}
plot(lmer.male.headbob.model)
qqnorm(resid(lmer.male.headbob.model))
qqline(resid(lmer.male.headbob.model))
shapiro.test(resid(lmer.male.headbob.model))
```
####log transform the data to meet model assumptions of normality
```{r, message=FALSE, warning=FALSE}
log.lmer.male.headbob.model= lmer(log1p(male_headbob)~female_lineage+male_lineage+female_lineage*male_lineage+(1|female_ID)+(1|male_ID)+(1|pairing_number),data=female_male_data, REML = FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(log.lmer.male.headbob.model)
anova(log.lmer.male.headbob.model, type = 3)
confint(log.lmer.male.headbob.model)
lsmeansLT(log.lmer.male.headbob.model, test.effs = NULL)
difflsmeans(log.lmer.male.headbob.model, test.effs = NULL)
r.squaredGLMM(log.lmer.male.headbob.model)
```

&nbsp;

###*Model 3: Number of female head-bobs*
####Generalised linear mixed model with female lineage, male lineage, female lineage * male lineage as fixed effects and female ID, male ID and pairing number (female's first or second trial) as random effects.

```{r, message=FALSE, warning=FALSE}
lmer.female.headbob.model= lmer(female_headbob~female_lineage+male_lineage+female_lineage*male_lineage+(1|female_ID)+(1|male_ID)+(1|pairing_number),data=female_male_data, REML = FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```
####Check the normality of the residuals
```{r}
plot(lmer.female.headbob.model)
qqnorm(resid(lmer.female.headbob.model))
qqline(resid(lmer.female.headbob.model))
shapiro.test(resid(lmer.female.headbob.model))
```
####log transform the data to meet model assumptions of normality
```{r, message=FALSE, warning=FALSE}
log.lmer.female.headbob.model= lmer(log1p(female_headbob)~female_lineage+male_lineage+female_lineage*male_lineage+(1|female_ID)+(1|male_ID)+(1|pairing_number),data=female_male_data, REML = FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(log.lmer.female.headbob.model)
anova(log.lmer.female.headbob.model, type = 3)
confint(log.lmer.female.headbob.model)
lsmeansLT(log.lmer.female.headbob.model, test.effs = NULL)
difflsmeans(log.lmer.female.headbob.model, test.effs = NULL)
r.squaredGLMM(log.lmer.female.headbob.model)
```

&nbsp;

###*Model 4: Number of female rejection behaviours*
####Generalised linear mixed model with female lineage, male lineage, female lineage * male lineage as fixed effects and female ID, male ID and pairing number (female's first or second trial) as random effects.

```{r, message=FALSE, warning=FALSE}
lmer.female.rejection.model= lmer(female_reject~female_lineage+male_lineage+female_lineage*male_lineage+(1|female_ID)+(1|male_ID)+(1|pairing_number),data=female_male_data, REML = FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```
####Check the normality of the residuals
```{r}
plot(lmer.female.rejection.model)
qqnorm(resid(lmer.female.rejection.model))
qqline(resid(lmer.female.rejection.model))
shapiro.test(resid(lmer.female.rejection.model))
```
####log transform the data to meet model assumptions of normality
```{r, message=FALSE, warning=FALSE}
log.lmer.female.rejection.model= lmer(log1p(female_reject)~female_lineage+male_lineage+female_lineage*male_lineage+(1|female_ID)+(1|male_ID)+(1|pairing_number),data=female_male_data, REML = FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(log.lmer.female.rejection.model)
anova(log.lmer.female.rejection.model, type = 3)
confint(log.lmer.female.rejection.model)
lsmeansLT(log.lmer.female.rejection.model, test.effs = NULL)
difflsmeans(log.lmer.female.rejection.model, test.effs = NULL)
r.squaredGLMM(log.lmer.female.rejection.model)
```
&nbsp;

##**Male-male behavioural trials**

```{r, echo=FALSE}
male_male_data <- read.csv("male_male_data.csv", header=TRUE, sep=",", dec=".", strip.white=TRUE)
male_male_data$fm_order<-factor(male_male_data$fm_order)
```
####The dataset consists of `r nrow(male_male_data)` behavioural trials involving 26 male tawny dragons (*Ctenophorus decresii*).

```{r, echo=FALSE}
options(contrasts = c("contr.sum","contr.poly"))
```

&nbsp;

###*Model 1: Focal male latency to emerge*
####Generalised linear mixed model with focal male behavioural group, opponent male behavioural group, focal male behavioural group * opponent male behavioural group, focal male condition and opponent male condition as fixed effects and focal male ID and focal male trial number as random effects.

```{r, message=FALSE, warning=FALSE}
lmer.latency.model= lmer(latency_min~fm_behav+om_behav+fm_behav*om_behav+fm_condition+om_condition+(1|fm_ID)+(1|fm_order),data=male_male_data, REML = FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```
####Check the normality of the residuals
```{r}
plot(lmer.latency.model)
qqnorm(resid(lmer.latency.model))
qqline(resid(lmer.latency.model))
shapiro.test(resid(lmer.latency.model))
```
```{r, message=FALSE, warning=FALSE}
summary(lmer.latency.model)
anova(lmer.latency.model, type = 3)
lsmeansLT(lmer.latency.model, test.effs = NULL)
difflsmeans(lmer.latency.model, test.effs = NULL)
r.squaredGLMM(lmer.latency.model)
```

&nbsp;

###*Model 2: Focal male display behaviour*
####Generalised linear mixed model with focal male behavioural group, opponent male behavioural group, focal male behavioural group * opponent male behavioural group, focal male condition and opponent male condition as fixed effects and focal male ID and focal male trial number as random effects.

```{r, message=FALSE, warning=FALSE}
lmer.display.model= lmer(display_num_prop~fm_behav+om_behav+fm_behav*om_behav+fm_condition+om_condition+(1|fm_ID)+(1|fm_order),data=male_male_data, REML = FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```
####Check the normality of the residuals
```{r}
plot(lmer.display.model)
qqnorm(resid(lmer.display.model))
qqline(resid(lmer.display.model))
shapiro.test(resid(lmer.display.model))
```
####log transform the data to meet model assumptions of normality
```{r, message=FALSE, warning=FALSE}
log.lmer.display.model= lmer(log1p(display_num_prop)~fm_behav+om_behav+fm_behav*om_behav+fm_condition+om_condition+(1|fm_ID)+(1|fm_order),data=male_male_data, REML = FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(log.lmer.display.model)
anova(log.lmer.display.model, type = 3)
lsmeansLT(log.lmer.display.model, test.effs = NULL)
difflsmeans(log.lmer.display.model, test.effs = NULL)
r.squaredGLMM(log.lmer.display.model)
```

&nbsp;

###*Model 3: Focal male physical aggression*
####Generalised linear mixed model with focal male behavioural group, opponent male behavioural group, focal male behavioural group * opponent male behavioural group, focal male condition and opponent male condition as fixed effects and focal male ID and focal male trial number as random effects.

```{r, message=FALSE, warning=FALSE}
lmer.aggression.model= lmer(agg_dur_prop~fm_behav+om_behav+fm_behav*om_behav+fm_condition+om_condition+(1|fm_ID)+(1|fm_order),data=male_male_data, REML = FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
```
####Check the normality of the residuals
```{r}
plot(lmer.aggression.model)
qqnorm(resid(lmer.aggression.model))
qqline(resid(lmer.aggression.model))
shapiro.test(resid(lmer.aggression.model))
```
####log transform the data to meet model assumptions of normality
```{r, message=FALSE, warning=FALSE}
log.lmer.aggression.model= lmer(log1p(agg_dur_prop)~fm_behav+om_behav+fm_behav*om_behav+fm_condition+om_condition+(1|fm_ID)+(1|fm_order),data=male_male_data, REML = FALSE, control = lmerControl(optimizer ="Nelder_Mead"))
summary(log.lmer.aggression.model)
anova(log.lmer.aggression.model, type = 3)
lsmeansLT(log.lmer.aggression.model, test.effs = NULL)
difflsmeans(log.lmer.aggression.model, test.effs = NULL)
r.squaredGLMM(log.lmer.aggression.model)
```