---
title: "Habituation Discrimination data analysis"
author: "James Tumulty"
date: "4/23/2021"
output: html_document
---

Analysis of aggressive responses of male golden rocket frogs (Anomaloglossus beebei) and male Kai rocket frogs (Anomaloglossus kaiei) to habituation-discrimination acoustic playback experiments in which half of the frogs experienced a "Different frog" experimental treatment in the discrimination phase, and the other half experienced a "Same frog" control treatment. For paper: 

Tumulty, Lange, Bee. Identity signals, identity reception, and the evolution of social recognition in a Neotropical frog


```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```


<br>
<br>

# Load and set up data

```{r, message=FALSE}

library(ggplot2)
library(lme4)
library(car)
library(lmtest)

#Load data
data <- read.csv("Habituation_discrimination_data.csv")


#Set up correct order for factors Phase, Block, and Treatment
data$Phase <- as.character(data$Phase)
data$Phase<- factor(data$Phase, levels=unique(data$Phase))

data$Block <- as.character(data$Block)
data$Block <- factor(data$Block, levels=unique(data$Block))

data$Treatment <- factor(data$Treatment, levels=c("Same_frog", "Different_frog"))


```

<br>
<br>

# Aggression Index (PCA)

Compute aggression index for each species as first principal component from PCA of 4 measures of aggression. 


### beebei aggression index

```{r}

beebeiHD <- data[data$Species=="beebei",]

beebeiPC <- prcomp(~ Agg_Calls + Pseudo_Agg + Approach + Prop_Dist, data=beebeiHD, center=TRUE, scale=TRUE)

summary(beebeiPC)
beebeiPC$rotation

biplot(beebeiPC)

#Need to invert PC1 so that higher values = more aggressive

#Add PC1 to beebei data frame
beebeiHD$PC1 <- 1-beebeiPC$x[,1]

```

<br>

### kaiei aggression index

```{r}
kaieiHD <- data[data$Species=="kaiei",]

kaieiPC <- prcomp(~ Agg_Calls + Pseudo_Agg + Approach + Prop_Dist, data=kaieiHD, center=TRUE, scale=TRUE)

summary(kaieiPC)
kaieiPC$rotation


biplot(kaieiPC)

#Add PC1 to kaiei data frame
kaieiHD$PC1 <- kaieiPC$x[,1]

```


<br>

### beebei Habituation-Discrimination plot

```{r}
ggplot(beebeiHD, aes(x=Block, y=PC1, fill=Treatment))+
  geom_boxplot(outlier.shape=NA)+
  scale_fill_manual(values=c("white", "gray"), labels=c("Same frog", "Different frog"))+
  coord_cartesian(ylim = c(-1,4))+
  scale_x_discrete(labels=c("First \n Habituation \n Block", "Last \n Habituation \n Block", 
                            "Discrimination \n Phase"))+
  ylab("Aggression index")+
  theme_bw()+
  theme(panel.grid=element_blank(), legend.title = element_blank(), 
        axis.title.x = element_blank())+
  theme(strip.background = element_blank())

```

<br>

### kaiei Habituation-Discrimination plot

```{r}

ggplot(kaieiHD, aes(x=Block, y=PC1, fill=Treatment))+
  geom_boxplot(outlier.shape=NA)+
  scale_fill_manual(values=c("white", "gray"), labels=c("Same frog", "Different frog"))+
  scale_x_discrete(labels=c("First \n Habituation \n Block", "Last \n Habituation \n Block", 
                            "Discrimination \n Phase"))+
  ylab("Aggression index")+
  coord_cartesian(ylim = c(-1.5,6))+
  theme_bw()+
  theme(panel.grid=element_blank(), legend.title = element_blank(), 
        axis.title.x = element_blank())+
  theme(strip.background = element_blank())
```

<br>

### Within-individual recovery of aggression plot

```{r}
#Remake data dataframe, now with PC1s
data <- rbind(beebeiHD, kaieiHD)

#New wide dataframe to compute within-individual differences
data_wide <- reshape(data[,names(data) %in% c("ID", "Treatment", "Block", "Species", "PC1")],
                         v.names=c("PC1"), timevar="Block", idvar="ID", direction="wide")

data_wide$zAgg <- data_wide$PC1.Discrimination - data_wide$PC1.Last_H

#Plot difference as boxplots
ggplot(data_wide, aes(x=Treatment, y=zAgg, fill=Treatment))+
  geom_boxplot(outlier.shape=NA)+
  facet_grid(~Species)+
  geom_hline(yintercept=0, lty=2)+
  scale_fill_manual(values=c("white", "gray"), labels=c("Same frog", "Different frog"))+
  scale_x_discrete(labels=c("Control", "Treatment"))+
  ylab("Difference in Aggression Index")+
  theme_bw()+
  theme(panel.grid = element_blank(), legend.title = element_blank(), 
        axis.title.x = element_blank())


```


<br>
<br>


## Analysis of decrease in aggression during Habituation Phase

This compares  the first block of habituation phase with the last block. Treatments are pooled together because treatment differences haven't been applied yet (they are applied during the discrimination phase). Non-parametric paired tests are used because data do not meet assumptions of parametric tests. 

### beebei

```{r}
beebei_wide <- data_wide[data_wide$Species=="beebei",]

wilcox.test(beebei_wide$PC1.First_H, beebei_wide$PC1.Last_H, paired=TRUE)

```

