[med-svn] [r-cran-calibrate] 01/02: Imported Upstream version 1.7.2
Dylan Aïssi
bob.dybian-guest at moszumanska.debian.org
Sun Mar 26 21:33:55 UTC 2017
This is an automated email from the git hooks/post-receive script.
bob.dybian-guest pushed a commit to branch master
in repository r-cran-calibrate.
commit eabfeb206ba72e6dd9c249f515738403306f954b
Author: Dylan Aïssi <bob.dybian at gmail.com>
Date: Sun Mar 26 22:48:22 2017 +0200
Imported Upstream version 1.7.2
---
DESCRIPTION | 14 +
MD5 | 38 ++
NAMESPACE | 1 +
R/calibrate.R | 116 ++++++
R/canocor.R | 50 +++
R/circle.R | 13 +
R/dlines.R | 14 +
R/half.R | 7 +
R/mhalf.R | 7 +
R/ones.R | 5 +
R/origin.R | 10 +
R/rad2degree.R | 6 +
R/rda.R | 100 ++++++
R/shiftvector.R | 48 +++
R/textxy.R | 19 +
data/calves.rda | Bin 0 -> 183 bytes
data/goblets.rda | Bin 0 -> 364 bytes
data/heads.rda | Bin 0 -> 365 bytes
data/linnerud.rda | Bin 0 -> 425 bytes
data/storks.rda | Bin 0 -> 479 bytes
inst/doc/CalibrationGuide.R | 367 +++++++++++++++++++
inst/doc/CalibrationGuide.Rnw | 777 +++++++++++++++++++++++++++++++++++++++++
inst/doc/CalibrationGuide.pdf | Bin 0 -> 285812 bytes
man/calibrate.Rd | 96 +++++
man/calves.Rd | 16 +
man/canocor.Rd | 56 +++
man/circle.Rd | 24 ++
man/dlines.Rd | 31 ++
man/goblets.Rd | 16 +
man/heads.Rd | 23 ++
man/linnerud.Rd | 18 +
man/ones.Rd | 27 ++
man/origin.Rd | 19 +
man/rad2degree.Rd | 25 ++
man/rda.Rd | 54 +++
man/shiftvector.Rd | 60 ++++
man/storks.Rd | 17 +
man/textxy.Rd | 39 +++
vignettes/CalibrationGuide.Rnw | 777 +++++++++++++++++++++++++++++++++++++++++
39 files changed, 2890 insertions(+)
diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..0a526b7
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,14 @@
+Package: calibrate
+Version: 1.7.2
+Date: 2013-09-09
+Title: Calibration of Scatterplot and Biplot Axes
+Author: Jan Graffelman <jan.graffelman at upc.edu>
+Maintainer: Jan Graffelman <jan.graffelman at upc.edu>
+Depends: R (>= 1.8.0), MASS
+Description: Package for drawing calibrated scales with tick marks on (non-orthogonal)
+ variable vectors in scatterplots and biplots.
+License: GPL-2
+Packaged: 2013-09-09 15:33:11 UTC; Beagle
+NeedsCompilation: no
+Repository: CRAN
+Date/Publication: 2013-09-10 11:10:21
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..587aa81
--- /dev/null
+++ b/MD5
@@ -0,0 +1,38 @@
+44671a8fade2a53939ecfbbd23d0e80b *DESCRIPTION
+f4c3a89b0a423584934958144d60629d *NAMESPACE
+92c4b9e372e384795e55b44dde66b0df *R/calibrate.R
+d29089cb6db670a9da6c62aa10f9c9b9 *R/canocor.R
+f010fe0dbfc45554ccbdb349eb3000a0 *R/circle.R
+b64da62d2f549501a556e3c5eae15904 *R/dlines.R
+afa03a1e580af939f82db7e117f12e63 *R/half.R
+64dcadda820ecb14578fea99376eefa1 *R/mhalf.R
+a0f081cc71f9da04a8268e3f890bd131 *R/ones.R
+3b3fc980be0859a385adab93db5b4c0b *R/origin.R
+140e137806866cd825d7eba7caf92e02 *R/rad2degree.R
+317ddceb7fa4b58394a2a27b3fb2bd54 *R/rda.R
+5ea3adfb04c8deb61b33c08c4ccf2ad1 *R/shiftvector.R
+e9e222e963766157fa203bdb8893b44a *R/textxy.R
+f1bb16e6ec7c42f94217e03aa4cf8041 *data/calves.rda
+a7f950a2dc032bf42bfb71844c48b0f3 *data/goblets.rda
+42bbaecdc84a7cd478609dec035ce980 *data/heads.rda
+3bbb8680e32deb7e0343adefbd84e931 *data/linnerud.rda
+7853c4419d06e40aeb506d476636e48f *data/storks.rda
+85040338b77619c141db2636531dd1a3 *inst/doc/CalibrationGuide.R
+ddd8ece673041df6728243c93dae30f2 *inst/doc/CalibrationGuide.Rnw
+64a550cf97b0cda90ad5515efdd1b886 *inst/doc/CalibrationGuide.pdf
+24a4c36dbe44f734067ae38be21fb8e0 *man/calibrate.Rd
+0908691b309ec6e5de94c2303c078225 *man/calves.Rd
+1b64a664fff01decf17c7b3f1174609a *man/canocor.Rd
+dc0b8f20cb7811bfc3112a7a2c2ac8c1 *man/circle.Rd
+e44e58382308c0882e552759d4ea3b3e *man/dlines.Rd
+660307f4e93682b29001fd9fc0d79089 *man/goblets.Rd
+d5f91be7b56eec8e0ad526a36c871d7f *man/heads.Rd
+e34e8979ed1ce76405fe3156d6903b7e *man/linnerud.Rd
+f106332b43b4e9218330fa790e6ba67f *man/ones.Rd
+22c9b09911d48e09e429301481089ac0 *man/origin.Rd
+8df0821fb8a396ffee5a895dad3e3a4c *man/rad2degree.Rd
+df46ed53135cca9da7a0e24ace44ca1b *man/rda.Rd
+430fad46bf95b8884d885c21c7fa5999 *man/shiftvector.Rd
+70b0962a116fb08c3297e34a1564de0b *man/storks.Rd
+5a9e2c7e7590ebd5f43aa70cd386b174 *man/textxy.Rd
+ddd8ece673041df6728243c93dae30f2 *vignettes/CalibrationGuide.Rnw
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..f8bb463
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1 @@
+export (calibrate,canocor,circle,dlines,origin,ones,rad2degree,rda,shiftvector,textxy)
diff --git a/R/calibrate.R b/R/calibrate.R
new file mode 100644
index 0000000..6c3111b
--- /dev/null
+++ b/R/calibrate.R
@@ -0,0 +1,116 @@
+calibrate <-
+function (g, y, tm, Fr, tmlab = tm, tl = 0.05, dt = TRUE, dp = FALSE,
+ lm = TRUE, verb = TRUE, axislab = "", reverse = FALSE,
+ alpha = NULL, labpos = 1, weights = diag(rep(1, length(y))),
+ axiscol = "blue", cex.axislab = 0.75, graphics = TRUE, where = 3,
+ laboffset = c(0, 0), m = matrix(c(0, 0), nrow = 1), markerpos = 3,
+ showlabel = TRUE, lwd = 1, shiftvec = c(0,0), shiftdir = "none", shiftfactor = 1.05 )
+{
+ g <- as.vector(g)
+ if (is.matrix(weights))
+ Dw <- weights
+ else if (is.vector(weights))
+ Dw <- diag(weights)
+ else stop("calibrate: weights is not a vector or matrix")
+ if (verb > 1)
+ print(Dw)
+ if(all(!as.logical(shiftvec))) { # c(0,0) vector
+ shiftvec <- switch(shiftdir,
+ right = shiftvector(g,Fr)$dr,
+ left = shiftvector(g,Fr)$dl,
+ none = c(0,0),
+ stop("invalid value for parameter shiftdir"))
+ }
+ shiftvec <- shiftfactor*shiftvec
+ optalpha <- t(g) %*% t(Fr) %*% Dw %*% Fr %*% g/((t(y) %*%
+ Dw %*% Fr %*% g) * (t(g) %*% g))
+ optalpha <- optalpha[1, 1]
+ useralpha <- NULL
+ if (is.null(alpha)) {
+ alpha <- optalpha
+ }
+ else {
+ useralpha <- alpha
+ }
+ Mmean <- matrix(rep(1, length(tm)), ncol = 1) %*% m
+ M <- alpha * (tm) %*% t(g) + Mmean
+ nrM <- nrow(M)
+ M <- M + matrix(rep(1, nrM), ncol = 1) %*% shiftvec
+ di <- 1/(alpha * t(g) %*% g)
+ yt <- di[1, 1] * Fr %*% g
+ e <- y - yt
+ Q <- t(e) %*% Dw %*% e
+ gos <- 1 - Q/(t(y) %*% Dw %*% y)
+ odi <- 1/(optalpha * t(g) %*% g)
+ oyt <- odi[1, 1] * Fr %*% g
+ oe <- y - oyt
+ oQ <- t(oe) %*% Dw %*% oe
+ gof <- 1 - oQ/(t(y) %*% Dw %*% y)
+ ang <- atan(g[2]/g[1]) * 180/pi
+ lengthoneunit <- alpha * sqrt(t(g) %*% g)
+ if (verb) {
+ cat("---------- Calibration Results for ", axislab, " ")
+ for (i in 1:(60 - (38 + nchar(axislab)))) cat("-")
+ cat("\n")
+ cat("Length of 1 unit of the original variable = ", round(lengthoneunit,
+ digits = 4), " \n")
+ cat("Angle = ", round(ang,
+ digits = 2), "degrees\n")
+ cat("Optimal calibration factor = ", round(optalpha,
+ digits = 4), " \n")
+ cat("Used calibration factor = ", round(alpha,
+ digits = 4), " \n")
+ cat("Goodness-of-fit = ", round(gof,
+ digits = 4), " \n")
+ cat("Goodness-of-scale = ", round(gos,
+ digits = 4), " \n")
+ cat("------------------------------------------------------------\n")
+ }
+ Fr2 <- Fr[, 1:2]
+ nn <- t(g) %*% g
+ scal <- (Fr2 %*% g)/nn[1, 1]
+ Dscal <- diag(as.vector(scal))
+ Fpr <- Dscal %*% matrix(rep(1, nrow(Fr)), ncol = 1) %*% t(g)
+ deltax <- tl * sin(ang * pi/180)
+ deltay <- tl * cos(ang * pi/180)
+ if (reverse == TRUE)
+ Mn <- cbind(M[, 1] - deltax, M[, 2] + deltay)
+ else Mn <- cbind(M[, 1] + deltax, M[, 2] - deltay)
+ if (graphics) {
+ lines(rbind(M[1, ], M[nrM, ]), col = axiscol, lwd = lwd)
+ if (lm) {
+ if (reverse == TRUE)
+ text(Mn[, 1], Mn[, 2], tmlab, pos = markerpos,
+ offset = 0.2, cex = cex.axislab, srt = ang)
+ else if (markerpos > 2)
+ text(Mn[, 1], Mn[, 2], tmlab, pos = markerpos -
+ 2, offset = 0.2, cex = cex.axislab, srt = ang)
+ else text(Mn[, 1], Mn[, 2], tmlab, pos = markerpos +
+ 2, offset = 0.2, cex = cex.axislab, srt = ang)
+ }
+ nm <- nrow(M)
+ if (dt == TRUE) {
+ for (i in 1:nm) lines(rbind(M[i, 1:2], Mn[i, 1:2]),
+ col = axiscol, lwd = lwd)
+ }
+ if (dp) {
+ nrFpr <- nrow(Fpr)
+ dlines(Fr2 + matrix(rep(1, nrFpr), ncol = 1) %*%
+ m, Fpr + matrix(rep(1, nrFpr), ncol = 1) %*%
+ m + matrix(rep(1, nrFpr), ncol = 1) %*% shiftvec)
+ }
+ if (showlabel) {
+ switch(where, text(M[1, 1] + laboffset[1], M[1, 2] +
+ laboffset[2], axislab, cex = cex.axislab, srt = ang, pos = labpos, adj = c(0.5,0.5)), text(M[round(nrow(M)/2),
+ 1] + laboffset[1], M[round(nrow(M)/2), 2] + laboffset[2],
+ axislab, cex = cex.axislab, srt = ang, pos = labpos,
+ adj = c(0.5,0.5) ), text(M[nrow(M), 1] + laboffset[1],
+ M[nrow(M), 2] + laboffset[2], axislab, cex = cex.axislab,
+ srt = ang, pos = labpos, adj = c(0.5,0.5) ))
+ }
+ }
+ return(list(useralpha = useralpha, gos = gos, optalpha = optalpha,
+ gof = gof, lengthoneunit = lengthoneunit, M = M, Q = Q,
+ ang = ang, shiftvec = shiftvec, yt = yt, e = e, Fpr = Fpr, Mn = Mn))
+}
+
diff --git a/R/canocor.R b/R/canocor.R
new file mode 100644
index 0000000..6ec160d
--- /dev/null
+++ b/R/canocor.R
@@ -0,0 +1,50 @@
+"canocor" <-
+function (X, Y)
+{
+ Xs <- scale(X)
+ Ys <- scale(Y)
+ Rxx <- cor(X)
+ Ryy <- cor(Y)
+ Rxy <- cor(X, Y)
+ d <- diag(eigen(Rxx)$values)
+ v <- eigen(Rxx)$vectors
+ Rxxmh <- v %*% sqrt(solve(d)) %*% t(v)
+ d <- diag(eigen(Ryy)$values)
+ v <- eigen(Ryy)$vectors
+ Ryymh <- v %*% sqrt(solve(d)) %*% t(v)
+ K <- Rxxmh %*% Rxy %*% Ryymh
+ D <- diag(svd(K)$d)
+ Ah <- svd(K)$u
+ Bh <- svd(K)$v
+ A <- Rxxmh %*% Ah
+ B <- Ryymh %*% Bh
+ U <- Xs %*% A
+ V <- Ys %*% B
+ Fs <- Rxx %*% A
+ Fp <- Rxx %*% A %*% D
+ Gs <- Ryy %*% B
+ Gp <- Ryy %*% B %*% D
+ Rxu <- Fs
+ Rxv <- Fp
+ Ryv <- Gs
+ Ryu <- Gp
+ lamb <- diag(D^2)
+ frac <- lamb/sum(lamb)
+ cumu <- cumsum(frac)
+ fitRxy <- rbind(lamb, frac, cumu)
+ AdeX <- apply(Rxu * Rxu, 2, mean)
+ AdeY <- apply(Ryv * Ryv, 2, mean)
+ RedX <- apply(Rxv * Rxv, 2, mean)
+ RedY <- apply(Ryu * Ryu, 2, mean)
+ cAdeX <- cumsum(AdeX)
+ cAdeY <- cumsum(AdeY)
+ cRedX <- cumsum(RedX)
+ cRedY <- cumsum(RedY)
+ fitXs <- rbind(AdeX, cAdeX)
+ fitXp <- rbind(RedX, cRedX)
+ fitYs <- rbind(AdeY, cAdeY)
+ fitYp <- rbind(RedY, cRedY)
+ return(list(ccor = D, A = A, B = B, U = U, V = V, Fs = Fs, Gs = Gs,
+ Fp = Fp, Gp = Gp, fitRxy = fitRxy, fitXs = fitXs, fitXp = fitXp,
+ fitYs = fitYs, fitYp = fitYp))
+}
diff --git a/R/circle.R b/R/circle.R
new file mode 100644
index 0000000..7ff3db8
--- /dev/null
+++ b/R/circle.R
@@ -0,0 +1,13 @@
+"circle" <- function(radius = 1, origin = c(0,0))
+{
+ t <- seq(-pi,pi,by=0.01)
+ a <- origin[1]
+ b <- origin[2]
+ r <- radius
+ x <- a + r*cos(t)
+ y <- b + r*sin(t)
+ points(x,y,type="l")
+ return(NULL)
+}
+
+
diff --git a/R/dlines.R b/R/dlines.R
new file mode 100644
index 0000000..2f60def
--- /dev/null
+++ b/R/dlines.R
@@ -0,0 +1,14 @@
+"dlines" <-
+function(SetA,SetB,lin="dotted") {
+#
+# Function DLINES connects the rows of SetA to the rows of SetB with lines
+#
+# Jan Graffelman
+# Universitat Politecnica de Catalunya
+# January 2004
+#
+np<-nrow(SetA)
+for(i in 1:np) lines(rbind(SetA[i,1:2],SetB[i,1:2]),lty=lin)
+return(NULL)
+}
+
diff --git a/R/half.R b/R/half.R
new file mode 100644
index 0000000..8e9f258
--- /dev/null
+++ b/R/half.R
@@ -0,0 +1,7 @@
+"half" <-
+function(A) {
+ d <- diag(eigen(A)$values)
+ v <- eigen(A)$vectors
+ Ah <- v%*%sqrt(d)%*%t(v)
+}
+
diff --git a/R/mhalf.R b/R/mhalf.R
new file mode 100644
index 0000000..03ad2b4
--- /dev/null
+++ b/R/mhalf.R
@@ -0,0 +1,7 @@
+"mhalf" <-
+function(A) {
+ d <- diag(eigen(A)$values)
+ v <- eigen(A)$vectors
+ Ah <- v%*%solve(sqrt(d))%*%t(v)
+}
+
diff --git a/R/ones.R b/R/ones.R
new file mode 100644
index 0000000..696db4a
--- /dev/null
+++ b/R/ones.R
@@ -0,0 +1,5 @@
+"ones" <-
+function(n,p=n) {
+ matrix(rep(1,n*p),nrow=n,ncol=p)
+}
+
diff --git a/R/origin.R b/R/origin.R
new file mode 100644
index 0000000..092246d
--- /dev/null
+++ b/R/origin.R
@@ -0,0 +1,10 @@
+"origin" <-
+function (m = c(0, 0), ...)
+{
+ vx <- c(par("usr")[1],par("usr")[2])
+ vy <- c(m[2],m[2])
+ lines(vx,vy,...)
+ vy <- c(par("usr")[3],par("usr")[4])
+ vx <- c(m[1],m[1])
+ lines(vx,vy,...)
+}
diff --git a/R/rad2degree.R b/R/rad2degree.R
new file mode 100644
index 0000000..1aa8493
--- /dev/null
+++ b/R/rad2degree.R
@@ -0,0 +1,6 @@
+rad2degree <-
+function(x) {
+ ang <- x*180/pi
+ return(ang)
+}
+
diff --git a/R/rda.R b/R/rda.R
new file mode 100644
index 0000000..cec8d75
--- /dev/null
+++ b/R/rda.R
@@ -0,0 +1,100 @@
+"rda" <-
+function(X,Y,scaling=1) {
+#
+# Function RDA performs a redundancy analysis of the data in X
+# and Y.
+#
+# scaling = 0 : use centred variables (X and Y)
+# scaling = 1 : use centred and standardized variables (X and Y)
+#
+# Jan Graffelman
+# Universitat Politecnica de Catalunya
+# January 2004
+#
+
+n<-nrow(X) # determine # of cases
+p<-ncol(X) # determine # of variables
+
+nY<-nrow(Y)
+q<-ncol(Y)
+
+if (scaling==0) {
+ Xa <- scale(X,scale=FALSE)
+ Ya <- scale(Y,scale=FALSE)
+}
+else {
+ if (scaling==1) {
+ Xa<-scale(X)
+ Ya<-scale(Y)
+ }
+ else
+ stop("rda: improper scaling parameter")
+}
+
+Rxx <- cor(X)
+Ryy <- cor(Y)
+
+B<-solve(t(Xa)%*%Xa)%*%t(Xa)%*%Ya
+
+Yh<-Xa%*%B
+
+pca.results <- princomp(Yh,cor=FALSE)
+Fp <- pca.results$scores
+Ds <- diag(sqrt(diag(var(Fp))))
+Fs <- Fp%*%solve(Ds)
+Gs <- pca.results$loadings
+Gp <- Gs%*%Ds
+
+la <- diag(var(Fp))
+laf <- la/sum(la)
+lac <- cumsum(laf)
+
+decom <- rbind(la,laf,lac)
+
+# It is important not to standardize Yh, and not to use cor=T. This will inflate the
+# variance of Yh, and give different eigenvalues.
+
+# Fs, Gp' biplot of fitted values
+
+Fs <- Fs[,1:min(p,q)]
+Fp <- Fp[,1:min(p,q)]
+Gs <- Gs[,1:min(p,q)]
+Gp <- Gp[,1:min(p,q)]
+
+Gxs <- solve(t(Xa)%*%Xa)%*%t(Xa)%*%Fs
+Gxp <- solve(t(Xa)%*%Xa)%*%t(Xa)%*%Fp
+
+Gyp <- Gp
+Gys <- Gs
+
+# Gxs Gyp', Gxp Gys' biplots of B (regression coefficients)
+
+# goodness of fit of regression coefficients
+
+decom <- decom[,1:min(p,q)]
+
+# alternative computations doing SVD of B.
+
+#W<-t(Xa)%*%Xa
+
+#result <- svd(1/(sqrt(n-1))*half(W)%*%B)
+#Gxxs <- sqrt(n-1)*mhalf(W)%*%result$u # = Gxs
+#Gyyp <- result$v%*%diag(result$d) # = Gyp
+#Gxxp <- sqrt(n-1)*mhalf(W)%*%result$u%*%diag(result$d) # = Gxp
+#Gyys <- result$v # = Gys
+
+#dd <- result$d*result$d
+#dds <- cumsum(dd)
+#ddf <- dd/sum(dd)
+#ddc <- cumsum(ddf)
+
+#decB <- rbind(dd,dds,ddf,ddc)
+
+#res<-t(Gyyp)%*%Gyyp # = D^2
+#res<-t(Gyys)%*%Gyys # = I
+#res<-t(Gxxp)%*%Rxx%*%Gxxp # = D^2
+#res<-t(Gxxs)%*%Rxx%*%Gxxs # = I
+
+return(list(Yh=Yh,B=B,decom=decom,Fs=Fs,Gyp=Gyp,Fp=Fp,Gys=Gys,Gxs=Gxs,Gxp=Gxp))
+}
+
diff --git a/R/shiftvector.R b/R/shiftvector.R
new file mode 100644
index 0000000..52d2baa
--- /dev/null
+++ b/R/shiftvector.R
@@ -0,0 +1,48 @@
+shiftvector <-
+function(g,X,x=c(1,0),verbose=FALSE) {
+ #
+ # compute the optimal shift vector for a calibrated axis.
+ #
+ costheta <- ((t(g)%*%x))/(sqrt(t(g)%*%g))
+ theta <- acos(costheta)
+ sintheta <- sin(theta)
+ ang <- rad2degree(theta)
+ if(verbose) cat("angle vector g -- x-axis:",theta," rad",ang," degrees\n")
+ if(g[2]>=0)
+ H <- matrix(c(costheta, sintheta, -sintheta, costheta), ncol = 2) # clockwise
+ else
+ H <- matrix(c(costheta, -sintheta, sintheta, costheta), ncol = 2) # counterclockwise
+ Xr <- X%*%H
+ above <- Xr[,2] > 0
+ below <- Xr[,2] < 0
+
+ mg <- matrix(rep(1, nrow(X)), ncol = 1)%*%g
+ Dp <- diag(as.vector((X%*%g)/(sum(g*g))))
+ dmat <- X-Dp%*%mg
+ vecdn <- diag(dmat%*%t(dmat) )
+
+ if(sum(above)>0) {
+ vecdna <- vecdn[above]
+ dmata <- dmat[above,]
+ if(is.vector(dmata)) dl <- dmata
+ else {
+ i <- which.max(vecdna)
+ dl <- dmata[i,]
+ }
+ } else dl <- c(0,0)
+ if(sum(below)>0) {
+ vecdnb <- vecdn[below]
+ dmatb <- dmat[below,]
+ if(is.vector(dmatb)) dr <- dmatb
+ else {
+ i <- which.max(vecdnb)
+ dr <- dmatb[i,]
+ }
+ } else dr <- c(0,0)
+ if(verbose) {
+ cat("Left shift vector: ",dl,"\n")
+ cat("Right shift vector: ",dr,"\n")
+ }
+ return(list(dr=dr,dl=dl))
+}
+
diff --git a/R/textxy.R b/R/textxy.R
new file mode 100644
index 0000000..f90cecd
--- /dev/null
+++ b/R/textxy.R
@@ -0,0 +1,19 @@
+"textxy" <- function (X, Y, labs, m = c(0, 0), cex = 0.5, offset = 0.8, ...)
+{
+ posXposY <- ((X >= m[1]) & ((Y >= m[2])))
+ posXnegY <- ((X >= m[1]) & ((Y < m[2])))
+ negXposY <- ((X < m[1]) & ((Y >= m[2])))
+ negXnegY <- ((X < m[1]) & ((Y < m[2])))
+ if (sum(posXposY) > 0)
+ text(X[posXposY], Y[posXposY], labs[posXposY], adj = c(0.5-offset,
+ 0.5-offset), cex = cex, ...)
+ if (sum(posXnegY) > 0)
+ text(X[posXnegY], Y[posXnegY], labs[posXnegY], adj = c(0.5-offset,
+ 0.5+offset), cex = cex, ...)
+ if (sum(negXposY) > 0)
+ text(X[negXposY], Y[negXposY], labs[negXposY], adj = c(0.5+offset,
+ 0.5-offset), cex = cex, ...)
+ if (sum(negXnegY) > 0)
+ text(X[negXnegY], Y[negXnegY], labs[negXnegY], adj = c(0.5+offset,
+ 0.5+offset), cex = cex, ...)
+}
diff --git a/data/calves.rda b/data/calves.rda
new file mode 100644
index 0000000..64e05e3
Binary files /dev/null and b/data/calves.rda differ
diff --git a/data/goblets.rda b/data/goblets.rda
new file mode 100644
index 0000000..501136e
Binary files /dev/null and b/data/goblets.rda differ
diff --git a/data/heads.rda b/data/heads.rda
new file mode 100644
index 0000000..9671876
Binary files /dev/null and b/data/heads.rda differ
diff --git a/data/linnerud.rda b/data/linnerud.rda
new file mode 100644
index 0000000..4dd04aa
Binary files /dev/null and b/data/linnerud.rda differ
diff --git a/data/storks.rda b/data/storks.rda
new file mode 100644
index 0000000..79cfe4e
Binary files /dev/null and b/data/storks.rda differ
diff --git a/inst/doc/CalibrationGuide.R b/inst/doc/CalibrationGuide.R
new file mode 100644
index 0000000..f851622
--- /dev/null
+++ b/inst/doc/CalibrationGuide.R
@@ -0,0 +1,367 @@
+### R code from vignette source 'CalibrationGuide.Rnw'
+### Encoding: ISO8859-1
+
+###################################################
+### code chunk number 1: CalibrationGuide.Rnw:79-92
+###################################################
+prefig <- function(){
+data(goblets)
+X <- goblets
+plot(X[,1],X[,2],pch=19,cex=0.5,xlab=expression(X[1]),
+ylab=expression(X[2]),xlim=c(5,25),ylim=c(5,25),asp=1)
+m <- apply(cbind(X[,1],X[,2]),2,mean)
+textxy(X[,1],X[,2],1:25,m=m,cex=0.75)
+points(m[1],m[2],col="red",pch=19,cex=0.5)
+abline(h = m[2])
+abline(v = m[1])
+}
+options("SweaveHooks"=list(aap=prefig))
+options("width"=60)
+
+
+###################################################
+### code chunk number 2: noot
+###################################################
+library(calibrate)
+
+
+###################################################
+### code chunk number 3: CalibrationGuide.Rnw:113-115
+###################################################
+data(goblets)
+X <- goblets
+
+
+###################################################
+### code chunk number 4: CalibrationGuide.Rnw:125-130
+###################################################
+plot(X[,1],X[,2],pch=19,cex=0.5,xlab=expression(X[1]),ylab=expression(X[2]),
+ xlim=c(5,25),ylim=c(5,25),asp=1)
+m <- apply(X[,1:2],2,mean)
+textxy(X[,1],X[,2],1:25,m=m,cex=0.75)
+origin(m)
+
+
+###################################################
+### code chunk number 5: CalibrationGuide.Rnw:137-143
+###################################################
+getOption("SweaveHooks")[["aap"]]()
+Xc <- scale(X,center=TRUE,scale=FALSE)
+b <- solve(t(Xc[,1:2])%*%Xc[,1:2])%*%t(Xc[,1:2])%*%Xc[,5]
+print(b)
+bscaled <- 20*b
+arrows(m[1],m[2],m[1]+bscaled[1],m[2]+bscaled[2],col="blue",length=0.1)
+arrows(m[1],m[2],m[1]-bscaled[1],m[2]-bscaled[2],length=0,lty="dashed",col="blue")
+
+
+###################################################
+### code chunk number 6: CalibrationGuide.Rnw:155-161
+###################################################
+getOption("SweaveHooks")[["aap"]]()
+print(range(X[,5]))
+yc <- scale(X[,5],scale=FALSE)
+tm <- seq(2,10,by=1)
+tmc <- tm - mean(X[,5])
+Calibrate.X5<-calibrate(b,yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.3,axislab="X_5",
+ labpos=4,cex.axislab=1)
+
+
+###################################################
+### code chunk number 7: CalibrationGuide.Rnw:198-203
+###################################################
+getOption("SweaveHooks")[["aap"]]()
+yc <- scale(X[,5],scale=FALSE)
+tm <- seq(2,10,by=1)
+tmc <- tm - mean(X[,5])
+Calibrate.X5<-calibrate(b,yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.3,axislab="X_5",labpos=4,
+ cex.axislab=1,shiftdir="right")
+
+
+###################################################
+### code chunk number 8: CalibrationGuide.Rnw:220-229
+###################################################
+getOption("SweaveHooks")[["aap"]]()
+opar <- par('xpd'=TRUE)
+tm <- seq(3,8,by=1)
+tmc <- (tm - mean(X[,5]))
+Calibrate.rightmargin.X5 <- calibrate(c(0,1),yc,tmc,Xc[,1:2],tmlab=tm,m=m,
+ axislab="X_5",tl=0.5,
+ shiftvec=c(par('usr')[2]-mean(X[,1]),0),
+ shiftfactor=1,where=2,
+ laboffset=c(1.5,1.5),cex.axislab=1)
+par(opar)
+
+
+###################################################
+### code chunk number 9: CalibrationGuide.Rnw:252-260
+###################################################
+getOption("SweaveHooks")[["aap"]]()
+tm <- seq(2,10,by=1)
+tmc <- (tm - mean(X[,5]))
+Calibrate.X5 <- calibrate(b,yc,tmc,Xc[,1:2],tmlab=tm,m=m,axislab="X_5",tl=0.5,
+ dp=TRUE,labpos=4)
+tm <- seq(2,10,by=0.1)
+tmc <- (tm - mean(X[,5]))
+Calibrate.X5 <- calibrate(b,yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.25,verb=FALSE,
+ lm=FALSE)
+
+
+###################################################
+### code chunk number 10: CalibrationGuide.Rnw:276-298
+###################################################
+getOption("SweaveHooks")[["aap"]]()
+opar <- par('xpd'=TRUE)
+tm <- seq(5,25,by=5)
+tmc <- (tm - mean(X[,1]))
+yc <- scale(X[,1],scale=FALSE)
+Calibrate.X1 <- calibrate(c(1,0),yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.5,
+ axislab="X_1",cex.axislab=1,showlabel=FALSE,
+shiftvec=c(0,-(m[2]-par("usr")[3])),shiftfactor=1,reverse=TRUE)
+tm <- seq(5,25,by=1); tmc <- (tm - mean(X[,1]))
+Calibrate.X1 <- calibrate(c(1,0),yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.25,
+ axislab="X_1",cex.axislab=1,showlabel=FALSE,
+shiftvec=c(0,-(m[2]-par("usr")[3])),shiftfactor=1,reverse=TRUE,
+ verb=FALSE,lm=FALSE)
+yc <- scale(X[,2],scale=TRUE)
+tm <- seq(-3,1,by=1)
+Calibrate.X2 <- calibrate(c(0,1),yc,tm,Xc[,1:2],tmlab=tm,m=m,tl=0.6,
+ axislab="X_2",cex.axislab=1,showlabel=FALSE,
+shiftvec=c(-(mean(X[,1])-par('usr')[1]),0),shiftfactor=1,verb=TRUE,lm=TRUE)
+tm <- seq(-3,1.5,by=0.1)
+Calibrate.X2 <- calibrate(c(0,1),yc,tm,Xc[,1:2],tmlab=tm,m=m,tl=0.3,
+ axislab="X_2",cex.axislab=1,showlabel=FALSE,
+shiftvec=c(-(mean(X[,1])-par('usr')[1]),0),shiftfactor=1,verb=FALSE,lm=FALSE)
+par(opar)
+
+
+###################################################
+### code chunk number 11: CalibrationGuide.Rnw:327-346
+###################################################
+# PCA and Biplot construction
+pca.results <- princomp(X,cor=FALSE)
+Fp <- pca.results$scores
+Gs <- pca.results$loadings
+plot(Fp[,1],Fp[,2],pch=16,asp=1,xlab="PC 1",ylab="PC 2",cex=0.5)
+textxy(Fp[,1],Fp[,2],rownames(X),cex=0.75)
+arrows(0,0,15*Gs[,1],15*Gs[,2],length=0.1)
+textxy(15*Gs[,1],15*Gs[,2],colnames(X),cex=0.75)
+# Calibration of X_3
+ticklab <- seq(5,30,by=5)
+ticklabc <- ticklab-mean(X[,3])
+yc <- (X[,3]-mean(X[,3]))
+g <- Gs[3,1:2]
+Calibrate.X3 <- calibrate(g,yc,ticklabc,Fp[,1:2],ticklab,tl=0.5,
+ axislab="X3",cex.axislab=0.75,where=1,labpos=4)
+ticklab <- seq(5,30,by=1)
+ticklabc <- ticklab-mean(X[,3])
+Calibrate.X3.fine <- calibrate(g,yc,ticklabc,Fp[,1:2],ticklab,lm=FALSE,tl=0.25,
+ verb=FALSE,cex.axislab=0.75,where=1,labpos=4)
+
+
+###################################################
+### code chunk number 12: CalibrationGuide.Rnw:353-376
+###################################################
+# PCA and Biplot construction
+pca.results <- princomp(X,cor=TRUE)
+Fp <- pca.results$scores
+Ds <- diag(pca.results$sdev)
+Fs <- Fp%*%solve(Ds)
+Gs <- pca.results$loadings
+Gp <- Gs%*%Ds
+#plot(Fs[,1],Fs[,2],pch=16,asp=1,xlab="PC 1",ylab="PC 2",cex=0.5)
+#textxy(Fs[,1],Fs[,2],rownames(X))
+plot(Gp[,1],Gp[,2],pch=16,cex=0.5,xlim=c(-1,1),ylim=c(-1,1),asp=1,
+ xlab="1st principal axis",ylab="2nd principal axis")
+arrows(0,0,Gp[,1],Gp[,2],length=0.1)
+textxy(Gp[,1],Gp[,2],colnames(X),cex=0.75)
+ticklab <- c(seq(-1,-0.2,by=0.2),seq(0.2,1.0,by=0.2))
+R <- cor(X)
+y <- R[,5]
+g <- Gp[5,1:2]
+Calibrate.X5 <- calibrate(g,y,ticklab,Gp[,1:2],ticklab,lm=TRUE,tl=0.05,dp=TRUE,
+ labpos=2,cex.axislab=0.75,axislab="X_5")
+ticklab <- seq(-1,1,by=0.1)
+Calibrate.X5 <- calibrate(g,y,ticklab,Gp[,1:2],ticklab,lm=FALSE,tl=0.05,verb=FALSE)
+ticklab <- seq(-1,1,by=0.01)
+Calibrate.X5 <- calibrate(g,y,ticklab,Gp[,1:2],ticklab,lm=FALSE,tl=0.025,verb=FALSE)
+
+
+###################################################
+### code chunk number 13: CalibrationGuide.Rnw:389-391
+###################################################
+print(R)
+print(cbind(R[,5],Calibrate.X5$yt))
+
+
+###################################################
+### code chunk number 14: CalibrationGuide.Rnw:430-472
+###################################################
+library(MASS)
+data(calves)
+ca.results <- corresp(calves,nf=2)
+Fs <- ca.results$rscore
+Gs <- ca.results$cscore
+Ds <- diag(ca.results$cor)
+Fp <- Fs%*%Ds
+Gp <- Gs%*%Ds
+plot(Gs[,1],Gs[,2],pch=16,asp=1,cex=0.5,xlab="1st principal axis",
+ ylab="2nd principal axis")
+textxy(Gs[,1],Gs[,2],colnames(calves),cex=0.75)
+points(Fp[,1],Fp[,2],pch=16,cex=0.5)
+textxy(Fp[,1],Fp[,2],rownames(calves),cex=0.75)
+origin()
+arrows(0,0,Gs[,1],Gs[,2])
+
+P <- as.matrix(calves/sum(calves))
+r <- apply(P,1,sum)
+k <- apply(P,2,sum)
+Dc <- diag(k)
+Dr <- diag(r)
+
+RP <- solve(Dr)%*%P
+print(RP)
+CRP <- RP - ones(nrow(RP), 1) %*% t(k)
+TCRP <- CRP%*%solve(Dc)
+
+y <- TCRP[,3]
+g <- Gs[3,1:2]
+
+ticklab <- c(0,seq(0,1,by=0.2))
+ticklabs <- (ticklab - k[3])/k[3]
+Calibrate.AI <- calibrate(g,y,ticklabs,Fp[,1:2],ticklab,lm=TRUE,tl=0.10,
+ weights=Dr,axislab="AI",labpos=4,dp=TRUE)
+ticklab <- c(0,seq(0,1,by=0.1))
+ticklabs <- (ticklab - k[3])/k[3]
+Calibrate.AI <- calibrate(g,y,ticklabs,Fp[,1:2],ticklab,lm=FALSE,tl=0.10,
+ weights=Dr,verb=FALSE)
+ticklab <- c(0,seq(0,1,by=0.01))
+ticklabs <- (ticklab - k[3])/k[3]
+Calibrate.AI <- calibrate(g,y,ticklabs,Fp[,1:2],ticklab,lm=FALSE,tl=0.05,
+ weights=Dr,verb=FALSE)
+
+
+###################################################
+### code chunk number 15: CalibrationGuide.Rnw:490-521
+###################################################
+data(heads)
+X <- cbind(heads$X1,heads$X2)
+Y <- cbind(heads$Y1,heads$Y2)
+
+Rxy<- cor(X,Y)
+Ryx<- t(Rxy)
+Rxx<- cor(X)
+Ryy<- cor(Y)
+
+cca.results <-canocor(X,Y)
+
+plot(cca.results$Gs[,1],cca.results$Gs[,2],pch=16,asp=1,xlim=c(-1,1),ylim=c(-1,1),
+ xlab=expression(V[1]),ylab=expression(V[2]))
+arrows(0,0,cca.results$Fp[,1],cca.results$Fp[,2],length=0.1)
+arrows(0,0,cca.results$Gs[,1],cca.results$Gs[,2],length=0.1)
+
+textxy(cca.results$Fp[1,1],cca.results$Fp[1,2],expression(X[1]),cex=0.75)
+textxy(cca.results$Fp[2,1],cca.results$Fp[2,2],expression(X[2]),cex=0.75)
+
+textxy(cca.results$Gs[1,1],cca.results$Gs[1,2],expression(Y[1]),cex=0.75)
+textxy(cca.results$Gs[2,1],cca.results$Gs[2,2],expression(Y[2]),cex=0.75)
+
+circle(1)
+ticklab <- seq(-1,1,by=0.2)
+
+y <- Rxy[,2]
+g <- cca.results$Gs[2,1:2]
+
+Cal.Cor.Y2 <- calibrate(g,y,ticklab,cca.results$Fp[,1:2],ticklab,lm=TRUE,tl=0.05,
+ dp=TRUE,reverse=TRUE,weights=solve(Rxx),
+axislab="Y_2",cex.axislab=0.75,showlabel=FALSE)
+
+
+###################################################
+### code chunk number 16: CalibrationGuide.Rnw:524-557
+###################################################
+
+plot(cca.results$Gs[,1],cca.results$Gs[,2],pch=16,asp=1,xlim=c(-2,2),ylim=c(-2,2),
+ xlab=expression(V[1]),ylab=expression(V[2]))
+#arrows(0,0,cca.results$Fp[,1],cca.results$Fp[,2],length=0.1)
+#arrows(0,0,cca.results$Gs[,1],cca.results$Gs[,2],length=0.1)
+
+textxy(cca.results$Fp[1,1],cca.results$Fp[1,2],expression(X[1]))
+textxy(cca.results$Fp[2,1],cca.results$Fp[2,2],expression(X[2]))
+
+textxy(cca.results$Gs[1,1],cca.results$Gs[1,2],expression(Y[1]))
+textxy(cca.results$Gs[2,1],cca.results$Gs[2,2],expression(Y[2]))
+
+points(cca.results$V[,1],cca.results$V[,2],pch=16,cex=0.5)
+textxy(cca.results$V[,1],cca.results$V[,2],1:nrow(X),cex=0.75)
+
+
+ticklab <- seq(135,160,by=5)
+ticklabc <- ticklab-mean(Y[,2])
+ticklabs <- (ticklab-mean(Y[,2]))/sqrt(var(Y[,2]))
+
+y <- (Y[,2]-mean(Y[,2]))/sqrt(var(Y[,2]))
+Fr <- cca.results$V[,1:2]
+g <- cca.results$Gs[2,1:2]
+
+#points(cca.results$V[,1],cca.results$V[,2],cex=0.5,pch=19,col="red")
+#textxy(cca.results$V[,1],cca.results$V[,2],rownames(Xn))
+
+Cal.Data.Y2 <- calibrate(g,y,ticklabs,Fr,ticklab,lm=TRUE,tl=0.1,dp=TRUE,
+ reverse=TRUE,verb=TRUE,axislab="Y_2",
+ cex.axislab=0.75,showlabel=FALSE)
+
+#cca.results<-lm.gls(Rxy[,5]~-1+Fr,W=solve(Rxx))
+
+
+
+###################################################
+### code chunk number 17: CalibrationGuide.Rnw:577-594
+###################################################
+data(linnerud)
+X <- linnerud[,1:3]
+Y <- linnerud[,4:6]
+rda.results <- rda(X,Y)
+plot(rda.results$Fs[,1],rda.results$Fs[,2],pch=16,asp=1,xlim=c(-2,2),ylim=c(-2,2),
+ cex=0.5,xlab="1st principal axis",ylab="2nd principal axis")
+arrows(0,0,2*rda.results$Gyp[,1],2*rda.results$Gyp[,2],length=0.1)
+textxy(rda.results$Fs[,1],rda.results$Fs[,2],rownames(X),cex=0.75)
+textxy(2*rda.results$Gyp[,1],2*rda.results$Gyp[,2],colnames(Y),cex=0.75)
+
+y <- rda.results$Yh[,3]
+g <- rda.results$Gyp[3,1:2]
+Fr <- rda.results$Fs[,1:2]
+
+ticklab <- c(seq(-0.6,-0.1,by=0.1),seq(0.1,0.6,by=0.1))
+Calibrate.Yhat3 <- calibrate(g,y,ticklab,Fr,ticklab,lm=TRUE,dp=TRUE,tl=0.1,
+ axislab="Sauts",showlabel=FALSE)
+
+
+###################################################
+### code chunk number 18: CalibrationGuide.Rnw:597-621
+###################################################
+plot(rda.results$Gxs[,1],rda.results$Gxs[,2],pch=16,asp=1,xlim=c(-2,2),
+ ylim=c(-2,2),cex=0.5,xlab="1st principal axis",
+ylab="2nd principal axis")
+arrows(0,0,rda.results$Gxs[,1],rda.results$Gxs[,2],length=0.1)
+arrows(0,0,rda.results$Gyp[,1],rda.results$Gyp[,2],length=0.1)
+textxy(rda.results$Gxs[,1],rda.results$Gxs[,2],colnames(X),cex=0.75)
+textxy(rda.results$Gyp[,1],rda.results$Gyp[,2],colnames(Y),cex=0.75)
+
+y <- rda.results$B[,3]
+g <- rda.results$Gyp[3,1:2]
+Fr <- rda.results$Gxs[,1:2]
+
+ticklab <- seq(-0.4,0.4,0.2)
+
+W <-cor(X)
+
+Calibrate.Y3 <- calibrate(g,y,ticklab,Fr,ticklab,lm=TRUE,dp=TRUE,tl=0.1,
+ weights=W,axislab="Sauts",showlabel=FALSE)
+ticklab <- seq(-0.4,0.4,0.1)
+Calibrate.Y3 <- calibrate(g,y,ticklab,Fr,ticklab,lm=FALSE,tl=0.05,verb=FALSE,
+ weights=W)
+ticklab <- seq(-0.4,0.4,0.01)
+Calibrate.Y3 <- calibrate(g,y,ticklab,Fr,ticklab,lm=FALSE,tl=0.025,verb=FALSE,
+ weights=W)
+
+
diff --git a/inst/doc/CalibrationGuide.Rnw b/inst/doc/CalibrationGuide.Rnw
new file mode 100644
index 0000000..2d6cc85
--- /dev/null
+++ b/inst/doc/CalibrationGuide.Rnw
@@ -0,0 +1,777 @@
+\documentclass[a4paper]{article}
+
+% \VignetteIndexEntry{A R package for calibration of biplot and scatterplot axis.}
+% \VignetteDepends{graphics,stats}
+% \VignetteKeyword{multivariate}
+
+% Documentation for the Calibrate package
+
+\usepackage[english]{babel}
+\usepackage{Sweave}
+%\usepackage{Rd}
+\usepackage{url}
+\usepackage{wrapfig}
+\usepackage{hyperref}
+\usepackage[utf8]{inputenc}
+
+\setlength{\parindent}{0cm}
+
+\newcommand{\bF}{\ensuremath{\mathbf F}}
+\newcommand{\bG}{\ensuremath{\mathbf G}}
+
+\begin{document}
+
+\begin{center}
+\sf
+{\sf \bf \Large A Guide to Scatterplot and Biplot Calibration}\\
+\vspace{4mm}
+{\sf \normalsize {\tt version 1.7.2}}\\%\normalsize
+\vspace{4mm}
+{\bf \large Jan Graffelman}\\
+\vspace{4mm} \rm \large
+Department of Statistics and Operations Research\\
+Universitat Polit\`ecnica de Catalunya\\
+Avinguda Diagonal 647, 08028 Barcelona, Spain.\\
+{\it email:} jan.graffelman at upc.edu\\
+\vspace{4mm}
+{\sc September 2012}
+\end{center}
+
+\section{Introduction}
+
+This guide gives detailed instructions on how to calibrate axes in scatterplots
+and biplots obtained in the statistical environment R~\cite{RRR} by using the package
+{\tt calibrate}. By calibration we
+refer to the procedure of drawing a (linear) scale along an axis in a plot with
+tick marks and numeric labels. In an ordinary scatter plot of two variables $x$ and $y$
+two calibrated perpendicular scales are typically automatically produced by the
+routine used for plotting the two variables. However, scatter plots can be
+extended with additional variables that are represented on oblique additional
+axes. The software described in this guide can be used to create calibrated
+scales on these oblique additional axes. Moreover, in a multivariate setting with more
+than two variables, raw data matrices, correlation matrices, contingency tables,
+regression coefficients, etc. are often represented graphically by means of biplots~\cite{Gabriel}. Biplots
+also contain oblique axes representing variables. The described software can also be
+used to construct scales on biplot axes.\\
+
+The outline of this guide is as follows. In Section~\ref{sec:install} we indicate how the
+R package {\tt calibrate} can be installed. Section~\ref{sec:scatter} describes in detail
+how to calibrate additional axes in scatter plots. Section~\ref{sec:biplot} treats the
+calibration of biplot axes. Several subsections follow with detailed
+instructions of how to calibrate biplot axis in principal component analysis
+(PCA, Section~\ref{sec:pca}), correspondence analysis (CA. Section~\ref{sec:ca}),
+canonical correlation analysis (CCA, Section~\ref{sec:cca}) and redundancy analysis (RDA, Section~\ref{sec:rda}).
+The online documentation of the main routine for calibration {\tt calibrate} is referenced in
+Section~\ref{sec:online}.\\
+
+This guide does not provide the theory for the construction of scales on scatterplot
+and biplot axes. For a theoretical account of biplot calibration, we refer to Graffelman
+\& van Eeuwijk~(2005) and to Gower and Hand~(1996). If you appreciate
+this software then please cite the following paper in your work:\\
+
+Graffelman, J. \& van Eeuwijk, F.A. (2005) Calibration of multivariate scatter plots for
+exploratory analysis of relations within and between sets of variables in genomic research
+{\it Biometrical Journal}, {\bf 47}(6) pp. 863-879. \href{http://dx.doi.org/10.1002/bimj.200510177}{(clic here to access the paper)}
+
+\section{Installation}
+\label{sec:install}
+
+<<fig=FALSE,echo=FALSE>>=
+prefig <- function(){
+data(goblets)
+X <- goblets
+plot(X[,1],X[,2],pch=19,cex=0.5,xlab=expression(X[1]),
+ylab=expression(X[2]),xlim=c(5,25),ylim=c(5,25),asp=1)
+m <- apply(cbind(X[,1],X[,2]),2,mean)
+textxy(X[,1],X[,2],1:25,m=m,cex=0.75)
+points(m[1],m[2],col="red",pch=19,cex=0.5)
+abline(h = m[2])
+abline(v = m[1])
+}
+options("SweaveHooks"=list(aap=prefig))
+options("width"=60)
+@
+
+Packages in R can be installed inside the program with the option "Packages"
+in the main menu and then choosing "Install package" and picking the package
+"calibrate". Typing:
+
+<<noot>>=
+library(calibrate)
+@
+
+will, among others, make the function {\tt calibrate, canocor} and {\tt rda} available. Several
+small data sets, also the ones used in this document, are included in the package ({\tt calves, goblets,
+heads, linnerud} and {\tt storks}).
+
+\section{Calibration of Scatterplot axes}
+\label{sec:scatter}
+
+We consider a archaeological data set concerning 6 size measurements ($X_1, \ldots, X_6$) on 25
+goblets. This data was published by Manly~(1989). The data can be loaded with
+the {\tt data} instruction.
+<<fig=FALSE,echo=TRUE>>=
+data(goblets)
+X <- goblets
+@
+
+\subsection*{Oblique additional axes in a scatterplot}
+
+We construct a scatterplot of $X_1$ versus $X_2$ and center a set of coordinate
+axes on the point $(\overline{x}_1,\overline{x}_2)$ with the function {\tt origin}.
+
+\setkeys{Gin}{width=\textwidth}
+
+<<fig=TRUE,echo=TRUE>>=
+plot(X[,1],X[,2],pch=19,cex=0.5,xlab=expression(X[1]),ylab=expression(X[2]),
+ xlim=c(5,25),ylim=c(5,25),asp=1)
+m <- apply(X[,1:2],2,mean)
+textxy(X[,1],X[,2],1:25,m=m,cex=0.75)
+origin(m)
+@
+
+Next, we perform the regression of $X_5$ onto $X_1$ and $X_2$ (all variables being centered) in order to
+obtain an additional axis for $X_5$. We represent $X_5$ in the plot as a simple
+arrow whose coordinates are given by the regression coefficients:
+
+<<fig=TRUE,echo=TRUE,aap=TRUE>>=
+Xc <- scale(X,center=TRUE,scale=FALSE)
+b <- solve(t(Xc[,1:2])%*%Xc[,1:2])%*%t(Xc[,1:2])%*%Xc[,5]
+print(b)
+bscaled <- 20*b
+arrows(m[1],m[2],m[1]+bscaled[1],m[2]+bscaled[2],col="blue",length=0.1)
+arrows(m[1],m[2],m[1]-bscaled[1],m[2]-bscaled[2],length=0,lty="dashed",col="blue")
+@
+
+A direction that is optimal in the least squares sense for $X_5$ is given by the vector of regression
+coefficients~\cite{Graffel13}. To make this direction more visible, we multiplied it by a constant (20).
+It is clear that the direction of increase for $X_5$ runs approximately North-East across the scatterplot.
+We now proceed to calibrate this direction with a scale for $X_5$. In order to choose sensible values for
+the scale of $X_5$, we first inspect the range of variation of $X_5$, and then choose a set of values we want
+to mark off on the scale ({\tt tm}) and also compute the deviations of these values from the mean
+({\tt tmc}). We specify a tick length of 0.3 ({\tt tl=0.3}). Depending on the data, some values of {\tt tl}
+typically have to be tried to see how to obtain a nice scale.
+
+<<fig=TRUE,echo=TRUE,aap=TRUE>>=
+print(range(X[,5]))
+yc <- scale(X[,5],scale=FALSE)
+tm <- seq(2,10,by=1)
+tmc <- tm - mean(X[,5])
+Calibrate.X5<-calibrate(b,yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.3,axislab="X_5",
+ labpos=4,cex.axislab=1)
+@
+
+The numerical output from routine {\tt calibrate} shows that one unit along the axis for $X_5$ occupies 2.47
+units in the plotting frame. The axis for $X_5$ makes an angle of 17.65 degrees with the positive x-axis.
+The calibration factor is 6.12. Multiplying the vector of regressions coefficients by this
+factor yields a vector that represents a unit change in the scale of $X_5$. E.g. for this data we have that
+the vector $6.12 \cdot (0.385, 0.123) = (2.358, 0.751)$ represents a unit change. This vector has
+norm $\sqrt{2.358^2 + 0.751^2} = 2.47$. Other calibration factors may be specified by using parameter
+{\tt alpha}. If {\tt alpha} is left unspecified the optimal value computed by least
+squares will be used. The goodness-of-fit of $X_5$ is 0.513. This means that 51.3\% of the variance
+of $X_5$ can be explained by a regression onto $X_1$ and $X_2$ ($R^2 = 0.513$). The goodness-of-scale has
+the same value. The goodness-of-scale is only relevant if we modify parameter {\tt alpha}. {\tt Calibrate.X5}
+is a list object containing all calibration results (calibration factor, fitted values according to the
+scale used, tick marker positions, etc.).
+
+\subsubsection*{Shifting a calibrated axis}
+
+Using many calibrated axes in a plot, all passing through the origin, leads to dense plots that become
+unreadable. It is therefore a good idea to shift calibrated axes towards the margins of the plot. This keeps
+the central cloud of data points clear and relegates all information on scales to the margins of the graph.
+There are two natural positions for a shifted axis: just above the largest data point in a direction perpendicular
+to the axis being calibrated, or just below the smallest data point in the perpendicular direction. The arguments
+{\tt shiftdir, shiftfactor} and {\tt shiftvec} can be used to control the shifting of a calibrated axis. {\tt
+shiftvec} allows the user to specify the shift vector manually. This is normally not needed, and good positions
+for an axis can be found by using only {\tt shiftdir} and {\tt shiftfactor}. Argument {\tt shiftdir} can be set
+to {\tt 'right'} or {\tt 'left'} and indicates in which direction the axis is to be shifted, with respect to the direction of
+increase of the calibrated axis. Setting {\tt shiftdir} shifts the axis automatically just above or below
+the most outlying data point in the direction perpendicular to the vector being calibrated. In order to
+move the calibrated axis farther out or to pull it more in, {\tt shiftfactor} can be used.
+Argument {\tt shiftfactor} stretches or shrinks the shift vector for the axis. A {\tt shiftfactor} larger
+than 1 moves the axis outwards, and a {\tt shiftfactor} smaller than 1 pulls the axis towards the origin of
+the plot. If set to 1 exactly, the shifted axis will cut through the most outlying data point.
+The default {\tt shiftfactor} is 1.05. We redo
+the previous plot, shifting the calibrated axis below the cloud of points, which is to the right w.r.t. the
+direction of increase of the variable.
+
+<<fig=TRUE,echo=TRUE,aap=TRUE>>=
+yc <- scale(X[,5],scale=FALSE)
+tm <- seq(2,10,by=1)
+tmc <- tm - mean(X[,5])
+Calibrate.X5<-calibrate(b,yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.3,axislab="X_5",labpos=4,
+ cex.axislab=1,shiftdir="right")
+@
+
+The shift of the axis does not affect the interpretation of the plot, because the projections of the points
+onto the axis remain the same.
+
+\subsubsection*{Second vertical axis in a scatterplot}
+
+The oblique direction in the previous section is the preferred direction for $X_5$, as this direction
+is optimal in the least squares sense. However, if desired, additional variables can also be represented
+as a second vertical axis on the right of the plot, or as a second horizontal axis on the top of the
+plot. We now proceed to construct a second vertical axis on the right hand of the scatter plot for
+$X_5$. This can be done by setting the vector to be calibrated (first argument of routine {\tt calibrate})
+to the (0,1) vector. By specifying a shiftvector explicitly ({\tt shiftvec}), the axis can be shifted. For
+this data, setting {\tt shiftvec} to {\tt c(par('usr')[2]-mean(X[,1]),0)} and {\tt shiftfactor = 1},
+makes the axis coincide with the right vertical borderline of the graph.
+
+<<fig=TRUE,echo=TRUE,aap=TRUE>>=
+opar <- par('xpd'=TRUE)
+tm <- seq(3,8,by=1)
+tmc <- (tm - mean(X[,5]))
+Calibrate.rightmargin.X5 <- calibrate(c(0,1),yc,tmc,Xc[,1:2],tmlab=tm,m=m,
+ axislab="X_5",tl=0.5,
+ shiftvec=c(par('usr')[2]-mean(X[,1]),0),
+ shiftfactor=1,where=2,
+ laboffset=c(1.5,1.5),cex.axislab=1)
+par(opar)
+@
+
+The second vertical axis has calibration factor 3.46, and a goodness of fit of 0.34. The fit of the
+variable is worse in comparison with the previous oblique direction given by the regression coefficients. Note
+that graphical clipping in temporarily turned off ({\tt par('xpd'=TRUE)}) to allow the calibration
+routine to draw ticks and labels outside the figure region, and that the range of the tick marks was
+shortened in order not surpass the figure region.
+
+\subsubsection*{Subscales and double calibrations}
+
+Scales with tick marks can be refined by drawing subscales with smaller tick marks.
+E.g. larger labelled
+tickmarks can be used to represent multiples of 10, and small unlabelled tick marks can be used to
+represent units. The subscale allows a more precise recovery of the data values. This can simply be
+achieved by calling the calibration routine twice, once with a coarse sequence and once with a finer
+sequence. For the second call one can specify {\tt verb=FALSE} in order to suppress the numerical output
+of the routine, and {\tt lm=FALSE} to supress the tick mark labels under the smaller ticks. The tickmarks
+for the finer scale are made smaller by modifying the tick length (e.g. {\tt tl=0.1}). Depending on the data,
+some trial and error with different values for {\tt tl} may be necessary before nice scales are obtained. This
+may be automatized in the future. Finally, reading off the (approximate) data values can further be enhanced
+by drawing perpendiculars from the points to the calibrated axis by setting {\tt dp=TRUE}.
+
+<<fig=TRUE,echo=TRUE,aap=TRUE>>=
+tm <- seq(2,10,by=1)
+tmc <- (tm - mean(X[,5]))
+Calibrate.X5 <- calibrate(b,yc,tmc,Xc[,1:2],tmlab=tm,m=m,axislab="X_5",tl=0.5,
+ dp=TRUE,labpos=4)
+tm <- seq(2,10,by=0.1)
+tmc <- (tm - mean(X[,5]))
+Calibrate.X5 <- calibrate(b,yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.25,verb=FALSE,
+ lm=FALSE)
+@
+
+A {\it double calibration} can be created by drawing two scales, one on each side of the axis. Double
+calibrations can be useful. For instance, one scale can be used for recovery of the original
+data values of the variable, whereas the second scale can be used for recovery of standardized values or
+of correlations with other variables. Double calibrations can also be used to
+graphically verify if two different calibration procedures give the same result or not.
+
+\subsubsection*{Recalibrating the original scatterplot axes}
+
+By calibrating the (0,1) and (1,0) vectors the original axes of the scatter plot can be
+redesigned. We illustrate the recalibration of the original axes by creating a second scale on
+the other side of the axes, a refined scale for $X_1$, and a scale for the standardized data for
+$X_2$. For the latter calibration one unit equals one standard deviation.
+
+<<fig=TRUE,echo=TRUE,aap=TRUE>>=
+opar <- par('xpd'=TRUE)
+tm <- seq(5,25,by=5)
+tmc <- (tm - mean(X[,1]))
+yc <- scale(X[,1],scale=FALSE)
+Calibrate.X1 <- calibrate(c(1,0),yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.5,
+ axislab="X_1",cex.axislab=1,showlabel=FALSE,
+shiftvec=c(0,-(m[2]-par("usr")[3])),shiftfactor=1,reverse=TRUE)
+tm <- seq(5,25,by=1); tmc <- (tm - mean(X[,1]))
+Calibrate.X1 <- calibrate(c(1,0),yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.25,
+ axislab="X_1",cex.axislab=1,showlabel=FALSE,
+shiftvec=c(0,-(m[2]-par("usr")[3])),shiftfactor=1,reverse=TRUE,
+ verb=FALSE,lm=FALSE)
+yc <- scale(X[,2],scale=TRUE)
+tm <- seq(-3,1,by=1)
+Calibrate.X2 <- calibrate(c(0,1),yc,tm,Xc[,1:2],tmlab=tm,m=m,tl=0.6,
+ axislab="X_2",cex.axislab=1,showlabel=FALSE,
+shiftvec=c(-(mean(X[,1])-par('usr')[1]),0),shiftfactor=1,verb=TRUE,lm=TRUE)
+tm <- seq(-3,1.5,by=0.1)
+Calibrate.X2 <- calibrate(c(0,1),yc,tm,Xc[,1:2],tmlab=tm,m=m,tl=0.3,
+ axislab="X_2",cex.axislab=1,showlabel=FALSE,
+shiftvec=c(-(mean(X[,1])-par('usr')[1]),0),shiftfactor=1,verb=FALSE,lm=FALSE)
+par(opar)
+@
+
+\section{Calibration of Biplot axes}
+\label{sec:biplot}
+
+In this section we give detailed instructions on how to calibrate biplot axes. We will consider biplots
+of raw data matrices and correlation matrices obtained by PCA, biplots of profiles obtained in CA,
+biplots of data matrices and correlation matrices (in particular the between-set correlation matrix)
+in CCA and biplots of fitted values and regression coefficients obtained by RDA. In principle, calibration
+of biplot axes has little additional complication in comparison with the calibration of additional
+axes in scatterplots explained above. The main issue is that, prior to calling the calibration routine,
+one needs to take care of the proper centring and standardisation of the tick marks.
+
+\subsection{Principal component analysis}
+\label{sec:pca}
+
+Principal component analysis can be performed by using routine {\tt princomp} from the {\tt stats }
+library. We use again Manly's goblets data to create a biplot of the data based on a
+PCA of the covariance matrix. We use {\tt princomp} to compute the scores for the rows and the columns
+of the data matrix. The first principal component is seen to be a size component, separating the
+smaller goblets on the right from the larger goblets on the left. The variable vectors are
+multiplied by a factor of 15 to facilitate interpretation. Next we
+calibrate the vector for $X_3$,
+using labelled tickmarks for multiples of 5 units, and shorter unlabelled tickmarks for the units. The
+goodness of fit of $X_3$ is very high (0.99), which means that $X_3$ is close to perfectly
+represented. {\tt Calibrate.X3} is a list object containing the numerical results of the
+calibration.
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+# PCA and Biplot construction
+pca.results <- princomp(X,cor=FALSE)
+Fp <- pca.results$scores
+Gs <- pca.results$loadings
+plot(Fp[,1],Fp[,2],pch=16,asp=1,xlab="PC 1",ylab="PC 2",cex=0.5)
+textxy(Fp[,1],Fp[,2],rownames(X),cex=0.75)
+arrows(0,0,15*Gs[,1],15*Gs[,2],length=0.1)
+textxy(15*Gs[,1],15*Gs[,2],colnames(X),cex=0.75)
+# Calibration of X_3
+ticklab <- seq(5,30,by=5)
+ticklabc <- ticklab-mean(X[,3])
+yc <- (X[,3]-mean(X[,3]))
+g <- Gs[3,1:2]
+Calibrate.X3 <- calibrate(g,yc,ticklabc,Fp[,1:2],ticklab,tl=0.5,
+ axislab="X3",cex.axislab=0.75,where=1,labpos=4)
+ticklab <- seq(5,30,by=1)
+ticklabc <- ticklab-mean(X[,3])
+Calibrate.X3.fine <- calibrate(g,yc,ticklabc,Fp[,1:2],ticklab,lm=FALSE,tl=0.25,
+ verb=FALSE,cex.axislab=0.75,where=1,labpos=4)
+@
+
+We do a PCA based on the correlation matrix, and proceed to construct a biplot of the correlation matrix. The
+correlations of $X_5$ with the other variables are computed, and the biplot axis for $X_5$ is calibrated with a
+correlation scale. Routine {\tt calibrate} is repeatedly called to create finer subscales.
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+# PCA and Biplot construction
+pca.results <- princomp(X,cor=TRUE)
+Fp <- pca.results$scores
+Ds <- diag(pca.results$sdev)
+Fs <- Fp%*%solve(Ds)
+Gs <- pca.results$loadings
+Gp <- Gs%*%Ds
+#plot(Fs[,1],Fs[,2],pch=16,asp=1,xlab="PC 1",ylab="PC 2",cex=0.5)
+#textxy(Fs[,1],Fs[,2],rownames(X))
+plot(Gp[,1],Gp[,2],pch=16,cex=0.5,xlim=c(-1,1),ylim=c(-1,1),asp=1,
+ xlab="1st principal axis",ylab="2nd principal axis")
+arrows(0,0,Gp[,1],Gp[,2],length=0.1)
+textxy(Gp[,1],Gp[,2],colnames(X),cex=0.75)
+ticklab <- c(seq(-1,-0.2,by=0.2),seq(0.2,1.0,by=0.2))
+R <- cor(X)
+y <- R[,5]
+g <- Gp[5,1:2]
+Calibrate.X5 <- calibrate(g,y,ticklab,Gp[,1:2],ticklab,lm=TRUE,tl=0.05,dp=TRUE,
+ labpos=2,cex.axislab=0.75,axislab="X_5")
+ticklab <- seq(-1,1,by=0.1)
+Calibrate.X5 <- calibrate(g,y,ticklab,Gp[,1:2],ticklab,lm=FALSE,tl=0.05,verb=FALSE)
+ticklab <- seq(-1,1,by=0.01)
+Calibrate.X5 <- calibrate(g,y,ticklab,Gp[,1:2],ticklab,lm=FALSE,tl=0.025,verb=FALSE)
+@
+
+The goodness of fit of the representation of the correlations of $X_5$ with the
+other variables is 0.98, the 6 correlations being close to perfectly represented.
+We compute the sample correlation matrix and compare the observed correlations of
+$X_5$ with those estimated from the calibrated biplot axis ({\tt yt}). Note that
+PCA also tries to approximate the correlation of a variable with itself, and that
+the arrow on representing $X_5$ falls short of the value 1 on its own calibrated scale.
+The refined subscale allows very precise graphical representation of the correlations
+as estimated by the biplot.
+
+\begin{center}
+<<fig=FALSE,echo=TRUE,aap=FALSE>>=
+print(R)
+print(cbind(R[,5],Calibrate.X5$yt))
+@
+\end{center}
+
+\subsection{Correspondence analysis}
+\label{sec:ca}
+
+We consider a contingency table of a sample of Dutch calves born in the late nineties,
+shown in Table~\ref{tab:calves2}. A total of 7257 calves were classified according
+to two categorical variables: the
+method of production (ET = Embryo Transfer, IVP = In Vitro Production, AI = Artificial
+Insemination) and the ease of delivery, scored on a scale from 1 (normal) to 6 (very heavy).
+The data in Table~\ref{tab:calves2} were provided by Holland Genetics.
+\begin{table}[htb]
+\centering
+\begin{tabular}{c|rrrr}
+ & \multicolumn{3}{c}{Type of calf}\\
+Ease of delivery & ET & IVP & AI\\
+\hline
+1 & 97 & 150 & 1686\\
+2 & 152 & 183 & 1339\\
+3 & 377 & 249 & 1209\\
+4 & 335 & 227 & 656\\
+5 & 42 & 136 & 277\\
+6 & 9 & 71 & 62\\
+\hline
+\end{tabular}
+\caption{Calves data from Holland Genetics.}
+\label{tab:calves2}
+\end{table}
+
+For this contingency table we obtain $\chi^2_{10} = 833.16$ with $p < 0.001$ and the null
+hypothesis of no association between ease of delivery and type of calf has to be rejected. However, what is
+the precise nature of this association? Correspondence analysis can be used to gain insight in the
+nature of this association. We use routine {\tt corresp} form the {\tt MASS} library~\cite{Venables}
+to perform correspondence analysis and to obtain the coordinates for a biplot of the row profiles.
+We compute the row profiles and then repeatedly call the calibration routine, each time with a
+different set of {\tt ticklabs}.
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+library(MASS)
+data(calves)
+ca.results <- corresp(calves,nf=2)
+Fs <- ca.results$rscore
+Gs <- ca.results$cscore
+Ds <- diag(ca.results$cor)
+Fp <- Fs%*%Ds
+Gp <- Gs%*%Ds
+plot(Gs[,1],Gs[,2],pch=16,asp=1,cex=0.5,xlab="1st principal axis",
+ ylab="2nd principal axis")
+textxy(Gs[,1],Gs[,2],colnames(calves),cex=0.75)
+points(Fp[,1],Fp[,2],pch=16,cex=0.5)
+textxy(Fp[,1],Fp[,2],rownames(calves),cex=0.75)
+origin()
+arrows(0,0,Gs[,1],Gs[,2])
+
+P <- as.matrix(calves/sum(calves))
+r <- apply(P,1,sum)
+k <- apply(P,2,sum)
+Dc <- diag(k)
+Dr <- diag(r)
+
+RP <- solve(Dr)%*%P
+print(RP)
+CRP <- RP - ones(nrow(RP), 1) %*% t(k)
+TCRP <- CRP%*%solve(Dc)
+
+y <- TCRP[,3]
+g <- Gs[3,1:2]
+
+ticklab <- c(0,seq(0,1,by=0.2))
+ticklabs <- (ticklab - k[3])/k[3]
+Calibrate.AI <- calibrate(g,y,ticklabs,Fp[,1:2],ticklab,lm=TRUE,tl=0.10,
+ weights=Dr,axislab="AI",labpos=4,dp=TRUE)
+ticklab <- c(0,seq(0,1,by=0.1))
+ticklabs <- (ticklab - k[3])/k[3]
+Calibrate.AI <- calibrate(g,y,ticklabs,Fp[,1:2],ticklab,lm=FALSE,tl=0.10,
+ weights=Dr,verb=FALSE)
+ticklab <- c(0,seq(0,1,by=0.01))
+ticklabs <- (ticklab - k[3])/k[3]
+Calibrate.AI <- calibrate(g,y,ticklabs,Fp[,1:2],ticklab,lm=FALSE,tl=0.05,
+ weights=Dr,verb=FALSE)
+@
+
+Because the calibration is done by weighted least squares, a diagonal matrix of weights ({\tt weights=Dr})
+is supplied as a parameter to the calibration routine
+Note that the calibrated axis for the row profiles with respect to AI has goodness of fit 1. This
+is due to the fact that the rank of the matrix of centred profiles is two, and that therefore
+all profiles can be perfectly represented in two dimensional space.
+
+\subsection{Canonical correlation analysis}
+\label{sec:cca}
+
+We consider a classical data set on the head sizes of the first and the second son of 25 families~\cite{Frets}.
+These data have been analysed by several authors~\cite{Anderson,Mardia,Graffel16}
+We first load the data and perform
+a canonical correlation analysis, using supplied function {\tt canocor} (a more fully
+fledged program for canonical correlation analysis in comparison with {\tt cancor} from the {\tt stats} package).
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+data(heads)
+X <- cbind(heads$X1,heads$X2)
+Y <- cbind(heads$Y1,heads$Y2)
+
+Rxy<- cor(X,Y)
+Ryx<- t(Rxy)
+Rxx<- cor(X)
+Ryy<- cor(Y)
+
+cca.results <-canocor(X,Y)
+
+plot(cca.results$Gs[,1],cca.results$Gs[,2],pch=16,asp=1,xlim=c(-1,1),ylim=c(-1,1),
+ xlab=expression(V[1]),ylab=expression(V[2]))
+arrows(0,0,cca.results$Fp[,1],cca.results$Fp[,2],length=0.1)
+arrows(0,0,cca.results$Gs[,1],cca.results$Gs[,2],length=0.1)
+
+textxy(cca.results$Fp[1,1],cca.results$Fp[1,2],expression(X[1]),cex=0.75)
+textxy(cca.results$Fp[2,1],cca.results$Fp[2,2],expression(X[2]),cex=0.75)
+
+textxy(cca.results$Gs[1,1],cca.results$Gs[1,2],expression(Y[1]),cex=0.75)
+textxy(cca.results$Gs[2,1],cca.results$Gs[2,2],expression(Y[2]),cex=0.75)
+
+circle(1)
+ticklab <- seq(-1,1,by=0.2)
+
+y <- Rxy[,2]
+g <- cca.results$Gs[2,1:2]
+
+Cal.Cor.Y2 <- calibrate(g,y,ticklab,cca.results$Fp[,1:2],ticklab,lm=TRUE,tl=0.05,
+ dp=TRUE,reverse=TRUE,weights=solve(Rxx),
+axislab="Y_2",cex.axislab=0.75,showlabel=FALSE)
+@
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+
+plot(cca.results$Gs[,1],cca.results$Gs[,2],pch=16,asp=1,xlim=c(-2,2),ylim=c(-2,2),
+ xlab=expression(V[1]),ylab=expression(V[2]))
+#arrows(0,0,cca.results$Fp[,1],cca.results$Fp[,2],length=0.1)
+#arrows(0,0,cca.results$Gs[,1],cca.results$Gs[,2],length=0.1)
+
+textxy(cca.results$Fp[1,1],cca.results$Fp[1,2],expression(X[1]))
+textxy(cca.results$Fp[2,1],cca.results$Fp[2,2],expression(X[2]))
+
+textxy(cca.results$Gs[1,1],cca.results$Gs[1,2],expression(Y[1]))
+textxy(cca.results$Gs[2,1],cca.results$Gs[2,2],expression(Y[2]))
+
+points(cca.results$V[,1],cca.results$V[,2],pch=16,cex=0.5)
+textxy(cca.results$V[,1],cca.results$V[,2],1:nrow(X),cex=0.75)
+
+
+ticklab <- seq(135,160,by=5)
+ticklabc <- ticklab-mean(Y[,2])
+ticklabs <- (ticklab-mean(Y[,2]))/sqrt(var(Y[,2]))
+
+y <- (Y[,2]-mean(Y[,2]))/sqrt(var(Y[,2]))
+Fr <- cca.results$V[,1:2]
+g <- cca.results$Gs[2,1:2]
+
+#points(cca.results$V[,1],cca.results$V[,2],cex=0.5,pch=19,col="red")
+#textxy(cca.results$V[,1],cca.results$V[,2],rownames(Xn))
+
+Cal.Data.Y2 <- calibrate(g,y,ticklabs,Fr,ticklab,lm=TRUE,tl=0.1,dp=TRUE,
+ reverse=TRUE,verb=TRUE,axislab="Y_2",
+ cex.axislab=0.75,showlabel=FALSE)
+
+#cca.results<-lm.gls(Rxy[,5]~-1+Fr,W=solve(Rxx))
+
+@
+
+We construct the biplot of the between-set correlation matrix (the joint
+plot of ${\bF}_p$ and ${\bG}_s$).
+Firstly we calibrate the biplot axis for $Y_2$ with a correlation scale.
+This calibration is done by generalised least squares with the inverse of the correlation matrix of
+the X-variables as a weight matrix ({\tt weights=solve(Rxx)}). Secondly, we calibrate the biplot axis
+for $Y_2$ with a scale for the original values. This second calibration has no weight
+matrix and is obtained by ordinary least squares. Both calibrations have a goodness of fit of 1 and
+allow perfect recovery of correlations and original data values.
+
+\subsection{Redundancy analysis}
+\label{sec:rda}
+
+Redundancy analysis can be seen as a constrained PCA. It allows two biplots, the biplot of the fitted
+values and a biplot of regression coefficients. Function {\tt rda} of the package provides a routine
+for redundancy analysis. We use Linnerud's data on physical exercise and body measurement
+variables~\cite{Tenenhaus} to illustrate calibrated biplots in redundancy analysis.
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+data(linnerud)
+X <- linnerud[,1:3]
+Y <- linnerud[,4:6]
+rda.results <- rda(X,Y)
+plot(rda.results$Fs[,1],rda.results$Fs[,2],pch=16,asp=1,xlim=c(-2,2),ylim=c(-2,2),
+ cex=0.5,xlab="1st principal axis",ylab="2nd principal axis")
+arrows(0,0,2*rda.results$Gyp[,1],2*rda.results$Gyp[,2],length=0.1)
+textxy(rda.results$Fs[,1],rda.results$Fs[,2],rownames(X),cex=0.75)
+textxy(2*rda.results$Gyp[,1],2*rda.results$Gyp[,2],colnames(Y),cex=0.75)
+
+y <- rda.results$Yh[,3]
+g <- rda.results$Gyp[3,1:2]
+Fr <- rda.results$Fs[,1:2]
+
+ticklab <- c(seq(-0.6,-0.1,by=0.1),seq(0.1,0.6,by=0.1))
+Calibrate.Yhat3 <- calibrate(g,y,ticklab,Fr,ticklab,lm=TRUE,dp=TRUE,tl=0.1,
+ axislab="Sauts",showlabel=FALSE)
+@
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+plot(rda.results$Gxs[,1],rda.results$Gxs[,2],pch=16,asp=1,xlim=c(-2,2),
+ ylim=c(-2,2),cex=0.5,xlab="1st principal axis",
+ylab="2nd principal axis")
+arrows(0,0,rda.results$Gxs[,1],rda.results$Gxs[,2],length=0.1)
+arrows(0,0,rda.results$Gyp[,1],rda.results$Gyp[,2],length=0.1)
+textxy(rda.results$Gxs[,1],rda.results$Gxs[,2],colnames(X),cex=0.75)
+textxy(rda.results$Gyp[,1],rda.results$Gyp[,2],colnames(Y),cex=0.75)
+
+y <- rda.results$B[,3]
+g <- rda.results$Gyp[3,1:2]
+Fr <- rda.results$Gxs[,1:2]
+
+ticklab <- seq(-0.4,0.4,0.2)
+
+W <-cor(X)
+
+Calibrate.Y3 <- calibrate(g,y,ticklab,Fr,ticklab,lm=TRUE,dp=TRUE,tl=0.1,
+ weights=W,axislab="Sauts",showlabel=FALSE)
+ticklab <- seq(-0.4,0.4,0.1)
+Calibrate.Y3 <- calibrate(g,y,ticklab,Fr,ticklab,lm=FALSE,tl=0.05,verb=FALSE,
+ weights=W)
+ticklab <- seq(-0.4,0.4,0.01)
+Calibrate.Y3 <- calibrate(g,y,ticklab,Fr,ticklab,lm=FALSE,tl=0.025,verb=FALSE,
+ weights=W)
+@
+
+The first biplot shown is a biplot of the fitted values (obtained
+from the regression of Y onto X). Vectors for the response variables are multiplied by a factor of 3 to increase
+readability. The fitted values of the regression of Sauts onto the body measurements have
+a goodness of fit of 0.9984 and can very well be recovered by projection onto the calibrated
+axis. The second biplot is a biplot of the matrix of regression coefficients. We
+calibrated the biplot axis for "Sauts", such that the regression coefficients of the
+explanory variables with respect to "Sauts" can be recovered. The goodness of fit for
+"Sauts" is over 0.99, which means that the regression coefficients are close to
+perfectly displayed. Note that the calibration for Sauts for the regression coefficients
+is done by GLS with weight matrix equal to the correlation matrix of the X variables
+({\tt weights=W}).
+
+\section{Online documentation}
+\label{sec:online}
+
+Online documentation for the package can be obtained by typing
+{\tt vignette("CalibrationGuide"} or by accessing the file {\tt CalibrationGuide.pdf} in the {\tt doc}
+directory of the installed package.
+
+\section{Version history}
+
+Version 1.6:\\
+
+\begin{itemize}
+
+\item Function {\tt rad2degree} and {\tt shiftvector} have been added.
+
+\item Function calibrate has changed. Argument {\tt shift} from previous versions is obsolete,
+ and replaced by {\tt shiftdir, shiftfactor} and {\tt shiftvec}.
+
+\end{itemize}
+
+Version 1.7.2:\\
+
+\begin{itemize}
+\item Function {\tt textxy} has been modified and improved. Arguments {\tt dcol} and {\tt cx} no longer work, and their role has been taken over by {\tt col} and {\tt cex}. A new argument {\tt offset} controls the distance between point and label.
+\end{itemize}
+
+
+\section*{Acknowledgements}
+
+This work was partially supported by the Spanish grant BEC2000-0983. I thank Holland Genetics
+({\tt http://www.hg.nl/}), Janneke van Wagtendonk and Sander de Roos for making the calves data
+available. This document was generated by Sweave~\cite{Leisch}.
+
+\bibliographystyle{humanbio}
+\begin{thebibliography}{10}
+
+\bibitem[Anderson (1984)]{Anderson}
+Anderson, T. W.
+(1984)
+{A}n {I}ntroduction to {M}ultivariate {S}tatistical {A}nalysis
+John Wiley,
+Second edition,
+New York.
+
+\bibitem[Frets (1921)]{Frets}
+Frets, G. P.
+(1921)
+Heredity of head form in man,
+Genetica,
+3,
+pp. 193-384.
+
+\bibitem[Gabriel, 1971]{Gabriel}
+Gabriel, K. R.
+(1971)
+The biplot graphic display of matrices with application to principal component analysis.
+Biometrika
+58(3)
+pp. 453-467.
+
+\bibitem[Gower and Hand (1996)]{Gower4}
+Gower, J. C. and Hand, D. J.
+(1996)
+Biplots
+Chapman \& Hall,
+London.
+
+\bibitem[Graffelman (2005)]{Graffel16}
+Graffelman, J.
+(2005)
+Enriched biplots for canonical correlation analysis
+Journal of Applied Statistics
+32(2)
+pp. 173-188.
+
+\bibitem[Graffelman and Aluja-Banet (2003)]{Graffel13}
+Graffelman, J. and Aluja-Banet, T.
+(2003)
+Optimal Representation of Supplementary Variables in Biplots from Principal Component Analysis and Correspondence Analysis
+Biometrical Journal,
+45(4)
+pp. 491-509.
+
+\bibitem[Graffelman and van Eeuwijk (2005)]{Graffel17}
+Graffelman, J. and van Eeuwijk, F. A.,
+(2005)
+Calibration of multivariate scatter plots for exploratory analysis of relations within and between sets of variables
+in genomic research,
+Biometrical Journal,
+47,
+6,
+863-879.
+
+\bibitem[Leisch (2002)]{Leisch}
+Leisch, F.
+(2002)
+Sweave: Dynamic generation of statistical reports using literate data analysis
+Compstat 2002, Proceedings in Computational Statistics
+pp. 575-580,
+Physica Verlag, Heidelberg,
+ISBN 3-7908-1517-9
+URL http:/www.ci.tuwien.ac.at/~leisch/Sweave.
+
+\bibitem[Manly (1989)]{Manly}
+Manly, B. F. J.
+(1989)
+Multivariate statistical methods: a primer
+Chapman and Hall, London.
+
+\bibitem[Mardia et al.(1979)]{Mardia}
+Mardia, K. V. and Kent, J. T. and Bibby, J. M.
+(1979)
+Multivariate Analysis
+Academic Press London.
+
+\bibitem[R Development Core Team (2004)]{RRR}
+R Development Core Team
+(2004)
+R: A language and environment forstatistical computing.
+R Foundation for Statistical Computing,
+Vienna, Austria,
+ISBN 3-900051-00-3,
+http://www.R-project.org.
+
+\bibitem[Tenenhaus (1998)]{Tenenhaus}
+Tenenhaus, M.
+(1998)
+La R\'{e}gression PLS
+Paris, \'Editions Technip.
+
+\bibitem[Venables and Ripley (2002)]{Venables}
+Venables, W. N. and Ripley, B. D.
+(2002)
+{M}odern {A}pplied {S}tatistics with {S}-{P}lus
+New York,
+Fourth edition,
+Springer.
+
+
+\end{thebibliography}
+
+\end{document}
diff --git a/inst/doc/CalibrationGuide.pdf b/inst/doc/CalibrationGuide.pdf
new file mode 100644
index 0000000..525241b
Binary files /dev/null and b/inst/doc/CalibrationGuide.pdf differ
diff --git a/man/calibrate.Rd b/man/calibrate.Rd
new file mode 100644
index 0000000..42ad652
--- /dev/null
+++ b/man/calibrate.Rd
@@ -0,0 +1,96 @@
+\name{calibrate}
+\alias{calibrate}
+\title{Calibration of Biplot and Scatterplot Axis}
+\description{Routine for the calibration of any axis (variable vector) in a biplot or a scatterplot}
+\usage{
+calibrate(g,y,tm,Fr,tmlab=tm,tl=0.05,dt=TRUE,dp=FALSE,
+ lm=TRUE,verb=TRUE,axislab="",reverse=FALSE,
+ alpha=NULL,labpos=1,weights=diag(rep(1,length(y))),
+ axiscol="blue",cex.axislab=0.75,graphics=TRUE,where=3,
+ laboffset=c(0,0),m=matrix(c(0,0),nrow=1),markerpos=3,
+ showlabel=TRUE,lwd=1,shiftvec=c(0,0),shiftdir="none",shiftfactor=1.05)
+}
+\arguments{
+ \item{g}{ the vector to be calibrated (2 x 1). }
+ \item{y}{ the data vector corresponding to \code{g}, appropriately centred and/or standardized. }
+ \item{tm}{ the vector of tick marks, appropiately centred and/or scaled. }
+ \item{Fr}{ the coordinates of the rows markers in the biplot. }
+ \item{tmlab}{ a list or vector of tick mark labels. }
+ \item{tl}{ the tick length. By default, the tick markers have length 0.05.}
+ \item{dt}{ draw ticks. By default, ticks markers are drawn. Set dt=F in order to compute
+calibration results without actually drawing the calibrated scale. }
+ \item{dp}{ drop perpendiculars. With dp=T perpendicular lines will be drawn from the row markers
+ specified by Fr onto the calibrated axis. This is a graphical aid to read off the values in the
+ corresponding scale. }
+ \item{lm}{ label markers. By default, all tick marks are labelled. Setting lm=F turns off the labelling
+of the tick marks. This allows for creating tick marks without labels. It is particularly useful for
+creating finer scales of tickmarks without labels.}
+ \item{verb}{ verbose parameter (F=be quiet, T=show results). }
+ \item{axislab}{ a label for the calibrated axis. }
+ \item{reverse}{ puts the tick marks and tick mark labels on the other side of the axis. }
+ \item{alpha}{ a value for the calibration factor. This parameter should only be specified if
+ a calibration is required that is different from the one that is optimal for data recovery.}
+ \item{labpos}{ position of the label for the calibrated axis (1,2,3 or 4). }
+ \item{laboffset}{ offset vector for the axis label. If specified, shifts the label by the specified amounts with respect to the current position. }
+ \item{weights}{ a matrix of weights (optional). }
+ \item{axiscol}{ color of the calibrated axis. }
+ \item{cex.axislab}{ character expansion factor for axis label and tick mark labels. }
+ \item{graphics}{ do graphics or not (F=no graphical output, T=draws calibrated scale). }
+ \item{where}{ label placement (1=beginning,2=middle,3=end). }
+ \item{m}{ vector of means. }
+ \item{markerpos}{ position specifier for the tick mark labels (1,2,3 or 4). }
+ \item{showlabel}{ show axis label in graph (T) or not (F). }
+ \item{lwd}{ line with for the calibrated axis }
+ \item{shiftvec}{ a shift vector for the calibrated axis ((0,0) by
+ default) }
+ \item{shiftdir}{ indicates in which direction the axis should be
+ shifted ("left","right" or "none"). This direction is w.r.t. vector \code{g}}
+ \item{shiftfactor}{ scalar by which the shift vector is stretched (or
+ shrunken). By default, the length of the shift vector is stretched
+ by 5 percent (\code{shiftfactor} = 1.05)}
+}
+\details{
+ This program calibrates variable vectors in biplots and scatterplots, by drawing tick marks along
+ a given the vector and labelling the tick marks with specified values. The optimal calibration is
+ found by (generalized) least squares. Non-optimal calibrations are possible by specifying a
+ calibration factor (alpha).
+}
+\value{
+ Returns a list with calibration results
+ \item{useralpha}{calibration factor specified by the user}
+ \item{optalpha}{optimal calibration factor}
+ \item{lengthoneunit}{length in the plot of one unit in the scale of the calibrated variable}
+ \item{gof}{goodness of fit (as in regression)}
+ \item{gos}{goodness of scale}
+ \item{M}{coordinates of the tick markers}
+ \item{ang}{angle in degrees of the biplot axis with the positive
+ x-axis}
+ \item{shiftvec}{the supplied or computed shift vector}
+ \item{yt}{fitted values for the variable according to the calibration}
+ \item{e}{errors according to the calibration}
+ \item{Fpr}{coordinates of the projections of the row markers onto the calibrated axis}
+ \item{Mn}{coordinates of the tick marker end points}
+}
+\references{
+Gower, J.C. and Hand, D.J., (1996) Biplots. Chapman & Hall, London
+
+Graffelman, J. and van Eeuwijk, F.A. (2005) Calibration of multivariate scatter plots for
+exploratory analysis of relations within and between sets of variables in genomic research
+Biometrical Journal, 47(6) pp. 863-879.
+
+Graffelman, J. (2006) A guide to biplot calibration.
+}
+\author{ Jan Graffelman \email{jan.graffelman at upc.edu} }
+\seealso{\code{\link{biplot}}}
+\examples{
+x <- rnorm(20,1)
+y <- rnorm(20,1)
+x <- x - mean(x)
+y <- y - mean(y)
+z <- x + y
+b <- c(1,1)
+plot(x,y,asp=1,pch=19)
+tm<-seq(-2,2,by=0.5)
+Calibrate.z <- calibrate(b,z,tm,cbind(x,y),axislab="Z",graphics=TRUE)
+}
+\keyword{multivariate}
diff --git a/man/calves.Rd b/man/calves.Rd
new file mode 100644
index 0000000..4395d10
--- /dev/null
+++ b/man/calves.Rd
@@ -0,0 +1,16 @@
+\name{calves}
+\docType{data}
+\alias{calves}
+\title{Delivery of Dutch Calves}
+\description{
+ This data set gives a cross classification of 7275
+ calves born in the late nineties according to type
+ of production and type of delivery.
+}
+\usage{data(calves)}
+\format{A data frame containing a contingency table of 7275 observations.}
+\source{Holland Genetics. \url{http://www.hg.nl} }
+\references{
+ Graffelman, J. (2005) \emph{A guide to scatterplot and biplot calibration}.
+}
+\keyword{datasets}
\ No newline at end of file
diff --git a/man/canocor.Rd b/man/canocor.Rd
new file mode 100644
index 0000000..6581f3b
--- /dev/null
+++ b/man/canocor.Rd
@@ -0,0 +1,56 @@
+\name{canocor}
+\alias{canocor}
+\title{Canonical correlation analysis}
+\description{
+ \code{canocor} performs canonical correlation analysis on the
+ basis of the standardized variables and stores extensive output
+ in a list object.
+}
+\usage{
+canocor(X, Y)
+}
+\arguments{
+ \item{X}{ a matrix containing the X variables }
+ \item{Y}{ a matrix containing the Y variables }
+}
+\details{
+ \code{canocor} computes the solution by a singular value
+ decomposition of the transformed between set correlation matrix.
+}
+\value{
+Returns a list with the following results
+\item{ccor }{ the canonical correlations }
+\item{A }{ canonical weights of the x variables }
+\item{B }{ canonical weights of the y variables }
+\item{U }{ canonical x variates }
+\item{V }{ canonical y variates }
+\item{Fs }{ biplot markers for x variables (standard coordinates) }
+\item{Gs }{ biplot markers for y variables (standard coordinates) }
+\item{Fp }{ biplot markers for x variables (principal coordinates) }
+\item{Gp }{ biplot markers for y variables (principal coordinates) }
+\item{fitRxy }{ goodness of fit of the between-set correlation matrix }
+\item{fitXs }{ adequacy coefficients of x variables }
+\item{fitXp }{ redundancy coefficients of x variables }
+\item{fitYs }{ adequacy coefficients of y variables }
+\item{fitYp }{ redundancy coefficients of y variables}
+}
+\references{
+ Hotelling, H. (1935) The most predictable criterion. Journal of Educational
+ Psychology (26) pp. 139-142.
+
+ Hotelling, H. (1936) Relations between two sets of variates. Biometrika
+ (28) pp. 321-377.
+
+ Johnson, R. A. and Wichern, D. W. (2002) Applied Multivariate Statistical Analysis.
+ New Jersey: Prentice Hall.
+}
+\author{ Jan Graffelman \email{jan.graffelman at upc.edu} }
+\seealso{\code{\link{cancor}}}
+\examples{
+set.seed(123)
+X <- matrix(runif(75),ncol=3)
+Y <- matrix(runif(75),ncol=3)
+cca.results <- canocor(X,Y)
+}
+\keyword{multivariate}
+
diff --git a/man/circle.Rd b/man/circle.Rd
new file mode 100644
index 0000000..717df85
--- /dev/null
+++ b/man/circle.Rd
@@ -0,0 +1,24 @@
+\name{circle}
+\alias{circle}
+\title{Draw a circle}
+\description{
+ \code{circle} draws a circle in an existing plot.
+}
+\usage{
+circle(radius,origin)
+}
+\arguments{
+ \item{radius}{ the radius of the circle }
+ \item{origin}{ the origin of the circle }
+}
+\value{
+ NULL
+}
+\author{ Jan Graffelman \email{jan.graffelman at upc.edu} }
+\examples{
+ set.seed(123)
+ X <- matrix(rnorm(20),ncol=2)
+ plot(X[,1],X[,2])
+ circle(1,c(0,0))
+}
+\keyword{aplot}
diff --git a/man/dlines.Rd b/man/dlines.Rd
new file mode 100644
index 0000000..be544a3
--- /dev/null
+++ b/man/dlines.Rd
@@ -0,0 +1,31 @@
+\name{dlines}
+\alias{dlines}
+\title{Connect two sets of points by lines}
+\description{
+ \code{dlines} connects two sets of points by lines in
+ a rowwise manner.
+}
+\usage{
+dlines(SetA, SetB, lin = "dotted")
+}
+%- maybe also 'usage' for other objects documented here.
+\arguments{
+ \item{SetA}{matrix with the first set of points}
+ \item{SetB}{matrix with teh second set of points}
+ \item{lin}{linestyle for the connecting lines}
+}
+\value{
+ NULL
+}
+\author{ Jan Graffelman (jan.graffelman at upc.edu) }
+\seealso{\code{\link{lines}}}
+\examples{
+X <- matrix(runif(20),ncol=2)
+Y <- matrix(runif(20),ncol=2)
+plot(rbind(X,Y))
+text(X[,1],X[,2],paste("X",1:10,sep=""))
+text(Y[,1],Y[,2],paste("Y",1:10,sep=""))
+dlines(X,Y)
+}
+\keyword{aplot}
+
diff --git a/man/goblets.Rd b/man/goblets.Rd
new file mode 100644
index 0000000..b5f5d1c
--- /dev/null
+++ b/man/goblets.Rd
@@ -0,0 +1,16 @@
+\name{goblets}
+\docType{data}
+\alias{goblets}
+\title{Size measurements of archeological goblets}
+\description{
+ This data set gives 6 different size measurements of 25 goblets
+}
+\usage{data(goblets)}
+\format{A data frame containing 25 observations.}
+\source{Manly, 1989}
+\references{
+ Manly, B. F. J. (1989) \emph{Multivariate statistical methods: a primer}.
+ London: Chapman and Hall, London
+}
+\keyword{datasets}
+
\ No newline at end of file
diff --git a/man/heads.Rd b/man/heads.Rd
new file mode 100644
index 0000000..ff785ab
--- /dev/null
+++ b/man/heads.Rd
@@ -0,0 +1,23 @@
+\name{heads}
+\docType{data}
+\alias{heads}
+\title{Dimensions of heads of first and second sons for 25 families}
+\description{
+ Variables X1 and X2 are the head length and head breadth of the first son and
+ Y1 and Y2 are the same variables for the second son.
+}
+\usage{data(heads)}
+\format{A data frame containing 25 observations.}
+\source{Mardia, 1979, p. 121}
+\references{
+ Frets, G. P. (1921) \emph{Heredity of head form in man},
+ Genetica 3, pp. 193-384.
+
+ Mardia, K. V. and Kent, J. T. and Bibby, J. M. (1979) \emph{Multivariate Analysis}.
+ Academic Press London.
+
+ Anderson, T. W. (1984) \emph{An Introduction to Multivariate Statistical Analysis}.
+ New York: John Wiley, Second edition.
+}
+\keyword{datasets}
+
\ No newline at end of file
diff --git a/man/linnerud.Rd b/man/linnerud.Rd
new file mode 100644
index 0000000..7e70f3f
--- /dev/null
+++ b/man/linnerud.Rd
@@ -0,0 +1,18 @@
+\name{linnerud}
+\docType{data}
+\alias{linnerud}
+\title{Linnerud's exercise and body measurements}
+\description{
+ The data set consist of 3 exercise variables (Tractions a la barre
+ fixe, Flexions, Sauts) and 3 body measurements (Poids, Tour de talle,
+ Pouls) of 20 individuals.
+}
+\usage{data(linnerud)}
+\format{A data frame containing 20 observations.}
+\source{Tenenhaus, 1998, table 1, page 15}
+\references{
+ Tenenhaus, M. (1998) \emph{La Regression PLS}.
+ Paris: Editions Technip.
+}
+\keyword{datasets}
+
\ No newline at end of file
diff --git a/man/ones.Rd b/man/ones.Rd
new file mode 100644
index 0000000..59080bc
--- /dev/null
+++ b/man/ones.Rd
@@ -0,0 +1,27 @@
+\name{ones}
+\alias{ones}
+\title{Generates a matrix of ones}
+\description{
+ \code{ones} generates a matrix of ones.
+}
+\usage{
+ones(n, p = n)
+}
+\arguments{
+ \item{n}{ number of rows }
+ \item{p}{ number of columns }
+}
+\details{
+ if only n is specified, the resulting matrix will be square.
+}
+\value{
+ a matrix filled with ones.
+}
+\author{ Jan Graffelman (jan.graffelman at upc.edu) }
+\seealso{\code{\link{matrix}}}
+\examples{
+Id <- ones(3)
+print(Id)
+}
+\keyword{multivariate}
+
diff --git a/man/origin.Rd b/man/origin.Rd
new file mode 100644
index 0000000..1590798
--- /dev/null
+++ b/man/origin.Rd
@@ -0,0 +1,19 @@
+\name{origin}
+\alias{origin}
+\title{Origin}
+\description{
+ Draws coordinate axes in a plot.
+}
+\usage{origin(m=c(0,0), ...)}
+\arguments{
+ \item{m}{ the coordinates of the means (2 x 1). }
+ \item{\dots}{ other arguments passed on to the \code{lines} function}
+}
+\author{ Jan Graffelman (jan.graffelman at upc.edu) }
+\seealso{\code{\link{lines}}}
+\examples{
+ X <- matrix(runif(40),ncol=2)
+ plot(X[,1],X[,2])
+ origin(m=c(mean(X[,1]),mean(X[,2])))
+}
+\keyword{multivariate}
diff --git a/man/rad2degree.Rd b/man/rad2degree.Rd
new file mode 100644
index 0000000..b0f26e8
--- /dev/null
+++ b/man/rad2degree.Rd
@@ -0,0 +1,25 @@
+\name{rad2degree}
+\alias{rad2degree}
+\title{
+ Convert radians to degrees.
+}
+\description{
+\code{rad2degree converts radians to degrees.}
+}
+\usage{
+rad2degree(x)
+}
+\arguments{
+ \item{x}{an angle in radians}
+}
+\value{
+ the angle with the positive x-axis in degrees.
+}
+\author{ Jan Graffelman (jan.graffelman at upc.edu) }
+\examples{
+ x <- pi/2
+ a <- rad2degree(x)
+ cat("angle is",a,"degrees\n")
+}
+\keyword{arith}
+
diff --git a/man/rda.Rd b/man/rda.Rd
new file mode 100644
index 0000000..180aa1f
--- /dev/null
+++ b/man/rda.Rd
@@ -0,0 +1,54 @@
+\name{rda}
+\alias{rda}
+\title{Redundancy analysis}
+\description{
+ \code{rda} performs redundancy analysis and stores extensive output
+ in a list object.
+}
+\usage{
+rda(X, Y, scaling = 1)
+}
+\arguments{
+ \item{X}{a matrix of x variables}
+ \item{Y}{a matrix of y variables}
+ \item{scaling}{scaling used for x and y variables. 0: x and y only centered. 1: x and y standardized}
+}
+\details{
+ Results are computed by doing a principal component analyis of the
+ fitted values of the regression of y on x.
+
+ Plotting the first two columns of Gxs and Gyp, or of Gxp and Gys provides a biplots of
+ the matrix of regression coefficients.
+
+ Plotting the first two columns of Fs and Gp or of Fp and Gs provides a biplot of the
+ matrix of fitted values.
+}
+\value{
+Returns a list with the following results
+\item{Yh}{ fitted values of the regression of y on x }
+\item{B}{ regression coefficients of the regresson of y on x }
+\item{decom}{ variance decomposition/goodness of fit of the fitted values AND of
+the regression coefficients }
+\item{Fs}{ biplot markers of the rows of Yh (standard coordinates) }
+\item{Fp}{ biplot markers of the rows of Yh (principal coordinates) }
+\item{Gys}{ biplot markers for the y variables (standard coordinates) }
+\item{Gyp}{ biplot markers for the y variables (principal coordinates) }
+\item{Gxs}{ biplot markers for the x variables (standard coordinates) }
+\item{Gxp}{ biplot markers for the x variables (principal coordinates) }
+}
+\references{
+ Van den Wollenberg, A.L. (1977) Redundancy Analysis, an alternative for canonical
+ correlation analysis. Psychometrika 42(2): pp. 207-219.
+
+ Ter Braak, C. J. F. and Looman, C. W. N. (1994) Biplots in Reduced-Rank Regression.
+ Biometrical Journal 36(8): pp. 983-1003.
+}
+\author{ Jan Graffelman (jan.graffelman at upc.edu) }
+\seealso{\code{\link{princomp}},\code{\link{canocor}},\code{\link{biplot}}}
+\examples{
+X <- matrix(rnorm(75),ncol=3)
+Y <- matrix(rnorm(75),ncol=3)
+rda.results <- rda(X,Y)
+}
+\keyword{multivariate}
+
diff --git a/man/shiftvector.Rd b/man/shiftvector.Rd
new file mode 100644
index 0000000..25268e7
--- /dev/null
+++ b/man/shiftvector.Rd
@@ -0,0 +1,60 @@
+\name{shiftvector}
+\alias{shiftvector}
+\title{
+ Compute a shift vector for a calibrated axis.
+}
+\description{
+ \code{shiftvector} computes two shift vectors perpendicular to the
+ supplied biplot or scatterplot axis \code{g}. The vector norm is
+ computed from the two most extreme data points.
+}
+\usage{
+shiftvector(g, X, x = c(1, 0), verbose = FALSE)
+}
+\arguments{
+ \item{g}{a biplot or scatterplot axis}
+ \item{X}{a n by 2 matrix of scatterplot or biplot coordinates}
+ \item{x}{reference axis, (1,0) by default}
+ \item{verbose}{print information or not}
+}
+\details{
+ \code{shiftvector} locates the tow most extreme datapoints in the
+ direction perpendicular to axis \code{g}.
+}
+\value{
+ \item{dr}{the right (w.r.t. the direction of \code{g}) shift vector}
+ \item{dl}{the left (w.r.t. the direction of \code{g}) shift vector}
+}
+\references{
+Graffelman, J. and van Eeuwijk, F.A. (2005) Calibration of multivariate scatter plots for
+exploratory analysis of relations within and between sets of variables in genomic research
+Biometrical Journal, 47(6) pp. 863-879.
+
+Graffelman, J. (2006) A guide to biplot calibration.
+}
+\author{ Jan Graffelman (jan.graffelman at upc.edu) }
+\seealso{\code{\link{calibrate}}}
+
+\examples{
+X <- matrix(rnorm(100),ncol=2)
+Xs <- scale(X)
+
+g <- c(1,1)
+
+plot(Xs[,1],Xs[,2],asp=1,pch=19)
+textxy(Xs[,1],Xs[,2],1:nrow(X))
+
+arrows(0,0,g[1],g[2])
+text(g[1],g[2],"g",pos=1)
+
+out <- shiftvector(g,X,verbose=TRUE)
+dr <- out$dr
+dl <- out$dl
+
+arrows(0,0,dl[1],dl[2])
+text(dl[1],dl[2],"dl",pos=1)
+
+arrows(0,0,dr[1],dr[2])
+text(dr[1],dr[2],"dr",pos=1)
+}
+\keyword{multivariate}
diff --git a/man/storks.Rd b/man/storks.Rd
new file mode 100644
index 0000000..45c0419
--- /dev/null
+++ b/man/storks.Rd
@@ -0,0 +1,17 @@
+\name{storks}
+\docType{data}
+\alias{storks}
+\title{Frequencies of nesting storks in Denmark}
+\description{
+ Danish data from 1953-1977 giving the frequency of nesting storks,
+ the human birth rate and the per capita electricity consumption.
+}
+\usage{data(storks)}
+\format{A data frame containing 25 observations.}
+\source{Gabriel and Odoroff, Table 1.}
+\references{
+ Gabriel, K. R. and Odoroff, C. L. (1990) Biplots in biomedical research.
+ \emph{Statistics in Medicine} 9(5): pp. 469-485.
+}
+\keyword{datasets}
+
\ No newline at end of file
diff --git a/man/textxy.Rd b/man/textxy.Rd
new file mode 100644
index 0000000..182faa2
--- /dev/null
+++ b/man/textxy.Rd
@@ -0,0 +1,39 @@
+\name{textxy}
+\alias{textxy}
+\title{ Nice placement of labels in a plot}
+\description{
+ Function \code{textxy} calls function \code{text} in order to add text
+ to points in a graph. \code{textxy} chooses a different position
+ for the text depending on the quadrant. This tends to
+ produces better readable plots, with labels fanning away from the origin.
+}
+\usage{
+textxy(X, Y, labs, m = c(0, 0), cex = 0.5, offset = 0.8, ...)
+}
+\arguments{
+ \item{X}{x coordinates of a set of points}
+ \item{Y}{y coordinates of a set of points}
+ \item{labs}{labels to be placed next to the points}
+ \item{m}{coordinates of the origin of the plot (default (0,0))}
+ \item{cex}{character expansion factor}
+ \item{offset}{controls the distance between the label and the point. A
+ value of 0 will plot labels on top of the point. Larger values give
+ larger separation between point and label. The default value is 0.8}
+ \item{\dots}{additiona arguments for function \code{text}.}
+}
+\value{
+ NULL
+}
+\references{ Graffelman, J. (2006) A guide to biplot calibration. }
+\author{ Jan Graffelman (jan.graffelman at upc.edu) }
+\seealso{ \code{\link{text}} }
+\examples{
+x <- rnorm(50)
+y <- rnorm(50)
+plot(x,y,asp=1)
+textxy(x,y,1:50,m=c(mean(x),mean(y)))
+}
+\keyword{aplot}
+\keyword{misc}
+
+
diff --git a/vignettes/CalibrationGuide.Rnw b/vignettes/CalibrationGuide.Rnw
new file mode 100644
index 0000000..2d6cc85
--- /dev/null
+++ b/vignettes/CalibrationGuide.Rnw
@@ -0,0 +1,777 @@
+\documentclass[a4paper]{article}
+
+% \VignetteIndexEntry{A R package for calibration of biplot and scatterplot axis.}
+% \VignetteDepends{graphics,stats}
+% \VignetteKeyword{multivariate}
+
+% Documentation for the Calibrate package
+
+\usepackage[english]{babel}
+\usepackage{Sweave}
+%\usepackage{Rd}
+\usepackage{url}
+\usepackage{wrapfig}
+\usepackage{hyperref}
+\usepackage[utf8]{inputenc}
+
+\setlength{\parindent}{0cm}
+
+\newcommand{\bF}{\ensuremath{\mathbf F}}
+\newcommand{\bG}{\ensuremath{\mathbf G}}
+
+\begin{document}
+
+\begin{center}
+\sf
+{\sf \bf \Large A Guide to Scatterplot and Biplot Calibration}\\
+\vspace{4mm}
+{\sf \normalsize {\tt version 1.7.2}}\\%\normalsize
+\vspace{4mm}
+{\bf \large Jan Graffelman}\\
+\vspace{4mm} \rm \large
+Department of Statistics and Operations Research\\
+Universitat Polit\`ecnica de Catalunya\\
+Avinguda Diagonal 647, 08028 Barcelona, Spain.\\
+{\it email:} jan.graffelman at upc.edu\\
+\vspace{4mm}
+{\sc September 2012}
+\end{center}
+
+\section{Introduction}
+
+This guide gives detailed instructions on how to calibrate axes in scatterplots
+and biplots obtained in the statistical environment R~\cite{RRR} by using the package
+{\tt calibrate}. By calibration we
+refer to the procedure of drawing a (linear) scale along an axis in a plot with
+tick marks and numeric labels. In an ordinary scatter plot of two variables $x$ and $y$
+two calibrated perpendicular scales are typically automatically produced by the
+routine used for plotting the two variables. However, scatter plots can be
+extended with additional variables that are represented on oblique additional
+axes. The software described in this guide can be used to create calibrated
+scales on these oblique additional axes. Moreover, in a multivariate setting with more
+than two variables, raw data matrices, correlation matrices, contingency tables,
+regression coefficients, etc. are often represented graphically by means of biplots~\cite{Gabriel}. Biplots
+also contain oblique axes representing variables. The described software can also be
+used to construct scales on biplot axes.\\
+
+The outline of this guide is as follows. In Section~\ref{sec:install} we indicate how the
+R package {\tt calibrate} can be installed. Section~\ref{sec:scatter} describes in detail
+how to calibrate additional axes in scatter plots. Section~\ref{sec:biplot} treats the
+calibration of biplot axes. Several subsections follow with detailed
+instructions of how to calibrate biplot axis in principal component analysis
+(PCA, Section~\ref{sec:pca}), correspondence analysis (CA. Section~\ref{sec:ca}),
+canonical correlation analysis (CCA, Section~\ref{sec:cca}) and redundancy analysis (RDA, Section~\ref{sec:rda}).
+The online documentation of the main routine for calibration {\tt calibrate} is referenced in
+Section~\ref{sec:online}.\\
+
+This guide does not provide the theory for the construction of scales on scatterplot
+and biplot axes. For a theoretical account of biplot calibration, we refer to Graffelman
+\& van Eeuwijk~(2005) and to Gower and Hand~(1996). If you appreciate
+this software then please cite the following paper in your work:\\
+
+Graffelman, J. \& van Eeuwijk, F.A. (2005) Calibration of multivariate scatter plots for
+exploratory analysis of relations within and between sets of variables in genomic research
+{\it Biometrical Journal}, {\bf 47}(6) pp. 863-879. \href{http://dx.doi.org/10.1002/bimj.200510177}{(clic here to access the paper)}
+
+\section{Installation}
+\label{sec:install}
+
+<<fig=FALSE,echo=FALSE>>=
+prefig <- function(){
+data(goblets)
+X <- goblets
+plot(X[,1],X[,2],pch=19,cex=0.5,xlab=expression(X[1]),
+ylab=expression(X[2]),xlim=c(5,25),ylim=c(5,25),asp=1)
+m <- apply(cbind(X[,1],X[,2]),2,mean)
+textxy(X[,1],X[,2],1:25,m=m,cex=0.75)
+points(m[1],m[2],col="red",pch=19,cex=0.5)
+abline(h = m[2])
+abline(v = m[1])
+}
+options("SweaveHooks"=list(aap=prefig))
+options("width"=60)
+@
+
+Packages in R can be installed inside the program with the option "Packages"
+in the main menu and then choosing "Install package" and picking the package
+"calibrate". Typing:
+
+<<noot>>=
+library(calibrate)
+@
+
+will, among others, make the function {\tt calibrate, canocor} and {\tt rda} available. Several
+small data sets, also the ones used in this document, are included in the package ({\tt calves, goblets,
+heads, linnerud} and {\tt storks}).
+
+\section{Calibration of Scatterplot axes}
+\label{sec:scatter}
+
+We consider a archaeological data set concerning 6 size measurements ($X_1, \ldots, X_6$) on 25
+goblets. This data was published by Manly~(1989). The data can be loaded with
+the {\tt data} instruction.
+<<fig=FALSE,echo=TRUE>>=
+data(goblets)
+X <- goblets
+@
+
+\subsection*{Oblique additional axes in a scatterplot}
+
+We construct a scatterplot of $X_1$ versus $X_2$ and center a set of coordinate
+axes on the point $(\overline{x}_1,\overline{x}_2)$ with the function {\tt origin}.
+
+\setkeys{Gin}{width=\textwidth}
+
+<<fig=TRUE,echo=TRUE>>=
+plot(X[,1],X[,2],pch=19,cex=0.5,xlab=expression(X[1]),ylab=expression(X[2]),
+ xlim=c(5,25),ylim=c(5,25),asp=1)
+m <- apply(X[,1:2],2,mean)
+textxy(X[,1],X[,2],1:25,m=m,cex=0.75)
+origin(m)
+@
+
+Next, we perform the regression of $X_5$ onto $X_1$ and $X_2$ (all variables being centered) in order to
+obtain an additional axis for $X_5$. We represent $X_5$ in the plot as a simple
+arrow whose coordinates are given by the regression coefficients:
+
+<<fig=TRUE,echo=TRUE,aap=TRUE>>=
+Xc <- scale(X,center=TRUE,scale=FALSE)
+b <- solve(t(Xc[,1:2])%*%Xc[,1:2])%*%t(Xc[,1:2])%*%Xc[,5]
+print(b)
+bscaled <- 20*b
+arrows(m[1],m[2],m[1]+bscaled[1],m[2]+bscaled[2],col="blue",length=0.1)
+arrows(m[1],m[2],m[1]-bscaled[1],m[2]-bscaled[2],length=0,lty="dashed",col="blue")
+@
+
+A direction that is optimal in the least squares sense for $X_5$ is given by the vector of regression
+coefficients~\cite{Graffel13}. To make this direction more visible, we multiplied it by a constant (20).
+It is clear that the direction of increase for $X_5$ runs approximately North-East across the scatterplot.
+We now proceed to calibrate this direction with a scale for $X_5$. In order to choose sensible values for
+the scale of $X_5$, we first inspect the range of variation of $X_5$, and then choose a set of values we want
+to mark off on the scale ({\tt tm}) and also compute the deviations of these values from the mean
+({\tt tmc}). We specify a tick length of 0.3 ({\tt tl=0.3}). Depending on the data, some values of {\tt tl}
+typically have to be tried to see how to obtain a nice scale.
+
+<<fig=TRUE,echo=TRUE,aap=TRUE>>=
+print(range(X[,5]))
+yc <- scale(X[,5],scale=FALSE)
+tm <- seq(2,10,by=1)
+tmc <- tm - mean(X[,5])
+Calibrate.X5<-calibrate(b,yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.3,axislab="X_5",
+ labpos=4,cex.axislab=1)
+@
+
+The numerical output from routine {\tt calibrate} shows that one unit along the axis for $X_5$ occupies 2.47
+units in the plotting frame. The axis for $X_5$ makes an angle of 17.65 degrees with the positive x-axis.
+The calibration factor is 6.12. Multiplying the vector of regressions coefficients by this
+factor yields a vector that represents a unit change in the scale of $X_5$. E.g. for this data we have that
+the vector $6.12 \cdot (0.385, 0.123) = (2.358, 0.751)$ represents a unit change. This vector has
+norm $\sqrt{2.358^2 + 0.751^2} = 2.47$. Other calibration factors may be specified by using parameter
+{\tt alpha}. If {\tt alpha} is left unspecified the optimal value computed by least
+squares will be used. The goodness-of-fit of $X_5$ is 0.513. This means that 51.3\% of the variance
+of $X_5$ can be explained by a regression onto $X_1$ and $X_2$ ($R^2 = 0.513$). The goodness-of-scale has
+the same value. The goodness-of-scale is only relevant if we modify parameter {\tt alpha}. {\tt Calibrate.X5}
+is a list object containing all calibration results (calibration factor, fitted values according to the
+scale used, tick marker positions, etc.).
+
+\subsubsection*{Shifting a calibrated axis}
+
+Using many calibrated axes in a plot, all passing through the origin, leads to dense plots that become
+unreadable. It is therefore a good idea to shift calibrated axes towards the margins of the plot. This keeps
+the central cloud of data points clear and relegates all information on scales to the margins of the graph.
+There are two natural positions for a shifted axis: just above the largest data point in a direction perpendicular
+to the axis being calibrated, or just below the smallest data point in the perpendicular direction. The arguments
+{\tt shiftdir, shiftfactor} and {\tt shiftvec} can be used to control the shifting of a calibrated axis. {\tt
+shiftvec} allows the user to specify the shift vector manually. This is normally not needed, and good positions
+for an axis can be found by using only {\tt shiftdir} and {\tt shiftfactor}. Argument {\tt shiftdir} can be set
+to {\tt 'right'} or {\tt 'left'} and indicates in which direction the axis is to be shifted, with respect to the direction of
+increase of the calibrated axis. Setting {\tt shiftdir} shifts the axis automatically just above or below
+the most outlying data point in the direction perpendicular to the vector being calibrated. In order to
+move the calibrated axis farther out or to pull it more in, {\tt shiftfactor} can be used.
+Argument {\tt shiftfactor} stretches or shrinks the shift vector for the axis. A {\tt shiftfactor} larger
+than 1 moves the axis outwards, and a {\tt shiftfactor} smaller than 1 pulls the axis towards the origin of
+the plot. If set to 1 exactly, the shifted axis will cut through the most outlying data point.
+The default {\tt shiftfactor} is 1.05. We redo
+the previous plot, shifting the calibrated axis below the cloud of points, which is to the right w.r.t. the
+direction of increase of the variable.
+
+<<fig=TRUE,echo=TRUE,aap=TRUE>>=
+yc <- scale(X[,5],scale=FALSE)
+tm <- seq(2,10,by=1)
+tmc <- tm - mean(X[,5])
+Calibrate.X5<-calibrate(b,yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.3,axislab="X_5",labpos=4,
+ cex.axislab=1,shiftdir="right")
+@
+
+The shift of the axis does not affect the interpretation of the plot, because the projections of the points
+onto the axis remain the same.
+
+\subsubsection*{Second vertical axis in a scatterplot}
+
+The oblique direction in the previous section is the preferred direction for $X_5$, as this direction
+is optimal in the least squares sense. However, if desired, additional variables can also be represented
+as a second vertical axis on the right of the plot, or as a second horizontal axis on the top of the
+plot. We now proceed to construct a second vertical axis on the right hand of the scatter plot for
+$X_5$. This can be done by setting the vector to be calibrated (first argument of routine {\tt calibrate})
+to the (0,1) vector. By specifying a shiftvector explicitly ({\tt shiftvec}), the axis can be shifted. For
+this data, setting {\tt shiftvec} to {\tt c(par('usr')[2]-mean(X[,1]),0)} and {\tt shiftfactor = 1},
+makes the axis coincide with the right vertical borderline of the graph.
+
+<<fig=TRUE,echo=TRUE,aap=TRUE>>=
+opar <- par('xpd'=TRUE)
+tm <- seq(3,8,by=1)
+tmc <- (tm - mean(X[,5]))
+Calibrate.rightmargin.X5 <- calibrate(c(0,1),yc,tmc,Xc[,1:2],tmlab=tm,m=m,
+ axislab="X_5",tl=0.5,
+ shiftvec=c(par('usr')[2]-mean(X[,1]),0),
+ shiftfactor=1,where=2,
+ laboffset=c(1.5,1.5),cex.axislab=1)
+par(opar)
+@
+
+The second vertical axis has calibration factor 3.46, and a goodness of fit of 0.34. The fit of the
+variable is worse in comparison with the previous oblique direction given by the regression coefficients. Note
+that graphical clipping in temporarily turned off ({\tt par('xpd'=TRUE)}) to allow the calibration
+routine to draw ticks and labels outside the figure region, and that the range of the tick marks was
+shortened in order not surpass the figure region.
+
+\subsubsection*{Subscales and double calibrations}
+
+Scales with tick marks can be refined by drawing subscales with smaller tick marks.
+E.g. larger labelled
+tickmarks can be used to represent multiples of 10, and small unlabelled tick marks can be used to
+represent units. The subscale allows a more precise recovery of the data values. This can simply be
+achieved by calling the calibration routine twice, once with a coarse sequence and once with a finer
+sequence. For the second call one can specify {\tt verb=FALSE} in order to suppress the numerical output
+of the routine, and {\tt lm=FALSE} to supress the tick mark labels under the smaller ticks. The tickmarks
+for the finer scale are made smaller by modifying the tick length (e.g. {\tt tl=0.1}). Depending on the data,
+some trial and error with different values for {\tt tl} may be necessary before nice scales are obtained. This
+may be automatized in the future. Finally, reading off the (approximate) data values can further be enhanced
+by drawing perpendiculars from the points to the calibrated axis by setting {\tt dp=TRUE}.
+
+<<fig=TRUE,echo=TRUE,aap=TRUE>>=
+tm <- seq(2,10,by=1)
+tmc <- (tm - mean(X[,5]))
+Calibrate.X5 <- calibrate(b,yc,tmc,Xc[,1:2],tmlab=tm,m=m,axislab="X_5",tl=0.5,
+ dp=TRUE,labpos=4)
+tm <- seq(2,10,by=0.1)
+tmc <- (tm - mean(X[,5]))
+Calibrate.X5 <- calibrate(b,yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.25,verb=FALSE,
+ lm=FALSE)
+@
+
+A {\it double calibration} can be created by drawing two scales, one on each side of the axis. Double
+calibrations can be useful. For instance, one scale can be used for recovery of the original
+data values of the variable, whereas the second scale can be used for recovery of standardized values or
+of correlations with other variables. Double calibrations can also be used to
+graphically verify if two different calibration procedures give the same result or not.
+
+\subsubsection*{Recalibrating the original scatterplot axes}
+
+By calibrating the (0,1) and (1,0) vectors the original axes of the scatter plot can be
+redesigned. We illustrate the recalibration of the original axes by creating a second scale on
+the other side of the axes, a refined scale for $X_1$, and a scale for the standardized data for
+$X_2$. For the latter calibration one unit equals one standard deviation.
+
+<<fig=TRUE,echo=TRUE,aap=TRUE>>=
+opar <- par('xpd'=TRUE)
+tm <- seq(5,25,by=5)
+tmc <- (tm - mean(X[,1]))
+yc <- scale(X[,1],scale=FALSE)
+Calibrate.X1 <- calibrate(c(1,0),yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.5,
+ axislab="X_1",cex.axislab=1,showlabel=FALSE,
+shiftvec=c(0,-(m[2]-par("usr")[3])),shiftfactor=1,reverse=TRUE)
+tm <- seq(5,25,by=1); tmc <- (tm - mean(X[,1]))
+Calibrate.X1 <- calibrate(c(1,0),yc,tmc,Xc[,1:2],tmlab=tm,m=m,tl=0.25,
+ axislab="X_1",cex.axislab=1,showlabel=FALSE,
+shiftvec=c(0,-(m[2]-par("usr")[3])),shiftfactor=1,reverse=TRUE,
+ verb=FALSE,lm=FALSE)
+yc <- scale(X[,2],scale=TRUE)
+tm <- seq(-3,1,by=1)
+Calibrate.X2 <- calibrate(c(0,1),yc,tm,Xc[,1:2],tmlab=tm,m=m,tl=0.6,
+ axislab="X_2",cex.axislab=1,showlabel=FALSE,
+shiftvec=c(-(mean(X[,1])-par('usr')[1]),0),shiftfactor=1,verb=TRUE,lm=TRUE)
+tm <- seq(-3,1.5,by=0.1)
+Calibrate.X2 <- calibrate(c(0,1),yc,tm,Xc[,1:2],tmlab=tm,m=m,tl=0.3,
+ axislab="X_2",cex.axislab=1,showlabel=FALSE,
+shiftvec=c(-(mean(X[,1])-par('usr')[1]),0),shiftfactor=1,verb=FALSE,lm=FALSE)
+par(opar)
+@
+
+\section{Calibration of Biplot axes}
+\label{sec:biplot}
+
+In this section we give detailed instructions on how to calibrate biplot axes. We will consider biplots
+of raw data matrices and correlation matrices obtained by PCA, biplots of profiles obtained in CA,
+biplots of data matrices and correlation matrices (in particular the between-set correlation matrix)
+in CCA and biplots of fitted values and regression coefficients obtained by RDA. In principle, calibration
+of biplot axes has little additional complication in comparison with the calibration of additional
+axes in scatterplots explained above. The main issue is that, prior to calling the calibration routine,
+one needs to take care of the proper centring and standardisation of the tick marks.
+
+\subsection{Principal component analysis}
+\label{sec:pca}
+
+Principal component analysis can be performed by using routine {\tt princomp} from the {\tt stats }
+library. We use again Manly's goblets data to create a biplot of the data based on a
+PCA of the covariance matrix. We use {\tt princomp} to compute the scores for the rows and the columns
+of the data matrix. The first principal component is seen to be a size component, separating the
+smaller goblets on the right from the larger goblets on the left. The variable vectors are
+multiplied by a factor of 15 to facilitate interpretation. Next we
+calibrate the vector for $X_3$,
+using labelled tickmarks for multiples of 5 units, and shorter unlabelled tickmarks for the units. The
+goodness of fit of $X_3$ is very high (0.99), which means that $X_3$ is close to perfectly
+represented. {\tt Calibrate.X3} is a list object containing the numerical results of the
+calibration.
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+# PCA and Biplot construction
+pca.results <- princomp(X,cor=FALSE)
+Fp <- pca.results$scores
+Gs <- pca.results$loadings
+plot(Fp[,1],Fp[,2],pch=16,asp=1,xlab="PC 1",ylab="PC 2",cex=0.5)
+textxy(Fp[,1],Fp[,2],rownames(X),cex=0.75)
+arrows(0,0,15*Gs[,1],15*Gs[,2],length=0.1)
+textxy(15*Gs[,1],15*Gs[,2],colnames(X),cex=0.75)
+# Calibration of X_3
+ticklab <- seq(5,30,by=5)
+ticklabc <- ticklab-mean(X[,3])
+yc <- (X[,3]-mean(X[,3]))
+g <- Gs[3,1:2]
+Calibrate.X3 <- calibrate(g,yc,ticklabc,Fp[,1:2],ticklab,tl=0.5,
+ axislab="X3",cex.axislab=0.75,where=1,labpos=4)
+ticklab <- seq(5,30,by=1)
+ticklabc <- ticklab-mean(X[,3])
+Calibrate.X3.fine <- calibrate(g,yc,ticklabc,Fp[,1:2],ticklab,lm=FALSE,tl=0.25,
+ verb=FALSE,cex.axislab=0.75,where=1,labpos=4)
+@
+
+We do a PCA based on the correlation matrix, and proceed to construct a biplot of the correlation matrix. The
+correlations of $X_5$ with the other variables are computed, and the biplot axis for $X_5$ is calibrated with a
+correlation scale. Routine {\tt calibrate} is repeatedly called to create finer subscales.
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+# PCA and Biplot construction
+pca.results <- princomp(X,cor=TRUE)
+Fp <- pca.results$scores
+Ds <- diag(pca.results$sdev)
+Fs <- Fp%*%solve(Ds)
+Gs <- pca.results$loadings
+Gp <- Gs%*%Ds
+#plot(Fs[,1],Fs[,2],pch=16,asp=1,xlab="PC 1",ylab="PC 2",cex=0.5)
+#textxy(Fs[,1],Fs[,2],rownames(X))
+plot(Gp[,1],Gp[,2],pch=16,cex=0.5,xlim=c(-1,1),ylim=c(-1,1),asp=1,
+ xlab="1st principal axis",ylab="2nd principal axis")
+arrows(0,0,Gp[,1],Gp[,2],length=0.1)
+textxy(Gp[,1],Gp[,2],colnames(X),cex=0.75)
+ticklab <- c(seq(-1,-0.2,by=0.2),seq(0.2,1.0,by=0.2))
+R <- cor(X)
+y <- R[,5]
+g <- Gp[5,1:2]
+Calibrate.X5 <- calibrate(g,y,ticklab,Gp[,1:2],ticklab,lm=TRUE,tl=0.05,dp=TRUE,
+ labpos=2,cex.axislab=0.75,axislab="X_5")
+ticklab <- seq(-1,1,by=0.1)
+Calibrate.X5 <- calibrate(g,y,ticklab,Gp[,1:2],ticklab,lm=FALSE,tl=0.05,verb=FALSE)
+ticklab <- seq(-1,1,by=0.01)
+Calibrate.X5 <- calibrate(g,y,ticklab,Gp[,1:2],ticklab,lm=FALSE,tl=0.025,verb=FALSE)
+@
+
+The goodness of fit of the representation of the correlations of $X_5$ with the
+other variables is 0.98, the 6 correlations being close to perfectly represented.
+We compute the sample correlation matrix and compare the observed correlations of
+$X_5$ with those estimated from the calibrated biplot axis ({\tt yt}). Note that
+PCA also tries to approximate the correlation of a variable with itself, and that
+the arrow on representing $X_5$ falls short of the value 1 on its own calibrated scale.
+The refined subscale allows very precise graphical representation of the correlations
+as estimated by the biplot.
+
+\begin{center}
+<<fig=FALSE,echo=TRUE,aap=FALSE>>=
+print(R)
+print(cbind(R[,5],Calibrate.X5$yt))
+@
+\end{center}
+
+\subsection{Correspondence analysis}
+\label{sec:ca}
+
+We consider a contingency table of a sample of Dutch calves born in the late nineties,
+shown in Table~\ref{tab:calves2}. A total of 7257 calves were classified according
+to two categorical variables: the
+method of production (ET = Embryo Transfer, IVP = In Vitro Production, AI = Artificial
+Insemination) and the ease of delivery, scored on a scale from 1 (normal) to 6 (very heavy).
+The data in Table~\ref{tab:calves2} were provided by Holland Genetics.
+\begin{table}[htb]
+\centering
+\begin{tabular}{c|rrrr}
+ & \multicolumn{3}{c}{Type of calf}\\
+Ease of delivery & ET & IVP & AI\\
+\hline
+1 & 97 & 150 & 1686\\
+2 & 152 & 183 & 1339\\
+3 & 377 & 249 & 1209\\
+4 & 335 & 227 & 656\\
+5 & 42 & 136 & 277\\
+6 & 9 & 71 & 62\\
+\hline
+\end{tabular}
+\caption{Calves data from Holland Genetics.}
+\label{tab:calves2}
+\end{table}
+
+For this contingency table we obtain $\chi^2_{10} = 833.16$ with $p < 0.001$ and the null
+hypothesis of no association between ease of delivery and type of calf has to be rejected. However, what is
+the precise nature of this association? Correspondence analysis can be used to gain insight in the
+nature of this association. We use routine {\tt corresp} form the {\tt MASS} library~\cite{Venables}
+to perform correspondence analysis and to obtain the coordinates for a biplot of the row profiles.
+We compute the row profiles and then repeatedly call the calibration routine, each time with a
+different set of {\tt ticklabs}.
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+library(MASS)
+data(calves)
+ca.results <- corresp(calves,nf=2)
+Fs <- ca.results$rscore
+Gs <- ca.results$cscore
+Ds <- diag(ca.results$cor)
+Fp <- Fs%*%Ds
+Gp <- Gs%*%Ds
+plot(Gs[,1],Gs[,2],pch=16,asp=1,cex=0.5,xlab="1st principal axis",
+ ylab="2nd principal axis")
+textxy(Gs[,1],Gs[,2],colnames(calves),cex=0.75)
+points(Fp[,1],Fp[,2],pch=16,cex=0.5)
+textxy(Fp[,1],Fp[,2],rownames(calves),cex=0.75)
+origin()
+arrows(0,0,Gs[,1],Gs[,2])
+
+P <- as.matrix(calves/sum(calves))
+r <- apply(P,1,sum)
+k <- apply(P,2,sum)
+Dc <- diag(k)
+Dr <- diag(r)
+
+RP <- solve(Dr)%*%P
+print(RP)
+CRP <- RP - ones(nrow(RP), 1) %*% t(k)
+TCRP <- CRP%*%solve(Dc)
+
+y <- TCRP[,3]
+g <- Gs[3,1:2]
+
+ticklab <- c(0,seq(0,1,by=0.2))
+ticklabs <- (ticklab - k[3])/k[3]
+Calibrate.AI <- calibrate(g,y,ticklabs,Fp[,1:2],ticklab,lm=TRUE,tl=0.10,
+ weights=Dr,axislab="AI",labpos=4,dp=TRUE)
+ticklab <- c(0,seq(0,1,by=0.1))
+ticklabs <- (ticklab - k[3])/k[3]
+Calibrate.AI <- calibrate(g,y,ticklabs,Fp[,1:2],ticklab,lm=FALSE,tl=0.10,
+ weights=Dr,verb=FALSE)
+ticklab <- c(0,seq(0,1,by=0.01))
+ticklabs <- (ticklab - k[3])/k[3]
+Calibrate.AI <- calibrate(g,y,ticklabs,Fp[,1:2],ticklab,lm=FALSE,tl=0.05,
+ weights=Dr,verb=FALSE)
+@
+
+Because the calibration is done by weighted least squares, a diagonal matrix of weights ({\tt weights=Dr})
+is supplied as a parameter to the calibration routine
+Note that the calibrated axis for the row profiles with respect to AI has goodness of fit 1. This
+is due to the fact that the rank of the matrix of centred profiles is two, and that therefore
+all profiles can be perfectly represented in two dimensional space.
+
+\subsection{Canonical correlation analysis}
+\label{sec:cca}
+
+We consider a classical data set on the head sizes of the first and the second son of 25 families~\cite{Frets}.
+These data have been analysed by several authors~\cite{Anderson,Mardia,Graffel16}
+We first load the data and perform
+a canonical correlation analysis, using supplied function {\tt canocor} (a more fully
+fledged program for canonical correlation analysis in comparison with {\tt cancor} from the {\tt stats} package).
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+data(heads)
+X <- cbind(heads$X1,heads$X2)
+Y <- cbind(heads$Y1,heads$Y2)
+
+Rxy<- cor(X,Y)
+Ryx<- t(Rxy)
+Rxx<- cor(X)
+Ryy<- cor(Y)
+
+cca.results <-canocor(X,Y)
+
+plot(cca.results$Gs[,1],cca.results$Gs[,2],pch=16,asp=1,xlim=c(-1,1),ylim=c(-1,1),
+ xlab=expression(V[1]),ylab=expression(V[2]))
+arrows(0,0,cca.results$Fp[,1],cca.results$Fp[,2],length=0.1)
+arrows(0,0,cca.results$Gs[,1],cca.results$Gs[,2],length=0.1)
+
+textxy(cca.results$Fp[1,1],cca.results$Fp[1,2],expression(X[1]),cex=0.75)
+textxy(cca.results$Fp[2,1],cca.results$Fp[2,2],expression(X[2]),cex=0.75)
+
+textxy(cca.results$Gs[1,1],cca.results$Gs[1,2],expression(Y[1]),cex=0.75)
+textxy(cca.results$Gs[2,1],cca.results$Gs[2,2],expression(Y[2]),cex=0.75)
+
+circle(1)
+ticklab <- seq(-1,1,by=0.2)
+
+y <- Rxy[,2]
+g <- cca.results$Gs[2,1:2]
+
+Cal.Cor.Y2 <- calibrate(g,y,ticklab,cca.results$Fp[,1:2],ticklab,lm=TRUE,tl=0.05,
+ dp=TRUE,reverse=TRUE,weights=solve(Rxx),
+axislab="Y_2",cex.axislab=0.75,showlabel=FALSE)
+@
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+
+plot(cca.results$Gs[,1],cca.results$Gs[,2],pch=16,asp=1,xlim=c(-2,2),ylim=c(-2,2),
+ xlab=expression(V[1]),ylab=expression(V[2]))
+#arrows(0,0,cca.results$Fp[,1],cca.results$Fp[,2],length=0.1)
+#arrows(0,0,cca.results$Gs[,1],cca.results$Gs[,2],length=0.1)
+
+textxy(cca.results$Fp[1,1],cca.results$Fp[1,2],expression(X[1]))
+textxy(cca.results$Fp[2,1],cca.results$Fp[2,2],expression(X[2]))
+
+textxy(cca.results$Gs[1,1],cca.results$Gs[1,2],expression(Y[1]))
+textxy(cca.results$Gs[2,1],cca.results$Gs[2,2],expression(Y[2]))
+
+points(cca.results$V[,1],cca.results$V[,2],pch=16,cex=0.5)
+textxy(cca.results$V[,1],cca.results$V[,2],1:nrow(X),cex=0.75)
+
+
+ticklab <- seq(135,160,by=5)
+ticklabc <- ticklab-mean(Y[,2])
+ticklabs <- (ticklab-mean(Y[,2]))/sqrt(var(Y[,2]))
+
+y <- (Y[,2]-mean(Y[,2]))/sqrt(var(Y[,2]))
+Fr <- cca.results$V[,1:2]
+g <- cca.results$Gs[2,1:2]
+
+#points(cca.results$V[,1],cca.results$V[,2],cex=0.5,pch=19,col="red")
+#textxy(cca.results$V[,1],cca.results$V[,2],rownames(Xn))
+
+Cal.Data.Y2 <- calibrate(g,y,ticklabs,Fr,ticklab,lm=TRUE,tl=0.1,dp=TRUE,
+ reverse=TRUE,verb=TRUE,axislab="Y_2",
+ cex.axislab=0.75,showlabel=FALSE)
+
+#cca.results<-lm.gls(Rxy[,5]~-1+Fr,W=solve(Rxx))
+
+@
+
+We construct the biplot of the between-set correlation matrix (the joint
+plot of ${\bF}_p$ and ${\bG}_s$).
+Firstly we calibrate the biplot axis for $Y_2$ with a correlation scale.
+This calibration is done by generalised least squares with the inverse of the correlation matrix of
+the X-variables as a weight matrix ({\tt weights=solve(Rxx)}). Secondly, we calibrate the biplot axis
+for $Y_2$ with a scale for the original values. This second calibration has no weight
+matrix and is obtained by ordinary least squares. Both calibrations have a goodness of fit of 1 and
+allow perfect recovery of correlations and original data values.
+
+\subsection{Redundancy analysis}
+\label{sec:rda}
+
+Redundancy analysis can be seen as a constrained PCA. It allows two biplots, the biplot of the fitted
+values and a biplot of regression coefficients. Function {\tt rda} of the package provides a routine
+for redundancy analysis. We use Linnerud's data on physical exercise and body measurement
+variables~\cite{Tenenhaus} to illustrate calibrated biplots in redundancy analysis.
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+data(linnerud)
+X <- linnerud[,1:3]
+Y <- linnerud[,4:6]
+rda.results <- rda(X,Y)
+plot(rda.results$Fs[,1],rda.results$Fs[,2],pch=16,asp=1,xlim=c(-2,2),ylim=c(-2,2),
+ cex=0.5,xlab="1st principal axis",ylab="2nd principal axis")
+arrows(0,0,2*rda.results$Gyp[,1],2*rda.results$Gyp[,2],length=0.1)
+textxy(rda.results$Fs[,1],rda.results$Fs[,2],rownames(X),cex=0.75)
+textxy(2*rda.results$Gyp[,1],2*rda.results$Gyp[,2],colnames(Y),cex=0.75)
+
+y <- rda.results$Yh[,3]
+g <- rda.results$Gyp[3,1:2]
+Fr <- rda.results$Fs[,1:2]
+
+ticklab <- c(seq(-0.6,-0.1,by=0.1),seq(0.1,0.6,by=0.1))
+Calibrate.Yhat3 <- calibrate(g,y,ticklab,Fr,ticklab,lm=TRUE,dp=TRUE,tl=0.1,
+ axislab="Sauts",showlabel=FALSE)
+@
+
+<<fig=TRUE,echo=TRUE,aap=FALSE>>=
+plot(rda.results$Gxs[,1],rda.results$Gxs[,2],pch=16,asp=1,xlim=c(-2,2),
+ ylim=c(-2,2),cex=0.5,xlab="1st principal axis",
+ylab="2nd principal axis")
+arrows(0,0,rda.results$Gxs[,1],rda.results$Gxs[,2],length=0.1)
+arrows(0,0,rda.results$Gyp[,1],rda.results$Gyp[,2],length=0.1)
+textxy(rda.results$Gxs[,1],rda.results$Gxs[,2],colnames(X),cex=0.75)
+textxy(rda.results$Gyp[,1],rda.results$Gyp[,2],colnames(Y),cex=0.75)
+
+y <- rda.results$B[,3]
+g <- rda.results$Gyp[3,1:2]
+Fr <- rda.results$Gxs[,1:2]
+
+ticklab <- seq(-0.4,0.4,0.2)
+
+W <-cor(X)
+
+Calibrate.Y3 <- calibrate(g,y,ticklab,Fr,ticklab,lm=TRUE,dp=TRUE,tl=0.1,
+ weights=W,axislab="Sauts",showlabel=FALSE)
+ticklab <- seq(-0.4,0.4,0.1)
+Calibrate.Y3 <- calibrate(g,y,ticklab,Fr,ticklab,lm=FALSE,tl=0.05,verb=FALSE,
+ weights=W)
+ticklab <- seq(-0.4,0.4,0.01)
+Calibrate.Y3 <- calibrate(g,y,ticklab,Fr,ticklab,lm=FALSE,tl=0.025,verb=FALSE,
+ weights=W)
+@
+
+The first biplot shown is a biplot of the fitted values (obtained
+from the regression of Y onto X). Vectors for the response variables are multiplied by a factor of 3 to increase
+readability. The fitted values of the regression of Sauts onto the body measurements have
+a goodness of fit of 0.9984 and can very well be recovered by projection onto the calibrated
+axis. The second biplot is a biplot of the matrix of regression coefficients. We
+calibrated the biplot axis for "Sauts", such that the regression coefficients of the
+explanory variables with respect to "Sauts" can be recovered. The goodness of fit for
+"Sauts" is over 0.99, which means that the regression coefficients are close to
+perfectly displayed. Note that the calibration for Sauts for the regression coefficients
+is done by GLS with weight matrix equal to the correlation matrix of the X variables
+({\tt weights=W}).
+
+\section{Online documentation}
+\label{sec:online}
+
+Online documentation for the package can be obtained by typing
+{\tt vignette("CalibrationGuide"} or by accessing the file {\tt CalibrationGuide.pdf} in the {\tt doc}
+directory of the installed package.
+
+\section{Version history}
+
+Version 1.6:\\
+
+\begin{itemize}
+
+\item Function {\tt rad2degree} and {\tt shiftvector} have been added.
+
+\item Function calibrate has changed. Argument {\tt shift} from previous versions is obsolete,
+ and replaced by {\tt shiftdir, shiftfactor} and {\tt shiftvec}.
+
+\end{itemize}
+
+Version 1.7.2:\\
+
+\begin{itemize}
+\item Function {\tt textxy} has been modified and improved. Arguments {\tt dcol} and {\tt cx} no longer work, and their role has been taken over by {\tt col} and {\tt cex}. A new argument {\tt offset} controls the distance between point and label.
+\end{itemize}
+
+
+\section*{Acknowledgements}
+
+This work was partially supported by the Spanish grant BEC2000-0983. I thank Holland Genetics
+({\tt http://www.hg.nl/}), Janneke van Wagtendonk and Sander de Roos for making the calves data
+available. This document was generated by Sweave~\cite{Leisch}.
+
+\bibliographystyle{humanbio}
+\begin{thebibliography}{10}
+
+\bibitem[Anderson (1984)]{Anderson}
+Anderson, T. W.
+(1984)
+{A}n {I}ntroduction to {M}ultivariate {S}tatistical {A}nalysis
+John Wiley,
+Second edition,
+New York.
+
+\bibitem[Frets (1921)]{Frets}
+Frets, G. P.
+(1921)
+Heredity of head form in man,
+Genetica,
+3,
+pp. 193-384.
+
+\bibitem[Gabriel, 1971]{Gabriel}
+Gabriel, K. R.
+(1971)
+The biplot graphic display of matrices with application to principal component analysis.
+Biometrika
+58(3)
+pp. 453-467.
+
+\bibitem[Gower and Hand (1996)]{Gower4}
+Gower, J. C. and Hand, D. J.
+(1996)
+Biplots
+Chapman \& Hall,
+London.
+
+\bibitem[Graffelman (2005)]{Graffel16}
+Graffelman, J.
+(2005)
+Enriched biplots for canonical correlation analysis
+Journal of Applied Statistics
+32(2)
+pp. 173-188.
+
+\bibitem[Graffelman and Aluja-Banet (2003)]{Graffel13}
+Graffelman, J. and Aluja-Banet, T.
+(2003)
+Optimal Representation of Supplementary Variables in Biplots from Principal Component Analysis and Correspondence Analysis
+Biometrical Journal,
+45(4)
+pp. 491-509.
+
+\bibitem[Graffelman and van Eeuwijk (2005)]{Graffel17}
+Graffelman, J. and van Eeuwijk, F. A.,
+(2005)
+Calibration of multivariate scatter plots for exploratory analysis of relations within and between sets of variables
+in genomic research,
+Biometrical Journal,
+47,
+6,
+863-879.
+
+\bibitem[Leisch (2002)]{Leisch}
+Leisch, F.
+(2002)
+Sweave: Dynamic generation of statistical reports using literate data analysis
+Compstat 2002, Proceedings in Computational Statistics
+pp. 575-580,
+Physica Verlag, Heidelberg,
+ISBN 3-7908-1517-9
+URL http:/www.ci.tuwien.ac.at/~leisch/Sweave.
+
+\bibitem[Manly (1989)]{Manly}
+Manly, B. F. J.
+(1989)
+Multivariate statistical methods: a primer
+Chapman and Hall, London.
+
+\bibitem[Mardia et al.(1979)]{Mardia}
+Mardia, K. V. and Kent, J. T. and Bibby, J. M.
+(1979)
+Multivariate Analysis
+Academic Press London.
+
+\bibitem[R Development Core Team (2004)]{RRR}
+R Development Core Team
+(2004)
+R: A language and environment forstatistical computing.
+R Foundation for Statistical Computing,
+Vienna, Austria,
+ISBN 3-900051-00-3,
+http://www.R-project.org.
+
+\bibitem[Tenenhaus (1998)]{Tenenhaus}
+Tenenhaus, M.
+(1998)
+La R\'{e}gression PLS
+Paris, \'Editions Technip.
+
+\bibitem[Venables and Ripley (2002)]{Venables}
+Venables, W. N. and Ripley, B. D.
+(2002)
+{M}odern {A}pplied {S}tatistics with {S}-{P}lus
+New York,
+Fourth edition,
+Springer.
+
+
+\end{thebibliography}
+
+\end{document}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-med/r-cran-calibrate.git
More information about the debian-med-commit
mailing list