---
title: Struthiolariid case-study for "Going round the twist – an empirical analysis
  of shell coiling in helicospiral gastropods"
author: "Katie S. Collins, Roman Klapaukh, James S. Crampton, Michael F. Gazley, C. Ian Schipper, A. Maksimenko, and Benjamin R. Hines"
date: "January 2021"
output:
  html_document:
    df_print: paged
  pdf_document: default
  word_document: default
---

To demonstrate the uses of the coiling and shape parameters, we subset out the Struthiolariidae from the main dataset.

This workflow is identical to that performed in **SOM_process_snails.Rmd** - that file should be used to run new datasets. You also need to run that file first to provide inputs for running this file.

```{r loadPackages, cache=F, warning=FALSE, message=FALSE}
library(magrittr)
library(dplyr)
library(stringr)
library(devtools)
library(data.table)
library(ggplot2)
library(reshape2)
library(scatterplot3d)
library(DT,quietly=T)
library(threejs)
library(gridExtra)
library(ggpubr)
#library(DiagrammeR)
#library(metricsgraphics)
#library(scatterD3)
library(jpeg)
library(pander)
library(magick)
library(plot3D)
library(tidyverse)
```


```{r readData}
snails <- read.csv(file="raw_snail_parameters.csv", header=T)

snails <- snails %>%
  filter(species %in% c("Struthiolaria_papulosa, Recent", "Conchothyra_australis, Wangaloan", 
                        "Conchothyra_parasitica, Haumurian", "Pelicaria_vermis, Castlecliffian", 
                        "Pelicaria_vermis, Nukumaruan", "Perissodonta_ellioti, Eocene-Oligocene",
                        "Perissodonta_nordenskjoldi, Eocene", "Struthiolaria_calcar, Waitakian", 
                        "Tylospira_coronata, Kalimnan", "Tylospira_coronata, Miocene", "Tylospira_coronata, Pliocene", 
                        "Tylospira_gilli, Pliocene", "Tylospira_glomerata, Miocene", "Monalaria_concinna, Bortonian",
                        "Tylospira_scutulata, Recent"))

```

```{r AIC_function}
AICc <- function(mod){
 n <- nrow(mod$model)
 k <- extractAIC(mod)[1]
 AIC <- extractAIC(mod)[2]
 AICc <- AIC + 2*k*(k+1)/(n-k-1)
 return(AICc)
}
```


```{r processSnailsT}

modelFittingT = snails %>%
  group_by(id, species) %>%
  mutate(num = n()) %>%
  filter(num >=3 ) %>%
  group_by(id,species)%>%
  do(linearModel = lm(y ~ theta, data = .),
     quadModel   = lm(y ~ poly(theta,2), data = . ),
     logModel    = lm(y ~ log(theta+0.0001), data = .), 
     powerModel  = lm(log(y+0.00001) ~ log(theta + 0.00001), data = . ),
     expModel    = lm(log(y+0.0001) ~ theta, data = (. ))) %>%
  apply(1, function(row){
           linear = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Linear",
                     Residual = residuals(row$linearModel),
                     Predicted = predict(row$linearModel),
                     Rsq = summary(row$linearModel)$adj.r.squared,
                     Rate = row$linearModel$coefficients[2],
                     AICc = AICc(row$linearModel))

           quad = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Polynomial",
                     Residual = residuals(row$quadModel),
                     Predicted = predict(row$quadModel),
                     Rsq = summary(row$quadModel)$adj.r.squared,
                     Rate = row$quadModel$coefficients[2],
                     AICc = AICc(row$quadModel))
           
           log = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Log",
                     Residual = residuals(row$logModel),
                     Predicted = predict(row$logModel),
                     Rsq = summary(row$logModel)$adj.r.squared,
                     Rate = row$logModel$coefficients[2],
                     AICc = AICc(row$logModel))

           power = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Power",
                     Residual = residuals(row$powerModel),
                     Predicted = predict(row$powerModel),
                     Rsq = summary(row$powerModel)$adj.r.squared,
                     Rate = row$powerModel$coefficients[2],
                     AICc = AICc(row$powerModel))

           exp = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Exponential",
                     Residual = residuals(row$expModel),
                     Predicted = predict(row$expModel),
                     Rsq = summary(row$expModel)$adj.r.squared,
                     Rate = row$expModel$coefficients[2],
                     AICc = AICc(row$expModel))

           return(rbind(linear,log,power,exp, quad))
       }) %>%
  rbindlist%>%
  group_by(ID, Species) %>%
  mutate(relativeLikelihood = exp((min(AICc) - AICc)/2)) %>%
  ungroup()

castTModels <- select(modelFittingT, Model, Species, ID, relativeLikelihood)%>%
  mutate(relativeLikelihood=round(relativeLikelihood*1000, 0))%>%
  unique()%>%
  dcast(Species+ID ~ Model, value.var="relativeLikelihood")

Tresiduals <- modelFittingT %>%
  ggplot(aes(Predicted, Residual, group=ID)) +
    geom_hline(yintercept=0,colour="purple") +
    #geom_point(aes(colour=Species)) + 
    geom_line(alpha=0.3, aes(colour=Species)) +
    facet_wrap(~ Model, nrow=1) +
    theme(legend.position="none")+ 
    theme(aspect.ratio = 1)

Tr2histo <- modelFittingT %>%
   ggplot(aes(Rsq)) +
   geom_histogram(binwidth=0.05) +
   geom_vline(aes(xintercept = meanR2), data = modelFittingT %>%
              group_by(Model) %>% summarise(meanR2 = mean(Rsq)),
            colour = "blue") +
   geom_text(aes(x=0, y=0, label=sprintf("%0.2f", meanR2)), data = modelFittingT %>%
              group_by(Model) %>% summarise(meanR2 = mean(Rsq)),
            colour = "blue",vjust="bottom", hjust="left") +
    facet_grid(~Model,scales="free_y")+ 
    theme(aspect.ratio = 1)

grid.arrange(Tresiduals, Tr2histo, nrow=2)

modelFittingT %>%
  group_by(Model) %>%
  summarise(`$Mean R^2$` = mean(Rsq)) %>%
  arrange(`$Mean R^2$`) %>%
  as.data.frame %>%
  pander

modelFittingT %>%
   ggplot(aes(relativeLikelihood, fill=Model)) +
   geom_histogram(position="identity")+
  theme(legend.position="bottom")

## The POWER-LAW model is best for this set of Struthiolariidae
TModel = modelFittingT 

```

