-
Notifications
You must be signed in to change notification settings - Fork 15
/
gsorth.R
74 lines (67 loc) · 2.5 KB
/
gsorth.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
##################################
#' Gram-Schmidt Orthogonalization of a Matrix
#'
#' Calculates a matrix with uncorrelated columns using the Gram-Schmidt process
#'
#' This function, originally from the \pkg{heplots} package has now been deprecated in \pkg{matlib}. Use
#' \code{\link{GramSchmidt}} instead.
#'
#' @param y a numeric matrix or data frame
#' @param order if specified, a permutation of the column indices of \code{y}
#' @param recenter logical; if \code{TRUE}, the result has same means as the original \code{y}, else means = 0 for cols 2:p
#' @param rescale logical; if \code{TRUE}, the result has same sd as original, else, sd = residual sd
#' @param adjnames logical; if \code{TRUE}, colnames are adjusted to Y1, Y2.1, Y3.12, ...
#' @return a matrix/data frame with uncorrelated columns
#' @export
#' @examples
#' \dontrun{
#' set.seed(1234)
#' A <- matrix(c(1:60 + rnorm(60)), 20, 3)
#' cor(A)
#' G <- gsorth(A)
#' zapsmall(cor(G))
#' }
# Return a matrix/data frame with uncorrelated columns
# recenter=TRUE -> result has same means as original, else means = 0 for cols 2:p
# rescale=TRUE -> result has same sd as original, else, sd = residual sd
# adjnames=TRUE -> colnames are adjusted to Y1, Y2.1, Y3.12, ...
# 12-5-2010: Fixed buglet when matrix has no row/col names
# 10-14-2011: Made sd() a local function to avoid deprecated warnings
gsorth <- function(y, order, recenter=TRUE, rescale=TRUE, adjnames=TRUE) {
# local function sd(), since sd(<matrix>) and sd(<data.frame>) now deprecated
sd <- function (x, na.rm = FALSE)
{
if (is.matrix(x))
apply(x, 2, sd, na.rm = na.rm)
else if (is.vector(x))
sqrt(var(x, na.rm = na.rm))
else if (is.data.frame(x))
sapply(x, sd, na.rm = na.rm)
else sqrt(var(as.vector(x), na.rm = na.rm))
}
.Deprecated("GramSchmidt")
n <- nrow(y)
if (missing(order)) order <- 1:ncol(y)
y <- y[,order]
p <- ncol(y)
if (is.data.frame(y)) {
numeric <- unlist(lapply(y, is.numeric))
if (!all(numeric)) stop("all columns of y must be numeric")
}
ybar <- colMeans(y)
ysd <- sd(y)
z <- scale(y, center=TRUE, scale=FALSE)
z <- qr.Q(qr(z))
zsd <- sd(z)
if (rescale) z <- z %*% diag( ysd/zsd )
if (recenter) z <- z + matrix(rep(ybar,times=n), ncol=p, byrow=TRUE)
rownames(z) <- rownames(y, do.NULL=FALSE)
colnames(z) <- colnames(y, do.NULL=FALSE)
if (adjnames) {
for (j in 2:p) {
colnames(z)[j] <- paste(colnames(z)[j], '.', sep="",
paste( 1:(j-1), collapse=""))
}
}
z
}