#'PC Score Correlation Test
#'
#'Given a set of covariance matrices and means for terminals, test the hypothesis
#'that observed divergence is larger/smaller than expected by drift alone using the correlation on
#'principal component scores.
#'
#'@param means list or array of species means being compared. array must have means in the rows.
#'@param cov.matrix ancestral covariance matrix for all populations
#'@param taxons names of taxons being compared. Must be in the same order of the means.
#'@param show.plots Logical. If TRUE, plot of eigenvalues of ancestral matrix by between group variance is showed.
#'@return list of results containing:
#'@return correlation matrix of principal component scores and p.values for each correlation. Lower triangle of output are correlations, and upper triangle are p.values.
#'@return if show.plots is TRUE, also returns a list of plots of all projections of the nth PCs, where n is the number of taxons.
#'@references Marroig, G., and Cheverud, J. M. (2004). Did natural selection or genetic drift
#'produce the cranial diversification of neotropical monkeys? The American Naturalist, 163(3), 417-428. doi:10.1086/381693
#'@author Ana Paula Assis, Diogo Melo
#'@export
#'@import plyr
#'@importFrom ggplot2 ggplot geom_text geom_smooth labs theme_bw aes_string
#'@examples
#' #Input can be an array with means in each row or a list of mean vectors
#'means = array(rnorm(40*10), c(10, 40))
#'cov.matrix = RandomMatrix(40, 1, 1, 10)
#'taxons = LETTERS[1:10]
#'PCScoreCorrelation(means, cov.matrix, taxons)
#'
#'\dontrun{
#'##Plots list can be displayed using plot_grid()
#'library(cowplot)
#'pc.score.output <- PCScoreCorrelation(means, cov.matrix, taxons, TRUE)
#'plot_grid(plotlist = pc.score.output$plots)
#'}
PCScoreCorrelation <- function(means, cov.matrix, taxons = names(means), show.plots = FALSE){
if(is.data.frame(means) | (!is.array(means) & !is.list(means)))
stop("means must be in a list or an array.")
if(!isSymmetric(cov.matrix)) stop("covariance matrix must be symmetric.")
if(is.list(means)){
mean.array <- laply(means, identity)
} else{
mean.array <- means
if(is.null(taxons)) taxons <- rownames(means)
}
if(is.null(taxons)) taxons <- paste0("taxon", 1:nrow(mean.array))
n.taxon <- length(taxons)-1
if(n.taxon > dim(cov.matrix)[1]) n.taxon <- dim(cov.matrix)[1]
W.pc <- eigen(cov.matrix)
proj.med <- as.matrix(mean.array) %*% W.pc$vectors
colnames(proj.med) <- paste0("PC", 1:ncol(proj.med))
mat.pcs.deriva <- matrix(1 , ncol = n.taxon, nrow = n.taxon)
if(show.plots){
plots <- vector("list", (n.taxon^2 - n.taxon)/2)
current.plot <- 1
}
x = NULL
y = NULL
for(i in 1:n.taxon){
for (j in 1:i){
if(j != i){
test = cor.test(proj.med[,i], proj.med[,j])
mat.pcs.deriva[i,j] <- test$estimate
mat.pcs.deriva [j,i] <- test$p.value
if(show.plots){
plots[[current.plot]] <- ggplot(data.frame(x = proj.med[,i],
y = proj.med[,j],
taxons = taxons),
aes(x, y)) +
geom_text(aes(label = taxons)) +
geom_smooth(method = "lm", color = 'black') +
labs(x = paste0("PC", i), y = paste0("PC", j)) +
theme_bw()
current.plot <- current.plot + 1
}
}
}
}
if(show.plots) return(list("correlation_p.value" = mat.pcs.deriva, "plots" = plots))
else return(mat.pcs.deriva) ## a diagonal inferior sao as correlacoes e a superior o valor de p
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.