```{r processSnailsW}

modelFittingW = snails %>%
  group_by(id, species) %>%
  mutate(num = n()) %>%
  filter(num >=3 ) %>%
  group_by(id,species)%>%
  do(linearModel = lm(ApertureArea ~ theta, data = .),
     quadModel   = lm(ApertureArea ~ poly(theta,2), data = . ),
     logModel    = lm(ApertureArea ~ log(theta+0.0001), data = .), 
     powerModel  = lm(log(ApertureArea+0.00001) ~ log(theta + 0.00001), data = . ),
     expModel    = lm(log(ApertureArea+0.0001) ~ theta, data = (. ))) %>%
  apply(1, function(row){
           linear = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Linear",
                     Residual = residuals(row$linearModel),
                     Predicted = predict(row$linearModel),
                     Rsq = summary(row$linearModel)$adj.r.squared,
                     Rate = row$linearModel$coefficients[2],
                     AICc = AICc(row$linearModel))

           quad = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Polynomial",
                     Residual = residuals(row$quadModel),
                     Predicted = predict(row$quadModel),
                     Rsq = summary(row$quadModel)$adj.r.squared,
                     Rate = row$quadModel$coefficients[2],
                     AICc = AICc(row$quadModel))
           
           log = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Log",
                     Residual = residuals(row$logModel),
                     Predicted = predict(row$logModel),
                     Rsq = summary(row$logModel)$adj.r.squared,
                     Rate = row$logModel$coefficients[2],
                     AICc = AICc(row$logModel))

           power = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Power",
                     Residual = residuals(row$powerModel),
                     Predicted = predict(row$powerModel),
                     Rsq = summary(row$powerModel)$adj.r.squared,
                     Rate = row$powerModel$coefficients[2],
                     AICc = AICc(row$powerModel))

           exp = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Exponential",
                     Residual = residuals(row$expModel),
                     Predicted = predict(row$expModel),
                     Rsq = summary(row$expModel)$adj.r.squared,
                     Rate = row$expModel$coefficients[2],
                     AICc = AICc(row$expModel))

           return(rbind(linear,log,power,exp, quad))
       }) %>%
  rbindlist%>%
  group_by(ID, Species) %>%
  mutate(relativeLikelihood = exp((min(AICc) - AICc)/2)) %>%
  ungroup()

castWModels <- select(modelFittingW, Model, Species, ID, relativeLikelihood)%>%
  mutate(relativeLikelihood=round(relativeLikelihood*1000, 0))%>%
  unique()%>%
  dcast(Species+ID ~ Model, value.var="relativeLikelihood")

Wresiduals <- modelFittingW %>%
  ggplot(aes(Predicted, Residual, group=ID)) +
    geom_hline(yintercept=0,colour="purple") +
    #geom_point(aes(colour=Species)) + 
    geom_line(alpha=0.3, aes(colour=Species)) +
    facet_wrap(~ Model, nrow=1) +
    theme(legend.position="none")+ 
    theme(aspect.ratio = 1)

Wr2histo <- modelFittingW %>%
   ggplot(aes(Rsq)) +
   geom_histogram(binwidth=0.05) +
   geom_vline(aes(xintercept = meanR2), data = modelFittingW %>%
              group_by(Model) %>% summarise(meanR2 = mean(Rsq)),
            colour = "blue") +
   geom_text(aes(x=0, y=0, label=sprintf("%0.2f", meanR2)), data = modelFittingW %>%
              group_by(Model) %>% summarise(meanR2 = mean(Rsq)),
            colour = "blue",vjust="bottom", hjust="left") +
    facet_grid(~Model,scales="free_y")+ 
    theme(aspect.ratio = 1)

grid.arrange(Wresiduals, Wr2histo, nrow=2)

modelFittingW %>%
  group_by(Model) %>%
  summarise(`$Mean R^2$` = mean(Rsq)) %>%
  arrange(`$Mean R^2$`) %>%
  as.data.frame %>%
  pander

modelFittingW %>%
   ggplot(aes(relativeLikelihood, fill=Model)) +
   geom_histogram(position="identity")+
  theme(legend.position="bottom")


## the best W model is EXPONENTIAL
WModel = modelFittingW
```

