### 3 Oct 2017: Created to analyze methodology study with data from Viewbox
### 11 Jan 2018: Cleaned up
### 8-16 Feb 2018: Incorporated comments & additions by Dr. Carmelo Fruciano & Dr. Vera Weisbecker to perform analyses on the symmetric component of shape only (including the Procrustes ANOVAs with nested terms) as well as export certain datasets into MorphoJ
### 27 Feb 2018: Improved readability of script for peer review

##############################################################
"########################  Data Set up  ######################"
##############################################################

library(geomorph)

## Read in data file and get it into format that geomorph recognizes
my_data <- read.csv(file.choose(), header = T) # choose Raw_Coordinates.csv

## Remove problematic fixed landmarks (these were identified during landmarking as unreliable and removed before any analyses were run)
temp <- my_data[-(907:912),] # A pt of ext auditory meatus
temp <- temp[-(889:894),] # A pt squamosal / Alisphenoid
temp <- temp[-(877:882),] # AL premax&maxillar intersect
temp <- temp[-(847:852),] # Bullae&basioccciotal/basisphenoid
temp <- temp[-(769:774),] # Deepest pt aud bullae
temp <- temp[-(652:657),] # L pt of occipital condyle
temp <- temp[-(508:513),] # Lat pt foramen under infraorbital
my_data <- temp[-(283:288),] # P pt of foramen ovale
# This leaves 58 landmarks, 145 semilandmarks, and 86 patches for a total of 289 points

data <- my_data[, 2:(dim(my_data)[2])] # gets rid of LM names for geomorph formatting reasons
data <- t(data)
A <- arrayspecs(data, 289, 3) # convert to 3D array with 289 points and 3 dimensions

## Create tables of metadata
library(stringr)

# Create CatNum and scan type metadata
foo <- my_data[0,]  # gets a list with mesh file name
sp_names <- names(foo[-1])  # creates a list of species names
categories <- strsplit(sp_names, "_") 
my_classifiers <- matrix(unlist(categories), ncol=3, byrow=T) 
colnames(my_classifiers) <- c("CatNum", "Type", "Rep")
sp_info <- as.data.frame(my_classifiers)
sp_info$Both <- with(sp_info, interaction(CatNum, Type, sep = "_"))
sp_info$All <- with(sp_info, interaction(CatNum, Type, Rep, sep = "_"))

## Remove problem specimen JM14056
# The meshes for this specimen were mirror-imaged before landmarking making the 3D and CT scans not comparable
good_data <- (sp_info$CatNum != "JM14056")
A <- A[ , ,good_data]
sp_info <- sp_info[good_data, ]

sp_info$Order <- c(1:nrow(sp_info)) # preserve landmarking order

# Create LM point name metadata
pt_names <-my_data[, 1]
foo2 <- str_sub(pt_names, 3, -1)
u_pt_names<- unique(foo2)

# Create vectors of patch points and semi-landmarks
patches <- str_detect(u_pt_names, "PAT")
pat_num <- which(patches == TRUE)
sli_matrix <- read.csv(file.choose(), header = T) # choose Smatrix.csv 

# Load in bilateral landmarks
land.pairs <- read.csv(file.choose(), header = F) # choose Bilateral_Landmarks.csv

##############################################################
"##########################  Analyses  #######################"
##############################################################

"### 1) MorphoJ Procrustes ANOVA ###"
"###################################"
### Table 1 ###

## Perform procrustes alignment in R with geomorph
Y <- gpagen(A, Proj = TRUE, ProcD = TRUE, curves = sli_matrix, surfaces = pat_num)

## Get data into MorphoJ, where error is designed as a nested factor in ANOVAs
library(Morpho)
r2morphoj(Y,"DataForTable1.txt")

## Analyze in MorphoJ:
# Make a new project with "DataForTable1.txt", tick 3D and object symmetry
# Perform: Preliminaries > New Procrustes Fit with alignment by principal axis
# If any pairings are inconsistent use the Bilateral_Landmark.csv doc to fix
# Perform: Preliminaries > Extract Identifiers from ID string, do so for CatNum (1, 7), Device (-4, -3), and Rep (-1, -1)
# Perform: Variation > Procrustes ANOVA with CatNum as Individual, Device as Error 1 and Rep as Error 2
# Scroll up to see results reported in Table 1

"### 2) Symmetric Procrustes ANOVA ###"
"#####################################"
### Table 2 ###