### kaiei

```{r}
kaiei_wide <- data_wide[data_wide$Species=="kaiei",]

wilcox.test(kaiei_wide$PC1.First_H, kaiei_wide$PC1.Last_H, paired=TRUE)
```

<br>
<br>


## Analysis of recovery of aggression during Discrimination Phase by Treatment - beebei

Plot what is being analyzed, connect points by individual

```{r}
ggplot(beebeiHD[beebeiHD$Block %in% c("Last_H", "Discrimination"),], 
       aes(x=Block, y=PC1))+
  facet_wrap(~Treatment)+
  geom_point(alpha=0.5)+
  geom_line(aes(group=ID), alpha=0.5)+
  scale_x_discrete(labels=c("Last \n Habituation \n Block", "Discrimination \n Phase"))+
  ylab("Aggression Index")+
  theme_bw()+
  theme(panel.grid=element_blank(), axis.title.x = element_blank(), 
        legend.title = element_blank())
```



### Linear mixed effects model with interaction term - beebei

Fit linear mixed effets model and look at model diagnostics

Here we are testing by a Treatment by Block interaction in aggression. Block refers to last block of habituation phase vs. Discrimination phase.

Model diagnostics look good. Residuals have pretty even variance across range of fitted values, qqplot is mostly linear indicating normally distributed residuals (other than perhaps one slight outlier in 21)

```{r}
#Linear mixed effects model
lm1 <- lmer(PC1 ~ Treatment*Block + (1|ID), 
            data=beebeiHD[beebeiHD$Block %in% c("Last_H", "Discrimination"),], REML=FALSE)

plot(lm1) #Looks ok
qqp(residuals(lm1)) #Also okay
```


Likelihood ratio test of significance of interaction term

```{r}
lrtest(lm1, update(lm1, .~. -Treatment:Block))
```


<br>

### Within treatment recovery of aggression - beebei Different frog Treatment

Fit linear mixed effects model of recovery of aggression within the different frog treatment for beebei

Evaluate significance with likelihood ratio test

```{r}

#Linear mixed effects model
lmchange <- lmer(PC1 ~ Block + (1|ID), 
                 data=beebeiHD[beebeiHD$Block %in% c("Last_H", "Discrimination") & 
                                 beebeiHD$Treatment=="Different_frog",], REML=FALSE)

lrtest(lmchange, update(lmchange, ~. -Block))
```

<br>

### Within treatment recovery of aggression - beebei Same frog control

Fit linear mixed effects model of recovery of aggression within the same frog control treatment for beebei

evaluate significance with likelihood ratio test


```{r}
#Linear mixed effects model
lmcont <- lmer(PC1 ~ Block + (1|ID), 
                 data=beebeiHD[beebeiHD$Block %in% c("Last_H", "Discrimination") & 
                                 beebeiHD$Treatment=="Same_frog",], REML=FALSE)

lrtest(lmcont, update(lmcont, ~. -Block))
```

<br>
<br>

## Analysis of recovery of aggression during Discrimination Phase by Treatment - kaiei

Plot what is being analyzed, points connected by individual

```{r}
ggplot(kaieiHD[kaieiHD$Block %in% c("Last_H", "Discrimination"),], 
       aes(x=Block, y=PC1))+
  facet_wrap(~Treatment)+
  geom_point(alpha=0.5)+
  geom_line(aes(group=ID), alpha=0.5)+
  scale_x_discrete(labels=c("Last \n Habituation \n Block", "Discrimination \n Phase"))+
  ylab("Aggression Index")+
  theme_bw()+
  theme(panel.grid=element_blank(), axis.title.x = element_blank(), 
        legend.title = element_blank())
```



Here again we are testing by a Treatment by Block interaction in aggression. 

Model diagnostics look good:

```{r}
#Linear mixed effects model
lm1 <- lmer(PC1 ~ Treatment*Block + (1|ID), 
            data=kaieiHD[kaieiHD$Block %in% c("Last_H", "Discrimination"),], REML=FALSE)

plot(lm1) #Looks good
qqp(residuals(lm1)) #Also good

```


Likelihood ratio test of interaction term

```{r}
lrtest(lm1, update(lm1, ~. -Treatment:Block))
```

<br>

### Within treatment recovery of aggression - kaiei Different frog Treatment

Fit linear mixed effects model of recovery of aggression within Different frog treatment for kaiei, evaluate significance with likelihood ratio test

```{r}
#Linear mixed effects model
lmchange <- lmer(PC1 ~ Block + (1|ID), 
                 data=kaieiHD[kaieiHD$Block %in% c("Last_H", "Discrimination") & 
                                 kaieiHD$Treatment=="Different_frog",], REML=FALSE)

lrtest(lmchange, update(lmchange, ~. -Block))
```


<br>

### Within treatment recovery of aggression - kaiei Same frog Control

Fit linear mixed effects model of recovery of aggression within Same frog control treatment for kaiei, evaluate significance with likelihood ratio test


```{r}
#Linear mixed effects model
lmcont <- lmer(PC1 ~ Block + (1|ID), 
                 data=kaieiHD[kaieiHD$Block %in% c("Last_H", "Discrimination") & 
                                 kaieiHD$Treatment=="Same_frog",], REML=FALSE)

lrtest(lmcont, update(lmcont, ~. -Block))
```





