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]]
Andrews, Chris
2013-Sep-23 11:22 UTC
[R] PLS1 NIPALS question: error with chemometrics package?
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:bpschn01 at gmail.com] Sent: Sunday, September 22, 2013 9:44 PM To: r-help at 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