Документ взят из кэша поисковой машины. Адрес оригинального документа : http://herba.msu.ru/shipunov/school/biol_240/en/r/pleiad.r Дата изменения: Sat Apr 2 10:44:26 2016 Дата индексирования: Mon Apr 11 03:52:54 2016 Кодировка:

# R functions to work with correlations
# Some of them are based on the works of Petr Terentjev's (Saint Petersburg) school
# Alexey Shipunov, dactylorhiza@gmail.com
# v. 20160402

#
# Stacks the matrix and selects values which are above the "level"
# Good for the analysis of correlation matrices
#
Topm <- function(X,
level=0.45, # treshold
N=10, # how many to output if values=FALSE
values=TRUE) # if FALSE, ignores "level" and outputs until reaches N
{
X.nam <- dimnames(X)[[1]]
X.col <- dimnames(X)[[2]]
X.rep.g <- rep(X.nam, length(X.col))
X.rep.e <- rep(X.col, each=length(X.nam))
X.vec <- as.vector(X)
X.df <- data.frame(Var1=X.rep.g, Var2=X.rep.e, Value=X.vec)
{if (values)
{X.df <- X.df[abs(X.df$Value) >= level, ]
X.df <- X.df[order(-abs(X.df$Value)), ]}
else
{X.df <- X.df[order(-abs(X.df$Value)), ]
X.df <- X.df[1:N, ]}}
row.names(X.df) <- seq(1, along=X.df$Value)
return(na.omit(X.df))
}

#
# Correlation matrix with p-values (accepts matrix or data frame as an argument)
#
Cor <- function(X, # matrix or data frame with values
stars=TRUE, # replaces p-values with stars if it not greater than "p.level"
p.level=.05, ...)
{
old.options <- options(warn=-1)
nc <- ncol(X)
cor.mat <- matrix(0, nc, nc)
p.mat <- matrix(0, nc, nc)
for (i in 1:nc)
{
for (j in 1:nc)
{
ifelse(i==j,
cor.res <- list(estimate=1, p.value=0),
cor.res <- cor.test(X[,i], X[,j], ...))
cor.mat[i,j] <- cor.res$estimate
p.mat[i,j] <- cor.res$p.value
}
}
cor.mat <- round(cor.mat, 4)
p.mat <- round(p.mat, 4)
if (stars)
{
p <- ifelse(p.mat <= p.level, "*", " ")
sep <- ""
}
else
{
p <- p.mat
sep <- "/"
}
result <- matrix(paste(cor.mat, p, sep=sep), ncol=nc)
result <- gsub("NANA", "NA", result)
dimnames(result) <- list(colnames(X), colnames(X))
return(data.frame(result))
options(old.options)}

#
# Next tree functions are from:
# N. S. Rostova. The variability of correlations systems between the
# morphological characters. 1. Natural populations of Leucanthemum
# vulgare (Asteraceae) // Bot. Journ. 1999. Vol. 84. N 11. P. 50--66
#

#
# Calculates correlation and converts results into the long vector
#
Cor.vec <- function(X, ...) {X.cor <- cor(X, ...)
X.nam <- row.names(X.cor)
X.tri <- as.vector(lower.tri(X.cor))
X.rep.g <- rep(X.nam, length(X.nam))
X.rep.e <- rep(X.nam, each=length(X.nam))
X.pas <- paste(X.rep.g, X.rep.e, sep=" & ")
X.vec <- as.vector(X.cor)[X.tri]
names(X.vec) <- X.pas[X.tri]
return(X.vec)}

#
# Calculates multiple correlation matrices (via GROUP) and stacks them together
# Suitable for PCA
#
Rostova.tbl <- function(X, GROUP, ...) {
old.options <- options(warn=-1)
r.table <- NULL
r.names <- unique(X[[GROUP]])
for (i in r.names)
r.table <- cbind (r.table, cor.vec(subset(X, X[[GROUP]]==i)[,-GROUP], ...))
dimnames(r.table)[[2]] <- r.names
r.table[is.na(r.table)] <- 0
r.table <- t(r.table)
return(r.table)
options(old.options)}

#
# Coefficient of determination for correlation matrix: allows to compare various correlation structures
#
Coeff.det <- function (X, ...) {X.cor <- cor(X, ...) # X is matrix or data frame with values
X.det <- NULL
X.dim <- dimnames(X.cor)[[2]]
for (i in 1:length(X.dim)) X.det <- c(X.det, mean(X.cor[,i]^2))
names(X.det) <- X.dim
return(X.det)}

