Dears,
I've been learning biplot (Gabriel, 1971) and I found the function
'biplot', inside of the package 'stats',
useful but, a bit limited.
So, I'm thinking to start a colaborative package to enhance this methods to
other multivariate methods. In this
way, I would like to start it, making public a new function (biplot.pca, still
in development, but running)
that make biplot more simple and power.
All users are free to modify and make it better.
Below the function and a small script to learn it.
#==============================================================================#
Name : biplot.pca
# Author : Jos? Cl?udio Faria (DCET/USC/BRAZIL)
# Date (dd/mm/yy): 4/6/2007 08:27:54
# Version : v3
# Aim : 2D and 3D (under scaterplot3d and rgl packages) PCA biplot
# Mail : joseclaudio.faria em terra.com.br
#==============================================================================#
Arguments:
# x Data (frame or matrix).
# center Either a logical value or a numeric vector of length equal
# to the number of columns of x (TRUE is the default).
# scale Either a logical value or a numeric vector of length equal
# to the number of columns of x (FALSE is the default).
# weight Way of factorize ('equal' is the default).
# plot Logical to produce or not a graphical representation of
# biplot (TRUE is the default).
# rgl.use If TRUE the 3D scatter will be under the rgl environment, in
# another way the scatterplot3d will be used.
# aspect3d Apparent ratios of the x, y, and z axes of the bounding box
# clear3d Logical to clear or not a 3D graphical representation of
# biplot before to make a new (TRUE is the default).
# simple.axes Whether to draw simple axes (TRUE or FALSE).
# box Whether to draw a box (the default is FALSE).
# spheres Logical to represent objects as spheres (the default is FALSE).
# sphere.factor Relative size factor of spheres representing points; the
# default size is dependent on the scale of observations.
# col.obj Color of spheres or labels of objects.
# col.var Color of lines and labels of variables.
# var.red Factor of reduction of the length of the lines of variables.
# graphical variables representation (<=1, 1 is the default).
# cex Character expansion.
biplot.pca = function (x,
n.values=2,
center=T,
scale=F,
weight=c('equal', 'samples',
'variables'),
plot=T,
rgl.use=T,
aspect3d=c(1, 1, 1),
clear3d=T,
simple.axes=T,
box=F,
spheres=T,
sphere.factor=1,
col.obj=1,
col.var=2,
var.red=1,
cex=.6 )
{
x = as.matrix(x)
x = scale(x, center=center, scale=scale)
if(is.null(rownames(x))) rownames = 1:nrow(x) else rownames = rownames(x)
if(is.null(colnames(x))) colnames = paste('V', 1:ncol(x),
sep='') else colnames = colnames(x)
s = svd(x)
s2 = diag(sqrt(s$d[1:n.values]))
#s2 = diag(s$d[1:n.values]) pca of pcurve is like this!?
switch(match.arg(weight),
equal = {
g = s$u[,1:n.values] %*% s2
h = s2 %*% t(s$v[,1:n.values])
hl = t(h)
},
samples = {
g = s$u[,1:n.values] %*% s2
h = t(s$v[,1:n.values])
hl = t(h)
},
variables = {
g = s$u[,1:n.values]
h = s2 %*% t(s$v[,1:n.values])
hl = t(h)
})
gencolnames = paste('PC', 1:n.values, sep='')
rownames(g) = rownames
colnames(g) = gencolnames
rownames(hl) = colnames
colnames(hl) = gencolnames
coo = rbind(g, hl)
rownames(coo) = c(rownames, colnames)
colnames(coo) = gencolnames
cooplot = rbind(g, hl*var.red)
cooplot = rbind(cooplot, rep(0, n.values)) # to correct visualization
if(plot) {
if(n.values == 2) {
plot(cooplot,
xlab='PC1', ylab='PC2',
type='n')
text(x=g[,1], y=g[,2],
labels=rownames(g),
cex=cex, col=col.obj)
arrows(x0=0, y0=0,
x1=hl[,1]*var.red, y1=hl[,2]*var.red,
length=0.1, angle=20,
col=col.var)
text(x=hl[,1]*var.red, y=hl[,2]*var.red,
labels = rownames(hl),
cex=cex, col=col.var)
}
if(n.values == 3) {
if (rgl.use) {
require(rgl)
require(mgcv)
size = max(g)/20 * sphere.factor
if (clear3d) clear3d()
if (spheres)
spheres3d(g, col=col.obj, radius=size, alpha=.5)
else
text3d(g, texts=rownames(g), col=col.obj, alpha=.5)
aspect3d(aspect3d)
for(i in 1:nrow(hl)) {
segments3d(rbind(matrix(0, nc=3),
hl[i,]*var.red),
col=col.var)
}
text3d(hl*var.red,
texts=rownames(hl),
col=col.var)
if(simple.axes) {
axes3d(c('x', 'y', 'z'))
}
else
decorate3d(xlab = 'PC1', ylab = 'PC2', zlab =
'PC3', box = box)
title3d(xlab = 'PC1', ylab = 'PC2', zlab =
'PC3')
} else {
require(scatterplot3d)
graph = scatterplot3d(cooplot,
type = if(spheres) 'p' else 'n',
xlab='PC1', ylab='PC2',
zlab='PC3',
grid=F,
box=box,
cex.symbols=cex,
color=col.obj,
pch=20)
if(!spheres)
text(graph$xyz.convert(g),
labels=rownames(g),
col=col.obj, cex=cex)
for(i in 1:nrow(hl)) {
graph$points3d(c(0, hl[i,1]*var.red),
c(0, hl[i,2]*var.red),
c(0, hl[i,3]*var.red),
type='l', col=col.var)
}
text(graph$xyz.convert(hl*var.red),
labels=rownames(hl),
col=col.var, cex=cex)
}
}
}
rlist = list(values=s$d,
objects=g,
variables=hl,
all=coo)
}
#==============================================================================#
Name : biplot.pca_test
# Author : Jos? Cl?udio Faria (DCET/USC/BRAZIL)
# Date (dd/mm/yy): 4/6/2007 08:27:54
# Version : v3
# Aim : to learn and to test the new 'biplot.pca' function
# Mail : joseclaudio.faria em terra.com.br
#==============================================================================
#mtrace(biplot.pca, T)
#mtrace(biplot.pca, F)
#???????????????????
# 2D with graphics package
#???????????????????
#x = matrix(c(42, 52, 48, 58, 4, 5, 4, 3), nc=2); x
#dimnames(x) = list(letters[1:nrow(x)], LETTERS[1:ncol(x)]); x
#x = stackloss; x
#x = cars; x
#x = longley; x
x = mtcars[,1:7]; x
#x = LifeCycleSavings; x
biplot.pca(x)
biplot.pca(x, scale=T)
biplot.pca(x, col.obj=3, col.var=4, var.red=.5)
biplot.pca(x, center=T, scale=F, weight='eq', cex=.5)
biplot.pca(x, center=T, scale=F, weight='eq', cex=.8)
biplot.pca(x, center=T, scale=F, weight='sa')
biplot.pca(x, center=T, scale=F, weight='va')
biplot.pca(x, center=T, scale=F, weight='va', var.red=.05)
#??????????????????????
# 3D with scatterplot3d package
#??????????????????????
x = stackloss; x
biplot.pca(x, n.values=3, rgl.use=F, cex=.5)
biplot.pca(x, n.values=3, rgl.use=F, spheres=F, simple.axes=F, box=T)
biplot.pca(x, n.values=3, rgl.use=F, spheres=F, col.obj=3, col.var=4,
var.red=.5)
biplot.pca(x, n.values=3, rgl.use=F, center=T, scale=F, weight='eq')
biplot.pca(x, n.values=3, rgl.use=F, center=T, scale=T, weight='eq')
biplot.pca(x, n.values=3, rgl.use=F, center=T, scale=T, weight='sa')
biplot.pca(x, n.values=3, rgl.use=F, spheres=F, center=T, scale=T,
weight='va')
biplot.pca(x, n.values=3, rgl.use=F, center=T, scale=T, weight='va',
var.red=.5)
#???????????????
# 3D with rgl package
#???????????????
x = iris[1:4]
#x = stackloss
x = LifeCycleSavings; x
clear3d()
rgl.bringtotop(stay=T)
biplot.pca(x, n.values=3, spheres=F)
rgl.bringtotop(stay=T)
biplot.pca(x, n.values=3, spheres=F, simple.axes=F, box=T, aspect3d=c(1, 1, 2))
rgl.bringtotop(stay=T)
biplot.pca(x, n.values=3, spheres=F, col.obj=3, col.var=4, var.red=.5)
rgl.bringtotop(stay=T)
biplot.pca(x, n.values=3, center=T, scale=F, weight='eq', plot=T)
rgl.bringtotop(stay=T)
biplot.pca(x, n.values=3, center=T, scale=T, weight='eq', plot=T)
rgl.bringtotop(stay=T)
biplot.pca(x, n.values=3, spheres=F, center=T, scale=T, weight='sa')
rgl.bringtotop(stay=T)
biplot.pca(x, n.values=3, center=T, scale=T, weight='va')
rgl.bringtotop(stay=T)
biplot.pca(x, n.values=3, center=T, scale=T, weight='va', var.red=.3)
Best regards,
Jose Claudio Faria
Estat?stica Experimental - Prof. Titular
Universidade Estadual de Santa Cruz - UESC
Departamento de Ciencias Exatas e Tecnologicas - DCET
Bahia - Brasil
Tels:
73-3634.2779 (fixo Ilheus)
19-9144.8979 (celular Piracicaba)
Hi Jose, Good idea. I haven't yet run your code, but it might be a good idea to take a look at the calibrate package (unfortunately not "upgraded" since its first release), and at what Chessel and his group have done in package ade4, as well as what Jari Oksanen & his co-authors have done in package vegan. There are also some interesting implementations of in package psy, and in the compositions package. Best Regards, Mark Difford. Jose Claudio Faria wrote:> > Dears, > > I've been learning biplot (Gabriel, 1971) and I found the function > 'biplot', inside of the package 'stats', > useful but, a bit limited. > > So, I'm thinking to start a colaborative package to enhance this methods > to other multivariate methods. In this > way, I would like to start it, making public a new function (biplot.pca, > still in development, but running) > that make biplot more simple and power. > > All users are free to modify and make it better. > Below the function and a small script to learn it. > > #==============================================================================> # Name : biplot.pca > # Author : Jos? Cl?udio Faria (DCET/USC/BRAZIL) > # Date (dd/mm/yy): 4/6/2007 08:27:54 > # Version : v3 > # Aim : 2D and 3D (under scaterplot3d and rgl packages) PCA > biplot > # Mail : joseclaudio.faria at terra.com.br > #==============================================================================> # Arguments: > # x Data (frame or matrix). > # center Either a logical value or a numeric vector of length equal > # to the number of columns of x (TRUE is the default). > # scale Either a logical value or a numeric vector of length equal > # to the number of columns of x (FALSE is the default). > # weight Way of factorize ('equal' is the default). > # plot Logical to produce or not a graphical representation of > # biplot (TRUE is the default). > # rgl.use If TRUE the 3D scatter will be under the rgl environment, > in > # another way the scatterplot3d will be used. > # aspect3d Apparent ratios of the x, y, and z axes of the bounding > box > # clear3d Logical to clear or not a 3D graphical representation of > # biplot before to make a new (TRUE is the default). > # simple.axes Whether to draw simple axes (TRUE or FALSE). > # box Whether to draw a box (the default is FALSE). > # spheres Logical to represent objects as spheres (the default is > FALSE). > # sphere.factor Relative size factor of spheres representing points; the > # default size is dependent on the scale of observations. > # col.obj Color of spheres or labels of objects. > # col.var Color of lines and labels of variables. > # var.red Factor of reduction of the length of the lines of > variables. > # graphical variables representation (<=1, 1 is the > default). > # cex Character expansion. > > biplot.pca = function (x, > n.values=2, > center=T, > scale=F, > weight=c('equal', 'samples', 'variables'), > plot=T, > rgl.use=T, > aspect3d=c(1, 1, 1), > clear3d=T, > simple.axes=T, > box=F, > spheres=T, > sphere.factor=1, > col.obj=1, > col.var=2, > var.red=1, > cex=.6 ) > { > x = as.matrix(x) > x = scale(x, center=center, scale=scale) > if(is.null(rownames(x))) rownames = 1:nrow(x) else rownames > rownames(x) > if(is.null(colnames(x))) colnames = paste('V', 1:ncol(x), sep='') else > colnames = colnames(x) > s = svd(x) > s2 = diag(sqrt(s$d[1:n.values])) > #s2 = diag(s$d[1:n.values]) pca of pcurve is like this!? > switch(match.arg(weight), > equal = { > g = s$u[,1:n.values] %*% s2 > h = s2 %*% t(s$v[,1:n.values]) > hl = t(h) > }, > samples = { > g = s$u[,1:n.values] %*% s2 > h = t(s$v[,1:n.values]) > hl = t(h) > }, > variables = { > g = s$u[,1:n.values] > h = s2 %*% t(s$v[,1:n.values]) > hl = t(h) > }) > gencolnames = paste('PC', 1:n.values, sep='') > rownames(g) = rownames > colnames(g) = gencolnames > rownames(hl) = colnames > colnames(hl) = gencolnames > coo = rbind(g, hl) > rownames(coo) = c(rownames, colnames) > colnames(coo) = gencolnames > cooplot = rbind(g, hl*var.red) > cooplot = rbind(cooplot, rep(0, n.values)) # to correct > visualization > if(plot) { > if(n.values == 2) { > plot(cooplot, > xlab='PC1', ylab='PC2', > type='n') > text(x=g[,1], y=g[,2], > labels=rownames(g), > cex=cex, col=col.obj) > arrows(x0=0, y0=0, > x1=hl[,1]*var.red, y1=hl[,2]*var.red, > length=0.1, angle=20, > col=col.var) > text(x=hl[,1]*var.red, y=hl[,2]*var.red, > labels = rownames(hl), > cex=cex, col=col.var) > } > if(n.values == 3) { > if (rgl.use) { > require(rgl) > require(mgcv) > size = max(g)/20 * sphere.factor > if (clear3d) clear3d() > if (spheres) > spheres3d(g, col=col.obj, radius=size, alpha=.5) > else > text3d(g, texts=rownames(g), col=col.obj, alpha=.5) > aspect3d(aspect3d) > for(i in 1:nrow(hl)) { > segments3d(rbind(matrix(0, nc=3), > hl[i,]*var.red), > col=col.var) > } > text3d(hl*var.red, > texts=rownames(hl), > col=col.var) > if(simple.axes) { > axes3d(c('x', 'y', 'z')) > } > else > decorate3d(xlab = 'PC1', ylab = 'PC2', zlab = 'PC3', box = box) > title3d(xlab = 'PC1', ylab = 'PC2', zlab = 'PC3') > } else { > require(scatterplot3d) > graph = scatterplot3d(cooplot, > type = if(spheres) 'p' else 'n', > xlab='PC1', ylab='PC2', zlab='PC3', > grid=F, > box=box, > cex.symbols=cex, > color=col.obj, > pch=20) > if(!spheres) > text(graph$xyz.convert(g), > labels=rownames(g), > col=col.obj, cex=cex) > for(i in 1:nrow(hl)) { > graph$points3d(c(0, hl[i,1]*var.red), > c(0, hl[i,2]*var.red), > c(0, hl[i,3]*var.red), > type='l', col=col.var) > } > text(graph$xyz.convert(hl*var.red), > labels=rownames(hl), > col=col.var, cex=cex) > } > } > } > rlist = list(values=s$d, > objects=g, > variables=hl, > all=coo) > } > > > #==============================================================================> # Name : biplot.pca_test > # Author : Jos? Cl?udio Faria (DCET/USC/BRAZIL) > # Date (dd/mm/yy): 4/6/2007 08:27:54 > # Version : v3 > # Aim : to learn and to test the new 'biplot.pca' function > # Mail : joseclaudio.faria at terra.com.br > #==============================================================================> > #mtrace(biplot.pca, T) > #mtrace(biplot.pca, F) > > #??????????????????? > # 2D with graphics package > #??????????????????? > #x = matrix(c(42, 52, 48, 58, 4, 5, 4, 3), nc=2); x > #dimnames(x) = list(letters[1:nrow(x)], LETTERS[1:ncol(x)]); x > #x = stackloss; x > #x = cars; x > #x = longley; x > x = mtcars[,1:7]; x > #x = LifeCycleSavings; x > biplot.pca(x) > biplot.pca(x, scale=T) > biplot.pca(x, col.obj=3, col.var=4, var.red=.5) > biplot.pca(x, center=T, scale=F, weight='eq', cex=.5) > biplot.pca(x, center=T, scale=F, weight='eq', cex=.8) > biplot.pca(x, center=T, scale=F, weight='sa') > biplot.pca(x, center=T, scale=F, weight='va') > biplot.pca(x, center=T, scale=F, weight='va', var.red=.05) > > #?????????????????????? > # 3D with scatterplot3d package > #?????????????????????? > x = stackloss; x > biplot.pca(x, n.values=3, rgl.use=F, cex=.5) > biplot.pca(x, n.values=3, rgl.use=F, spheres=F, simple.axes=F, box=T) > biplot.pca(x, n.values=3, rgl.use=F, spheres=F, col.obj=3, col.var=4, > var.red=.5) > biplot.pca(x, n.values=3, rgl.use=F, center=T, scale=F, weight='eq') > biplot.pca(x, n.values=3, rgl.use=F, center=T, scale=T, weight='eq') > biplot.pca(x, n.values=3, rgl.use=F, center=T, scale=T, weight='sa') > biplot.pca(x, n.values=3, rgl.use=F, spheres=F, center=T, scale=T, > weight='va') > biplot.pca(x, n.values=3, rgl.use=F, center=T, scale=T, weight='va', > var.red=.5) > > #??????????????? > # 3D with rgl package > #??????????????? > x = iris[1:4] > #x = stackloss > x = LifeCycleSavings; x > > clear3d() > rgl.bringtotop(stay=T) > biplot.pca(x, n.values=3, spheres=F) > rgl.bringtotop(stay=T) > biplot.pca(x, n.values=3, spheres=F, simple.axes=F, box=T, aspect3d=c(1, > 1, 2)) > rgl.bringtotop(stay=T) > biplot.pca(x, n.values=3, spheres=F, col.obj=3, col.var=4, var.red=.5) > rgl.bringtotop(stay=T) > biplot.pca(x, n.values=3, center=T, scale=F, weight='eq', plot=T) > rgl.bringtotop(stay=T) > biplot.pca(x, n.values=3, center=T, scale=T, weight='eq', plot=T) > rgl.bringtotop(stay=T) > biplot.pca(x, n.values=3, spheres=F, center=T, scale=T, weight='sa') > rgl.bringtotop(stay=T) > biplot.pca(x, n.values=3, center=T, scale=T, weight='va') > rgl.bringtotop(stay=T) > biplot.pca(x, n.values=3, center=T, scale=T, weight='va', var.red=.3) > > Best regards, > > Jose Claudio Faria > Estat?stica Experimental - Prof. Titular > Universidade Estadual de Santa Cruz - UESC > Departamento de Ciencias Exatas e Tecnologicas - DCET > Bahia - Brasil > Tels: > 73-3634.2779 (fixo Ilheus) > 19-9144.8979 (celular Piracicaba) > > ______________________________________________ > R-help at stat.math.ethz.ch mailing list > https://stat.ethz.ch/mailman/listinfo/r-help > PLEASE do read the posting guide > http://www.R-project.org/posting-guide.html > and provide commented, minimal, self-contained, reproducible code. > >-- View this message in context: http://www.nabble.com/biplot-package-tf3869013.html#a10965497 Sent from the R help mailing list archive at Nabble.com.
joseclaudio.faria <joseclaudio.faria <at> terra.com.br> writes:> > Dears, > > I've been learning biplot (Gabriel, 1971) and I found the function 'biplot',inside of the package 'stats',> useful but, a bit limited. > > So, I'm thinking to start a colaborative package to enhance this methods toother multivariate methods. In this> way, I would like to start it, making public a new function (biplot.pca, stillin development, but running)> that make biplot more simple and power. > > All users are free to modify and make it better. > Below the function and a small script to learn it. >Dear Jose Claudio Faria, Looks nice. However, I'd suggets you articulate this into R core and tradition. There the biplot function is now defined as: biplot <- function (x, ...) UseMethod("biplot") with methods("biplot") [1] biplot.default* biplot.prcomp* biplot.princomp* You function is now called biplot.pca which sounds like biplot method for the class "pca" (which, I think, exists in labdsv and perhaps somewhere else). What you do is to propose a biplot method for class-less function svd. Or perhaps this could be something like biplot.data.frame so that this function is launched when user supplies as a data.frame as the first argument for biplot. The good old R tradition is to separate plotting from calculation so that you can have these separately (which obviously is related to the unix toolbox thinking: do one thing well, and pipe the result to the next function). It might make sense to have a biplot method for data.frame (if I remember correctly, some of the core members fancied a spanning tree method for data.frame which is along similar lines). However, I would like to have enhanced biplot methods for other methods as well, such as analyses using prcomp or princomp, or multivariate analyses in packages. There is no sense to replicate all mv analyses within biplot functions, but you should just try cope with the outputs of various methods. Then you perhaps have to change the name so that you can have biplot.prcomp alongside with your new hyperbiplot.prcomp. (Namespace is an answer to no real problem, but a reason of many new problems.) Somebody already promised to write a biplot method for rda in vegan (howdy Gav), but I haven't heard of this for a long time. It would be nice to be able to have this kind of interface for your enhanced code, too. cheers, jari oksanen