So there is little explanation needed. There are probably implementations around there, but I got bored and decided to do mine as well. May contain errors.

pls <- function(y,X){ ## Standarize variables p <- length(X[1,]) n <- length(X[,1]) for(i in 1:p){ X[,i] <- (X[,i]-mean(X[,i]))/sd(X[,i]) } # Compute PLS z <- matrix(NA,n,p) t <- c() for(m in 1:p){ phi <- t(y)%*%X z[,m] <- phi%*%t(X) theta <- y%*%z[,m]/(z[,m]%*%z[,m]) t <- c(t,theta) X <- X - (as.numeric(t(X)%*%z[,m])/(z[,m]%*%z[,m]))*z[,m] } return(list(z=z,theta=theta)) }