```{r processSnailsD}
modelFittingD = snails %>%
  group_by(id, species) %>%
  mutate(num = n()) %>%
  filter(num >=3 ) %>%
  group_by(id,species)%>%
  do(linearModel = lm(center_gap ~ theta, data = .),
     quadModel   = lm(center_gap ~ poly(theta,2), data = . ),
     logModel    = lm(center_gap ~ log(theta+0.0001), data = .), 
     powerModel  = lm(log(center_gap+0.00001) ~ log(theta + 0.00001), data = . ),
     expModel    = lm(log(center_gap+0.0001) ~ theta, data = (. ))) %>%
  apply(1, function(row){
           linear = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Linear",
                     Residual = residuals(row$linearModel),
                     Predicted = predict(row$linearModel),
                     Rsq = summary(row$linearModel)$adj.r.squared,
                     Rate = row$linearModel$coefficients[2],
                     AICc = AICc(row$linearModel))

           quad = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Polynomial",
                     Residual = residuals(row$quadModel),
                     Predicted = predict(row$quadModel),
                     Rsq = summary(row$quadModel)$adj.r.squared,
                     Rate = row$quadModel$coefficients[2],
                     AICc = AICc(row$quadModel))
           
           log = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Log",
                     Residual = residuals(row$logModel),
                     Predicted = predict(row$logModel),
                     Rsq = summary(row$logModel)$adj.r.squared,
                     Rate = row$logModel$coefficients[2],
                     AICc = AICc(row$logModel))

           power = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Power",
                     Residual = residuals(row$powerModel),
                     Predicted = predict(row$powerModel),
                     Rsq = summary(row$powerModel)$adj.r.squared,
                     Rate = row$powerModel$coefficients[2],
                     AICc = AICc(row$powerModel))

           exp = data.table(ID= row$id,
                     Species = row$species,
                     Model = "Exponential",
                     Residual = residuals(row$expModel),
                     Predicted = predict(row$expModel),
                     Rsq = summary(row$expModel)$adj.r.squared,
                     Rate = row$expModel$coefficients[2],
                     AICc = AICc(row$expModel))

           return(rbind(linear,log,power,exp, quad))
       }) %>%
  rbindlist%>%
  group_by(ID, Species) %>%
  mutate(relativeLikelihood = exp((min(AICc) - AICc)/2)) %>%
  ungroup()

castDModels <- select(modelFittingD, Model, Species, ID, relativeLikelihood)%>%
  mutate(relativeLikelihood=round(relativeLikelihood*1000, 0))%>%
  unique()%>%
  dcast(Species+ID ~ Model, value.var="relativeLikelihood")

Dresiduals <- modelFittingD %>%
  ggplot(aes(Predicted, Residual, group=ID)) +
    geom_hline(yintercept=0,colour="purple") +
    #geom_point(aes(colour=Species)) + 
    geom_line(alpha=0.3, aes(colour=Species)) +
    facet_wrap(~ Model, nrow=1) +
    theme(legend.position="none")+ 
    theme(aspect.ratio = 1)

Dr2histo <- modelFittingD %>%
   ggplot(aes(Rsq)) +
   geom_histogram(binwidth=0.05) +
   geom_vline(aes(xintercept = meanR2), data = modelFittingD %>%
              group_by(Model) %>% summarise(meanR2 = mean(Rsq)),
            colour = "blue") +
   geom_text(aes(x=0, y=0, label=sprintf("%0.2f", meanR2)), data = modelFittingD %>%
              group_by(Model) %>% summarise(meanR2 = mean(Rsq)),
            colour = "blue",vjust="bottom", hjust="left") +
    facet_grid(~Model,scales="free_y")+ 
    theme(aspect.ratio = 1)

grid.arrange(Dresiduals, Dr2histo, nrow=2)

modelFittingD %>%
  group_by(Model) %>%
  summarise(`$Mean R^2$` = mean(Rsq)) %>%
  arrange(`$Mean R^2$`) %>%
  as.data.frame %>%
  pander

modelFittingD %>%
   ggplot(aes(relativeLikelihood, fill=Model)) +
   geom_histogram(position="identity")+
  theme(legend.position="bottom")

## best model for D is POLYNOMIAL
DModel = modelFittingD 

```

