## ================================================================ g.biplot <- function(object, cmp=FALSE, cor=NULL, plot=TRUE, choices=1:2, circle=cor, circlelwd=1, cextext=0.8, col=1:2, ...) { ## Purpose: biplot ## ---------------------------------------------------------------------- ## Arguments: ## object data or prcomp or princomp or factanal object ## cmp if TRUE, column metric preserving biplot is produced ## cor if TRUE, data will be scaled (princomp on correlation mx.) ## plot if FALSE, no plotting will be done, but a "biplot" object ## will be returned ## choices components to be plotted (do not change) ## circle if TRUE and biplot is done for scaled data, ## a circle (for CMP) or ellipse (for RMP biplot) will be drawn ## circlelwd line width for circle ## cextext character expansion for labels ## col colors for points and arrows ## ... further arguments to be passed to prcomp or princomp or ## biplot for factor analysis, like ## xlabs, ylabs labels for x and y elements ## xcol, ycol colors ... ## ---------------------------------------------------------------------- ## Author: Werner Stahel, Date: 23 Apr 2003 / Feb 2005 ## if object is data, call princomp if ((!is.null(dim(object)))||is.call(object)) {# is.formula does not exist if (is.null(cor)) cor <- TRUE object <- prcomp(object, scale=cor, scores=TRUE, ...) } else if (is.null(cor)) { lscl <- object$scale cor <- if (length(lscl)>0) { if (is.logical(lscl)) lscl else any(lscl!=1) } else FALSE } lsc <- object$scores if (length(lsc)==0) lsc <- object$x if (length(lsc)==0) stop(">g.biplot> no scores found") lld <- object$loadings if (length(lld)==0) lld <- object$rotation if (length(lld)==0) stop(">g.biplot> no loadings found") lsc <- lsc[,choices] lld <- lld[,choices] lsd <- object$sdev[choices] ## change signs and, if cmp, scale lfac <- 1-2*(apply(lld,2,mean)<0) if (cmp) { if (length(lsd)==0) warning(">g.biplot> if cmp=T, an sdev component is needed") else lfac <- lfac*lsd } object$loadings <- lld <- sweep(lld,2,lfac,"*") object$scores <- lsc <- sweep(lsc,2,lfac,"/") ## plotting col <- rep(col,length=2) cextext <- rep(cextext,length=2) lusr <- g.biplot.default(lsc,lld, expand=0.8,cex=cextext,col=col, ...) if (circle) { la <- seq(0,2*pi,length=73) lx <- cos(la) ly <- sin(la) if (!cmp) { if (length(lsd)<2) warning(">g.biplot> if circle=T, an sdev component is needed") else lines(lx/lsd[1],ly/lsd[2],lty=3, col = col[2], lwd=circlelwd) } else lines(lx,ly,lty=3, col = col[2], lwd=circlelwd) } object$usr <- lusr invisible(object) } ## ================================================================ g.biplot.default <- function (x, y, var.axes = TRUE, col, cex = rep(par("cex"), 2), xlabs = NULL, ylabs = NULL, xcol = col[1], ycol = col[2], expand = 1, xlim = NULL, ylim = NULL, arrow.len = 0.1, main = NULL, sub = NULL, xlab = NULL, ylab = NULL, labelexpand = 1.1, ...) { ## slightly modified version of biplot.default n <- nrow(x) p <- nrow(y) if (missing(xlabs)|is.null(xlabs)) xlabs <- dimnames(x)[[1]] if (is.null(xlabs)) xlabs <- as.character(1:n) xlabs <- rep(xlabs,length=n) dimnames(x) <- list(xlabs, dimnames(x)[[2]]) if (missing(ylabs)) { ylabs <- dimnames(y)[[1]] if (is.null(ylabs)) ylabs <- paste("Var", 1:p) } ylabs <- rep(as.character(ylabs),length=p) # !!! dimnames(y) <- list(ylabs, dimnames(y)[[2]]) if (length(cex) == 1) cex <- c(cex, cex) if (missing(col)) { col <- par("col") if (!is.numeric(col)) col <- match(col, palette()) col <- c(col, col + 1) } else if (length(col) == 1) col <- c(col, col) ## plot ranges unsigned.range <- function(x) c(-abs(min(x)), abs(max(x))) rangx1 <- unsigned.range(x[, 1]) rangx2 <- unsigned.range(x[, 2]) rangy1 <- unsigned.range(y[, 1]) rangy2 <- unsigned.range(y[, 2]) if (is.null(xlim) && is.null(ylim)) xlim <- ylim <- rangx1 <- rangx2 <- range(rangx1, rangx2) if (is.null(xlim)) xlim <- rangx1 if (is.null(ylim)) ylim <- rangx2 ratio <- max(rangy1/rangx1, rangy2/rangx2)/expand on.exit(par(op)) op <- par(pty = "s") if (!is.null(main)) op <- c(op, par(mar = par("mar") + c(0, 0, 1, 0))) ## do plot plot(x, axes = FALSE, type = "n", xlim = xlim, ylim = ylim, col = xcol, xlab = xlab, ylab = ylab, sub = sub, main = main, ...) axis(1, col = col[1]) axis(2, col = col[1]) if (is.character(xlabs)) text(x, labels = xlabs, cex = cex[1], col = xcol, ...) else points(x, pch = xlabs, cex = cex[1], col = xcol, ...) usrr <- par("usr") par(new = TRUE) plot(y, axes = FALSE, type = "n", xlim = xlim * ratio, ylim = ylim * ratio, xlab = "", ylab = "", col = ycol, ...) axis(3, col = col[2]) axis(4, col = col[2]) box() text(labelexpand*y[, 1],labelexpand*y[, 2], labels = ylabs, cex = cex[2], col = ycol, ...) if (var.axes) arrows(0, 0, y[, 1], y[, 2], col = ycol, length = arrow.len) usrc <- par("usr") invisible(cbind(usrr=usrr,usrc=usrc)) } ## ================================================================ g.ellipse <- function(sigma=NULL, mu=NULL, data=NULL, b=NULL, prob=0.5, n=128, plot=TRUE, lty=2, col=1, ...) { ## Purpose: draw an ellipse characterizing a normal distribution ## ---------------------------------------------------------------------- ## Arguments: ## sigma covariance matrix ## b a "square root" of sigma (e.g. cholesky) ## either data or sigma or b must be provided ## mu expectation vector ## prob probability included in ellipse ## n number of points on the ellipse ## plot should ellipse be drawn? ## lty line type ## ---------------------------------------------------------------------- ## Author: Werner Stahel, Date: 3 May 2003, 12:06 on.exit(browser()) if (!is.null(data)) { b <- chol(var(data,na.rm=TRUE)) mu <- apply(data,2,mean,na.rm=TRUE) } if (is.null(b)) { if (is.null(sigma)) stop ("!g.ellipse!: either sigma or b must be provided") else b <- chol(sigma) } if (is.null(mu)) mu <- rep(0,nrow(b)) lty <- rep(lty,length=length(prob)) col <- rep(col,length=length(prob)) lphi <- seq(0,2*pi,length=n+1) for (li in 1:length(prob)) { lzz <- cbind(cos(lphi),sin(lphi)) * sqrt(qchisq(prob[li],nrow(b))) lz <- matrix(mu,n+1,2,byrow=T)+lzz%*%b if (plot) lines(lz[,1],lz[,2], lty=lty[li], col=col[li]) } if (!Browse) on.exit() invisible(lz) }