[R] tricks

November 17, 2009

Quantitative link strength for APE cophyloplot

Filed under: Plot — Tags: , — Timothée Poisot @ 10:42 pm

Just add a third column with link strength to the association matrix

plotCophylo2 <- function (x, y, assoc = assoc, use.edge.length = use.edge.length,
space = space, length.line = length.line, gap = gap, type = type,
return = return, col = col, show.tip.label = show.tip.label,
font = font)
{
if(ncol(assoc)==2)
{
assoc <- cbind(assoc,rep(1,nrow(assoc)))
}
res <- list()
left <- max(nchar(x$tip.label, type = "width")) + length.line
right <- max(nchar(y$tip.label, type = "width")) + length.line
space.min <- left + right + gap * 2
if ((space <= 0) || (space < space.min))
space <- space.min
N.tip.x <- Ntip(x)
N.tip.y <- Ntip(y)
res$N.tip.x <- N.tip.x
res$N.tip.y <- N.tip.y
a <- plotPhyloCoor(x, use.edge.length = use.edge.length,
type = type)
res$a <- a
b <- plotPhyloCoor(y, use.edge.length = use.edge.length,
direction = "leftwards", type = type)
a[, 2] <- a[, 2] - min(a[, 2])
b[, 2] <- b[, 2] - min(b[, 2])
res$b <- b
b2 <- b
b2[, 1] <- b[1:nrow(b), 1] * (max(a[, 1])/max(b[, 1])) +
space + max(a[, 1])
b2[, 2] <- b[1:nrow(b), 2] * (max(a[, 2])/max(b[, 2]))
res$b2 <- b2
c <- matrix(ncol = 2, nrow = nrow(a) + nrow(b))
c[1:nrow(a), ] <- a[1:nrow(a), ]
c[nrow(a) + 1:nrow(b), 1] <- b2[, 1]
c[nrow(a) + 1:nrow(b), 2] <- b2[, 2]
res$c <- c
plot(c, type = "n", xlim = NULL, ylim = NULL, log = "", main = NULL,
sub = NULL, xlab = NULL, ylab = NULL, ann = FALSE, axes = FALSE,
frame.plot = FALSE)
if (type == "cladogram") {
for (i in 1:(nrow(a) - 1)) segments(a[x$edge[i, 1], 1],
a[x$edge[i, 1], 2], a[x$edge[i, 2], 1], a[x$edge[i,
2], 2])
for (i in 1:(nrow(b) - 1)) segments(b2[y$edge[i, 1],
1], b2[y$edge[i, 1], 2], b2[y$edge[i, 2], 1], b2[y$edge[i,
2], 2])
}
if (type == "phylogram") {
for (i in (N.tip.x + 1):nrow(a)) {
l <- length(x$edge[x$edge[, 1] == i, ][, 1])
for (j in 1:l) {
segments(
a[x$edge[x$edge[, 1] == i, ][1, 1], 1],
a[x$edge[x$edge[, 1] == i, 2], 2][1],
a[x$edge[x$edge[, 1] == i, ][1, 1], 1],
a[x$edge[x$edge[, 1] == i, 2], 2][j]
)
segments(
a[x$edge[x$edge[, 1] == i, ][1, 1], 1],
a[x$edge[x$edge[, 1] == i, 2], 2][j],
a[x$edge[x$edge[, 1] == i, 2], 1][j],
a[x$edge[x$edge[, 1] == i, 2], 2][j]
)
}
}
for (i in (N.tip.y + 1):nrow(b)) {
l <- length(y$edge[y$edge[, 1] == i, ][, 1])
for (j in 1:l) {
segments(b2[y$edge[y$edge[, 1] == i, ][1, 1],
1], b2[y$edge[y$edge[, 1] == i, 2], 2][1],
b2[y$edge[y$edge[, 1] == i, ][1, 1], 1], b2[y$edge[y$edge[,
1] == i, 2], 2][j])
segments(b2[y$edge[y$edge[, 1] == i, ][1, 1],
1], b2[y$edge[y$edge[, 1] == i, 2], 2][j],
b2[y$edge[y$edge[, 1] == i, 2], 1][j], b2[y$edge[y$edge[,
1] == i, 2], 2][j])
}
}
}
if (show.tip.label) {
text(a[1:N.tip.x, ], cex = 0, font = font, pos = 4, labels = x$tip.label)
text(b2[1:N.tip.y, ], cex = 1, font = font, pos = 2,
labels = y$tip.label)
}
lsa <- 1:N.tip.x
lsb <- 1:N.tip.y
decx <- array(nrow(assoc))
decy <- array(nrow(assoc))
for (i in 1:nrow(assoc)) {
if (show.tip.label) {
decx[i] <- strwidth(x$tip.label[lsa[x$tip.label ==
assoc[i, 1]]])
decy[i] <- strwidth(y$tip.label[lsb[y$tip.label ==
assoc[i, 2]]])
}
else {
decx[i] <- decy[i] <- 0
}
segments(a[lsa[x$tip.label == assoc[i, 1]], 1] + decx[i] +
gap, a[lsa[x$tip.label == assoc[i, 1]], 2], a[lsa[x$tip.label ==
assoc[i, 1]], 1] + gap + left, a[lsa[x$tip.label ==
assoc[i, 1]], 2], col = col , lwd = as.numeric(assoc[i,3]))
segments(b2[lsb[y$tip.label == assoc[i, 2]], 1] - (decy[i] +
gap), b2[lsb[y$tip.label == assoc[i, 2]], 2], b2[lsb[y$tip.label ==
assoc[i, 2]], 1] - (gap + right), b2[lsb[y$tip.label ==
assoc[i, 2]], 2], col = col , lwd = as.numeric(assoc[i,3]))
segments(a[lsa[x$tip.label == assoc[i, 1]], 1] + gap +
left, a[lsa[x$tip.label == assoc[i, 1]], 2], b2[lsb[y$tip.label ==
assoc[i, 2]], 1] - (gap + right), b2[lsb[y$tip.label ==
assoc[i, 2]], 2], col = col , lwd = as.numeric(assoc[i,3]))
}
if (return == TRUE)
return(res)
}
Advertisements