```{r hangleRead}
hangle <- read.csv(file="S_Table_3_hangle_output.csv", header=T)%>%
  select(-ID)%>%
  dplyr::rename(ID=Species)
```


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

## we need to filter the Hangle data to just Struthiolariidae, so first we need to find which IDs belong to which genera

IDs <- snails %>%
  select(id, species)%>%
  separate(species, c("Genus"), extra='drop')%>%
  filter(Genus %in% c("Struthiolaria", "Pelicaria", "Conchothyra", "Monalaria", "Tylospira", "Perissodonta"))%>%
  pull(id)

## ALL apertures
apertures_full <- hangle %>%
  filter(!grepl("P135697a_", ID))%>%
  filter(!grepl("060_", ID))%>%
  dplyr::group_by(ID)%>%
  dplyr::mutate(n_ap = row_number())%>%
  ungroup()%>%
  filter(ID %in% IDs)

pcaData <- apertures_full %>%
  select(-ID) %>%
  as.matrix()%>%
  na.omit()

pca <- princomp(na.omit(pcaData))

## how much variance does each variable explain?
summary(pca)

## extract the PC scores of interest (in this case, Comp.1)
scores <- as.data.frame(pca$scores)%>%
  select(Comp.1)

## what's the min value of Comp.1? Subtract that from every element of Comp.1 (because it's a negative number)
minPC1 <- min(scores$Comp.1)

scores <- scores - minPC1

snails2 <- snails %>%
  dplyr::group_by(id)%>%
  dplyr::mutate(Aperture = row_number())
  
scores$ID <- apertures_full$ID
scores$Aperture <- apertures_full$Aperture
scores <- scores %>%
  left_join(snails2, by=c("ID"="id", "Aperture"))%>%
  select(-X, -Y, -center_gap_x, -r, -gap_2, -Degrees)%>%
  ungroup()%>%
  drop_na()
  
## now fit the models

modelFittingS = scores %>%
  group_by(ID, species) %>%
  mutate(num = n()) %>%
  group_by(ID, species)%>%
  do(linearModel = lm(Comp.1 ~ theta, data = .),
     quadModel   = lm(Comp.1 ~ poly(theta,2), data = . ),
     logModel    = lm(Comp.1 ~ log(theta+0.0001), data = .), 
     powerModel  = lm(log(Comp.1+0.00001) ~ log(theta + 0.00001), data = . ),
     expModel    = lm(log(Comp.1+0.0001) ~ theta, data = (. ))) %>%
  apply(1, function(row){
           linear = data.table(ID= row$ID,
                     Species = row$species,
                     Model = "Linear",
                     Residual = residuals(row$linearModel),
                     Predicted = predict(row$linearModel),
                     Rsq = summary(row$linearModel)$adj.r.squared,
                     Rate = row$linearModel$coefficients[2],
                     AICc = AICc(row$linearModel))

           quad = data.table(ID= row$ID,
                     Species = row$species,
                     Model = "Polynomial",
                     Residual = residuals(row$quadModel),
                     Predicted = predict(row$quadModel),
                     Rsq = summary(row$quadModel)$adj.r.squared,
                     Rate = row$quadModel$coefficients[2],
                     AICc = AICc(row$quadModel))
           
           log = data.table(ID= row$ID,
                     Species = row$species,
                     Model = "Log",
                     Residual = residuals(row$logModel),
                     Predicted = predict(row$logModel),
                     Rsq = summary(row$logModel)$adj.r.squared,
                     Rate = row$logModel$coefficients[2],
                     AICc = AICc(row$logModel))

           power = data.table(ID= row$ID,
                     Species = row$species,
                     Model = "Power",
                     Residual = residuals(row$powerModel),
                     Predicted = predict(row$powerModel),
                     Rsq = summary(row$powerModel)$adj.r.squared,
                     Rate = row$powerModel$coefficients[2],
                     AICc = AICc(row$powerModel))

           exp = data.table(ID= row$ID,
                     Species = row$species,
                     Model = "Exponential",
                     Residual = residuals(row$expModel),
                     Predicted = predict(row$expModel),
                     Rsq = summary(row$expModel)$adj.r.squared,
                     Rate = row$expModel$coefficients[2],
                     AICc = AICc(row$expModel))

           return(rbind(linear,log,power,exp, quad))
       }) %>%
  rbindlist%>%
  group_by(ID, Species) %>%
  mutate(relativeLikelihood = exp((min(AICc) - AICc)/2)) %>%
  ungroup()

