Re: [R] biplot package

From: Mark Difford <mark_difford_at_yahoo.co.uk>
Date: Tue, 05 Jun 2007 01:38:58 -0700 (PDT)

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.

______________________________________________
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.
Received on Tue 05 Jun 2007 - 08:47:26 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Tue 05 Jun 2007 - 09:31:30 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help. Please read the posting guide before posting to the list.