library(FD)
library(cluster)
library(mclust)
library(MASS)
library(vegan)
library(plot3D)
library(vegan)
library(reshape2)
library(ape)
library(clustMD)
library(RColorBrewer)
library(car)
library(rgl)
library(colorspace)
library(ggplot2)
library(plyr)
library(reshape2)
library(gridExtra)
library(ggrepel)
library(scatterplot3d)


temp.data <- read.csv(file="RawData.csv")
pilosa.data <- temp.data
colnames(pilosa.data)
pilosa.data[,c(3:12,15,22:23)] <- log(pilosa.data[,c(3:12,15,22:23)]+1)  

pilosa.wide <- pilosa.data[,2:37]

pilosa.wide.reduced <- pilosa.wide[,c(2:11,14,21:28)] # remove component traits

pilosa.all.diss <- daisy(pilosa.wide.reduced[,1:19], metric='gower', type=list(symm=c(16,17,18)))

################################################################
#
# PCoA
#
################################################################
par(mfrow=c(1,1))
pilosa.all.pcoa <- pcoa(pilosa.all.diss, correction='cailliez')

# basic plot
biplot.pcoa(pilosa.all.pcoa)
plot(pilosa.all.pcoa$vectors.cor[,1], pilosa.all.pcoa$vectors.cor[,2],
     pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
text(pilosa.all.pcoa$vectors.cor[,1], pilosa.all.pcoa$vectors.cor[,2],
     label=pilosa.wide$Collector_Collection_number,cex=0.3, pos=4)

# examine all axes
pair1<- pairs(pilosa.all.pcoa$vectors.cor[,1:6],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair2 <- pairs(pilosa.all.pcoa$vectors.cor[,7:12],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair3 <- pairs(pilosa.all.pcoa$vectors.cor[,13:18],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair4 <- pairs(pilosa.all.pcoa$vectors.cor[,19:24],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair5 <- pairs(pilosa.all.pcoa$vectors.cor[,25:30],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair6 <- pairs(pilosa.all.pcoa$vectors.cor[,31:36],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair7 <- pairs(pilosa.all.pcoa$vectors.cor[,37:42],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair8 <- pairs(pilosa.all.pcoa$vectors.cor[,43:48],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair9 <- pairs(pilosa.all.pcoa$vectors.cor[,49:54],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair10 <- pairs(pilosa.all.pcoa$vectors.cor[,55:60],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair11 <- pairs(pilosa.all.pcoa$vectors.cor[,61:66],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair12 <- pairs(pilosa.all.pcoa$vectors.cor[,67:72],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair13 <- pairs(pilosa.all.pcoa$vectors.cor[,73:78],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair14 <- pairs(pilosa.all.pcoa$vectors.cor[,79:84],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair15 <- pairs(pilosa.all.pcoa$vectors.cor[,85:90],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair16 <- pairs(pilosa.all.pcoa$vectors.cor[,91:96],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair17 <- pairs(pilosa.all.pcoa$vectors.cor[,97:102],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])
pair18 <- pairs(pilosa.all.pcoa$vectors.cor[,103:106],
      pch=20, cex=1.5, col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4])


# ellipses
pilosa.ellipses1 <- dataEllipse(pilosa.all.pcoa$vectors.cor[,1], pilosa.all.pcoa$vectors.cor[,2],
                                groups=pilosa.data$Taxon_ID_4, plot.points=FALSE,
                                col=c("#fc8d59","#fee090","#4575b4","#d73027"),
                                lwd=1, group.labels=NULL)

par(mfrow=c(1,1))
# Basic scatterplot
plot(pilosa.all.pcoa$vectors.cor[,1:2], type='n',
     #main="colorbrewer2",
     xlab='PCoA Axis 1',
     ylab='PCoA Axis 2',
     cex.axis=1.5, cex.lab=1.5,
     xlim=rev(c(min(pilosa.all.pcoa$vectors.cor[,1] + (-.83)), max(pilosa.all.pcoa$vectors.cor[,1] + (.37)))),
     ylim=rev(c(min(pilosa.all.pcoa$vectors.cor[,2] + (-.15)), max(pilosa.all.pcoa$vectors.cor[,2] + (.05))))
)
rect(par('usr')[1],par('usr')[3],par('usr')[2],par('usr')[4],col='gray93')
grid(lty=1, lwd=1, col='white')
points(pilosa.all.pcoa$vectors.cor[,1], pilosa.all.pcoa$vectors.cor[,2],
       #longispica, nana, pilosa, steenensis
       bg=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4], col='darkgrey',
       pch=21, cex=3)

# Add 95% CI Ellipses
lines(pilosa.ellipses1$longispica$`0.95`, col='#fc8d59', lty=3, lwd=2)
lines(pilosa.ellipses1$pilosa$`0.95`, col='#4575b4', lty=3, lwd=2)
lines(pilosa.ellipses1$nana$`0.95`, col='#fee090', lty=3, lwd=4)
lines(pilosa.ellipses1$steenensis$`0.95`, col='#d73027', lty=3, lwd=2)



##3D, including axis 3 of PCoA
x <- pilosa.all.pcoa$vectors.cor[,1]
y <- pilosa.all.pcoa$vectors.cor[,2]
z <- pilosa.all.pcoa$vectors.cor[,3]

scatter3D_fancy <- function(x,y,z,...,colvar=NULL)
{
  panelfirst <- function(pmat) {
    XY <- trans3D(x,y,z = rep(min(z), length(z)), pmat=pmat)
    scatter2D(XY$x, XY$y, colvar=NULL, pch=20,
              #taxon names for colors below are longispica, pilosa, nana, steenensis
              col=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4],
              cex=0.6, alpha=0.7, add=TRUE, colkey=FALSE)
  }
  scatter3d(x,y,z, ..., colvar=NULL, panel.first=panelfirst,
            colkey=FALSE)
}
scatter3D_fancy(x,y,z, pch=21, 
                cex=2, 
                type='h', 
                lwd=0.3,
                bty='u',
                bg=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4],
                col='darkgrey',
                col.panel='gray93', col.grid='white',
                ticktype='detailed',
                theta=20,
                phi=5,
                d=2, 
                xlab="PCoA Axis 1",
                ylab="PCoA Axis 2",
                zlab="PCoA Axis 3",
                line=1.5)





################################################################
#
# FUZZY CLUSTERING
#
################################################################

#pilosa.wide -> pilosa.wide.reduced -> pilosa.diss, all in alphabetical order
#numbers in pilosa.diss, are used as the labels for fuzzy clustering
#so number 1 in pilosa.diss = the Atwood collection
#results of clustering are ordered by silhouette widths

############## 4 clusters ##############
clust.num <- 4
pilosa.all.k4.fuz <- fanny(pilosa.all.diss, k=clust.num, memb.exp=1.3, maxit=100000)
summary(pilosa.all.k4.fuz)
plot(pilosa.all.k4.fuz)
k4.pilosa.sil <- silhouette(pilosa.all.k4.fuz, pilosa.all.diss)

trial.k4 <- as.data.frame(cbind(pilosa.all.k4.fuz$silinfo$widths, pilosa.all.k4.fuz$membership[as.integer(rownames(k4.pilosa.sil)),], as.character(pilosa.data$Collector_Collection_number), 
                             as.character(pilosa.data$Taxon_ID_4[as.integer(rownames(k4.pilosa.sil))])))


############## 3 clusters ##############
clust.num <- 3
pilosa.all.k3.fuz <- fanny(pilosa.all.diss, k=clust.num, memb.exp=1.3, maxit=100000)
summary(pilosa.all.k3.fuz)
plot(pilosa.all.k3.fuz)
k3.pilosa.sil <- silhouette(pilosa.all.k3.fuz, pilosa.all.diss)

trial.k3 <- as.data.frame(cbind(pilosa.all.k3.fuz$silinfo$widths, pilosa.all.k3.fuz$membership[as.integer(rownames(k3.pilosa.sil)),], pilosa.data$Collector_Collection_number, 
                             as.character(pilosa.data$Taxon_ID_4[as.integer(rownames(k3.pilosa.sil))])))


############## 2 clusters ##############
clust.num <- 2
pilosa.all.k2.fuz <- fanny(pilosa.all.diss, k=clust.num, memb.exp=1.3, maxit=100000)
summary(pilosa.all.k2.fuz)
plot(pilosa.all.k2.fuz)
k2.pilosa.sil <- silhouette(pilosa.all.k2.fuz, pilosa.all.diss)

trial.k2 <- as.data.frame(cbind(pilosa.all.k2.fuz$silinfo$widths, #
                                pilosa.all.k2.fuz$membership[as.integer(rownames(k2.pilosa.sil)),],
                                pilosa.data$Collector_Collection_number, 
                                as.character(pilosa.data$Taxon_ID_4[as.integer(rownames(k2.pilosa.sil))])))



########## plot it all together
par(mar=c(3,3,3,3))
par(mfrow=c(1,3))

plot(k4.pilosa.sil, col=c("#fc8d59","#fee090","#4575b4","#d73027")[trial.k4$V9],
     main="", sub="")

plot(k3.pilosa.sil, col=c("#fc8d59","#fee090","#4575b4","#d73027")[trial.k3$V8],
     main="", sub="")


plot(k2.pilosa.sil, col=c("#fc8d59","#fee090","#4575b4","#d73027")[trial.k2$V7],
     main="", sub="")


### assignment probabilities for each individual to each cluster

par(mfrow=c(2,2))
par(mar=c(2,2,2,2))

plot(pilosa.all.pcoa$vectors.cor[,1:2], type='n',
     #main="colorbrewer2",
     xlab='PCoA Axis 1',
     ylab='PCoA Axis 2',
     cex.axis=1.5, cex.lab=1.5
     #xlim=c(min(pilosa.all.pcoa$vectors.cor[,1] + (-.83)), max(pilosa.all.pcoa$vectors.cor[,1] + (.37))),
     #ylim=c(min(pilosa.all.pcoa$vectors.cor[,2] + (-.15)), max(pilosa.all.pcoa$vectors.cor[,2] + (.05)))
)
points(pilosa.all.pcoa$vectors.cor[,1], pilosa.all.pcoa$vectors.cor[,2],
       #longispica, nana, pilosa, steenensis
       bg=c("#fc8d59","#fee090","#4575b4","#d73027")[pilosa.wide$Taxon_ID_4], col='darkgrey',
       pch=21, cex=2)
text(pilosa.all.pcoa$vectors.cor[,1], pilosa.all.pcoa$vectors.cor[,2],
     label=rownames(pilosa.wide),cex=0.3, pos=4)


plot(pilosa.all.pcoa$vectors.cor[,1:2], type='n',
     main="Cluster (k) = 4",
     xlab='PCoA Axis 1',
     ylab='PCoA Axis 2',
     cex.axis=1.5, cex.lab=1.5
     #xlim=c(min(pilosa.all.pcoa$vectors.cor[,1] + (-.83)), max(pilosa.all.pcoa$vectors.cor[,1] + (.37))),
     #ylim=c(min(pilosa.all.pcoa$vectors.cor[,2] + (-.15)), max(pilosa.all.pcoa$vectors.cor[,2] + (.05)))
)
stars(pilosa.all.k4.fuz$membership, locations = pilosa.all.pcoa$vectors.cor[,1:2],
      draw.segments=TRUE, add=TRUE, scale=FALSE,
      len=0.04, col.segments = c("#7b3294","#c2a5cf","#a6dba0","#008837"))


plot(pilosa.all.pcoa$vectors.cor[,1:2], type='n',
     main="Cluster (k) = 3",
     xlab='PCoA Axis 1',
     ylab='PCoA Axis 2',
     cex.axis=1.5, cex.lab=1.5
     #xlim=c(min(pilosa.all.pcoa$vectors.cor[,1] + (-.83)), max(pilosa.all.pcoa$vectors.cor[,1] + (.37))),
     #ylim=c(min(pilosa.all.pcoa$vectors.cor[,2] + (-.15)), max(pilosa.all.pcoa$vectors.cor[,2] + (.05)))
)
stars(pilosa.all.k3.fuz$membership, locations = pilosa.all.pcoa$vectors.cor[,1:2],
      draw.segments=TRUE, add=TRUE, scale=FALSE,
      len=0.04, col.segments = c("#7b3294","#a6dba0","#008837"))


plot(pilosa.all.pcoa$vectors.cor[,1:2], type='n',
     main="Cluster (k) = 2",
     xlab='PCoA Axis 1',
     ylab='PCoA Axis 2',
     cex.axis=1.5, cex.lab=1.5
     #xlim=c(min(pilosa.all.pcoa$vectors.cor[,1] + (-.83)), max(pilosa.all.pcoa$vectors.cor[,1] + (.37))),
     #ylim=c(min(pilosa.all.pcoa$vectors.cor[,2] + (-.15)), max(pilosa.all.pcoa$vectors.cor[,2] + (.05)))
)
stars(pilosa.all.k2.fuz$membership, locations = pilosa.all.pcoa$vectors.cor[,1:2],
      draw.segments=TRUE, add=TRUE, scale=FALSE,
      len=0.04, col.segments = c("#7b3294","#008837"))






###################################
# VIOLIN / dot plots for quantitative traits
# omit labels by commenting out 'geom_label_repel'
###################################

plant.height <- ggplot(temp.data, aes(x=temp.data$Taxon_ID_4, y=temp.data$Plant_height_cm)) +
  geom_violin(trim=FALSE) +
  geom_dotplot(aes(group=as.factor(temp.data$Taxon_ID_4), color=as.factor(temp.data$Taxon_ID_4), fill=as.factor(temp.data$Taxon_ID_4)),
               binaxis = 'y', stackdir = 'center', binwidth = 0.8) +
  scale_fill_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  scale_color_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  theme(legend.title=element_blank()) +
  theme(legend.position="none") +
  ylab("Plant height (cm)") +
  xlab(NULL) +
  geom_label_repel(label=rownames(temp.data))

bract.length <- ggplot(temp.data, aes(x=temp.data$Taxon_ID_4, y=temp.data$Bract_length_mm)) +
  geom_violin(trim=FALSE) +
  geom_dotplot(aes(group=as.factor(temp.data$Taxon_ID_4), color=as.factor(temp.data$Taxon_ID_4), fill=as.factor(temp.data$Taxon_ID_4)),
               binaxis = 'y', stackdir = 'center', binwidth = 0.8) +
  scale_fill_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  scale_color_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  theme(legend.title=element_blank()) +
  theme(legend.position="none") +
  ylab("Bract length (mm)") +
  xlab(NULL)+
  geom_label_repel(label=rownames(temp.data))

total.length <- ggplot(temp.data, aes(x=temp.data$Taxon_ID_4, y=temp.data$Total_length_mm)) +
  geom_violin(trim=FALSE) +
  geom_dotplot(aes(group=as.factor(temp.data$Taxon_ID_4), color=as.factor(temp.data$Taxon_ID_4), fill=as.factor(temp.data$Taxon_ID_4)),
               binaxis = 'y', stackdir = 'center', binwidth = 0.4) +
  scale_fill_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  scale_color_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  theme(legend.title=element_blank()) +
  theme(legend.position="none") +
  ylab("Total length (mm)") +
  xlab(NULL)+
  geom_label_repel(label=rownames(temp.data))

length.raceme <- ggplot(temp.data, aes(x=temp.data$Taxon_ID_4, y=temp.data$Length_of_raceme_mm)) +
  geom_violin(trim=FALSE) +
  geom_dotplot(aes(group=as.factor(temp.data$Taxon_ID_4), color=as.factor(temp.data$Taxon_ID_4), fill=as.factor(temp.data$Taxon_ID_4)),
               binaxis = 'y', stackdir = 'center', binwidth = 3) +
  scale_fill_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  scale_color_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  theme(legend.title=element_blank()) +
  theme(legend.position="none") +
  ylab("Length of raceme (mm)") +
  xlab(NULL)+
  geom_label_repel(label=rownames(temp.data))

bract.width <- ggplot(temp.data, aes(x=temp.data$Taxon_ID_4, y=temp.data$Bract_width_mm)) +
  geom_violin(trim=FALSE) +
  geom_dotplot(aes(group=as.factor(temp.data$Taxon_ID_4), color=as.factor(temp.data$Taxon_ID_4), fill=as.factor(temp.data$Taxon_ID_4)),
               binaxis = 'y', stackdir = 'center', binwidth = 0.15) +
  scale_fill_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  scale_color_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  theme(legend.title=element_blank()) +
  theme(legend.position="none") +
  ylab("Bract width (mm)") +
  xlab(NULL)+
  geom_label_repel(label=rownames(temp.data))

beak.tube.ratio <- ggplot(temp.data, aes(x=temp.data$Taxon_ID_4, y=temp.data$Beak_tube_ratio)) +
  geom_violin(trim=FALSE) +
  geom_dotplot(aes(group=as.factor(temp.data$Taxon_ID_4), color=as.factor(temp.data$Taxon_ID_4), fill=as.factor(temp.data$Taxon_ID_4)),
               binaxis = 'y', stackdir = 'center', binwidth = 0.08) +
  scale_fill_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  scale_color_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  theme(legend.title=element_blank()) +
  theme(legend.position="none") +
  ylab("Beak to tube ratio") +
  xlab(NULL)+
  geom_label_repel(label=rownames(temp.data))

leaf.length <- ggplot(temp.data, aes(x=temp.data$Taxon_ID_4, y=temp.data$Leaf_length_mm)) +
  geom_violin(trim=FALSE) +
  geom_dotplot(aes(group=as.factor(temp.data$Taxon_ID_4), color=as.factor(temp.data$Taxon_ID_4), fill=as.factor(temp.data$Taxon_ID_4)),
               binaxis = 'y', stackdir = 'center', binwidth = 1.5) +
  scale_fill_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  scale_color_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  theme(legend.title=element_blank()) +
  theme(legend.position="none") +
  ylab("Leaf length (mm)") +
  xlab(NULL)+
  geom_label_repel(label=rownames(temp.data))

point.lobe.attach <- ggplot(temp.data, aes(x=temp.data$Taxon_ID_4, y=temp.data$Point_of_lobe_attachment_mm)) +
  geom_violin(trim=FALSE) +
  geom_dotplot(aes(group=as.factor(temp.data$Taxon_ID_4), color=as.factor(temp.data$Taxon_ID_4), fill=as.factor(temp.data$Taxon_ID_4)),
               binaxis = 'y', stackdir = 'center', binwidth = .3) +
  scale_fill_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  scale_color_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  theme(legend.title=element_blank()) +
  theme(legend.position="none") +
  ylab("Point of lobe attachment (mm)") +
  xlab(NULL)+
  geom_label_repel(label=rownames(temp.data))

ammt.calyx.subequal <- ggplot(temp.data, aes(x=temp.data$Taxon_ID_4, y=temp.data$Absolute_value_between_calyx_sinus1_sinus2)) +
  geom_violin(trim=FALSE) +
  geom_dotplot(aes(group=as.factor(temp.data$Taxon_ID_4), color=as.factor(temp.data$Taxon_ID_4), fill=as.factor(temp.data$Taxon_ID_4)),
               binaxis = 'y', stackdir = 'center', binwidth = 0.08) +
  scale_fill_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  scale_color_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  theme(legend.title=element_blank()) +
  theme(legend.position="none") +
  ylab("Amount of calyx subequality") +
  xlab(NULL)+
  geom_label_repel(label=rownames(temp.data))

leaf.width <- ggplot(temp.data, aes(x=temp.data$Taxon_ID_4, y=temp.data$Leaf_width_mm)) +
  geom_violin(trim=FALSE) +
  geom_dotplot(aes(group=as.factor(temp.data$Taxon_ID_4), color=as.factor(temp.data$Taxon_ID_4), fill=as.factor(temp.data$Taxon_ID_4)),
               binaxis = 'y', stackdir = 'center', binwidth = 0.1) +
  scale_fill_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  scale_color_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  theme(legend.title=element_blank()) +
  theme(legend.position="none") +
  ylab("Leaf width (mm)") +
  xlab(NULL)+
  geom_label_repel(label=rownames(temp.data))

calyx.length <- ggplot(temp.data, aes(x=temp.data$Taxon_ID_4, y=temp.data$Calyx_length_mm)) +
  geom_violin(trim=FALSE) +
  geom_dotplot(aes(group=as.factor(temp.data$Taxon_ID_4), color=as.factor(temp.data$Taxon_ID_4), fill=as.factor(temp.data$Taxon_ID_4)),
               binaxis = 'y', stackdir = 'center', binwidth = .5) +
  scale_fill_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  scale_color_manual(values=c("#fc8d59","#fee090","#4575b4","#d73027")) +
  theme(legend.title=element_blank()) +
  theme(legend.position="none") +
  ylab("Calyx length (mm)") +
  xlab(NULL)+
  geom_label_repel(label=rownames(temp.data))

########## organized by 'type' of trait
grid.arrange(plant.height,length.raceme,leaf.length,leaf.width, ncol=1,nrow=4) # plant height; raceme length; length/width of leaves

grid.arrange(bract.length,bract.width,point.lobe.attach,calyx.length, ncol=1, nrow=4) # bract length/width; point of lobe attachment; calyx length

grid.arrange(total.length,beak.tube.ratio,ammt.calyx.subequal, ncol=1, nrow=4) # total corolla length; beak:tube ratio; calyx subequality





#############################################################################
#
# Box Plots of discrete traits
#
#############################################################################

decumbent4 <- table(pilosa.data$Taxon_ID_4, pilosa.data$Decumbent)
plot(decumbent4, main="Decumbent habit", col=c("grey","white","black"), cex=2)
decumbent4b <- table(pilosa.data$Decumbent, pilosa.data$Taxon_ID_4)
plot(decumbent4b, main="Decumbent habit", col=c("#fc8d59","#fee090","#4575b4","#d73027"), cex=2)
#dev.copy2pdf(file="decumbent4.barplot.pdf", family="sans", useDingbats=FALSE)

hair.length4 <- table(pilosa.data$Taxon_ID_4, pilosa.data$Length)
plot(hair.length4, main="Herbage pubescence length", col=c("grey","white","black"), cex=2)
hair.length4b <- table(pilosa.data$Length, pilosa.data$Taxon_ID_4)
plot(hair.length4b, main="Herbage pubescence length", col=c("#fc8d59","#fee090","#4575b4","#d73027"), cex=2)
#dev.copy2pdf(file="hair.length4.barplot.pdf", family="sans", useDingbats=FALSE)

recurved4 <- table(pilosa.data$Taxon_ID_4, pilosa.data$Recurved_presence)
plot(recurved4, main="Presence of recurved hairs", col=c("white", "black"), cex=2)
recurved4b <- table(pilosa.data$Recurved_presence, pilosa.data$Taxon_ID_4)
plot(recurved4b, main="Presence of recurved hairs", col=c("#fc8d59","#fee090","#4575b4","#d73027"), cex=2)
#dev.copy2pdf(file="recurved4.barplot.pdf", family="sans", useDingbats=FALSE)

glandular4 <- table(pilosa.data$Taxon_ID_4, pilosa.data$Glandular_present)
plot(glandular4, main="Presence of glandular hairs", col=c("white", "black"), cex=2)
glandular4b <- table(pilosa.data$Glandular_present, pilosa.data$Taxon_ID_4)
plot(glandular4b, main="Presence of glandular hairs", col=c("#fc8d59","#fee090","#4575b4","#d73027"), cex=2)
#dev.copy2pdf(file="glandular4.barplot.pdf", family="sans", useDingbats=FALSE)

leaf.lobing4 <- table(pilosa.data$Taxon_ID_4, pilosa.data$Leaf_lobing)
plot(leaf.lobing4, main="Lobes on leaves", col=c("white", "black"), cex=2)
leaf.lobing4b <- table(pilosa.data$Leaf_lobing, pilosa.data$Taxon_ID_4)
plot(leaf.lobing4b, main="Lobes on leaves", col=c("#fc8d59","#fee090","#4575b4","#d73027"), cex=2)
#dev.copy2pdf(file="lobing4.barplot.pdf", family="sans", useDingbats=FALSE)

seg.shape4 <- table(pilosa.data$Taxon_ID_4, pilosa.data$Shape_of_segment)
plot(seg.shape4, main="Shape of calyx lobes",
     col=c("white","darkgrey","grey","lightgrey","black"), cex=0.8)
seg.shape4b <- table(pilosa.data$Shape_of_segment, pilosa.data$Taxon_ID_4)
plot(seg.shape4b, main="Shape of calyx lobes",
     col=c("#fc8d59","#fee090","#4575b4","#d73027"), cex=0.8)
#dev.copy2pdf(file="seg.shape4.barplot.pdf", family="sans", useDingbats=FALSE)

par(mfrow=c(3,2))
plot(glandular4, main="Presence of glandular hairs", col=c("white", "black"), cex=2)
plot(hair.length4, main="Herbage pubescence length", col=c("grey","white","black"), cex=2)
plot(recurved4, main="Presence of recurved hairs", col=c("white", "black"), cex=2)
plot(decumbent4, main="Decumbent habit", col=c("grey","white","black"), cex=2)
plot(leaf.lobing4, main="Lobes on leaves", col=c("white", "black"), cex=2)
plot(seg.shape4, main="Shape of calyx lobes",
     col=c("white","darkgrey","grey","lightgrey","black"), cex=0.8)


plot(glandular4b, main="Presence of glandular hairs", col=c("#fc8d59","#fee090","#4575b4","#d73027"), cex=2)
plot(hair.length4b, main="Herbage pubescence length", col=c("#fc8d59","#fee090","#4575b4","#d73027"), cex=2)
plot(recurved4b, main="Presence of recurved hairs", col=c("#fc8d59","#fee090","#4575b4","#d73027"), cex=2)
plot(decumbent4b, main="Decumbent habit", col=c("#fc8d59","#fee090","#4575b4","#d73027"), cex=2)
plot(leaf.lobing4b, main="Lobes on leaves", col=c("#fc8d59","#fee090","#4575b4","#d73027"), cex=2)
plot(seg.shape4b, main="Shape of calyx lobes",
     col=c("#fc8d59","#fee090","#4575b4","#d73027"), cex=0.8)




#########################################################################################
#
# print mean and stdev of each continuous trait for each focal taxon
#
#########################################################################################

library(plyr)
library(reshape2)
nana.data <- subset(temp.data, temp.data$Taxon_ID_4=='nana')
pilosasub.data <- subset(temp.data, temp.data$Taxon_ID_4=='pilosa')
long.data <- subset(temp.data, temp.data$Taxon_ID_4=='longispica')
steen.data <- subset(temp.data, temp.data$Taxon_ID_4=='steenensis')
all.pilosa.data <- subset(temp.data, temp.data$Taxon_ID_2=='pilosa')

nana.means <- lapply(nana.data[,3:23], mean)
nana.sd <- lapply(nana.data[,3:23], sd)
pilosa.means <- lapply(pilosasub.data[,3:23], mean)
pilosa.sd <- lapply(pilosasub.data[,3:23], sd)
long.means <- lapply(long.data[,3:23], mean)
long.sd <- lapply(long.data[,3:23], sd)
steen.means <- lapply(steen.data[,3:23], mean)
steen.sd <- lapply(steen.data[,3:23], sd)
all.pilosa.means <- lapply(all.pilosa.data[,3:23], mean)
all.pilosa.sd <- lapply(all.pilosa.data[,3:23], sd)

all.means.and.sd <- cbind(nana.means, nana.sd, pilosa.means, pilosa.sd, long.means, long.sd, steen.means, steen.sd, all.pilosa.means, all.pilosa.sd)
write.csv(all.means.and.sd, file='all.means.and.sd.csv')

table(pilosa.data$Glandular_present, pilosa.data$Taxon_ID_4)