#Empirical Example: Data from Marquez 2008 (Sigmodontine mandible shape)

library(geomorph)
source('compare.CR.r')

######### Mouse Mandibles: Marquez 2008
myland <- readland.tps("Marquez2008.tps")  #NOTE data already GPA-aligned
covariates <- read.csv("Covariates.csv",header=TRUE)
  Csize <- covariates$Size
  gp <- covariates$Group
mod.hyp <- read.csv('ModularHypoth.csv', header=TRUE, row.names = 1)
links <- read.table('links.txt')
gdf <- geomorph.data.frame(coords=myland, Csize = Csize, gp = gp)

#Account for allometry and group differences
reg.res <- procD.lm(coords~Csize+gp, data = gdf)
  shape.reg <- reg.res$GM$residuals + replicate(length(Csize), mshape(myland))
plotAllSpecimens(shape.reg)

##### 1: Compare Modular Hypotheses for Marquez (2008) data: all specimens
AllSpec.CR <- lapply(1:ncol(mod.hyp), function(j) modularity.test(shape.reg, partition.gp = as.factor(mod.hyp[,j])))
res.hypoth <- compare.CR(CR.null = TRUE, AllSpec.CR)
res.hypoth$sample.z
res.hypoth$pairwise.P #6-module hypothesis is best supported
AllSpec.CR[[5]]$CR.mat

################## 2: Compare Strength of modular signal across species
spec.CR <- lapply(1:nlevels(gp), function(j) modularity.test(shape.reg[,,which(gp==levels(gp)[j])],
                                                             partition.gp = as.factor(mod.hyp[,1])))
res.spec <- compare.CR(CR.null = FALSE, spec.CR)
  #Levels of modularity consistent across species 
res.spec

#### PLOTS

library(gplots)
#Species
Z.sp<-res.spec$sample.z
CI.sp<-qnorm(res.spec$sample.r.sd)/2

#Hypotheses
Z.h<-res.hypoth$sample.z
CI.h<-qnorm(res.hypoth$sample.r.sd)/2; CI.h[1] <- 0

par(mfrow=c(1,2))
plotCI(Z.h,ui =(Z.h+CI.h), li = (Z.h-CI.h), ylab="Z-CR", xlab="Alternative Modular Hypotheses",pch=21,cex=2, pt.bg="black",lwd=3,lty=1)
plotCI(Z.sp,ui =(Z.sp+CI.sp), li = (Z.sp-CI.sp), ylab="Z-CR", xlab = "Species Modularity",pch=21,cex=2, pt.bg="black",lwd=3,lty=1)
par(mfrow=c(1,1))
