偏最小二乘迴歸分析建模步驟的R實現(康復俱樂部20名成員測試數據)+補充pls迴歸係數矩陣的算法實現

kf=read.csv('d:/kf.csv') # 讀取康復數據
kf
sl=as.matrix(kf[,1:3]) #生成生理指標矩陣
xl=as.matrix(kf[,4:6]) #生成訓練指標矩陣
x=sl
x
y=xl
y
x0=scale(x)
x0
y0=scale(y)
y0
m=t(x0)%*%y0%*%t(y0)%*%x0
m
eigen(m)
w1=eigen(m)$vectors[,1]
v1=t(y0)%*%x0%*%w1/sqrt(as.matrix(eigen(m)$values)[1,])
v1
t1=x0%*%w1 #第一對潛變量得分向量
t1 # 以上爲第一步(1)分別提取兩變量組的第一對成分，並使之相關性達最大。
u1=y0%*%v1
u1 #第一對潛變量得分向量
library("pracma")
α1=inv(t(t1)%*%t1)%*%t(t1)%*%x0 #也可由t(x0)%*%t1/norm(t1,'2')^2算得α1 #α1在pls中稱爲模型效應負荷量
β1=inv(t(t1)%*%t1)%*%t(t1)%*%y0 #也可由t(y0)%*%t1/norm(t1,'2')^2算得β1
t(x0)%*%t1/norm(t1,'2')^2 # norm(t1,'2')爲svd(t1)即t1的最大奇異值,也可用sqrt(t(t1)%*%t1)求得
t(y0)%*%t1/norm(t1,'2')^2 # 以上爲第二步(2)創建y對T1的迴歸及x對T1的迴歸。
α1
β1
lm(x0~t1) #驗證α1即爲x0作應變量，t1作自變量的最小二回歸的t1的迴歸係數(分別爲weight、waist和pulse的迴歸係數，共3個)
lm(y0~t1) #驗證β1即爲y0作應變量，t1作自變量的最小二回歸的t1的迴歸係數(分別爲chins、situps和jumps的迴歸係數，共3個)
B=t(x0)%*%u1%*%inv(t(t1)%*%x0%*%t(x0)%*%u1)%*%t(t1)%*%y0 #保留第一對潛變量對應的標準化自變量x和標準化應變量y的pls迴歸係數矩陣(該矩陣公式參見‘KernelPartialLeastSquaresRegressionin Reproducing KernelHilbert Space’p102)

B
library("pls")
pls1=plsr(y0~x0,ncomp=1,validation='LOO',jackknife=T)
coef(pls1) #上式中B的求解等價於R的pls包中保留一個主成分的結果,係數爲標準化迴歸係數，能夠經過逆標準化過程還原爲原始自變量x和應變量y的迴歸係數。如下保留2個主成分的結果中有具體逆標準化過程。
px0=t1%*%α1 #求x0的預測值矩陣
E1=x0-px0 #求x0的殘差矩陣
py0=t1%*%β1 #求y0的預測值矩陣
F1=y0-py0 # #求y0的殘差矩陣  
m2=t(E1)%*%F1%*%t(F1)%*%E1 #（3）用殘差陣E1和F1代替x0和y0重複以上步驟。
eigen(m2)
w2=eigen(m2)$vectors[,1]
w2
v2=t(F1)%*%E1%*%w2/sqrt(as.matrix(eigen(m2)$values)[1,])
v2
t2=E1%*%w2
t2
u2=F1%*%v2
u2
α2=inv(t(t2)%*%t2)%*%t(t2)%*%E1 #也可由t(E1)%*%t2/norm(t2,'2')^2算得α2
β2=inv(t(t2)%*%t2)%*%t(t2)%*%F1 #也可由t(F1)%*%t2/norm(t2,'2')^2算得β2
α2
β2
library("pls")
pls1=plsr(y0~x0,ncomp=2,validation='LOO',jackknife=T) #如下爲R中pls包運算結果，顯示迴歸結果(包括預測值偏差平方和PRESS與變異解釋度)，與上述純算法結果進行對比和補充，
summary(pls1) #其中對於解釋變量潛變量T1對應變量y的總變異解釋的比例爲chins(23.26%)、situps(35.06%)和jumps(4.14%)等價於SAS中對y的綜合結果20.9447≈mean(23.26%,35.06%,4.14%)四捨五入形成的。2 comps列顯示的爲引入第二解釋變量潛變量後的對應變量y的總變異解釋的比例。
coef(pls1) #以應變量situps爲例得situps關於各自變量的迴歸方程（*表示標準化）：situps*=-0.13846688weight*-0.52444579waist*-0.08542029pulse* 據此標準化迴歸方程能夠推導出原始變量y與x的迴歸方程:(situps-mean(situps))/sd(situps)=-0.13846688*(weight-mean(weight))/sd(weight)-0.52444579*(waist-mean(waist))/sd(waist)-0.08542029*(pulse-mean(pulse))/sd(pulse)——>situps=sd(situps)[-0.13846688*(weight-mean(weight))/sd(weight)-0.52444579*(waist-mean(waist))/sd(waist)-0.08542029*(pulse-mean(pulse))/sd(pulse)]+mean(waist)
sd(y[,2])*-0.1384668393/sd(x[,1]) #weight的迴歸係數
sd(y[,2])*-0.52444579/sd(x[,2]) #waist的迴歸係數
sd(y[,2])*-0.08542029/sd(x[,3]) #pulse的迴歸係數
sd(y[,2])*(-0.13846688*-mean(x[,1])/sd(x[,1])+-0.52444579*-mean(x[,2])/sd(x[,2])+-0.08542029*-mean(x[,3])/sd(x[,3]))+mean(y[,2]) #原始變量y與x的迴歸方程截距
model="SITUPS=612.56712-0.35088WEIGHT-10.24768WAIST-0.74122PULSE——耶！與SAS給出的結果徹底一致。"
model
jack.test(pls1) #即對coef(pls1)生成的係數進行假設檢驗
scores(pls1) #即求第一解釋潛變量的得分向量t1=x0%*%w1和第二解釋變量潛變量的得分向量t2=E1%*%w2
loadings(pls1) #即求α1
plot(pls1)
validationplot(pls1) #validationplot()函數能夠畫出PLS模型在不一樣主成分數下對應的RMSEP(由留一交叉驗證法算得的均方預測偏差根)
predict(pls1) #即求py0=t1%*%β1
#關於決定係數算法還需研究