---
title: "Acoustic property summary statistics"
author: "James Tumulty"
output: html_document
---

Analysis of individual variation in acoustic properties of golden rocket frogs and Kai rocket frogs 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)
```

## Load Data

```{r, message=FALSE, warning=FALSE}

set.seed(42)

### Load Packages ###
library(ggplot2)
library(MASS)
library(lme4)
library(car)
library(lmtest)



### beebei ###
beebei <- read.csv("beebei_call_properties.csv")

print("beebei counts of calls with different numbers of pulses")
table(beebei$Num_Pulses)



#Remove variables related to pulses 4, 5 and 6 because of so much missing data
p56 <- c("Pulse4_Interval", "Pulse4_Period","Pulse5_Duration", "Pulse5_Rise", "Pulse5_Fall",
              "Pulse5_HalfRise", "Pulse5_HalfFall", "Pulse5_Interval", "Pulse5_Period", "Pulse6_Duration",
              "Pulse6_Rise", "Pulse6_Fall", "Pulse6_HalfRise", "Pulse6_HalfFall",
              "Pulse5_Max_Freq", "Pulse6_Max_Freq")

beebei <- beebei[,!(names(beebei) %in% p56)]

print(paste("No. beebei individuals (full dataset) =", length(levels(beebei$ID))))
print(paste("No. beebei calls (full dataset) =", nrow(beebei)))



### kaiei ###
kaiei <- read.csv("kaiei_call_properties.csv")


print("kaiei counts of calls with different numbers of pulses")
table(kaiei$Num_Pulses)

print(paste("No. kaiei individuals (full dataset) =", length(levels(kaiei$ID))))
print(paste("No. kaiei calls (full dataset) =", nrow(kaiei)))


```



## Create temperature corrected values

Adjust to mean temperatures of 23.6 for beebei and 22.7 for kaiei

```{r}
### beebei ###

#New data frame to hold adjusted values
beebeitemp <- beebei[,1:4]

meanbeebeitemp <- mean(beebei$Temp)

print(paste("mean beebei temp =", round(meanbeebeitemp, 1)))

#Run the loop
for (i in 5:ncol(beebei)){ # for each variable,
  fit <- lm(beebei[,i] ~ Temp,  data=beebei, na.action=na.exclude) #Fit lm with temp
  intercept <- coefficients(fit)[[1]] #extract the intercept
  slope <- coefficients(fit)[[2]] #extract the slope
  pred <- intercept + slope*meanbeebeitemp #predict value mean temp
  resids <- resid(fit, na.action=na.exclude) #get model residuals, exclude NAs
  beebeitemp[,i] <- pred + resids #Adjust residuals by predicted value into new data frame
}

#Get the right column names
colnames(beebeitemp) <- colnames(beebei)


### kaiei ### 

#New data frame to hold adjusted values
kaitemp <- kaiei[,1:4]

meankaitemp <- mean(kaiei$Temp)

print(paste("mean kaiei temp =", round(meankaitemp, 1))) #Mean = 22.7 degrees

#Run the loop
for (i in 5:ncol(kaiei)){ # for each variable,
  fit <- lm(kaiei[,i] ~ Temp, data=kaiei, na.action=na.exclude) #Fit lm with temp
  intercept <- coefficients(fit)[[1]] #extract the intercept
  slope <- coefficients(fit)[[2]] #extract the slope
  pred <- intercept + slope*meankaitemp #predict value at mean temp
  resids <- resid(fit, na.action=na.exclude) #get model residuals, exclude NAs
  kaitemp[,i] <- pred + resids #Adjust residuals by predicted value into new data frame
}


#Get the right column names
colnames(kaitemp) <- colnames(kaiei)


```



## beebei Summary Statistics

Compute summary statistics on temperature corrected values for beebei data, output is rounded to 3 decimal places

```{r, warning=FALSE}

### beebei ###

#New data frame to store output
outmat <- array(NA, dim=c(8, ncol(beebeitemp)))
colnames(outmat) <- colnames(beebeitemp)
rownames(outmat) <- c("Grand Mean", "Grand SD", "CVa", "Mean CVw", "Min CVw",
                      "Max CVw", "PIC",  "etasq")
outmat <- as.data.frame(outmat)

datatemp <- beebeitemp

#Run loop!
for (i in 5:ncol(datatemp)){
  x <- tapply(datatemp[,i], datatemp$ID, mean, na.rm=TRUE)
  sd <- tapply(datatemp[,i], datatemp$ID, sd, na.rm=TRUE)
  CVw <- (sd/x)*100
  MeanCVw <- mean(CVw, na.rm=TRUE)
  GrandMean <- mean(x, na.rm=TRUE)
  GrandSD <- sd(x, na.rm=TRUE)
  CVa <- (GrandSD/GrandMean)*100
  outmat["Mean CVw", i] <- MeanCVw
  outmat["Min CVw", i] <- min(CVw, na.rm=TRUE)
  outmat["Max CVw", i] <- max(CVw, na.rm=TRUE)
  outmat["Grand Mean", i] <- GrandMean
  outmat["Grand SD", i] <- GrandSD
  outmat["CVa", i] <- CVa
  outmat["PIC", i] <- CVa/MeanCVw
  fit <- aov(datatemp[,i] ~ ID, data=datatemp)
  outmat["etasq", i] <- summary.lm(fit)$r.squared
}

#Remove 1st 4 columns and transpose
outmat <- outmat[, 5:ncol(outmat)]
outmat <- t(outmat)

#Write CSV file
#write.csv(outmat, file="beebei_summary_stats.csv")

round(outmat, 3)



```


## kaiei Summary Statistics

Compute summary statistics on temperature corrected values for kaiei data, output is rounded to 3 decimal places

```{r, warning=FALSE}
### kaiei ###

#New data frame to store output
outmat <- array(NA, dim=c(8, ncol(kaitemp)))
colnames(outmat) <- colnames(kaitemp)
rownames(outmat) <- c("Grand Mean", "Grand SD", "CVa", "Mean CVw", "Min CVw",
                      "Max CVw", "PIC",  "etasq")
outmat <- as.data.frame(outmat)

datatemp <- kaitemp

#Run loop!
for (i in 5:ncol(datatemp)){
  x <- tapply(datatemp[,i], datatemp$ID, mean, na.rm=TRUE)
  sd <- tapply(datatemp[,i], datatemp$ID, sd, na.rm=TRUE)
  CVw <- (sd/x)*100
  MeanCVw <- mean(CVw, na.rm=TRUE)
  GrandMean <- mean(x, na.rm=TRUE)
  GrandSD <- sd(x, na.rm=TRUE)
  CVa <- (GrandSD/GrandMean)*100
  outmat["Mean CVw", i] <- MeanCVw
  outmat["Min CVw", i] <- min(CVw, na.rm=TRUE)
  outmat["Max CVw", i] <- max(CVw, na.rm=TRUE)
  outmat["Grand Mean", i] <- GrandMean
  outmat["Grand SD", i] <- GrandSD
  outmat["CVa", i] <- CVa
  outmat["PIC", i] <- CVa/MeanCVw
  fit <- aov(datatemp[,i] ~ ID, data=datatemp)
  outmat["etasq", i] <- summary.lm(fit)$r.squared
}

#Remove 1st 4 columns and transpose
outmat <- outmat[, 5:ncol(outmat)]
outmat <- t(outmat)

#write.csv(outmat, file="kaiei summary stats.csv")

round(outmat, 3)

```


