---
title: "Instar Duration and Time to Molt"
author: "Ben Jahnes"
date: "Mar. 16, 2018"
output:
  pdf_document: default
  html_notebook: default
---

Data 
```{r}
#####################################################################
setwd("~/Manuscripts/Plots")
gutstats<- read.csv("gutlength.cumulative.GF.GT.WT.csv", header =T,row.names=1)
#Datav<- read.csv("GrowTestEnv.csv", header =T,row.names=1)

#########################################################################
gutstats$Trt<-as.factor(gutstats$Trt)
gutstats$coh<-as.factor(gutstats$coh)
gutstats$ind<-as.factor(gutstats$ind)
gutstats$ins<-as.factor(gutstats$ins)
gutstats$trt<-as.factor(gutstats$trt)
##gutlength$ins<-as.numeric(gutlength$ins)
Data=gutstats
#before n=173
Data.na <-Data[complete.cases(Data[ , 7:11]),]
#after n=138
```

Libraries
```{r}

library(vegan)
library(MASS)
library(vegan3d)
library(mgcv)
library(RColorBrewer)
library(scatterplot3d)
library(car)
citation()
citation("vegan")
```

SplitData
```{r}

insplit<-split(Data,Data$ins)
ins3<-insplit$'3'
ins4<-insplit$'4' 
ins5<-insplit$'5'
ins6<-insplit$'6'



#N<-c("8","14","14")

#envsplit<-split(Datav,Datav$age)
#Envins1<-envsplit$'1'
#Envins2<-envsplit$'2' 
#Envins3<-envsplit$'3' 
#Data=na.omit(Envins3)
#Data$ins<-as.numeric(Data$age)
```


```{r}
#variance explained by variables
Data.na.var<-Data[complete.cases(Data[ , 7:14]),]
varp <- varpart (Data.na.var[ ,7:14],  ~wg, ~ hg,  ~mg, data = Data.na.var[ ,7:14])
varp
plot(varp)
```

```{r}
#variance explained by variables
Data.na.var<-Data[complete.cases(Data[ , 7:14]),]
varp2 <- varpart (Data.na.var[ ,7:14], ~ fg,  ~gmass, ~ bmass,~ bw, data = Data.na.var[ ,7:14])
varp2
plot(varp2)
```


PCA5
```{r}
Data.na5 <-ins5[complete.cases(ins5[ , 7:11]),]
ord5 <- rda(Data.na5[ ,7:11]~ Trt+mg+hg+wg,data=Data.na5)
ord5
summary(ord5)
```

```{r}
na5.dist<-vegdist(Data.na5[,7:11], method="euclidean")
na5mr<-mrpp(na5.dist,Data.na5$Trt)
na5mr

na5ad<-adonis(na5.dist ~ Trt, data=Data.na5, perm=1e3,method="euclidean")
na5an <- anosim(Data.na5[,7:11],Data.na5$Trt )
na5an <- anosim(na5.dist,Data.na5$Trt )
na5ad
summary(na5an)
plot(na5an)

```




```{r}
ef <- envfit(ord5 ~ Trt+mg+hg+wg, data=Data.na5, permu=1000)
ef

```
Distance between centroids: instar 5
```{r}
ord5centroids<-ord5$CCA$centroids

ord5centroids<-ord5centroids[,1:3] 
distord5<-dist(ord5centroids,method = "euclidean", diag = T)
distord5
site.scr5 <- scores(ord5, display = "sites", scaling = "symmetric", choices = 1:3)
distord5sites<-dist(site.scr5,method = "euclidean", diag = T)
```