castSModels <- select(modelFittingS, Model, Species, ID, relativeLikelihood)%>%
  mutate(relativeLikelihood=round(relativeLikelihood*1000, 0))%>%
  unique()%>%
  dcast(Species+ID ~ Model, value.var="relativeLikelihood")%>%
  write.csv(file="S1ModelsAIC.csv", row.names=FALSE)

```

The residual plots:

```{r AllS1ResidualsSum, warning=FALSE, message=FALSE}
Sresiduals <- modelFittingS %>%
  ggplot(aes(Predicted, Residual, group=ID)) +
    geom_hline(yintercept=0,colour="purple") +
    #geom_point(aes(colour=Species)) + 
    geom_line(alpha=0.3, aes(colour=Species)) +
    facet_wrap(~ Model, nrow=1) +
    theme(legend.position="none")+ 
    theme(aspect.ratio = 1)
```

```{r S1r2sum, warning=FALSE, message=FALSE}
Sr2histo <- modelFittingS %>%
   ggplot(aes(Rsq)) +
   geom_histogram(binwidth=0.05) +
   geom_vline(aes(xintercept = meanR2), data = modelFittingS %>%
              group_by(Model) %>% summarise(meanR2 = mean(Rsq)),
            colour = "blue") +
   geom_text(aes(x=0, y=0, label=sprintf("%0.2f", meanR2)), data = modelFittingS%>%
              group_by(Model) %>% summarise(meanR2 = mean(Rsq)),
            colour = "blue",vjust="bottom", hjust="left") +
    facet_grid(~Model,scales="free_y")+ 
    theme(aspect.ratio = 1)
```


```{r S1R2residuals, warning=FALSE, message=FALSE}
grid.arrange(Sresiduals, Sr2histo, nrow=2)
```

```{r S1R2SummaryTable, warning=FALSE, message=FALSE}
modelFittingS %>%
  group_by(Model) %>%
  summarise(`$Mean R^2$` = mean(Rsq)) %>%
  arrange(`$Mean R^2$`) %>%
  as.data.frame %>%
  pander
```
 
The relative likelihoods:

```{r S1AICc, warning=FALSE, message=FALSE}
modelFittingS %>%
   ggplot(aes(relativeLikelihood, fill=Model)) +
   geom_histogram(position="identity")+
  theme(legend.position="bottom")
```

```{r acceptedSModels, warning=FALSE, message=FALSE}
SModel = modelFittingS 
```


# Combine the "growth" parameters

```{r generateCombinedOutput}
combined <- TModel %>%
  group_by(ID,Species,Model) %>%
  summarise(Rate = unique(Rate)) %>%
  mutate(Model = paste0("T", Model)) %>%
  dcast(ID + Species ~ Model,value.var = "Rate") %>%
  left_join(
    WModel%>%
      group_by(ID, Species, Model) %>%
      summarise(Rate = unique(Rate)) %>%
      mutate(Model = paste0("W", Model)) %>%
      dcast(ID + Species ~ Model,value.var = "Rate"),
    by = c("ID","Species")) %>%
  left_join(
    DModel%>%
      group_by(ID, Species, Model) %>%
      summarise(Rate = unique(Rate)) %>%
      mutate(Model = paste0("D", Model)) %>%
      dcast(ID + Species ~ Model,value.var = "Rate"),
    by = c("ID","Species")) %>%
  left_join(
    SModel%>%
      group_by(ID, Species, Model) %>%
      summarise(Rate = unique(Rate)) %>%
      mutate(Model = paste0("S", Model)) %>%
      dcast(ID + Species ~ Model,value.var = "Rate"),
    by = c("ID","Species"))

```
We can visualise the fits of each model to the raw data:
```{r ModelCurvePlots, warning=FALSE, message=FALSE}

## make dataframes to plot using 'scores' and the model dataframe

Tscores <- TModel %>%
  select(ID, Species, Model, Predicted)%>%
  group_by(ID, Species, Model)%>%
  mutate(Aperture=row_number())%>%
  left_join(scores)

