[R] PLS1 NIPALS question: error with chemometrics package?

[Andrews, Chris]

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