pca_biplot.r

#-----------------------------------------------------------------------------------------------------------------------------------------------
# Beroenden
require(ggplot2)
require(grid) # For arrow

# PCA Biplot
# for models x generated by princomp


pca_biplot <- function(x, choices = 1L:2L, main = 'PCA', plot.scores = TRUE, plot.loadings = TRUE,...) { 
	argList<-list(...)

	# Handle input
	if (length(choices) != 2L) 
	stop("length of choices must be 2")
	if (!length(scores <- x$scores)) 
	stop(gettextf("object '%s' has no scores", deparse(substitute(x))), 
	domain = NA)

	# Set up scores and loadings
	scores = x$scores[,choices]
	scores = as.data.frame(scores,)
	loads = loadings(x)[,choices]

	# scaling part
	lam <- x$sdev[choices]
	if (is.null(n <- x$n.obs)) 
	n <- 1
	lam <- lam * sqrt(n)
	scores = as.data.frame(t(t(scores)/lam))


	# as.data.frame(FAP_Functional.pc$loadings[])

	loads = as.data.frame(loads[]*lam*0.01)

	# Use load.scale to scale loads
	if( !is.null(argList$load.scale)) {
		load.scale = argList$load.scale
		# loads = as.data.frame(t(t(loads)*lam*load.scale))
		loads = as.data.frame(t(t(loads)*load.scale))
	}

	# Trick to make plot work
	names(loads) <- c('Comp.1','Comp.2')
	names(scores) <- c('Comp.1','Comp.2')

	# Add Rownames as a column
	loads$labels = rownames(loads)

	# Get the proportions of variace for labels and title
	prop_var <- x$sdev^2/sum(x$sdev^2)
	cum_var <- sum(prop_var[[choices[1]]],prop_var[[choices[2]]]) # The sum of variances from both components
	prop1 <- paste(',',format(prop_var[[choices[1]]], digits = 4))
	prop2 <- paste(',',format(prop_var[[choices[2]]], digits = 4))

	# paste('Canonical scores, FULL, cor = ',format(cc_full$cor[1],digits = 4))


	# The plotting part
	require(ggplot2); 
	# quartz()
	p <- ggplot(data = scores,aes(0,0))

	# # Add Loadings
	if(plot.loadings){
		p <- p + geom_segment(data = loads, aes(xend = Comp.1,yend = Comp.2),arrow=arrow(length=unit(0.2,"cm")),color='red',alpha = I(1/3),size = 1) +	# Plot loadings Arrows
		geom_text(data = loads, aes(x = Comp.1*1.1,y = Comp.2*1.1, label = labels ), angle=30)	# Text till loadings
	}

	# Add scores , pch and colouring of scores
	if(plot.scores){
		if( is.null(argList$colour)){
			p <- p + geom_point(aes(x = Comp.1,y = Comp.2),data = scores,alpha = I(2/3))
		}
		else{
			scores$color = argList$colour  
			p <- p + geom_point(aes(x = Comp.1,y = Comp.2, colour = color),data = scores,alpha = I(2/3))
		}
	}

	# shape
	if( !is.null(argList$shape)){
		p <- p + geom_shape

	}

	# textLabels and colouring
	if( !is.null(argList$text.labels)){
		if( !is.null(argList$colour)){
			scores$color = argList$colour
		} else {
			scores$color = 'grey80'
		}
		scores$text.labels = argList$text.labels
		p <- p + geom_text(aes(x = Comp.1,y = Comp.2, label = text.labels, colour = color),data = scores)
	}
	
	# Rug
	if( !is.null(argList$colour)){
		p <- p + geom_rug(data= scores, aes(x = Comp.1, y = Comp.2, colour = color))
	}

	# Ellipse
	if( !is.null(argList$ellipse)){
		require(geozoo)
		ellip <- ellipsoid(.1,.2,100)
		print(ellip)
		p <- p + geom_point(data <- ellip, mapping <- aes(x = x, y = y))

	}

	p <- p + 
	geom_hline(yintercept = 0,colour = 'grey50') +					# Axlar genom origo
	geom_vline(xintercept = 0,colour = 'grey50') +
	scale_x_continuous(paste("Comp.",choices[1],prop1, sep = '')) + # Set up axis labels
	scale_y_continuous(paste("Comp.",choices[2],prop2, sep = '')) + 	
	coord_equal( ratio = 1) +											# locks the aspect ratio so that you can't skew your data when adjusting the window size
	theme_bw() + 													# BW theme for 
	opts(title = paste(main, format(cum_var, digits = 4), sep = ', '))
	return(p)
	p
}