Plot instar 5
```{r}

col <- c("#66C2A5","#8DA0CB","#FC8D62")
pch <- 17:20

par(mfrow=c(1,1))
ordiplot(ord5, display="sites", xlab="PCA1", ylab="PCA2")
with(Data.na5,points(ord5, display = "sites", col = col[trt],pch = pch[ins],cex=1, 
font=2))
ordispider(ord5, show.groups=1,groups=Data.na5$trt,col=c(c("#66C2A5")), label=F)
ordispider(ord5, show.groups=2,groups=Data.na5$trt,col=c(c("#8DA0CB")), label=F)
ordispider(ord5, show.groups=3,groups=Data.na5$trt,col=c(c("#FC8D62")), label=F)
#ordiellipse(ord5,kind="se", show.groups=1,groups=Data.na5$trt,col=c(c("#66C2A5")), label=F)
#ordiellipse(ord5,kind="se", show.groups=2,groups=Data.na5$trt,col=c(c("#8DA0CB")), label=F)
#ordiellipse(ord5,kind="se", show.groups=3,groups=Data.na5$trt,col=c(c("#FC8D62")), label=F)
legend("topleft",title="Treatment", legend = c("WT","GT","GF"), bty = "n",
       col = col, pch = 19, pt.bg = col)
title(main="PCA: Roach Morphological Characteristics
5th Instar", cex.main=0.80, xlim=c(-.9,1.7), ylim=c(-1.5,1))
ef <- envfit(ord5 ~ mg+hg+wg, data=Data.na5, permu=1000)
plot(ef)
#ordisurf(ord5 ~ hg,data=Data.na5, add=TRUE, bubble=TRUE)
```
```{r}
ordirgl(ord5, display="sites",type="t")

orglspider (ord5, groups = Data.na5$trt)

#orglhull (ord, groups = Data.na$ins, col = 'tomato', alpha = 0.5)
```


```{r}
#na5.dist<-vegdist(Data.na5[,7:11], method="euclidean")
na5mr<-mrpp(distord5sites,Data.na5$Trt)
na5mr

na5ad<-adonis(distord5sites ~ Trt, data=Data.na5, perm=1e3,method="euclidean")
na5an <- anosim(Data.na5[,7:11],Data.na5$Trt )
na5an <- anosim(distord5sites,Data.na5$Trt )
na5ad
summary(na5an)
plot(na5an)

```


```{r}
Data.na.all <-Data[complete.cases(Data[ , 7:11]),]
ordall <- rda(Data.na.all[ ,7:11] ~Trt+mg+hg+wg,data=Data.na.all)
ordall
summary(ordall)
```


```{r}
ef <- envfit(ordall ~ Trt+mg+hg+wg, data=Data.na.all, permu=1000)
ef

```
```{r}
all.dist<-vegdist(Data.na.all[,7:11], method="euclidean")
allmr<-mrpp(all.dist,Data.na.all$Trt)
allmr

allad<-adonis(all.dist ~ Trt, data=Data.na.all, perm=1e3,method="euclidean")
allan <- anosim(Data.na.all[,7:11],Data.na.all$Trt )
allan <- anosim(all.dist,Data.na.all$Trt )
allad
summary(allan)
plot(allan)

```
Distance between centroids: All Instar, 4 variables

```{r}
ordallcentroids<-ordall$CCA$centroids

ordallcentroids<-ordallcentroids[,1:3] 
distordall<-dist(ordallcentroids,method = "euclidean", diag = T)
distordall
site.scrall <- scores(ordall, display = "sites", scaling = "symmetric", choices = 1:3)
distordallsites<-dist(site.scrall,method = "euclidean", diag = T)
```

Plot All
```{r}
Data.na.all$ins<-as.numeric(Data.na.all$ins)
col <- c("#66C2A5","#8DA0CB","#FC8D62")
pch <- 17:20

par(mfrow=c(1,1))
plot(ordall, display="sites", xlab="PCA1", ylab="PCA2")
ef <- envfit(ordall ~ mg+hg+wg, data=Data.na.all, permu=1000)
plot(ef)
with(Data.na.all,points(ordall, display = "sites", col = col[trt],pch = pch[ins],cex=1, 
font=2))

ordispider(ordall, show.groups=1,groups=Data.na.all$trt,col=c(c("#66C2A5")), label=F)
ordispider(ordall, show.groups=2,groups=Data.na.all$trt,col=c(c("#8DA0CB")), label=F)
ordispider(ordall, show.groups=3,groups=Data.na.all$trt,col=c(c("#FC8D62")), label=F)
#ordiellipse(ordall, kind="se", show.groups=1,groups=Data.na$trt,col=c(c("#66C2A5")), label=F)
#ordiellipse(ordall, kind="se", show.groups=2,groups=Data.na$trt,col=c(c("#8DA0CB")), label=F)
#ordiellipse(ordall, kind="se", show.groups=3,groups=Data.na$trt,col=c(c("#FC8D62")), label=F)
#ordisurf(ordall ~ ins,data=Data.na.all, add=TRUE, bubble=F)
legend("topleft",title="Treatment", legend = c("WT","GT","GF"), bty = "n",
       col = col, pch = 15, pt.bg = col)
legend("topright",title="instar", legend = c("3rd","4th","5th","6th"), bty = "n",
                     col = 1, pch = pch, pt.bg = col)
title(main="PCA: Roach Morphological Characteristics
      All Instars", cex.main=0.80, xlim=c(-.9,1.7), ylim=c(-1.5,1))
#ordisurf(ordall ~ wg,data=Data.na, add=TRUE)


```