Wscores <- WModel %>%
  select(ID, Species, Model, Predicted)%>%
  group_by(ID, Species, Model)%>%
  mutate(Aperture=row_number())%>%
  left_join(scores)

Dscores <- DModel %>%
  select(ID, Species, Model, Predicted)%>%
  group_by(ID, Species, Model)%>%
  mutate(Aperture=row_number())%>%
  left_join(scores)

Sscores <- SModel %>%
  select(ID, Species, Model, Predicted)%>%
  group_by(ID, Species, Model)%>%
  mutate(Aperture=row_number())%>%
  left_join(scores)


## Te

TeLinear <- Tscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Linear")%>%
  select(theta, y, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, y), colour="black")+
    geom_point(aes(theta, Predicted), colour="red")+
    ylab("Te Lin")
TeQuad <- Tscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Polynomial")%>%
  select(theta, y, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, y), colour="black")+
    geom_point(aes(theta, Predicted), colour="orange")+
    ylab("Te Quad")
TeLog <- Tscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Log")%>%
  select(theta, y, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, y), colour="black")+
    geom_point(aes(theta, Predicted), colour="yellow")+
    scale_y_log10()+
    ylab("Te Log")
TePower <- Tscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Power")%>%
  select(theta, y, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, y), colour="black")+
    geom_point(aes(theta, exp(Predicted)), colour="green")+
    #scale_y_log10()+
    ylab("Te Powr")
TeExp <- Tscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Exponential")%>%
  select(theta, y, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, y), colour="black")+
    geom_point(aes(theta, exp(Predicted)), colour="blue")+
    ylab("Te Exp")


## We

WeLinear <- Wscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Linear")%>%
  select(theta, ApertureArea, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, ApertureArea), colour="black")+
    geom_point(aes(theta, Predicted), colour="red")+
    ylab("We Lin")
WeQuad <- Wscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Polynomial")%>%
  select(theta, ApertureArea, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, ApertureArea), colour="black")+
    geom_point(aes(theta, Predicted), colour="orange")+
    ylab("We Quad")
WeLog <- Wscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Log")%>%
  select(theta, ApertureArea, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, ApertureArea), colour="black")+
    geom_point(aes(theta, Predicted), colour="yellow")+
    scale_y_log10()+
    ylab("We Log")
WePower <- Wscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Power")%>%
  select(theta, ApertureArea, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, ApertureArea), colour="black")+
    geom_point(aes(theta, exp(Predicted)), colour="green")+
    ylab("We Powr")
WeExp <- Wscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Exponential")%>%
  select(theta, ApertureArea, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, ApertureArea), colour="black")+
    geom_point(aes(theta, exp(Predicted)), colour="blue")+
    ylab("We Exp")

## De

DeLinear <- Dscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Linear")%>%
  select(theta, center_gap, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, center_gap), colour="black")+
    geom_point(aes(theta, Predicted), colour="red")+
    ylab("De Lin")
DeQuad <- Dscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Polynomial")%>%
  select(theta, center_gap, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, center_gap), colour="black")+
    geom_point(aes(theta, Predicted), colour="orange")+
    ylab("De Quad")
DeLog <- Dscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Log")%>%
  select(theta, center_gap, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, center_gap), colour="black")+
    geom_point(aes(theta, Predicted), colour="yellow")+
    scale_y_log10()+
    ylab("De Log")
DePower <- Dscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Power")%>%
  select(theta, center_gap, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, center_gap), colour="black")+
    geom_point(aes(theta, exp(Predicted)), colour="green")+
    ylab("De Powr")
DeExp <- Dscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Exponential")%>%
  select(theta, center_gap, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, center_gap), colour="black")+
    geom_point(aes(theta, exp(Predicted)), colour="blue")+
    ylab("De Exp")

## Se

SeLinear <- Sscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Linear")%>%
  select(theta, Comp.1, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, Comp.1), colour="black")+
    geom_point(aes(theta, Predicted), colour="red")+
    ylab("Se Lin")
SeQuad <- Sscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Polynomial")%>%
  select(theta, Comp.1, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, Comp.1), colour="black")+
    geom_point(aes(theta, Predicted), colour="orange")+
    ylab("Se Quad")
SeLog <- Sscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Log")%>%
  select(theta, Comp.1, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, Comp.1), colour="black")+
    geom_point(aes(theta, Predicted), colour="yellow")+
    scale_y_log10()+
    ylab("Se Log")
SePower <- Sscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Power")%>%
  select(theta, Comp.1, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, Comp.1), colour="black")+
    geom_point(aes(theta, exp(Predicted)), colour="green")+
    scale_y_log10()+
    ylab("Se Powr")