mybiplot <- function (x, choices = 1L:2L, scale = 1, pc.biplot = FALSE, ellipse = NULL, ellipse.col = 'black', ...) 
{
  if (length(choices) != 2L) 
    stop("length of choices must be 2")
  if (!length(scores <- x$x)) 
    stop(gettextf("object '%s' has no scores", deparse(substitute(x))), 
         domain = NA)
  if (is.complex(scores)) 
    stop("biplots are not defined for complex PCA")
  lam <- x$sdev[choices]
  n <- NROW(scores)
  lam <- lam * sqrt(n)
  if (scale < 0 || scale > 1) 
    warning("'scale' is outside [0, 1]")
  if (scale != 0) 
    lam <- lam^scale
  else lam <- 1
  if (pc.biplot) 
    lam <- lam/sqrt(n)
  
  if (missing(ellipse))
    ellipse <- NULL
  else ellipse_lam <- t(t(ellipse)/lam)
  
  mybiplot2(t(t(scores[, choices])/lam), t(t(x$rotation[,choices]) * lam), ellipse = ellipse_lam, ellipse.col = ellipse.col, ...)
  invisible()
}

mybiplot2 <- function (x, y, var.axes = TRUE, col, cex = rep(par("cex"), 2), 
                       xlabs = NULL, ylabs = NULL, expand = 1, xlim = NULL, ylim = NULL, 
                       arrow.len = 0.1, main = NULL, sub = NULL, xlab = NULL, ylab = NULL,
                       scorecol = NULL, ellipse = NULL, ellipse.col = 'black',
                       ...) 
{
  n <- nrow(x)
  p <- nrow(y)
  
  # labels
  if (missing(xlabs)) {
    xlabs <- dimnames(x)[[1L]]
    if (is.null(xlabs)) 
      xlabs <- 1L:n
  }
  xlabs <- as.character(xlabs)
  dimnames(x) <- list(xlabs, dimnames(x)[[2L]])
  if (missing(ylabs)) {
    ylabs <- dimnames(y)[[1L]]
    if (is.null(ylabs)) 
      ylabs <- paste("Var", 1L:p)
  }
  ylabs <- as.character(ylabs)
  dimnames(y) <- list(ylabs, dimnames(y)[[2L]])
  if (length(cex) == 1L) 
    cex <- c(cex, cex)
  if (missing(col)) {
    col <- par("col")
    if (!is.numeric(col)) 
      col <- match(col, palette(), nomatch = 1L)
    col <- c(col, col + 1L)
  }
  else if (length(col) == 1L) 
    col <- c(col, col)
  unsigned.range <- function(x) c(-abs(min(x, na.rm = TRUE)), 
                                  abs(max(x, na.rm = TRUE)))
  rangx1 <- unsigned.range(x[, 1L])
  rangx2 <- unsigned.range(x[, 2L])
  rangy1 <- unsigned.range(y[, 1L])
  rangy2 <- unsigned.range(y[, 2L])
  if (missing(xlim) && missing(ylim)) 
    xlim <- ylim <- rangx1 <- rangx2 <- range(rangx1, rangx2)
  else if (missing(xlim)) 
    xlim <- rangx1
  else if (missing(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)))
  
  
  if (missing(scorecol)) {
    scorecol <- 'black'
  }
  
  plot(x, type = "p", xlim = xlim, ylim = ylim, col = scorecol, 
       xlab = xlab, ylab = ylab, sub = sub, main = main, ...)
  
  #   Lägg till ellips
  if (!missing(ellipse)) {
    lines(ellipse[,1], ellipse[,2], lty = 2, col = ellipse.col)
  }
  
  #   text(x, xlabs, cex = cex[1L], col = col[1L], ...)
  par(new = TRUE)
  dev.hold()
  on.exit(dev.flush(), add = TRUE)
  plot(y, axes = FALSE, type = "n", xlim = xlim * ratio, ylim = ylim * 
         ratio, xlab = "", ylab = "", col = col[1L], ...) # tomplotta
  axis(3, col = col[2L], ...)
  axis(4, col = col[2L], ...)
  box(col = col[1L])
  text(y, labels = ylabs, cex = cex[2L], col = col[2L], ...)
  if (var.axes) 
    arrows(0, 0, y[, 1L] * 0.8, y[, 2L] * 0.8, col = col[2L], 
           length = arrow.len)
  invisible()
}