if(!require(pls)) {install.packages("pls"); library(pls)}
#[1]========== preparation ==========
ifname<- "DataAgonist.csv"
alldata <- read.delim(ifname, row.names=1,header=T, sep= ",", as.is=TRUE, 
                      strip.white=FALSE)
dataUsed<-alldata
yname<-colnames(alldata[1]); xname<-colnames(dataUsed[,-1])
nf<-ncol(dataUsed)-1
#[2] PLS---------------------------
result1.pls<-plsr(y~x1+x2+x3+x4+x5+x6+x7+x8,data=dataUsed,scale=FALSE,validation="LOO")
#result1.pls<-plsr(as.formula(paste0(yname,"~",paste0(xname,collapse="+"))),
#                  data=dataUsed,scale=FALSE,validation="LOO")
names(result1.pls); summary(result1.pls)
msepout <- MSEP(result1.pls); msepout
plot(msepout,legendpos="topright",main="Mean squared error of prediction")
selectNcomp(result1.pls,method="onesigma",plot=T)
#[3][sum((dataUsed$Y-Ypred1)^2) == plspress]--------------------
plspress<-result1.pls$validation$PRESS; plspress
ncomponent<-which.min(plspress); ncomponent
Ypred1<-result1.pls$validation$pred[,,ncomponent]; Ypred1
#[4]Q2---
tss<-sum((dataUsed[,1]-mean(dataUsed[,1]))^2)
Q2<-1-plspress/tss; Q2
barplot(Q2,main="Q2")
cor(dataUsed[,1],Ypred1)
s1 <- paste("PLS: r = ",sprintf("%5.3f",cor(dataUsed[,1],Ypred1)))
plot(dataUsed[,1],Ypred1,xlab="original",ylab="leave-one-out",main=s1)
abline(0,1,col="black")
#[5]R2---
ncomponent<-which.min(plspress)
Ypred2<-predict(result1.pls,newdata=dataUsed, ncomp = ncomponent)
Ypred2
rss<-sum((dataUsed[,1]-Ypred2)^2)
R2<-1-rss/tss; R2
#[6]--------------------------------------
cor(dataUsed[,1],Ypred2)
s2 <- paste("PLS: r = ",sprintf("%5.3f",cor(dataUsed[,1],Ypred2)))
plot(dataUsed[,1],Ypred2,xlab="original",ylab="prediction",main=s2)
abline(0,1,col="black")
barplot(coefficients(result1.pls,ncomp=ncomponent)[,1,1],horiz=TRUE,las=1,
        main=colnames(alldata[1]))
abline(v=0,lty=1,col="black")
#------------------------------------------------
coefficients(result1.pls,ncomp=ncomponent)[,1,1]
#[7]--- MRA: Multiregression Analysis
yregmodel<-lm(y~x1+x2+x3+x4+x5+x6+x7+x8,data=dataUsed)
#yregmodel<-lm(as.formula(paste0(yname,"~",paste0(xname,collapse="+"))),data=dataUsed)
summary(yregmodel)
ypred3<-predict(yregmodel)
cor(dataUsed[,1],ypred3)
s3 <- paste("MRA: r = ",sprintf("%5.3f",cor(dataUsed[,1],ypred3)))
plot(dataUsed[,1],ypred3,xlab="original",ylab="prediction",main=s3)
abline(0,1,col="black")