## Perform bilateral symmetry test and extract symmetric component of shape
bY <- bilat.symmetry(Y$coords, ind = sp_info$All, object.sym = T, replicate = NULL, side = NULL, land.pairs = land.pairs)

ind.coords <- bY$symm.shape # Analogous to Y$coords

# Create new table of specimen information since the bilat function rearranges the coordinate array
b_names <- dimnames(ind.coords)[[3]]
b_names <- str_replace(b_names, "ind", "")
categories <- strsplit(b_names, "_")
my_classifiers <- matrix(unlist(categories), ncol=3, byrow=T)
colnames(my_classifiers) <- c("CatNum", "Type", "Rep")
b_names <- as.data.frame(my_classifiers)
b_names$Both <- with(b_names, interaction(CatNum, Type, sep = "_"))
b_names$All <- with(b_names, interaction(CatNum, Type, Rep, sep = "_"))

## Run Procrustes ANOVA
DFforGeneralProcANOVA <- geomorph.data.frame(coords = ind.coords, ind = factor(b_names$CatNum), Dev = b_names$Type) 

GeneralProcANOVA <- procD.lm(coords ~ ind + Dev %in% ind, data = DFforGeneralProcANOVA, iter = 999)

GeneralProcANOVA # Results reported in Table 2

# Calculate variation explained among repetitions within scan type
GeneralProcANOVA$SS[3]/GeneralProcANOVA$SS[4]

"### 3) Symmetric Exploratory PCA ###"
"####################################"
### Figures 4 & 5 ###

## Graph the PCA with the symmetric component
gp <- as.factor(b_names$CatNum)
col.gp <- rainbow(20)
names(col.gp) <- levels(gp)
col.gp <- col.gp[match(gp, names(col.gp))] # colors by CatNum

pca <- plotTangentSpace(ind.coords, groups = col.gp, axis1 = 1, axis2 = 2, label = b_names$Type, verbose = T) # change axis1 to see warpplots for different PCs

# We changed colors in Illustrator for better visibility

"### 4) Repeatability ###"
"########################"
### Table 3 ###

## Perform separate GPAs & bilateral symmetry for scan subsets to get repeatability for individual+scan type
is_ct <- (sp_info$Type == "CT")
ct_set <- A[, ,is_ct]
ct_sp_info <- sp_info[is_ct, ]

ct_Y <- gpagen(ct_set, Proj = TRUE, ProcD = TRUE, curves = sli_matrix, surfaces = pat_num)

ct_bY <- bilat.symmetry(ct_Y, ind = ct_sp_info$All, object.sym = T, replicate = NULL, side = NULL, land.pairs = land.pairs)

ct_ind.coords <- ct_bY$symm.shape

# Get metadata table again
ct_b_names <- dimnames(ct_ind.coords)[[3]]
ct_b_names <- str_replace(ct_b_names, "ind", "")
categories <- strsplit(ct_b_names, "_")
my_classifiers <- matrix(unlist(categories), ncol=3, byrow=T)
colnames(my_classifiers) <- c("CatNum", "Type", "Rep")
ct_b_names <- as.data.frame(my_classifiers)
ct_b_names$All <- with(ct_b_names, interaction(CatNum, Type, Rep, sep = "_"))

## Calculate repeatability for CT scans only
ct.gdf <- geomorph.data.frame(coords = ct_ind.coords, ct_ind = factor(ct_b_names$CatNum))
rep.er <- procD.lm(coords ~ ct_ind, data = ct.gdf, iter = 999)
repeatability <- ((rep.er$aov.table$MS[1] - rep.er$aov.table$MS[2])/2) / (rep.er$aov.table$MS[2] + ((rep.er$aov.table$MS[1] - rep.er$aov.table$MS[2])/2))

# See results reported in Table 3a
rep.er$aov.table
repeatability

## Table 3b: Perform separate GPAs on repeatability on just 3D scans
is_td <- (sp_info$Type == "3D")
td_set <- A[, ,is_td]
td_sp_info <- sp_info[is_td, ]

td_Y <- gpagen(td_set, Proj = TRUE, ProcD = TRUE, curves = sli_matrix, surfaces = pat_num)

td_bY <- bilat.symmetry(td_Y, ind = td_sp_info$All, object.sym = T, replicate = NULL, side = NULL, land.pairs = land.pairs)

