# Load librarys
library(shapes)
library(car)
library(vegan)
library(heplots)
mshape<-function(A) { apply(A, c(1,2), mean) }
# Use the gorilla skull lankmark data in the shapes pacakge
# 30 females
gorf.dat
# 3-D array: 8 traits (landmarks), 2D coordinates, 30 individuals
dim(gorf.dat)

# 29 males
gorm.dat
dim(gorm.dat)

# let's combine them, first put all the data into one vector
gorilla_vector<-c(as.vector(gorf.dat),as.vector(gorm.dat))
# the new combined configuration
gor_mf<-array(gorilla_vector,dim=c(8,2,59))
gor_mf
# now we should have 8 landmarks for 59
dim(gor_mf)

# need vector for sex factor to do later comparisons
female<-rep("female",30)
male<-rep("male",29)
sex_factor<-as.factor(c(female,male))

# First, plot the raw data by sex
plot(gor_mf[1:8,1,1:30],gor_mf[1:8,2,1:30],pch=16,col="red", xlab="X", ylab="Y")
points(gor_mf[1:8,1,31:59],gor_mf[1:8,2,31:59],pch=16,col="blue")

# do a generalized procrustes analysis and PCA
GPA<-procGPA(gor_mf)
# plot all indiviuals post-rotation
plotshapes(GPA$rotated)
points(GPA$mshape,pch=16,col="red")
# plot post-rotation by sex
plot(GPA$rotated[1:8,1,1:30],GPA$rotated[1:8,2,1:30],pch=16,col="red",cex=.5, xlab="rotated X", ylab="rotated Y")
points(GPA$rotated[1:8,1,31:59],GPA$rotated[1:8,2,31:59],pch=16,col="blue",cex=.5)
# Add 95% confidence ellipses
for(l in 1:8)
  {
  dataEllipse(GPA$rotated[l,1,1:30],GPA$rotated[l,2,1:30],add=T, levels=0.95, col="red", center.pch=17, plot.points=FALSE)
  dataEllipse(GPA$rotated[l,1,31:59],GPA$rotated[l,2,31:59],add=T, levels=0.95, col="blue", center.pch=17, plot.points=FALSE)
  }

# plot the mean shape
plotshapes(GPA$mshape)
# plot individuals in PCA space
plot(GPA$stdscores[,1],GPA$stdscores[,2],pch=16,col=c("red","blue")[sex_factor])

# plot shape change vectors for first three PCA axes, magnifiy vectors 3x
shapepca(GPA,type="v",mag=3)

# To see how any individual or group shape changes, you can visualize with the tpsgrid function
# Compare average male or female to overall consensis
tpsgrid(GPA$mshape,mshape(GPA$rotated[,,1:30]),mag=6)
tpsgrid(GPA$mshape,mshape(GPA$rotated[,,31:59]),mag=6)
# Directly compare average male to average female
tpsgrid(mshape(GPA$rotated[,,1:30]),mshape(GPA$rotated[,,31:59]),mag=6)

# What about allometry?
# Look at male and female size
# centroid sizes
gor_size<-GPA$size
gor_size

# Plot size vs. PCA axis 1
plot(gor_size, GPA$stdscores[,1],pch=16,col=c("red","blue")[sex_factor])

# Mean size by sex
tapply(gor_size,sex_factor,mean)
# SD in size
tapply(gor_size,sex_factor,sd)
# Males are bigger
# In fact, the smallest male in the dataset is larger than the largest female
tapply(gor_size,sex_factor,min)
tapply(gor_size,sex_factor,max)

# The procGPA function is very "black box", it does the Procrustes rotation and subsequent PCA
# Lets extract the Procrustes part, then do a PCA as we did earlier in class
# make a matrix of 16 shape variables for the 59 individuals
shape_variables<-t(matrix(GPA$rotated,16,59))

# Do a PCA with the rda function
pca<-rda(shape_variables)
biplot(pca,type="text")

# Let's make sure these are the same, plot the first axis of the two approaches, they should be identical
plot(GPA$scores[,1],scores(pca)$sites[,1])
# They are...we can plot the PCA like we did before, should be identiacl to the earlier plot
plot(scores(pca)$sites[,1],scores(pca)$sites[,2],pch=16, col=c("red","blue")[sex_factor])

# Let's do a MANOVA!
# First, get the shape PCA scores. Note that we have 14 PCA axes (16-2). Two were dropped. There are 16 in the GPA$scores matrix, but the last two are essentially 0.
# How many axes should we include? The first 10 axes explain 98% of morphological variance
response_matrix<-cbind(scores(pca,choices = c(1:10))$sites[,1:10])
model<-manova(response_matrix~gor_size*sex_factor)
anova(model)
summary(model,test="Wilks")
etasq(model)