##2020.02.03 by jingshan lu wrote
rm(list =ls() ) #clear all
classiftime1 <- proc.time()
library("openxlsx")
jianmo=read.xlsx("C:\\Users\\20162\\Desktop\\建模.xlsx")# this is calibration dataset, first column is dependent variable, other columns are independent variables
yanzheng=read.xlsx("C:\\Users\\20162\\Desktop\\验证.xlsx")#this isvalidation dataset
jianmoy=as.matrix(jianmo[,1])
jianmox=as.data.frame(jianmo[,-1])
yanzhengy=as.matrix(yanzheng[,1])
yanzhengx=as.data.frame(yanzheng[,-1])
##
library(pls)
library(plsVarSel)#VIP of this package can be used to calculate VIP of variables
## set the parameter unchanged
set.seed(1)
y.plsr=plsr(jianmoy~.,data=jianmox,scale=T,validation="CV")
component=which.min(y.plsr$validation$PRESS)
jianmopre=as.matrix(as.data.frame(y.plsr$validation$pred)[,component])
jianmoshice=jianmoy
#calculate VIP
y.plsr.vip=as.matrix(VIP(y.plsr,component))
#
jianmomoxing=cbind(jianmoshice,jianmopre)
plot(jianmoshice,jianmopre)
write.xlsx(jianmomoxing,file="C:\\Users\\20162\\Desktop\\建模模型.xlsx")#export calibration model that two columns data including measured and predicted dependent variable
write.xlsx(y.plsr.vip,file="C:\\Users\\20162\\Desktop\\模型VIP.xlsx")#export VIP of plsr model
## model validation
yanzhengyuce=as.matrix(as.data.frame(predict(y.plsr,yanzhengx))[,component])
yanzhengmoxing=cbind(yanzhengy,yanzhengyuce)
plot(yanzhengy,yanzhengyuce)
write.xlsx(yanzhengmoxing,file="C:\\Users\\20162\\Desktop\\验证模型.xlsx")#export validation model

# validationplot(y.plsr,val.type="RMSEP") show the RMSEP of the model
##

classiftime <-  proc.time() - classiftime1