##########################################
##### Multivariate Ratio Analysis (MRA)
##### Shape PCA and Ratio Spectra, v1.02
##### Baur H., Leuenberger C. (2011) Systematic Biology 60, p. 818-819
##########################################

##### Arguments
 # U           a numeric matrix (or data frame) that provides the data for the principal component analysis
 # npc         number of shape components shown in numerical output, as usually only the first two or three components contain useful information. Default value is npc=3
 # rpc         number of decimal places used for printing components and loadings. Default value is rpc=4 

##### Values
 # shapePCA    returns the results of the shape PCA as an object of class "prcomp". Mainly used for displaying the screeplot
 # loadings    a matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors)
 # components  a matrix containing the variance, proportion of varaiance and cumulative proportion of variance
 # PCmatrix    a matrix of rotated data (i.e., the centered data multiplied by the loadings)
 # isosize     a matrix containing isometric size (i.e., the centered data multiplide by the isometric size vector)

ShapePCA <- function(U, npc=4, rpc=1){
  X <- log(U)
  X <- scale(X, center=TRUE, scale=FALSE)
  p <- dim(X)[2]
  I <- diag(1,p,p)
  a0 <- as.vector(rep(1,p)/p)
  P <- I-p*(a0%*%t(a0))
  Y <- X%*%P
  colnames(Y) <- colnames(X)
  isosize <- X%*%a0; colnames(isosize) <- "isosize"
  PCA <- prcomp(Y,center=FALSE, scale=FALSE)
  loadings <- PCA$rotation
  p <- dim(loadings)[2]
  PCmatrix <- PCA$x
  components.a <- PCA$sdev^2; components.b <- components.a/sum(components.a); components.c <- cumsum(components.b); components <- rbind(components.a, components.b, components.c)
  colnames(loadings) <- paste("shape.PC", 1:p, sep="")
  colnames(PCmatrix) <- colnames(loadings)
  colnames(components) <- colnames(loadings); rownames(components) <- c("Variance", "Proportion of Variance", "Cumulative Proportion")
  list(PCA=PCA, PCmatrix=PCmatrix[,1:npc], isosize=as.numeric(isosize), loadings=loadings[,1:npc], components=components[,1:npc], pc_var=round(components.b[1:npc]*100, rpc))
}








##########################################
##### Isometric size (to be calculated separately from function ShapePCA)
##########################################

isosize <- function(U) {
	X <- log(U)
	X <- scale(X, center=TRUE, scale=FALSE)
	p <- dim(X)[2]
	I <- diag(1,p,p)
	a0 <- as.vector(rep(1,p)/p)
	isosize <- X%*%a0
	isosize	
}







##########################################
##### Matrix of isometry free shape values (to be calculated separately from function ShapePCA)
##########################################

isofreeshapes <- function(U) {

X <- log(U)
X <- scale(X, center=TRUE, scale=FALSE)
p <- dim(X)[2]
I <- diag(1,p,p)
a0 <- as.vector(rep(1,p)/p)
P <- I-p*(a0%*%t(a0))
if.shapes <- X%*%P; colnames(if.shapes) <- colnames(X)
if.shapes
}






##########################################
##### PCA Ratio Spectrum v1.0
##########################################

##### Arguments
# x          a numeric matrix (or data frame) that provides the data for the PCA ratio spectrum
# pc         the shape component for which a PCA ratio spectrum should be drawn. Default value is pc=1
# bootrep    the number of bootstrap replicates on which the 68% confidence intervals are based. Default value is bootrep=500 
# barcol     colour of horizontal bars (confidence intervals). Default value is barcol="blue" 
# barlwd     line width of horizontal bars (confidence intervals). Default value is barlwd=2.5 
# linecol    colour of vertical line. Default value is linecol="black" 
# linelwd    line width of vertical line. Default value is linelwd=0.5 
# labelsize  size of variable labels. Default value is labelsize=0.9 
# labelfont  font of variable labels. Default value is labelfont=1 
# nosize     size of number labels. Default value is nosize=0.7 
# nofont     font of number labels. Default value is nofont=1
# maina      first part of main title. Default value is maina="PCA Ratio Spectrum for PC"
# mainb      second part of main title. Default value is mainb=pc (i.e., the number of the component, for which the PCA ratio spectrum is drawn)
# mainc      third part of main title. Default value is mainc=""
# suba       first part of sub title. Default value is suba="bars = 68% confidence intervals based on "
# subb       second part of sub title. Default value is subb=bootrep (i.e., the number of bootstrap replicates)
# subc       third part of sub title. Default value is subc=" bootstrap replicates"

