x1 <- as.numeric(round(runif(10,-40,40),2))
x2 <- as.numeric(round(x1*1.4+60,2))
x3 <- as.numeric(round(runif(10,20,60),2))
x4 <- as.numeric(round(x2*0.9+60,2))
x5 <- as.numeric(round(x2*x3*0.9+60,2))
x6 <- as.numeric(round(x2*x3*x4*x5/1000000,2))

y <-  as.numeric(round(runif(10,50,150),2))

df <- data.frame(y,x1,x2,x3,x4,x5,x6)


library(pls)
# plsr, RMSEP
mod.plsr <- plsr(y~x1+x2+x3+x4+x5+x6, data=df,
             ncomp=5, validation="CV")

## delta vector contains RMSEP differences
err.CV = c()
for (i in 1:10) {err.CV[i] = RMSEP(mod.plsr)$val[i*2+1]}
delta = err.CV[1:9] - err.CV[2:10]
comp.plsr = min(which(delta<0.05))
plot(RMSEP(mod.plsr),legendpos="topright", main="")

## mixed model regression coefficients
mod.plsr.opt = plsr(y~x1+x2+x3+x4+x5+x6, data=df,
                ncomp = comp.plsr)

coef(mod.plsr.opt)

, , 1 comps

           y
x1 4.324635e-05
x2 6.054166e-05
x3 3.218208e-05
x4 5.449111e-05
x5 4.142277e-03
x6 4.653091e-03




library(plsRglm)
# plsrglm, BIC
mod.plsrglm = plsRglm(y~x1+x2+x3+x4+x5+x6, data=df,
                  nt=5, model="pls")

# use BIC to determine optimal number of components
comp.plsrglm = which(mod.plsrglm$InfCrit[,2] ==    min(mod.plsrglm$InfCrit[,2]))-1

# refit model and extract beta coefficients from the optimal model
mod.plsrglm.opt = plsRglm(y~x1+x2+x3+x4+x5+x6, data=df,
                      nt=comp.plsr, model="pls")

mod.plsrglm.opt$Coeffs

               [,1]
Intercept -4.422569e+05
x1        -3.150225e+03
x2        -2.355536e+03
x3         4.523422e+00
x4         5.120661e+03
x5        -1.490321e-01
x6         7.920704e-02

coef(mod.plsr.opt, intercept = TRUE)