x, # output object from candisc

which=1:2, # canonical dimensions to plot

conf=0.95, # confidence coverage of circles for class means

col, # vector of unique colors used for plotting the canonical scores

pch, # vector of unique point symbols

scale, # scale factor for variable vectors in can space

asp=1, # aspect ratio, to ensure equal units

var.col="blue",

var.lwd=par("lwd"),

var.labels,

var.cex = 1,

var.pos,

rev.axes=c(FALSE, FALSE),

ellipse=FALSE, # draw data ellipses for canonical scores?

ellipse.prob = 0.68,

fill.alpha=0.1,

prefix = "Can", # prefix for labels of canonical dimensions

suffix = TRUE, # add label suffix with can % ?

titles.1d = c("Canonical scores", "Structure"),

points.1d = FALSE,

... # extra args passed to plot, e.g., notch = TRUE for boxplot

) {

term <- x$term

factors <- x$factors

rev.axes <- rep(rev.axes, length.out=2)

nlev <- if(length(dim(x$means))>1) nrow(x$means) else length(x$means) # number of groups

if (missing(col)) col <- rep(palette()[-1], length.out=nlev)

fill.col <- heplots::trans.colors(col, fill.alpha)

if (missing(pch)) pch <- rep(1:18, length.out=nlev)

# ----------- 1D plot -------------

if (x$ndim < 2 || length(which)==1) {

which <- which[1]

### is there a better way to show the 1D distributions of canonical scores?

### Show the canonical structure coeffs -- vectors from 0

op <- par(no.readonly = TRUE) # save default, for resetting...

ng <- length(x$means)

structure <- as.vector(x$structure[,which])

canvar <- paste('Can', which, sep="") # names of canonical variable to plot

scores <- x$scores[, canvar]

if(isTRUE(rev.axes[1])) {

x$scores[, canvar] <- -x$scores[, canvar]

structure <- -structure

}

ns <- length(structure)

wid <- if (ns < 2*ng) c(2,1) else c(1.2,1)

#cat("ng:", ng, "\tns:", ns, "\twid:", wid, "\n")

layout(matrix(c(1,2),1,2), widths=wid)

par(mar=c(5,4,4,0)+.1)

if (is.logical(suffix) & suffix)

suffix <- paste( " (", round(x$pct[which],1), "%)", sep="" ) else suffix <- NULL

canlab <- paste(prefix, which, suffix, sep="")

formule <- formula( paste0(canvar, " ~", term) )

boxplot(formule, data=x$scores, ylab=canlab, xlab=term, col=fill.col, main=titles.1d[1], ...)

# add points?

if (points.1d) points(jitter(as.numeric(x$scores[, term]), factor=0.5), x$scores[,canvar],

col=col[x$scores[, term]],

pch=pch[x$scores[, term]] )

xx <- 1:ns

par(mar=c(5,0,4,1)+.1)

# TODO: should set ylim = c(0,1) or c(-1,1) and maybe draw tick marks

ylim <- range( c(structure,

if (all(structure > 0)) +1 else 0,

if (all(structure < 0)) -1 else 0)

)

xlim <- c(0.5, ns+0.5)

plot(xx, structure, type="n", ylab="", xlab="",

xaxt="n", yaxt="n", xlim=xlim, ylim=ylim, main=titles.1d[2])

arrows(xx, 0, xx, structure, length=.1, angle=15,

col=var.col, lwd=var.lwd )

abline(h=0, lty=2, col="gray")

vars <- if(missing(var.labels)) rownames(x$structure) else var.labels

adj1 <- as.vector(ifelse (structure > 0, 1, 0))

adj2 <- rep(-0.3, ns)

adj2[1] <- 1.1

for (i in 1:ns) text(xx[i], structure[i], paste(" ", vars[i] ," "),

adj=c(adj1[i], adj2[i]), col=var.col, srt=90, cex=var.cex, xpd=TRUE)

par(op)

return(invisible())

}

# ----------- 2D plot -------------

canvar <- paste('Can', which, sep="") # names of canonical variables to plot

if (is.logical(suffix) & suffix)

suffix <- paste( " (", round(x$pct[which],1), "%)", sep="" ) else suffix <- NULL

canlab <- paste(prefix, which, suffix, sep="")

nlev <- nrow(x$means) # number of groups

# TODO: can we be more clever about assigning default col & pch by taking the

# structure of x$factors into account?

if (missing(col)) col <- rep(palette(), length.out=nlev)

scores <- x$scores[,canvar]

means <- x$means[,which]

labels <- rownames(x$means)

structure <- x$structure[,which]

if(isTRUE(rev.axes[1])) {

scores[, 1] <- -scores[, 1]

means[, 1] <- -means[, 1]

structure[, 1] <- -structure[, 1]

}

if(isTRUE(rev.axes[2])) {

scores[, 2] <- -scores[, 2]

means[, 2] <- -means[, 2]

structure[, 2] <- -structure[, 2]

}

# use asp=1 to make the plot equally scaled

Ind <- if(length(factors)==1)

match(x$scores[,term], levels(x$scores[,term]))

else dataIndex(x$scores, term)

plot(scores, asp=asp,

xlab=canlab[1], ylab=canlab[2],

col=col[Ind], pch=pch[Ind], ...)

abline(h=0, v=0,

lty=2, col="grey")

if(ellipse) {

fill.col <- heplots::trans.colors(col, fill.alpha)

radius <- sqrt(qchisq(ellipse.prob, df = 2))

angles <- (0:60)*2*pi/60

circle <- radius * cbind( cos(angles), sin(angles))

for(i in 1:nlev) {

sigma <- var(scores[Ind==i,])

mu <- as.numeric(means[i,])

Q <- chol(sigma, pivot = TRUE)

order <- order(attr(Q, "pivot"))

ell <- sweep(circle %*% Q[, order], 2,

mu, FUN = "+")

polygon(ell, col=fill.col[i], border=col[i], lty=1)

}

}

points(means[,1], means[,2], col=col, pch="+", cex=2)

pos <- ifelse(means[,2]>0, 3, 1)

text(means[,1], means[,2],

col = col,

labels=labels, pos=pos)

# plot variable vectors

# DONE: replaced previous scaling with vecscale()

if (missing(scale)) {

scale <- vecscale(structure)

message("Vector scale factor set to ", round(scale, 3))

}

# DONE: replaced with call to vectors()

cs <- scale * structure

if(!missing(var.labels)) rownames(cs) <- var.labels

vectors(cs, col=var.col, cex=var.cex, lwd=var.lwd, pos=var.pos, xpd=TRUE)

### why doesn't this work???

circle <- function( center, radius, segments=41, ...) {

angles <- (0:segments) * 2 * pi/segments

unit.circle <- cbind(cos(angles), sin(angles))

circle <- t(as.numeric(center) + radius*t(unit.circle))

lines(circle, col = col, ...)

}

# plot confidence circles for canonical means

if ( conf>0 ) {

n <- as.vector(table(factors))

radii <- sqrt(qchisq(conf, 2) / n)

symbols(means, circles=radii, inches=FALSE, add=TRUE, fg=col)

}