```{r}
cols <- bgs <- c("tomato","steelblue","palegreen4")
pchs <- 15:19
par(mfrow=c(1,1))
plot(ordall, display="sites", xlab="PCA1", ylab="PCA2",scaling=2)
ef <- envfit(ordall ~ bl+mg+hg+wg, data=Data.na.all, permu=1000)
plot(ef)
with(Data.na,points(site.scrall[84:107,1], site.scrall[84:107,2],
              col = cols[1], bg = bgs[1], pch = 19))
with(Data.na,points(site.scrall[108:128,1], site.scrall[108:128,2],
              col = cols[2], bg = bgs[2], pch = 19))
with(Data.na,points(site.scrall[129:149,1], site.scrall[129:149,2],
              col = cols[3], bg = bgs[3], pch = 19))

```





Distance between centroids: all Instar all variables

```{r}
ordall2centroids<-ordall2$CCA$centroids

ordall2centroids<-ordall2centroids[,1:3] 
distordall2<-dist(ordall2centroids,method = "euclidean", diag = T)
distordall2
site.scrall2 <- scores(ordall2, display = "sites", scaling = "symmetric", choices = 1:3)
distordall2sites<-dist(site.scrall2,method = "euclidean", diag = T)
```

Plot All2
```{r}

col <- c("#66C2A5","#8DA0CB","#FC8D62")
pch <- 17:20

par(mfrow=c(1,1))
plot(ordall2, display="sites", xlab="PCA1", ylab="PCA2")
ef <- envfit(ordall2 ~ bl+mg+hg+wg, data=Data.na.all2, permu=1000)
plot(ef)
with(Data.na.all2,points(ordall2, display = "sites", col = col[trt],pch = pch[ins],cex=1, 
font=2))

ordispider(ordall2, show.groups=1,groups=Data.na.all2$trt,col=c(c("#66C2A5")), label=F)
ordispider(ordall2, show.groups=2,groups=Data.na.all2$trt,col=c(c("#8DA0CB")), label=F)
ordispider(ordall2, show.groups=3,groups=Data.na.all2$trt,col=c(c("#FC8D62")), label=F)
ordiellipse(ordall2,kind="se", show.groups=1,groups=Data.na.all2$trt,col=c(c("#66C2A5")), label=F)
ordiellipse(ordall2,kind="se", show.groups=2,groups=Data.na.all2$trt,col=c(c("#8DA0CB")), label=F)
ordiellipse(ordall2,kind="se", show.groups=3,groups=Data.na.all2$trt,col=c(c("#FC8D62")), label=F)
legend("topleft",title="Treatment", legend = c("WT","GT","GF"), bty = "n",
       col = col, pch = 15, pt.bg = col)
legend("topright",title="instar", legend = c("3rd","4th","5th","6th"), bty = "n",
                     col = 1, pch = pch, pt.bg = col)
title(main="PCA: Roach Morphological Characteristics
      All Instars, All Parameters", cex.main=0.80, xlim=c(-.9,1.7), ylim=c(-1.5,1))
#ordisurf(ordall2 ~ bl,data=Data.na.all2, add=TRUE)

```

```{r}
Data.na.5.2 <-ins5[complete.cases(ins5[ , 8:12]),]
ord5.2 <- rda(Data.na.5.2[ ,8:12] ~ Trt+bl+mg+hg+wg+bmass,data=Data.na.5.2)
ord5.2
summary(ord5.2)
```
```{r}
#variance explained by variables

varp2 <- varpart (Data.na.5.2[ ,7:14], ~ wg, ~ hg,  ~mg, ~ bl , data = Data.na.5.2[ ,7:14])
varp2
plot(varp2)
```