#
# Correlation circles (correlation pleiads)
# Depends on Recode()
Recode <- function(var, from, to){
x <- as.vector(var)
x.tmp <- x
for (i in 1:length(from)){x <- replace(x, x.tmp == from[i], to[i])}
if(is.factor(var)) factor(x) else x}
# Validations are rudimentary so be careful playing with options
#
Pleiad <- function(table,
corr=FALSE, # if TRUE, uses absolute values istead of real -- good for correlation matrix
dist=FALSE, # if TRUE, converts distance matrix to the data frame -- good for dist() and vegdist() outputs
treshold=FALSE, # if this is (saying) =.5, selects for plotting (only one lty) those values which are >.5
circle=1, # line type for the cirle, set it to 0 for no cicrle
breaks=5, # how many groups you want, should be concerted with "lty"'s and "lwd"'s
lwd=c(1,1,1,2,3), # line widths to use, should be concerted with "breaks" and "lty"'s
lty=c(3,2,1,1,1), # line types to use, should be concerted with "breaks" and "lwd"'s
abbr=-1, # if =-1, no abbreviation; if =0, no labels; other values use abbreviate()
lbltext="internal", # if this is a vector starting from something else, will replace dimnames
off=1.09, # radial offset of labels, be careful here
legend=TRUE, # if FALSE, no legend
legtext=1, # if =1 then "weaker ... stronger"; if =2, shows cutting intervals; if =3, then 1:5; if >3, issues error
legpos="topright", # this is from legend()
leghoriz=FALSE, ...) # equal to horiz= from legend()
{
if (dist) table <- as.matrix(table)
table <- data.frame(table)
ddu <- unique(unlist(dimnames(table)))
#
ddc <- t(combn(ddu, 2))
ddn <- apply(ddc, 1, function(.x) {.y <- table[.x[1],.x[2]]; ifelse(is.null(.y),NA,.y)}) # fragile, depends on strange features of data.frame!
ddn[is.na(ddn)] <- 0
if (corr) ddn <- abs(ddn)
#
x <- sin(seq(0, 2*pi, length.out=length(ddu)+1))
y <- cos(seq(0, 2*pi, length.out=length(ddu)+1))
#
fromx <- as.numeric(Recode(ddc[ddn!=0,1], ddu, x))
fromy <- as.numeric(Recode(ddc[ddn!=0,1], ddu, y))
tox <- as.numeric(Recode(ddc[ddn!=0,2], ddu, x))
toy <- as.numeric(Recode(ddc[ddn!=0,2], ddu, y))
#
seglwd <- cbind(1:breaks, lwd)
seglty <- cbind(1:breaks, lty)
if (corr) segcut <- cut(ddn[ddn!=0], seq(0, 1, length.out=(breaks+1))) else segcut <- cut(ddn[ddn!=0], breaks, dig.lab=1)
segp <- as.numeric(segcut)
segpt <- Recode(segp, seglty[,1], seglty[,2])
segpw <- Recode(segp, seglwd[,1], seglwd[,2])
if (treshold) segpw <- segpt <- (ddn > treshold) * 1
#
if (abbr > -1) dda <- abbreviate(ddu, abbr) else dda <- ddu
seglev <- sub("(","", levels(segcut), fixed=TRUE)
seglev <- sub("]","", seglev, fixed=TRUE)
seglev <- sub(","," - ", seglev, fixed=TRUE)
legtxtable <- cbind(c("weaker","","","","stronger"), seglev, 1:5)
legtxt <- legtxtable[,legtext]
if (treshold)
{
lty <- lwd <- 1
legtxt <- paste(">",treshold)
}
oldpar <- par(mar=c(0,0,0,0))
plot(x, y, xlim=c(-1,1)*1.15, ylim=c(-1,1)*1.15, axes=FALSE, type="n", ...)
polygon(sin(seq(0, 2*pi, length.out=100)), cos(seq(0, 2*pi, length.out=100)), lty=circle)
segments(fromx, fromy, tox, toy, lwd=segpw, lty=segpt)
points(x, y, ...)
if (lbltext[1]=="internal") lbltxt <- dda else lbltxt <- lbltext
text(x[-length(x)]*off, y[-length(y)]*off, labels=lbltxt)
if (legend) legend(legpos, horiz=leghoriz, lty=lty, lwd=lwd, legend=legtxt, bty="n", seg.len=1.2)
par(oldpar)
invisible(data.frame(x=x[-length(x)], y=y[-length(y)], row.names=lbltxt))
} # retunts the data frame with position of points, helps in subsequent plot enhancing