[R] biplot package
joseclaudio.faria
joseclaudio.faria at terra.com.br
Tue Jun 5 03:25:06 CEST 2007
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)
More information about the R-help
mailing list