pcaRS <- function (U, pc=1, sequence=1, bootrep=999, barcol="#1F78B4", barlwd=4, linecol="black", linelwd=0.5, labelsize=0.9, labelfont=1, nosize=0.7, nofont=1, maina="PCA Ratio Spectrum for shape PC", mainb=pc, mainc="", suba="bars = 68% confidence intervals based on ", subb=bootrep, subc=" bootstrap replicates", main=paste(maina, mainb, mainc, sep=""), sub=paste(suba, subb, subc, sep=""), ...) {
  trait.names <- colnames(U)
  Y <- log(U)
  Y <- scale(Y, center=TRUE, scale=FALSE)
  Sigma <- cov(Y)
  p <- length(trait.names)
  n <- nrow(Y)
  I <- diag(1, p, p)
  a0 <- as.vector(rep(1,p)/p)
  P <- I-a0%*%solve(t(a0)%*%a0)%*%t(a0)
  Sigma.1 <- P%*%Sigma%*%P
  u <- (eigen(Sigma.1)$vectors[, pc])*sequence
  rep.boot=bootrep
  U.boot <- matrix(0, nrow=rep.boot, ncol=p)
  for (i in c(1:rep.boot)){
    index.boot=sample(c(1:n), n, replace=TRUE)
    Y.boot <- Y[index.boot, ]
    Sigma.boot <- P%*%cov(Y.boot)%*%P
    u.boot <- (eigen(Sigma.boot)$vectors[, pc])*sequence
    if (u%*%u.boot<0){
      u.boot <- -u.boot
    }
    U.boot[i, ] <- u.boot
  }
  mean.u <- apply(U.boot, 2, mean)
  sd.u  <- apply(U.boot, 2, sd)
  mean.minus.sd <- mean.u-sd.u
  mean.plus.sd <- mean.u+sd.u
  m <- min(mean.u); M <- max(mean.u)
  u.sorted <- sort(mean.u, decreasing=TRUE)
  index <- sort(mean.u, index.return=TRUE, decreasing=TRUE)$ix
  trait.names.sorted <- trait.names[index]
  mean.minus.sd <- mean.minus.sd[index]
  mean.plus.sd <- mean.plus.sd[index]
  plot(c(0, 0), c(m, M), xlab="", ylab="", type="n", cex=1, asp=1, main=main, sub=sub, cex.lab=1, axes=FALSE, ...)
  lines(c(0, 0), c(m, M), col=linecol, lwd=linelwd)
  text(0, u.sorted[1], labels=signif(u.sorted[1], digits=2), pos=3, srt=0, col="black", font=nofont, cex=nosize)
  text(0, u.sorted[p], labels=signif(u.sorted[p], digits=2), pos=1, srt=0, col="black", font=nofont, cex=nosize)
  for (k in 1:p){
    lines(c(mean.minus.sd[k]-u.sorted[k], mean.plus.sd[k]-u.sorted[k]), c(u.sorted[k], u.sorted[k]), col=barcol, lwd=barlwd)
    if(k%%2==0){
      text(max(sd.u), u.sorted[k], labels=trait.names.sorted[k], pos=4, srt=0, col="black", font=labelfont, cex=labelsize)
    }
    if(k%%2==1){
      text(-max(sd.u), u.sorted[k], labels=trait.names.sorted[k], pos=2, srt=0, col="black", font=labelfont, cex=labelsize)
    }}
}











##########################################
##### Allometry Ratio Spectrum v1.0
##########################################