SeExp <- Sscores %>%
  filter(grepl("WM862",ID)) %>%
  filter(Model == "Exponential")%>%
  select(theta, Comp.1, Predicted)%>%
  distinct()%>%
  ggplot()+
    geom_point(aes(theta, Comp.1), colour="black")+
    geom_point(aes(theta, exp(Predicted)), colour="blue")+
    ylab("Se Exp")


grid.arrange(TeLinear, TeQuad, TeLog, TePower, TeExp,
             WeLinear, WeQuad, WeLog, WePower, WeExp,
             DeLinear, DeQuad, DeLog, DePower, DeExp,
             SeLinear, SeQuad, SeLog, SePower, SeExp, ncol=5, nrow=4)



```

Moving forward, though, we only want one model per parameter. Using the R2, AICc and residual plots, we can choose which ones. 

```{r selectAllModels}
growthrates <- combined %>%
  tidyr::separate(Species, into=c("Species", "Age"), sep=",")%>%
  filter(!Species %in% c("Tylospira_coronata_replicate1", "Tylospira_coronata_replicate2", "Tylospira_coronata_replicate3"))%>%
  mutate(grp=paste("all"))%>%
  dplyr::select(ID, Species, Age, TPower, 
                WExponential, DPolynomial, SLinear) %>%
  dplyr::rename(Wpara=WExponential, 
                Tpara=TPower,
                Dpara=DPolynomial, 
                Spara=SLinear)
```


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

useparas <- growthrates%>%
  unique()%>%
  drop_na()

## for the sake of ease of plotting, let's pull Genus name out of the Species variable
useparas$Genus <- str_extract(useparas$Species, "[^_]+")

write.csv(useparas, file="final_parameters_for_use.csv", row.names=F)

```

Now we can plot the Raup cube!

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

## Convex hulls

chullsTD <- useparas %>%
  split(list(useparas$Species)) %>%
  lapply(. , function(x) {
    x <- as.data.frame(x)
    x[chull(x[ , "Tpara"], x[ , "Dpara"]), ]
  }) %>%
  bind_rows()

chullsTS <- useparas %>%
  split(list(useparas$Species)) %>%
  lapply(. , function(x) {
    x <- as.data.frame(x)
    x[chull(x[ , "Tpara"], x[ , "Spara"]), ]
  }) %>%
  bind_rows()

chullsTW <- useparas %>%
  split(list(useparas$Species)) %>%
  lapply(. , function(x) {
    x <- as.data.frame(x)
    x[chull(x[ , "Tpara"], x[ , "Wpara"]), ]
  }) %>%
  bind_rows()

chullsDS <- useparas %>%
  split(list(useparas$Species)) %>%
  lapply(. , function(x) {
    x <- as.data.frame(x)
    x[chull(x[ , "Dpara"], x[ , "Spara"]), ]
  }) %>%
  bind_rows()

chullsDW <- useparas %>%
  split(list(useparas$Species)) %>%
  lapply(. , function(x) {
    x <- as.data.frame(x)
    x[chull(x[ , "Wpara"], x[ , "Dpara"]), ]
  }) %>%
  bind_rows()

chullsSW <- useparas %>%
  split(list(useparas$Species)) %>%
  lapply(. , function(x) {
    x <- as.data.frame(x)
    x[chull(x[ , "Spara"], x[ , "Wpara"]), ]
  }) %>%
  bind_rows()

## all snails

## substituting in the commented out versions of lines will produce black and white versions that form the panels from Figure 5 in the paper

## T x W
P1 <- ggplot(useparas, aes(Wpara, Tpara, colour=Species))+
#P1 <- ggplot(useparas, aes(Wpara, Tpara, shape=Species))+
  #geom_polygon(data=chullsTW, alpha=0.2)+
  geom_polygon(data=chullsTW, alpha=0.2, aes(fill=Species))+
  geom_point(size=2)+
  #scale_shape_manual(values=c(19, 17, 15, 5, 4, 3, 2, 1, 0, 8, 10))+
  #scale_color_grey()+
  #scale_fill_grey()+
  xlab("W")+
  ylab("T")+
  theme_bw()

## W x D
P2 <- ggplot(useparas, aes(Dpara,Wpara, colour=Species))+
#P2 <- ggplot(useparas, aes(Dpara, Wpara, shape=Species))+
  #geom_polygon(data=chullsDW, alpha=0.2)+
  geom_polygon(data=chullsDW, alpha=0.2, aes(fill=Species))+
  geom_point(size=2)+
  #scale_shape_manual(values=c(19, 17, 15, 5, 4, 3, 2, 1, 0, 8, 10))+
  #scale_color_grey()+
  #scale_fill_grey()+
  xlab("D")+
  ylab("W")+
  theme_bw()