```{r}
ef <- envfit(ord5.2 ~ Trt+bl+mg+hg+wg+bmass, data=Data.na.5.2, permu=1000)
ef

```
Distance between centroids: instar 5 all parameters
```{r}
ord5.2centroids<-ord5.2$CCA$centroids

ord5.2centroids<-ord5.2centroids[,1:3] 
distord5.2<-dist(ord5.2centroids,method = "euclidean", diag = T)
distord5.2
site.scr5.2 <- scores(ord5.2, display = "sites", scaling = "symmetric", choices = 1:3)
distord5.2sites<-dist(site.scr5.2,method = "euclidean", diag = T)
```

```{r}
#na5.2.dist<-vegdist(Data.na5.2[,7:11], method="euclidean")
na5.2mr<-mrpp(distord5.2sites,Data.na.5.2$Trt)
na5.2mr

na5.2ad<-adonis(distord5.2sites ~ Trt, data=Data.na.5.2, perm=1e3,method="euclidean")
na5.2an <- anosim(Data.na.5.2[,7:14],Data.na.5.2$Trt )
na5.2an <- anosim(distord5.2sites,Data.na.5.2$Trt )
na5.2ad
summary(na5.2an)
plot(na5.2an)

```

Plot instar 5.2
```{r}

col <- c("#66C2A5","#8DA0CB","#FC8D62")
pch <- 17:20

par(mfrow=c(1,1))
ordiplot(ord5.2, display="sites", xlab="PCA1", ylab="PCA2")
with(Data.na.5.2, points(ord5.2, display = "sites", col = col[trt],pch = pch[ins],cex=1, 
font=2))
ordispider(ord5.2, show.groups=1,groups=Data.na.5.2$trt,col=c(c("#66C2A5")), label=F)
ordispider(ord5.2, show.groups=2,groups=Data.na.5.2$trt,col=c(c("#8DA0CB")), label=F)
ordispider(ord5.2, show.groups=3,groups=Data.na.5.2$trt,col=c(c("#FC8D62")), label=F)
ordiellipse(ord5.2,kind="se", show.groups=1,groups=Data.na.5.2$trt,col=c(c("#66C2A5")), label=F)
ordiellipse(ord5.2,kind="se", show.groups=2,groups=Data.na.5.2$trt,col=c(c("#8DA0CB")), label=F)
ordiellipse(ord5.2,kind="se", show.groups=3,groups=Data.na.5.2$trt,col=c(c("#FC8D62")), label=F)

legend("topleft",title="Treatment", legend = c("WT","GT","GF"), bty = "n",
       col = col, pch = 19, pt.bg = col)
title(main="PCA: Roach Morphological Characteristics
5th Instar, all parameters", cex.main=0.80, xlim=c(-.9,1.7), ylim=c(-1.5,1))
ef <- envfit(ord5.2 ~ bl+mg+hg+wg, data=Data.na.5.2, permu=1000)
plot(ef)
#ordisurf(ord5.2 ~ hg,data=Data.na.5.2, add=TRUE)

```



Plot instar 5
```{r}

col <- c("#66C2A5","#8DA0CB","#FC8D62")
pch <- 17:20

par(mfrow=c(1,1))
ordiplot(ord5, display="sites", xlab="PCA1", ylab="PCA2")
with(Data.na5,points(ord5, display = "sites", col = col[trt],pch = pch[ins],cex=1, 
font=2))
ordispider(ord5, show.groups=1,groups=Data.na5$trt,col=c(c("#66C2A5")), label=F)
ordispider(ord5, show.groups=2,groups=Data.na5$trt,col=c(c("#8DA0CB")), label=F)
ordispider(ord5, show.groups=3,groups=Data.na5$trt,col=c(c("#FC8D62")), label=F)
#ordiellipse(ord5,kind="se", show.groups=1,groups=Data.na5$trt,col=c(c("#66C2A5")), label=F)
#ordiellipse(ord5,kind="se", show.groups=2,groups=Data.na5$trt,col=c(c("#8DA0CB")), label=F)
#ordiellipse(ord5,kind="se", show.groups=3,groups=Data.na5$trt,col=c(c("#FC8D62")), label=F)
legend("topleft",title="Treatment", legend = c("WT","GT","GF"), bty = "n",
       col = col, pch = 19, pt.bg = col)
title(main="PCA: Roach Morphological Characteristics
5th Instar", cex.main=0.80, xlim=c(-.9,1.7), ylim=c(-1.5,1))
ef <- envfit(ord5 ~ bl+mg+hg+wg, data=Data.na5, permu=1000)
#plot(ef)
ordisurf(ord5 ~ hg,data=Data.na5, add=TRUE, col="black")
```