October 30, 2009

Decimal log scale on a plot

Filed under: Plot — Tags: , — Timothée Poisot @ 8:30 pm

R only does natural (neperian) log scales by default, and this is lame.

Here is a simple code to do decimal log scale, pretty much a requirement for scientists…

The force(x/y)lim options works for natural and log scales (for the later case, you need to specify the power of 10 that you want : c(-2,2) fixes the limit from 0.01 to 100)

drawlogaxis <- function(side,range)
{
	par(tck=0.02)
#	d <- log(range,10)
	d <- range
	mlog <- floor(min(d))
	Mlog <- ceiling(max(d))
	SeqLog <- c(mlog:Mlog)
	Nlog <- (Mlog-mlog)+1
	axis(side,at=SeqLog,labels=10^SeqLog)
	ats <- log(seq(from=2,to=9,by=1),10)
	mod <- NULL
	for(i in SeqLog)
	{
		mod <- c(mod,rep(i,length(ats)))
	}
	ats <- rep(ats,Nlog)
	ats <- ats+mod
	par(tck=0.02/3)
	axis(side,at=ats,labels=NA)
}

logplot <- function(x,y,log='xy',...,forceylim=c(0,0),forcexlim=c(0,0))
{
	par(tck=0.02)
	xlg <- FALSE
	ylg <- FALSE
	if('x'%in%strsplit(log,'')[[1]]){x <- log(x,10);xlg=TRUE}
	if('y'%in%strsplit(log,'')[[1]]){y <- log(y,10);ylg=TRUE}
	yl <- ifelse(forceylim==c(0,0),range(y),forceylim)
	xl <- ifelse(forcexlim==c(0,0),range(x),forcexlim)
	plot(x,y,...,axes=FALSE,ylim=yl,xlim=xl)
	if(xlg){drawlogaxis(1,xl)}else{axis(1,at=pretty(xl),labels=pretty(xl))}
	if(ylg){drawlogaxis(2,yl)}else{axis(2,at=pretty(yl),labels=pretty(yl))}
	box()
}

addlog <- function(x,y,log='xy',...)
{
	xlg <- FALSE
	ylg <- FALSE
	if('x'%in%strsplit(log,'')[[1]]){x <- log(x,10);xlg=TRUE}
	if('y'%in%strsplit(log,'')[[1]]){y <- log(y,10);ylg=TRUE}
	points(x,y,...)

}

October 26, 2009

MultBar : Advanced multiple barplot with SEM

Filed under: Plot — Tags: — Timothée Poisot @ 8:40 pm

Producing this kind of graphs (below) in R can be a pain in the a*s.

Rplot001

Here is a simple code that requires that data are presented in lists (see the example below).

