<4D F736F F D20B5DAC8FDCAAED5C22020C6ABD7EED0A1B6FEB3CBBBD8B9E92E646F63>

Size: px

Start display at page:

Download "<4D F736F F D20B5DAC8FDCAAED5C22020C6ABD7EED0A1B6FEB3CBBBD8B9E92E646F63>"

军融
5 years ago
Views:

1 第三十章偏最小二乘回归在实际问题中, 经常遇到需要研究两组多重相关变量间的相互依赖关系, 并研究用一组变量 ( 常称为自变量或预测变量 ) 去预测另一组变量 ( 常称为因变量或响应变量 ), 除了最小二乘准则下的经典多元线性回归分析 (MLR), 提取自变量组主成分的主成分回归分析 (PCR) 等方法外, 还有近年发展起来的偏最小二乘 (PLS) 回归方法偏最小二乘回归提供一种多对多线性回归建模的方法, 特别当两组变量的个数很多, 且都存在多重相关性, 而观测数据的数量 ( 样本量 ) 又较少时, 用偏最小二乘回归建立的模型具有传统的经典回归分析等方法所没有的优点偏最小二乘回归分析在建模过程中集中了主成分分析, 典型相关分析和线性回归分析方法的特点, 因此在分析结果中, 除了可以提供一个更为合理的回归模型外, 还可以同时完成一些类似于主成分分析和典型相关分析的研究内容, 提供更丰富深入的一些信息本章介绍偏最小二乘回归分析的建模方法 ; 通过例子从预测角度对所建立的回归模型进行比较偏最小二乘回归考虑个变量 y, y,, y 与个自变量 x, x,, x 的建模问题偏最小二乘回归的基本作法是首先在自变量集中提出第一成分 t ( t 是 x,, x 的线性组合, 且尽可能多地提取原自变量集中的变异信息 ); 同时在因变量集中也提取第一成分 u, 并要求 t 与 u 相关程度达到最大然后建立因变量 y,, y 与 t 的回归, 如果回归方程已达到满意的精度, 则算法中止否则继续第二对成分的提取, 直到能达到满意的精度为止若最终对自变量集提取个成分 t, t,, t, 偏最小二乘回归将通过建立 y,, y 与 t, t,, t 的回归式, 然后再表示为 y,, y 与原自变量的回归方程式, 即偏最小二乘回归方程式为了方便起见, 不妨假定个因变量 y,, y x,, x 均为标准与个自变量化变量因变量组和自变量组的 n 次标准化观测数据阵分别记为 F y yn y y n, 偏最小二乘回归分析建模的具体步骤如下 : E x xn x x n -5-

2 () 分别提取两变量组的第一对成分, 并使之相关性达最大假设从两组变量分别提出第一对成分为 t 和 u, t 是自变量集线性组合 : t 合 : u ( x,, x 的 X ) wx + + w x w X, u 是因变量集 Y ( y,, y ) 的线性组 v y + + v y v Y 为了回归分析的需要, 要求 : t 和 u 各自尽可能多地提取所在变量组的变异信息 ; t 和 u 的相关程度达到最大由两组变量集的标准化观测数据阵 E 和 F, 可以计算第一对成分的得分向量, 记为 ˆt 和 û : E w ˆ u F v x xn y yn x x n w w t t n y v u y n v un 第一对成分 t 和 u 的协方差 Cov( t, u ) 可用第一对成分的得分向量 ˆt 和 û 的内积来计算故而以上两个要求可化为数学上的条件极值问题 : < ˆ t ˆ, u >< Ew, Yv > w E F x ax w w w, v v v -5- 利用 Lagange 乘数法, 问题化为求单位向量 w 和 v, 使 θ w E Fv 最大问题的求解只须通过计算矩阵 M E F F E 的特征值和特征向量, 且 M 的最大特 F θ 征值为 θ, 相应的单位特征向量就是所求的解 w, 而 v 可由 w 计算得到 y,, y () 建立对 t 的回归及假定回归模型为 x,, x 对 t 的回归 v E w