##### Arguments
# x          a numeric matrix (or data frame) that provides the data for the allometry ratio spectrum
# bootrep    the number of bootstrap replicates on which the 68% confidence intervals are based. Default value is bootrep=500 
# barcol     colour of horizontal bars (confidence intervals). Default value is barcol="blue" 
# barlwd     line width of horizontal bars (confidence intervals). Default value is barlwd=2.5 
# linecol    colour of vertical line. Default value is linecol="black" 
# linelwd    line width of vertical line. Default value is linelwd=0.5 
# labelsize  size of variable labels. Default value is labelsize=0.9 
# labelfont  font of variable labels. Default value is labelfont=1 
# nosize     size of number labels. Default value is nosize=0.7 
# nofont     font of number labels. Default value is nofont=1
# maina      first part of main title. Default value is maina="Allometry Ratio Spectrum"
# mainb      second part of main title. Default value is mainb=""
# mainc      third part of main title. Default value is mainc=""
# suba       first part of sub title. Default value is suba="bars = 68% confidence intervals based on "
# subb       second part of sub title. Default value is subb=bootrep (i.e., the number of bootstrap replicates)
# subc       third part of sub title. Default value is subc=" bootstrap replicates"

allometryRS <- function (U, bootrep=999, sequence=1, barcol="#33A02C", barlwd=4, linecol="black", linelwd=0.5, labelsize=0.9, labelfont=1, nosize=0.7, nofont=1, maina="Allometry Ratio Spectrum", mainb="", mainc="", suba="bars = 68% confidence intervals based on ", subb=bootrep, subc=" bootstrap replicates", main=paste(maina, mainb, mainc, sep=""), sub=paste(suba, subb, subc, sep=""), ...) {
  trait.names <- colnames(U)
  Y <- log(U)
  Y <- scale(Y,center=TRUE,scale=FALSE)
  Sigma <- cov(Y)
  p <- length(trait.names)
  n <- nrow(Y)
  a0 <- as.vector(rep(1,p)/p)
  a <- a0 
  d <- as.vector(Sigma%*%a)*sequence
  rep.boot=bootrep
  D.boot <- matrix(0, nrow=rep.boot, ncol=p)
  for (i in c(1:rep.boot)){
    index.boot=sample(c(1:n), n, replace=TRUE)
    Y.boot <- Y[index.boot, ]
    Sigma.boot <- cov(Y.boot)
    d.boot <- as.vector(Sigma.boot%*%a)*sequence
    d.boot <- d.boot/as.vector(sqrt(d.boot%*%d.boot))
    if (d%*%d.boot<0){
      d.boot <- -d.boot
    }
    D.boot[i, ] <- d.boot
  }
  mean.d <- apply(D.boot, 2, mean); mean.d
  sd.d <- apply(D.boot, 2, sd)
  mean.minus.sd <- mean.d-sd.d
  mean.plus.sd <- mean.d+sd.d
  m <- min(mean.d); M <- max(mean.d)
  d.sorted <- sort(mean.d, decreasing=TRUE)
  index <- sort(mean.d, index.return=TRUE, decreasing=TRUE)$ix
  trait.names.sorted <- trait.names[index]
  mean.minus.sd <- mean.minus.sd[index]
  mean.plus.sd <- mean.plus.sd[index]
  plot(c(0, 0), c(m, M), xlab=NA, ylab=NA, type="n", axes=FALSE, asp=1, main=main, sub=sub, cex=1, cex.lab=1, ...)
  lines(c(0, 0), c(m, M), col=linecol, lwd=linelwd)
  text(0, d.sorted[1], labels=signif(d.sorted[1], digits=2), pos=3, srt=0, col="black", font=nofont, cex=nosize)
  text(0, d.sorted[p], labels=signif(d.sorted[p], digits=2), pos=1, srt=0, col="black", font=nofont, cex=nosize)
  for (k in 1:p){
    lines(c(mean.minus.sd[k]-d.sorted[k], mean.plus.sd[k]-d.sorted[k]), c(d.sorted[k], d.sorted[k]), col=barcol, lwd=barlwd)
    if(k%%2==0){
      text(max(sd.d), d.sorted[k], labels=trait.names.sorted[k], pos=4, srt=0, col="black", font=labelfont, cex=labelsize)
    }
    if(k%%2==1){
      text((-max(sd.d)), d.sorted[k], labels=trait.names.sorted[k], pos=2, srt=0, col="black", font=labelfont, cex=labelsize)}}
}