multbar <- function(list.of.lists,...,condnames=0,pal=colorRampPalette(c('grey','cornsilk')),seriesnames=0,legendpos='topleft',legh=TRUE,do.pty='s')
{
	par(pty=do.pty,mgp=c(1.9,0.8,0),oma=c(0,0,0,0),mar=c(4,3,2,1),bg='transparent',bty='o',tck=0.02,yaxs='i')
	NofList <- length(list.of.lists)
	NofSubList <- length(list.of.lists[[1]])
	if(condnames==0){condnames=c(1:NofList)}
	if(seriesnames==0){seriesnames=c(1:NofSubList)}
	dim.mat.treat <- (NofList+1)*NofSubList
	pos <- c(1:dim.mat.treat)
	pos <- matrix(pos,nrow=(NofList+1))
	nbreaks <- dim.mat.treat+1
	par(xaxs='i',yaxs='i',bty='o')
	plot(0,0,pch=NA,xlim=c(0,nbreaks),ylim=c(min(pretty(range(list.of.lists))),max(pretty(range(list.of.lists)))),xaxt='n',yaxt='n',...)
	calc.ylab <- pretty(range(list.of.lists))
	lab.ylab <- round(calc.ylab,2)
	abline(h=pretty(range(list.of.lists)),lty=3,col='grey',lwd=2)
	abline(h=0,col='darkgrey',lwd=1.4)
	palette <- pal(NofList)
	for(condition in 1:NofList)
	{
		for(treatment in 1:NofSubList)
		{
			CD <- list.of.lists[[condition]][[treatment]]
			CDm <- mean(CD)
			CDv <- sd(CD)
			x.coord <- pos[condition,treatment]
			rect(x.coord,0,x.coord+1,CDm,col=palette[condition])
			arrows(x.coord+0.5,CDm+CDv,x.coord+0.5,CDm-CDv,angle=90,code=3,length=0.05)
		}
	}
	axis(side=2,at=calc.ylab,labels=lab.ylab,tick=TRUE)
	#abline(h=calc.ylab[1])
	abline(h=pretty(range(list.of.lists))[1])
	fills <- c(1:NofList)
	legend(legendpos,legend=condnames,fill=palette[fills],bty='n',horiz=legh,cex=1)
	cl.tab <- (c(1:NofSubList)-1)*NofSubList + NofList/2 +1
	lab.cond <- NULL
	int.lab.cond <- pos[c(1:(nrow(pos)-1)),]
	for(COL in 1:ncol(int.lab.cond))
	{
		lab.cond[COL] <- mean(int.lab.cond[,COL])+0.5
	}
	axis(side=1,at=lab.cond,labels=seriesnames,tick=FALSE)
	par(xaxs='r')
	if(min(pretty(range(list.of.lists)))<0){abline(h=0,lwd=1.1)}
	box()
}

list.1 <- list(rnorm(10,mean=3.4),rnorm(10,mean=6.8))
list.2 <- list(rnorm(10,mean=4),rnorm(10,mean=8.8))
list.3 <- list(rnorm(10,mean=3.1),rnorm(10,mean=7.4))

data <- list(list.1,list.2,list.3)

png()
par(family='serif')
multbar(data,main='Faked data',seriesnames=c('Control','Stress'),condnames=c('Ind. 1','Ind. 2','Ind. 3'),xlab='Treatment',ylab='Response')
dev.off()

October 15, 2009

Filled contour with log-log scale

Filed under: Plot — Tags: , , , — Timothée Poisot @ 8:23 pm

A quick workaround to have a filled.contour plot with natural log10-log10 scale (instead of the default natural log scale)

plotmat <- function(mat,main='',factor='M',MeasuredResponse='Coexistence')
{
X <- as.numeric(rownames(mat))
Y <- as.numeric(colnames(mat))
if(factor=='C')
{
Y <- Y/0.16
}
rownames(mat) <- as.numeric(X)
colnames(mat) <- as.numeric(Y)
colorFun <- colorRampPalette(c("black","darkblue","blue","green",
"orange",'yellow',"red","darkred",'white'))

lX <- log(X, 10)
lY <- log(Y, 10)

pretty.X.at <- pretty(range(lX),n=6)
pretty.X.lab <- round(10^pretty.X.at,0)
pretty.Y.at <- pretty(lY,n=4)
pretty.Y.lab <- round(10^pretty.Y.at,2) pretty.Y.lab[pretty.Y.lab>1] <- round(pretty.Y.lab[pretty.Y.lab>1],0)
pretty.Y.lab[(pretty.Y.lab>0.1)&(pretty.Y.lab<1)] <- round(pretty.Y.lab[(pretty.Y.lab>0.1)&(pretty.Y.lab<1)],1)

filled.contour(lX,lY,mat,
axes=FALSE,
frame.plot=TRUE,
color=colorFun,
ylab= MeasuredResponse,
xlab='Time between perturbation events',
main=main,
key.title=title(main=""),
key.axes=axis(4,at=pretty(vmat)),
plot.axes={ axis(1,at=pretty.X.at,labels=pretty.X.lab)
axis(2,at=pretty.Y.at,labels=pretty.Y.lab) })
}

May 3, 2009

Colouring a 3D plot according to z-values

Filed under: Plot — Tags: , , — Timothée Poisot @ 11:29 pm

persp3dHere is a script that colors a 3D plot according to the z value (height) of each point.

persp.withcol <- function(x,y,z,pal,nb.col,...,xlg=TRUE,ylg=TRUE)
{
	colnames(z) <- y
	rownames(z) <- x

	nrz <- nrow(z)
	ncz <- ncol(z) 

	color <- pal(nb.col)
	zfacet <- z[-1, -1] + z[-1, -ncz] + z[-nrz, -1] + z[-nrz, -ncz]
	facetcol <- cut(zfacet, nb.col)
	par(xlog=xlg,ylog=ylg)
	persp(
		as.numeric(rownames(z)),
		as.numeric(colnames(z)),
		as.matrix(z),
		col=color[facetcol],
		...
		)
}

dd

Create a free website or blog at WordPress.com.