## T x D
P3 <- ggplot(useparas, aes(Dpara, Tpara, colour=Species))+
#P3 <- ggplot(useparas, aes(Dpara, Tpara, shape=Species))+  
  #geom_polygon(data=chullsTD, alpha=0.2)+
  geom_polygon(data=chullsTD, alpha=0.2, aes(fill=Species))+
  geom_point(size=2)+
  #scale_shape_manual(values=c(19, 17, 15, 5, 4, 3, 2, 1, 0, 8, 10))+
  #scale_color_grey()+
  #scale_fill_grey()+
  xlab("D")+
  ylab("T")+
  theme_bw()

## S x W
P4 <- ggplot(useparas, aes(Spara, Wpara, colour=Species))+
#P4 <- ggplot(useparas, aes(Spara, Wpara, shape=Species))+
    #geom_polygon(data=chullsSW, alpha=0.2)+
  geom_polygon(data=chullsSW, alpha=0.2, aes(fill=Species))+
  geom_point(size=2)+
  #scale_shape_manual(values=c(19, 17, 15, 5, 4, 3, 2, 1, 0, 8, 10))+
  #scale_color_grey()+
  #scale_fill_grey()+
  xlab("S")+
  ylab("W")+
  theme_bw()

## S x D
P5 <- ggplot(useparas, aes(Spara, Dpara, colour=Species))+
#P5 <- ggplot(useparas, aes(Spara, Dpara, shape=Species))+
  #geom_polygon(data=chullsDS, alpha=0.2)+
  geom_polygon(data=chullsDS, alpha=0.2, aes(fill=Species))+
  geom_point(size=2)+
  #scale_shape_manual(values=c(19, 17, 15, 5, 4, 3, 2, 1, 0, 8, 10))+
  #scale_color_grey()+
  #scale_fill_grey()+
  xlab("S")+
  ylab("D")+
  theme_bw()
 

## S x T
P6 <- ggplot(useparas, aes(Spara, Tpara, colour=Species))+
#P6 <- ggplot(useparas, aes(Spara, Tpara, shape=Species))+
  #geom_polygon(data=chullsTS, alpha=0.2)+
  geom_polygon(data=chullsTS, alpha=0.2, aes(fill=Species))+
  geom_point(size=2)+
  #scale_shape_manual(values=c(19, 17, 15, 5, 4, 3, 2, 1, 0, 8, 10))+
  #scale_color_grey()+
  #scale_fill_grey()+
  xlab("S")+
  ylab("T")+
  theme_bw()

## BLANK GROB IS KEY
P7 <- ggplot(useparas, aes(Spara, Tpara, colour=Species))+
  geom_blank()+
  theme(axis.text=element_blank(),
        axis.title=element_blank(),
        axis.ticks=element_blank(),
        axis.line=element_blank(),
        line=element_blank(),
        rect=element_blank())

figure1 <- ggarrange(P1, P3, P6, P7, P2, P4, P7, P7, P5, ncol=3, nrow=3, common.legend = TRUE, legend="bottom")

figure1
```


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

## Let's export the mean R2 and relative likelihood parameters as a table for reporting

Dsummary <- modelFittingD %>%
  group_by(Model)%>%
  summarise(meanR2 = mean(Rsq),
            meanRelLik = mean(relativeLikelihood))%>%
  mutate(Parameter=paste("De"))

Tsummary <- modelFittingT %>%
  group_by(Model)%>%
  summarise(meanR2 = mean(Rsq),
            meanRelLik = mean(relativeLikelihood))%>%
  mutate(Parameter=paste("Te"))

Ssummary <- modelFittingS %>%
  group_by(Model)%>%
  summarise(meanR2 = mean(Rsq),
            meanRelLik = mean(relativeLikelihood))%>%
  mutate(Parameter=paste("Se"))

Wsummary <- modelFittingW %>%
  group_by(Model)%>%
  summarise(meanR2 = mean(Rsq),
            meanRelLik = mean(relativeLikelihood))%>%
  mutate(Parameter=paste("We"))

modelSummary <- Dsummary %>%
  bind_rows(Tsummary)%>%
  bind_rows(Ssummary)%>%
  bind_rows(Wsummary)

write.csv(modelSummary, file="struth_likelihood_and_R2.csv", row.names=F)

```

And a 3d plot! 

```{r 3dCube}
## 3D interactive

library(plotly)

plot_ly(x=useparas$Spara, y=useparas$Wpara, z=useparas$Dpara, type="scatter3d",
        mode="markers", color=useparas$Species, colors=scales::hue_pal()(12))
```