3 其中 α E F ˆ tα + E uˆ β + F ( α,, α, ) β (,, ) β β E 和 F 是残差阵回归系数向量 α, β 的最小二乘估计为分别是多对一的回归模型中的参数向量, 称 α E β F α, β 为模型效应负荷量, () 用残差阵 E 和 F 代替 E 和 F 重复以上步骤记 ˆ ˆ, ˆ E tα ˆ, 则残差阵 E E Ê, F F ˆF 如果残差阵 F F tβ 中元素的绝对值近似为, 则认为用第一个成分建立的回归式精度已满足需要了, 可以停止抽取成分否则用残差阵 E 和 F 代替 E 和 F 重复以上步骤即得 : w ) ( w,, w ; v ) ( v,, v 分别为第二对成分的权数而 ˆ t Ew, u F v ˆ E α, 为第二对成分的得分向量 F β 分别为 X, Y 的第二对成分的负荷量这时有 E F α β + α + β + E + F (4) 设 n 数据阵 E 的秩为 in( n, ), 则存在个成分 t, t,, t, 使得 E F α β + + α + + β + E + F 把 t k wk x + + wkx ( k,,, ), 代入 Y t β + + t β, 即得个因变量的偏最小二乘回归方程式 -5-

4 x + a x,(,,, y a + ) (5) 交叉有效性检验一般情况下, 偏最小二乘法并不需要选用存在的个成分 t, t,, t 来建立回归式, 而像主成分分析一样, 只选用前 l 个成分 ( l ), 即可得到预测能力较好的回归模型对于建模所需提取的主成分个数 l, 可以通过交叉有效性检验来确定每次舍去第 i 个观测 ( i,,, n ), 用余下的 n 个观测值按偏最小二乘回归方法建模, 并考虑抽取个成分后拟合的回归式, 然后把舍去的第 i 个观测点代入所拟合的回归方程式, 得到 y (,,, ) 在第 i 个观测点上的预测值 y ( i ( ) 对 i,,, n 重复以上的验证, 即得抽取个成分时第个因变量 (,,, ) 的预测误差平方和为 Y ) n PRESS ( ) ( y yˆ ( ) ( )) (,,, ) i ( y,, y 的预测误差平方和为 i i y ˆ ) -54- PRESS ( ) PRESS ( ) i 另外, 再采用所有的样本点, 拟合含个成分的回归方程这时, 记第 i 个样本点的预测值为 yˆ i ( ), 则可以定义 y 的误差平方和为 SS ( ) 定义 Y 的误差平方和为 SS ( ) n i ( y i yˆ SS ( ) i ( )) 当 PRESS( ) 达到最小值时, 对应的即为所求的成分个数通常, 总有 PRESS( ) 大于 SS( ), 而 SS( ) 则小于 SS( ) 因此, 在提取成分时, 总希望比值 PRESS( ) SS( ) 越小越好 ; 一般可设定限制值为.5, 即当

5 PRESS ( ) SS( ) (.5).95 时, 增加成分 t 有利于模型精度的提高或者反过来说, 当 PRESS ( ) SS( ) >.95 时, 就认为增加新的成分 t, 对减少方程的预测误差无明显的改善作用为此, 定义交叉有效性为 PRESS( ) SS( ), 这样, 在建模的每一步 Q 计算结束前, 均进行交叉有效性检验, 如果在第步有 Q < , 则模型达到精度要求, 可停止提取成分 ; 若 Q. 975 贡献显著, 应继续第 + 步计算, 表示第步提取的 t 成分的边际一种更简洁的计算方法上节介绍的算法原则和推导过程的思路在目前的文献中是最为常见的然而, 还有一种更为简洁的计算方法, 即直接在 E, E 矩阵中提取成分 t,,t ( ) 要, 求 t 能尽可能多地携带 X 中的信息, 同时, t 对因变量系统 F 有最大的解释能力注意, 无需在 F 中提取成分得分 u, 这可以使计算过程大为简化, 并且对算法结论的解释也更为方便偏最小二乘法的简记算法的步骤如下 : () 求矩阵 E F F E 最大特征值所对应的特征向量 w, 求得成分 t w X, 计算成分得分向量 t ˆ Ew, 和残差矩阵 ˆα, 其中 α E E E t () 求矩阵 E F F E 最大特征值所对应的特征向量 w, 求得成分 t w X, 计算成分得分向量 t ˆ Ew, 和残差矩阵 E E t ˆ α, 其中 α E () 至第步, 求矩阵 E F F E 最大特征值所对应的特征向量 w, 求得成分 t w X, 计算成分得分向量 E w -55-

6 如果根据交叉有效性, 确定共抽取个成分 t,,t 可以得到一个满意的预测模型, 则求 F 在 t, ˆ, 上的普通最小二乘回归方程为 ˆ β ˆ F t + + tβ + F 把 t k w x w x ( k,,, ), 代入 Y t β + + t β, 即得个 k + + k 因变量的偏最小二乘回归方程式 x + a x,(,,, y a + ) 这里的 w 满足 ˆ t Ew, w ( I w α ) w 案例分析本节采用兰纳胡德 (Linneud) 给出的关于体能训练的数据进行偏最小二乘回归建模在这个数据系统中被测的样本点, 是某健身俱乐部的位中年男子被测变量分为两组第一组是身体特征指标 X, 包括 : 体重腰围脉搏第二组变量是训练结果指标 Y, 包括 : 单杠弯曲跳高原始数据见表表体能训练数据 No 体重 (x ) 腰围 (x ) 脉搏 (x ) 单杠 (y ) 弯曲 (y ) 跳高 (y )

7 均值标准差表给出了这 6 个变量的简单相关系数矩阵从相关系数矩阵可以看出, 体重与腰围是正相关的 ; 体重腰围与脉搏负相关 ; 而在单杠弯曲与跳高之间是正相关的从两组变量间的关系看, 单杠弯曲和跳高的训练成绩与体重腰围负相关, 与脉搏正相关表相关系数矩阵利用如下的 MALAB 程序 : clc,clea load z.txt % 原始数据存放在纯文本文件 z.txt 中 uean(z);sigstd(z); % 求均值和标准差 cocoef(z); % 求相关系数矩阵 datazscoe(z); % 数据标准化 n;; %n 是自变量的个数, 是因变量的个数 xz(:,:n);yz(:,n+:end); edata(:,:n);fdata(:,n+:end); nusize(e,);% 求样本点的个数 cgeye(n); %w 到 w 变换矩阵的初始化 fo i:n % 以下计算 w,w 和 t 的得分向量, atixe'ff'e; [vec,val]eig(atix); % 求特征值和特征向量 valdiag(val); % 提出对角线元素 [val,ind]sot(val,'descend'); w(:,i)vec(:,ind()); % 提出最大特征值对应的特征向量 -57-

8 w_sta(:,i)cgw(:,i); % 计算 w 的取值 t(:,i)ew(:,i); % 计算成分 ti 的得分 alae't(:,i)/(t(:,i)'t(:,i)); % 计算 ala_i cgcg(eye(n)-w(:,i)ala'); % 计算 w 到 w 的变换矩阵 ee-t(:,i)ala'; % 计算残差矩阵 ee; % 以下计算 ss(i) 的值 beta[t(:,:i),ones(nu,)]\f; % 求回归方程的系数 beta(end,:)[]; % 删除回归分析的常数项 cancaf-t(:,:i)beta; % 求残差矩阵 ss(i)su(su(canca.^)); % 求误差平方和 % 以下计算 ess(i) fo :nu tt(:,:i);ff; se_tt(,:);se_ff(,:); % 把舍去的第个样本点保存起来 t(,:)[];f(,:)[]; % 删除第个观测值 beta[t,ones(nu-,)]\f; % 求回归分析的系数 beta(end,:)[]; % 删除回归分析的常数项 cancase_f-se_tbeta; % 求残差向量 ess_i()su(canca.^); end ess(i)su(ess_i); if i> Q_(i)-ess(i)/ss(i-); else Q_(); end if Q_(i)<.975 fintf(' 提出的成分个数 %d',i); i; beak end end beta_z[t(:,:),ones(nu,)]\f; % 求 Y 关于 t 的回归系数 beta_z(end,:)[]; % 删除常数项 xisuw_sta(:,:)beta_z; % 求 Y 关于 X 的回归系数, 且是针对标准数据的回归系数, 每一列是一个回归方程 u_xu(:n);u_yu(n+:end); sig_xsig(:n);sig_ysig(n+:end); -58-

9 fo i: c(i)u_y(i)-u_x./sig_xsig_y(i)xisu(:,i); % 计算原始数据的回归方程的常数项 end fo i: xis(:,i)xisu(:,i)./sig_x'sig_y(i); % 计算原始数据的回归方程的系数, 每一列是一个回归方程 end sol[c;xis] % 显示回归方程的系数, 每一列是一个方程, 每一列的第一个数是常数项 save ydata x y nu xisu c xis 计算得只要提出两个成分 t,t 即可, 交叉有效性 Q.969 w 与见表, 成分 t 的得分见表 4 w 的取值表 w 与 w 的取值自变量 w w w x x x w 表 4 成分 t 的得分 No ˆt ˆt No ˆt ˆt 标准化变量 y~ k 关于成分 t 的回归模型如下 ~ y t, k,, k kt + k -59-

10 由于成分 t 可以写成原变量的标准化变量 x~ 的函数, 即有 t w ~ x w ~ x w ~ x + + 由此可得由成分 t 所建立的偏最小二乘回归模型为 ~ y k k ( w ( w k ~ x + w ~ x + w ~ x) + ~ ~ ~ k ( wx + wx + wx) + w ) ~ x + ( w + w ) ~ x + ( w + w k k k k k ) ~ x 有关,, ) 的计算结果见表 5 ( 表 5 回归系数 k 所以, 有 ~ y.778 ~ x.4989 ~ x. ~ x ~ y.85 ~ x.544 ~ x.854 ~ x ~ y.64 ~ x.559 ~ x.7 ~ x 将标准化变量 ~ y, ~ ( k xk k,, ) 分别还原成原始变量 yk, xk ( k,, ), 则回归方程为 y y y.97.67x.87x x x.477x x x.4969x x 为了更直观迅速地观察各个自变量在解释 ( k,,) 时的边际作用, 可以绘制回归系数图, 见图这个图是针对标准化数据的回归方程的从回归系数图中可以立刻观察到, 腰围变量在解释三个回归方程时起到了极为重要的作用然而, 与单杠及弯曲相比, 跳高成绩的回归方程显然不够理想, 三个自变量对它的解释能力均很低 -54- y k

11 为了考察这三个回归方程的模型精度, 我们以 y ˆ, y ) 为坐标值, 对所有的样本 ( ik ik 点绘制预测图 ŷ 是第 k 个变量, 第 i 个样本点 y ) 的预测值在这个预测图上, 如 ik 果所有点都能在图的对角线附近均匀分布, 则方程的拟合值与原值差异很小, 这个方程的拟合效果就是满意的体能训练的预测图见图 ( ik 单杠弯曲跳高 -.7 图回归系数的直方图单杠成绩预测图弯曲成绩预测图 5 5 跳高成绩预测图图体能训练预测图画直方图的 MALAB 程序为 :ba(xisu') 画体能训练的预测图的 MALAB 程序如下 : load ydata nu ceat(c,nu,); yatc+xxis; % 计算 y 的预测值 yaxax(yat); yaxax(y); -54-

12 yaxax([yax;yax]) cancayat-y; % 计算残差 sublot(,,) lot(:yax(),:yax(),yat(:,),y(:,),'') sublot(,,) lot(:yax(),:yax(),yat(:,),y(:,),'o') sublot(,,) lot(:yax(),:yax(),yat(:,),y(:,),'h') -54-

<4D F736F F F696E74202D20CAB5CFB0C1F920CAFDD6B5B1E4C1BFD7CAC1CFB5C4CDB3BCC6CDC6B6CF E707074>

<4D F736F F F696E74202D20CAB5CFB0C1F920CAFDD6B5B1E4C1BFD7CAC1CFB5C4CDB3BCC6CDC6B6CF E707074> 卫生统计学实习何平平北京大学公共卫生学院流行病与卫生统计学系 Tel: 82801619 实习六数值变量资料的统计推断 ( 三 ) 第 237~249 页一直线回归 (linear regression) ( 一 ) 定义 : 用直线方程表达 X( 自变量,independent variable) 和 Y ( 应变量,dependent variable) 之间的数量关系 Yˆ = a+