pls.cox=
function (X, Y, ncomp = 2, mode = c("regression", "canonical",
"invariant", "classic"), max.iter = 500, tol = 1e-06, scale.X = TRUE,
scale.Y = TRUE, ...)
{
force(ncomp)
if (length(dim(X)) != 2)
stop("'X' must be a numeric matrix.")
X = as.matrix(X)
Y = as.matrix(Y)
if (!is.numeric(X) || !is.numeric(Y))
stop("'X' and/or 'Y' must be a numeric matrix.")
n = nrow(X)
q = ncol(Y)
if ((n != nrow(Y)))
stop("unequal number of rows in 'X' and 'Y'.")
if (is.null(ncomp) || !is.numeric(ncomp) || ncomp <= 0)
stop("invalid number of variates, 'ncomp'.")
nzv = mixOmics::nearZeroVar(X, ...)
if (length(nzv$Position > 0)) {
warning("Zero- or near-zero variance predictors. \n Reset predictors matrix to not near-zero variance predictors.\n See $nzv for problematic predictors.")
X = X[, -nzv$Position]
}
p = ncol(X)
ncomp = round(ncomp)
if (ncomp > p) {
warning("Reset maximum number of variates 'ncomp' to ncol(X) = ",
p, ".")
ncomp = p
}
mode = match.arg(mode)
X.names = dimnames(X)[[2]]
if (is.null(X.names))
X.names = paste("X", 1:p, sep = "")
if (dim(Y)[2] == 1)
Y.names = "Y"
else {
Y.names = dimnames(Y)[[2]]
if (is.null(Y.names))
Y.names = paste("Y", 1:q, sep = "")
}
ind.names = dimnames(X)[[1]]
if (is.null(ind.names)) {
ind.names = dimnames(Y)[[1]]
rownames(X) = ind.names
}
if (is.null(ind.names)) {
ind.names = 1:n
rownames(X) = rownames(Y) = ind.names
}
if (scale.X) {
X = scale(X, center = TRUE, scale = TRUE)
}
if (scale.Y) {
Y = scale(Y, center = TRUE, scale = TRUE)
}
X.temp = X
Y.temp = Y
mat.t = matrix(nrow = n, ncol = ncomp)
mat.u = matrix(nrow = n, ncol = ncomp)
mat.a = matrix(nrow = p, ncol = ncomp)
mat.b = matrix(nrow = q, ncol = ncomp)
mat.c = matrix(nrow = p, ncol = ncomp)
mat.d = matrix(nrow = q, ncol = ncomp)
mat.e = matrix(nrow = q, ncol = ncomp)
n.ones = rep(1, n)
p.ones = rep(1, p)
q.ones = rep(1, q)
na.X = FALSE
na.Y = FALSE
is.na.X = is.na(X)
is.na.Y = is.na(Y)
if (any(is.na.X))
na.X = TRUE
if (any(is.na.Y))
na.Y = TRUE
for (h in 1:ncomp) {
u = Y.temp[, 1]
if (any(is.na(u)))
u[is.na(u)] = 0
a.old = 0
b.old = 0
iter = 1
if (na.X) {
X.aux = X.temp
X.aux[is.na.X] = 0
}
if (na.Y) {
Y.aux = Y.temp
Y.aux[is.na.Y] = 0
}
repeat {
if (na.X) {
a = crossprod(X.aux, u)
U = drop(u) %o% p.ones
U[is.na.X] = 0
u.norm = crossprod(U)
a = a/diag(u.norm)
a = a/drop(sqrt(crossprod(a)))
t = X.aux %*% a
A = drop(a) %o% n.ones
A[t(is.na.X)] = 0
a.norm = crossprod(A)
t = t/diag(a.norm)
}
else {
a = crossprod(X.temp, u)/drop(crossprod(u))
a = a/drop(sqrt(crossprod(a)))
t = X.temp %*% a/drop(crossprod(a))
}
if (na.Y) {
b = crossprod(Y.aux, t)
T = drop(t) %o% q.ones
T[is.na.Y] = 0
t.norm = crossprod(T)
b = b/diag(t.norm)
u = Y.aux %*% b
B = drop(b) %o% n.ones
B[t(is.na.Y)] = 0
b.norm = crossprod(B)
u = u/diag(b.norm)
}
else {
b = crossprod(Y.temp, t)/drop(crossprod(t))
u = Y.temp %*% b/drop(crossprod(b))
}
if (crossprod(a - a.old) < tol)
break
if (iter == max.iter) {
warning(paste("Maximum number of iterations reached for dimension",
h), call. = FALSE)
break
}
a.old = a
b.old = b
iter = iter + 1
}
if (na.X) {
X.aux = X.temp
X.aux[is.na.X] = 0
c = crossprod(X.aux, t)
T = drop(t) %o% p.ones
T[is.na.X] = 0
t.norm = crossprod(T)
c = c/diag(t.norm)
}
else {
c = crossprod(X.temp, t)/drop(crossprod(t))
}
X.temp = X.temp - t %*% t(c)
if (mode == "canonical") {
if (na.Y) {
Y.aux = Y.temp
Y.aux[is.na.Y] = 0
e = crossprod(Y.aux, u)
U = drop(u) %o% q.ones
U[is.na.Y] = 0
u.norm = crossprod(U)
e = e/diag(u.norm)
}
else {
e = crossprod(Y.temp, u)/drop(crossprod(u))
}
Y.temp = Y.temp - u %*% t(e)
}
if (mode == "classic")
Y.temp = Y.temp - t %*% t(b)
if (mode == "regression") {
if (na.Y) {
Y.aux = Y.temp
Y.aux[is.na.Y] = 0
d = crossprod(Y.aux, t)
T = drop(t) %o% q.ones
T[is.na.Y] = 0
t.norm = crossprod(T)
d = d/diag(t.norm)
}
else {
d = crossprod(Y.temp, t)/drop(crossprod(t))
}
Y.temp = Y.temp - t %*% t(d)
}
if (mode == "invariant")
Y.temp = Y
mat.t[, h] = t
mat.u[, h] = u
mat.a[, h] = a
mat.b[, h] = b
mat.c[, h] = c
if (mode == "regression")
mat.d[, h] = d
if (mode == "canonical")
mat.e[, h] = e
}
rownames(mat.a) = rownames(mat.c) = X.names
rownames(mat.b) = Y.names
rownames(mat.t) = rownames(mat.u) = ind.names
comp = paste("comp", 1:ncomp)
colnames(mat.t) = colnames(mat.u) = comp
colnames(mat.a) = colnames(mat.b) = colnames(mat.c) = comp
cl = match.call()
cl[[1]] = as.name("pls")
result = list(call = cl, X = X, Y = Y, ncomp = ncomp, mode = mode,
mat.c = mat.c, mat.d = mat.d, mat.e = mat.e, variates = list(X=mat.t,
Y = mat.u), loadings = list(X = mat.a, Y = mat.b),
names = list(X = X.names, Y = Y.names, indiv = ind.names))
if (length(nzv$Position > 0))
result$nzv = nzv
class(result) = "pls"
return(invisible(result))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.