td_ind.coords <- td_bY$symm.shape

# Get metadata table again
td_b_names <- dimnames(td_ind.coords)[[3]]
td_b_names <- str_replace(td_b_names, "ind", "")
categories <- strsplit(td_b_names, "_")
my_classifiers <- matrix(unlist(categories), ncol=3, byrow=T)
colnames(my_classifiers) <- c("CatNum", "Type", "Rep")
td_b_names <- as.data.frame(my_classifiers)
td_b_names$All <- with(td_b_names, interaction(CatNum, Type, Rep, sep = "_"))

# Check if names and coordinates are matched up
gp <- as.factor(td_b_names$CatNum)
col.gp <- rainbow(20)
names(col.gp) <- levels(gp)
col.gp <- col.gp[match(gp, names(col.gp))] #colors by CatNum

pca <- plotTangentSpace(td_bY$symm.shape, groups = col.gp, axis1 = 1, axis2 = 2, label = td_b_names$CatNum, verbose = T)

## Calculate repeatability for 3D scans only
td.gdf <- geomorph.data.frame(coords = td_ind.coords, td_ind = factor(td_b_names$CatNum))
rep.er <- procD.lm(coords ~ td_ind, data = td.gdf, iter = 999)
repeatability <- ((rep.er$aov.table$MS[1] - rep.er$aov.table$MS[2])/2) / (rep.er$aov.table$MS[2] + ((rep.er$aov.table$MS[1] - rep.er$aov.table$MS[2])/2))

# See results reported in Table 3b
rep.er$aov.table
repeatability

"### 5) Analysis of Disparity ###"
"################################"
### Figure 6 and in-text reports of p-value and mean variance ###

## Make individual/type combos a factor to use the model's residuals to assess intra-observer error
Cat <- b_names$Both
Dev <- b_names$Type
gdf_CF <- geomorph.data.frame(ind.coords, Cat = Cat, Dev= Dev)
Disparity <- morphol.disparity(ind.coords ~Cat, groups= NULL, iter = 999) 

## Create box plot of disparity that shows how variances within each Device/Individual combination compare - i.e. how operator error within replicates compares
library(stringi)

DisparityFactor <- as.factor(unlist(lapply(names(Disparity$Procrustes.var), function(x) stri_sub(x,-2))))

ProcVarCF <- data.frame(ProcVar = Disparity$Procrustes.var, Device = DisparityFactor)

library(ggplot2)
BoxPlot <- ggplot(ProcVarCF, aes(x=Device, y=ProcVar, fill=Device)) + 
        geom_boxplot()+
        theme_classic()

BoxPlot # Results shown in Figure 6

dev.off()

# Test for differences in variance and get P-values for those differences
Disparity2 <- morphol.disparity(ind.coords ~Cat, groups= ~Dev, iter = 999)

Disparity2 # Results reported in text

"### 6) Procrustes ANOVA w Sex in MorphoJ ###"
"############################################"
### Table 4 ###

library(dplyr)
pse_sex <- read.csv(file.choose(), header = T) # choose Pse_sex.csv
sex <- dplyr::filter(pse_sex, Sex != "u")
sp_info_sex <- merge(sp_info, sex, by = "CatNum", all.x = FALSE)
good_set <- A[,,sp_info_sex$Order] # subsets the coordinate data AND re-orders it so they match

# Redo gpa & bilateral symmetry with sexed-only dataset
Y_sex <- gpagen(good_set, Proj = TRUE, ProcD = TRUE, curves = sli_matrix, surfaces = pat_num)

bY_sex <- bilat.symmetry(Y_sex, ind = sp_info_sex$All, object.sym = T, replicate = NULL, side = NULL, land.pairs = land.pairs)

ind.coords_sex <- bY_sex$symm.shape

# Create new table of specimen information since the bilat function rearranges the coordinate array
b_names_sex <- dimnames(ind.coords_sex)[[3]]
b_names_sex <- str_replace(b_names_sex, "ind", "")
categories <- strsplit(b_names_sex, "_")
my_classifiers <- matrix(unlist(categories), ncol=3, byrow=T)
colnames(my_classifiers) <- c("CatNum", "Type", "Rep")
b_names_sex <- as.data.frame(my_classifiers)
b_names_sex$All <- with(b_names_sex, interaction(CatNum, Type, Rep))
b_names_sex$Order <- c(1:nrow(b_names_sex))

