I think you need to divide by sqrt(sum(th^2)) rather than sum(th^2)
ch <- as.numeric(t(yh) %*% th)/sqrt(sum(th^2)) #
modified: / sqrt(SS)
#ch <- as.numeric(t(yh) %*% th)/sum(th^2) # modified:
/ SS
# ch <- ch/as.vector(sqrt(t(th) %*% th)) # modified:
removed normalization of ch
Chris
-----Original Message-----
From: Brad P [mailto:bpsch...@gmail.com]
Sent: Sunday, September 22, 2013 9:44 PM
To: r-help@r-project.org
Subject: [R] PLS1 NIPALS question: error with chemometrics package?
I am doing a self study, trying to understand PLS.
I have run across the following and hope that someone here can clarify as
to what is going on.
These are the data from Wold et al. (1984)
dat <- structure(list(t = 1:15, y = c(4.39, 4.42, 5, 5.85, 4.35, 4.51,
6.33, 6.37, 4.68, 5.04, 7.1, 5.04, 6, 5.48, 7.1), x1 = c(4.55,
4.74, 5.07, 5.77, 4.62, 4.41, 6.17, 6.17, 4.33, 4.62, 7.22, 4.64,
5.62, 6.19, 7.85), x2 = c(8.93, 8.93, 9.29, 9.9, 9.9, 9.93, 9.19,
9.19, 10.03, 10.29, 9.29, 10.22, 9.94, 9.77, 9.29), x3 = c(1.14,
1.14, 1.14, 1.14, 1.14, 1.14, 1.14, 1.14, 1.14, 1.14, 1.14, 1.14,
-0.07, -0.07, -0.07), x4 = c(0.7, 1.23, 0.19, 0.19, 1.23, 1.23,
0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.19, 0.19), x5 = c(0.19,
0.19, 0.7, 1.64, 1.64, 2.35, 2.83, 2.56, 2.42, 3.36, 2.43, 2.95,
1.64, 1.64, 3.8), x6 = c(0.49, 0.49, 0, -0.1, -0.1, -0.2, -0.13,
-0.13, -0.08, -0.13, -0.3, -0.08, -0.19, -0.19, -0.3), x7 = c(1.24,
1.24, 0, -0.47, -0.47, -0.51, -0.93, -0.93, -0.38, -0.93, -1.6,
-0.38, -0.47, -0.47, -1.6), x8 = c(1L, 1L, 1L, 1L, 1L, 1L, 1L,
1L, 0L, 0L, 1L, 0L, 1L, 1L, 1L)), .Names = c("t", "y", "x1",
"x2", "x3", "x4", "x5", "x6", "x7", "x8"), class = "data.frame", row.names
= c(NA,
-15L))
In Wold et al. (1984) (pg 741, Table 2) beta coefficients are given for PLS
(1) and PLS (2) where (1) and (2) correspond to the number of PLS
components retained. The coefficients are not the same as those when I run
the following:
library("chemometrics")
# retaining 1 component:
pls1_nipals(X=dat[,c(3:10)], y=dat[2], a=1, scale=TRUE)$b
# retaining 2 components:
pls1_nipals(X=dat[,c(3:10)], y=dat[2], a=2, scale=TRUE)$b
However, if I modify the source code like this:
pls1_nipals_mod <- function(X, y, a, it = 50, tol = 1e-08, scale = FALSE)
{
Xh <- scale(X, center = TRUE, scale = scale)
yh <- scale(y, center = TRUE, scale = scale)
T <- NULL
P <- NULL
C <- NULL
W <- NULL
for (h in 1:a) {
wh <- t(Xh) %*% yh/sum(yh^2) # modified: / SS
wh <- wh/as.vector(sqrt(t(wh) %*% wh))
th <- Xh %*% wh
ch <- as.numeric(t(yh) %*% th)/sum(th^2) # modified: / SS
# ch <- ch/as.vector(sqrt(t(th) %*% th)) # modified: removed
normalization of ch
ph <- t(Xh) %*% th/as.vector(t(th) %*% th)
Xh <- Xh - th %*% t(ph)
yh <- yh - th * ch
T <- cbind(T, th)
P <- cbind(P, ph)
C <- c(C, ch)
W <- cbind(W, wh)
}
b <- W %*% solve(t(P) %*% W) %*% C
list(P = P, T = T, W = W, C = C, b = b)
}
pls1_nipals_mod(X=dat[,c(3:10)], y=dat[2], a=1, scale=TRUE)$b
pls1_nipals_mod(X=dat[,c(3:10)], y=dat[2], a=2, scale=TRUE)$b
These beta coefficients are exactly the same as in Wold et al. (1984)
Furthermore, if I do a leave-one-out CV, my modified version has a PRESS =
1.27, whereas the original pls1_nipals() function has a PRESS = 18.11!
That's not good, right? Here is my LOOCV code, 1:1 lines added in plots:
### LOOCV for original function
out.j <- vector("list", length=nrow(dat))
for(j in c(2:nrow(dat), 1)){
b <- pls1_nipals(X=dat[-j,c(3:10)], y=dat[-j,2], a=2, scale=TRUE)$b
dats <- scale(dat)
y.est <- dats[j,c(3:10)] %*% b
y.obs <- dats[j,2]
out.j[[j]] <- data.frame(y.obs, y.est)
}
out <- do.call(rbind, out.j)
sqrt(sum((out[,1]-out[,2])^2) )
plot(out[,2]~ out[,1], ylab="pred", xlab="obs")
abline(0,1, col="grey")
### LOOCV for modified function
out.j <- vector("list", length=nrow(dat))
for(j in c(2:nrow(dat), 1)){
b <- pls1_nipals_mod(X=dat[-j,c(3:10)], y=dat[-j,2], a=2, scale=TRUE)$b
dats <- scale(dat)
y.est <- dats[j,c(3:10)] %*% b
y.obs <- dats[j,2]
out.j[[j]] <- data.frame(y.obs, y.est)
}
out <- do.call(rbind, out.j)
sqrt(sum((out[,1]-out[,2])^2) )
plot(out[,2]~ out[,1], ylab="pred", xlab="obs")
abline(0,1, col="grey")
Is this an error with the chemometrics function; or am I simply
understanding something incorrectly?
Thank you.
Citation: Wold, S., A. Ruhe, H. Wold, W. Dunn. (1984) The collinearity
problem in linear regression. The partial least squares (PLS) approach to
generalized inverses*" SIAM J. Sci. Stat. Comput.
[[alternative HTML version deleted]]
**********************************************************
Electronic Mail is not secure, may not be read every day, and should not be
used for urgent or sensitive issues
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