b_names_sex <- merge(b_names_sex, sex, by = "CatNum", all.x = FALSE)
b_names_sex$Class <- with(b_names_sex, interaction(All, Sex, sep = "."))

ind.coords_sex <- ind.coords_sex[,,b_names_sex$Order]

dimnames(ind.coords_sex)[[3]] <- b_names_sex$Class

# Export as a dataset into MorphoJ
r2morphoj(ind.coords_sex,"Sex_dataset.txt")

## In MorphoJ:
# Make new project with "Sex_dataset.txt" and tick 3D & No object symmetry (since its already been taken into account)
# Perform: Preliminaries > Extract Identifiers from ID string, do so for CatNum (1, 7), Device (-6, -5), Rep (-3, -3), and Sex (-1, -1)
# Perform: Preliminaries > Procrustes fit with alignment by principal axis
# Perform: Variation > Procrustes ANOVA with Device as error 1 and Sex as Error 2
# Scroll up to see results reported in Table 4

"### 7) Cross Validation Classification Tests ###"
"################################################"
### Table 5 ###

## Subset sexed-only dataset into 3D and CT scan datasets ##
bcoordsSex_CT <- ind.coords_sex[,,which(b_names_sex$Type == "CT")]
b_nameSex_CT <- b_names_sex[which(b_names_sex$Type == "CT"),]

## Calculate the mean coordinates for each specimen
MeanCoords_CTsex <- aggregate(vecx(bcoordsSex_CT), by = list(b_nameSex_CT$CatNum), FUN = mean) 

# Update metadata to reflect just one of each specimen and take only the CatNum and Sex columns 
mb_nameSex_CT <- unique(b_nameSex_CT[,c(1,6)])

rownames(mb_nameSex_CT) <- mb_nameSex_CT$CatNum
rownames(MeanCoords_CTsex) <- MeanCoords_CTsex$Group.1
mb_nameSex_CT <- mb_nameSex_CT[rownames(MeanCoords_CTsex),]

## Do the same for 3D data
bcoordsSex_3D <- ind.coords_sex[,,which(b_names_sex$Type=="3D")]
b_nameSex_3D <- b_names_sex[which(b_names_sex$Type=="3D"),]

## Calculate the mean coordinates for each specimen
MeanCoords_3Dsex <- aggregate(vecx(bcoordsSex_3D), by = list(b_nameSex_3D$CatNum), FUN = mean) 

# Update metadata to reflect just one of each specimen and take only the CatNum and Sex columns
mb_nameSex_3D <- unique(b_nameSex_3D[,c(1,6)])

rownames(mb_nameSex_3D) <- mb_nameSex_3D$CatNum
rownames(MeanCoords_3Dsex) <- MeanCoords_3Dsex$Group.1
mb_nameSex_3D <- mb_nameSex_3D[rownames(MeanCoords_3Dsex),]

## Create the between group PCA
s_bwgPCA_CT <- groupPCA(MeanCoords_CTsex[,2:ncol(MeanCoords_CTsex)], mb_nameSex_CT$Sex, rounds = 10000, cv = T, weighting = F)

s_bwgPCA_3D <- groupPCA(MeanCoords_3Dsex[,2:ncol(MeanCoords_3Dsex)], mb_nameSex_3D$Sex, rounds = 10000, cv = T, weighting = F)

# Results reported in Table 5
s_bwgPCA_CT
s_bwgPCA_3D

"### 8) bwgPCAs of Sex Data ###"
"##############################"
### Figure 7 ### 

## Make PCA of sexed-only CT data
# Color points by sex
gp <- as.factor(b_nameSex_CT$Sex)
col.gp <- c("red","blue", "grey")
names(col.gp) <- levels(gp)
col.gp <- col.gp[match(gp, names(col.gp))] 

pca <- plotTangentSpace(bcoordsSex_CT, groups = col.gp, axis1 = 1, axis2 = 2, label = NULL, verbose = T)

## Repeat for PCA of sexed-only 3D data
gp <- as.factor(b_nameSex_3D$Sex)
col.gp <- c("red","blue", "grey")
names(col.gp) <- levels(gp)
col.gp <- col.gp[match(gp, names(col.gp))] 

pca <- plotTangentSpace(bcoordsSex_3D, groups = col.gp, axis1 = 1, axis2 = 2, label = NULL, verbose = T)

# PCAs seen in Figure 7