典范相关分析

实用多变数分析

4、典范相关分析

CanonicalCorrelationAnalysis1、典范变数与典范相关

若对一个个体观察了一组(p+q)个、可分成两种不同类型（或不同性质的）性状：x’=(x1,x2,…,xp),y’=(y1,y2,...,yq)如：对小麦品系（单株）考察了株高、茎粗、（剑）叶长、叶宽、穗下节间长、单株成穗数、主穗小穗数、每穗粒数、千粒重、单株产量等性状，可将前面的5个性状看成是株型性状，后面6个性状看成是穗部或产量性状，它们分别以x，y表示。

这样的性状分类事实上很常见，如株型与产量性状，产量与品质性状，淀粉与蛋白质性状，RVA特征值与淀粉蛋白质性状，长度与重量性状，价格与消费量性状等等。当考察了n个个体以后，我们往往要了解两组变数在整体上有无关系？有多大的关系？用典范相关的语言，指的是一组变数主要方向上的变异能否由另一组变数主要方向上的变异所说明？及其这种说明的程度？设x变数的线性组合为：

ξ=a’x

=a1x1+a2x2+…+apxp;y变数的线性组合为：

η=b’y=b1y1+b2y2+…+bqyqa’=(a1,a2,…,ap),b’=(b1,b2,…,bq),能否有a,b，使ξ与η之间有一个最大的相关？即，ρξη=max在求

a,b时，须满足的条件有：E(ξ)=E(η)=0,

V(ξ)=V(η)=1求最大的相关系数ρξη=max，即是根据条件：V(ξ)=V(a’x)=a’V(x)a=a’Σxxa=1V(η)=V(b’y)=b’V(y)b=b’Σyyb=1构造一个函数(G):这里,λ1,λ2称为拉格朗日乘子(Lagrangemultiplier)，这在求条件极值时经常采用。使a’Σxyb最大，亦即使G最大，常采用：因此，假定Σxx,Σyy

有逆，解上述方程组：由可见，λ1，这一拉格朗日乘子，就是所求的ξ与η的相关系数。再由：前式可写成：将（3）代入（2）式得：以及：上式可理解为：相对于的特征值和特征向量;相对于其λ2是共有的特征值。

λ2

和

b

是λ2和

a

是的特征值和特征向量。及上式中λ2和a是的特征值和特征向量;的特征值和特征向量。λ2和b是其中的λ2是共有的特征值。若令：代入式得到：这里，λ2和γ是B=T’*T=T*T’的特征值和特征向量，其中具有这种形式的矩阵B是非负定的，其特征值一定大于0，即：(λ21,λ22,

…,λ2r)λ2≥0,r=min(p,q)。

通过B=以及(B-λ2

I)=0，求解λ2的r次方程，得到r个根，将其由大到小顺序排列由此再求解r个特征向量（γ1,γ2,...,γr）。(λ21

≥λ22

≥

…≥

λ2r≥0),这ξ1，η1称为第1对典范变数，λ1为第1典范相关系数（注：不是

)。第1典范相关系数具有最大值。ξ2，η2称为第2对典范变数，λ2为第2典范相关系数，依此类推。典范变数的几点特征：1）ξi与ξj是相互独立的，ηi与ηj也是独立的，这可从特征值、特征向量（亦即主成分）的特性可知；2）ξi与ηj是相互独立的，它们间的相关系数为0；由于到此，似乎解决了所有的问题，但中的怎么算？对于若有特征值θi与特征向量li，则：该式称为矩阵Σ的谱分解（spectrumdecomposition）以上的分析是在原始数据的基础上进行的，典范相关更多地是在数据标准化的基础上进行，即每一个变数具有同等的权重。因而，典范相关分析从相关阵开始。若从相关系数矩阵开始：将用R取代Σ，用RXX取代ΣXX，用RXY取代ΣXY，用RYY取代ΣYY，与上述过程进行相同的运算，求出a,b,以及λ。2、典范变数的应用第1典范变数ξ1，η1代表了具有最大相关的两组变数最大变异度方向的线性组合：

典范变数的系数向量a,b的含义。aki和bli大小反映了xk和yl变数对第i对典范变数的作用。它们的正负代表了作用的方向，其大小代表了作用程度。绝对值越小，说明这些变数在第i对典范变数中作用微弱，aki和bli（绝对值）越大，说明这些变数在第i对典范变数中有较大的作用。根据其正负、大小，可为揭示两组变数的关系（的解释）起到一定作用。

由于两组多维变数的关系退化成两个1维变数之间的关系，这种关系可用图形的形式表示出来，使这种关系一目了然，便于对结果的分析与解释。3、典范相关系数的测验

Bartlett(1941)提出了一个测验方法，对于第1典范相关系数λ1：对于第2典范相关系数λ2：对于第j个典范相关系数λj：4、一些例子Matlab中的命令：[a,b,r,u,v,stats]=canoncorr(x,y)helpcanoncorrCANONCORRCanonicalcorrelationanalysis.[A,B]=CANONCORR(X,Y)computesthesamplecanonicalcoefficientsfortheN-by-P1andN-by-P2datamatricesXandY.XandYmusthavethesamenumberofobservations(rows)butcanhavedifferentnumbersofvariables(cols).AandBareP1-by-DandP2-by-Dmatrices,whereD=min(rank(X),rank(Y)).ThejthcolumnsofAandBcontainthecanonicalcoefficients,i.e.thelinearcombinationofvariablesmakingupthejthcanonicalvariableforXandY,respectively.ColumnsofAandBarescaledtomakeCOV(U)andCOV(V)(seebelow)theidentitymatrix.IfXorYarelessthanfullrank,CANONCORRgivesawarningandreturnszerosintherowsofAorBcorrespondingtodependentcolumnsofXorY.[A,B,R]=CANONCORR(X,Y)returnsthe1-by-DvectorRcontainingthesamplecanonicalcorrelations.ThejthelementofRisthecorrelationbetweenthejthcolumnsofUandV(seebelow).[A,B,R,U,V]=CANONCORR(X,Y)returnsthecanonicalvariables,alsoknownasscores,intheN-by-DmatricesUandV.UandVarecomputedasU=(X-repmat(mean(X),N,1))*AandV=(Y-repmat(mean(Y),N,1))*B.[A,B,R,U,V,STATS]=CANONCORR(X,Y)returnsastructurecontaininginformationrelatingtothesequenceofDnullhypothesesH0_K,thatthe(K+1)stthroughDthcorrelationsareallzero,forK=0:(D-1).STATScontainseightfields,eacha1-by-DvectorwithelementscorrespondingtovaluesofK:

Wilks:Wilks'lambda(likelihoodratio)statistic

chisq:Bartlett'sapproximatechi-squaredstatisticforH0_K,withLawley'smodification

pChisq:theright-tailsignificancelevelforCHISQF:Rao'sapproximateFstatisticforH0_KpF:theright-tailsignificancelevelforFdf1:thedegreesoffreedomforthechi-squaredstatistic,alsothenumeratordegreesoffreedomfortheFstatisticdf2:thedenominatordegreesoffreedomfortheFstatisticExample:loadcarbig;X=[DisplacementHorsepowerWeightAccelerationMPG];nans=sum(isnan(X),2)>0;[ABrUV]=canoncorr(X(~nans,1:3),X(~nans,4:5));plot(U(:,1),V(:,1),'.');xlabel('0.0025*Disp+0.020*HP-0.000025*Wgt');ylabel('-0.17*Accel+-0.092*MPG')实用多变数分析

5、因子分析

FactorAnalysis1、导言

20世纪初，KarlPearson&CharlesSpearman就开始用此方法分析一些问题，后者曾对一些学生的考试成绩进行分析。根据这些变数（成绩）及其它们之间的相互关系，用少量但本质的成分如：理解力、逻辑思维能力、记忆力等来说明学生的考试成绩。一开始由于涉及智力和精神分析方面的解释，带有一点神秘感，另加上计算繁琐的阻滞，此法一经提出，未见有多大进展。计算机的出现和广泛使用，因子分析计算量大的问题基本得到解决，再加上软件的普遍使用，因子分析的应用有所抬头，经常看到有关因子分析的文献。因子分析的基本目的在于通过变数之间的相互关系，可能的话，用少量不可观察的(unobservable)、但起根本性作用(underlying)的因子(factors)来描述。

田径运动中运动员在各项目的表现，如短跑、跳高、跳远、中长跑、铅球、铁饼、标枪的成绩，能否可以分解为几个基本因子如速度、爆发力、弹跳力和耐力所组成，即这些项目的成绩可以用速度、爆发力等因子的线性组合所表述。某些运动员若在某些能力上较弱，可以针对性地加以训练，以便收取较好效果。同样，学生的考试成绩可以用理解力、逻辑思维能力、记忆力等因子的线性组合所表达。也能对某些能力缺陷进行重点训练、培养。或，组织人事部门能掌握（洞察）有关人员的不可观察的但起根本作用的一些能力如记忆力、语言表达能力、逻辑思维能力，就能因才施教、知人善任、人尽其才。2、因子分析模式设x是一个p维变数，具有平均数μ和方差Σ，x变数可由一些未知的、不可观察的、具有本质作用的m个公共因子(F1,F2,...,Fm,m<p)再加上p个特殊因子的线性可加模型所表示：系数lij称为第i个变数在第j个因子上的负荷量(loading)，L称为因子负荷量矩阵(matrixoffactorloadings)。F1,F2,...,Fm对所有xi都有作用（如果有作用的话），称为公共因子（commonfactors），而特异因子εi只与第i个变数有关，称为特殊（异）因子（specialfactor）。

p个变数的离差x1-μ1,x2-μ2,...,xp-μp,是以p+m个变数F1,F2,...,Fm+ε1,ε2,...,εp所表达，而这m+p个起根本作用的因子又是不可观察的，因此，解出这L,F,ε是比较困难的，必须加一些条件：ⅰF与ε是独立的：cov(F,ε)=0m×pE(F)=0,V(F)=IE(ε)=0,V(ε)=diag(e1,e2,...,ep)=EⅱV(x-μ)=V(LF+ε)=LV(x-μ)L’+V(ε)+2cov(LF,ε)=LIL’+E+2Lcov(F,ε)=LL’+Eⅲcov(x,F)=cov(x-μ,F)=cov(LF+ε,F)=Lcov(F,F)+cov(ε,F)=LI+0=LⅳV(x)=Σ=LL’+E

上述模式是p个变数的由m个公共因子的线性表达。如果p个变数实际上与基本因子是有关的，但不一定成线性，则上述表述不一定正确。但线性模型是最简单的，也能说明较多的问题，因此这一线性模型应用最广泛，可视为经典模型。再回头看上式：x变数的方差可分解成两部分，其与公共因子有关，称为公因子方差(communality,

)由特异因子所引起的方差为特异因子方差(specialvariance，ei）。对一组p维x数据，能否用尽可能少的公共因子F来说明，即求出L,E来，其中m<p。当m=p时，Σ能精确地分解为LL’，但此时的E是个0矩阵，这一般也没有多大用处。只有当m<p时，因子分析才较为有用。且m越小，以越少的公共因子说明原x变数的信息，因子分析才越有用。

一般而言，Σ不能精确地分解为LL’+E，尤其是因子数m远小于p时，这样的分解难于实现，这就给因子分析的应用带来问题。另外，当m>1时，有时因子模式也会产生不确定的结果。若T是m×m的正交矩阵:F与F*具有相同的统计学性质。虽然L与L*不同，但它们也可产生相同的协方差：

其中正交矩阵T可称为旋转矩阵，它给我们带来了因子负荷量矩阵不确定的问题，但也给出了因子旋转原理，由于正交矩阵对应于x坐标轴系统的旋转，找到合适的旋转矩阵，能帮助我们对因子进行旋转，这对因子的合理解释有较好地作用。3、估计方法给出p个相关的变数：x1,x2,...,xp以及由此得到的方差协方差阵Σ，因子分析寻求解答这样的问题，是否这些因子模式（具有少数几个因子），正确地表达这些数据。一是要估计合适的因子数m，二是确定合适的因子负荷量矩阵L，三是估计特殊方差E。样本方差协方差阵S是Σ的估计，一般因子分析是从S或相关阵R开始的。如果S的非对角线元素很小，或相关矩阵R的非对角线元素接近于0，变数是不相关的，则因子分析不会有用。L=0,

在这种情况下，公共因子没有起作用，只有特异因子在起作用，每一个x变数皆是特异的。

因子分析的目的是要鉴别出少量起公共作用的因子，它们对所有x变数皆有作用。反过来，若x变数存在相关，也即说明有这样的公共因子在起作用，且相关程度越大，公共因子的作用就越明显。因子分析的参数估计有两种方法：主成分(主因子)方法和极大似然法。不管哪种方法，都能对因子旋转，便于因子的解释。有时可采用两种方法进行估计，选择其中之一。一般说来，两种方法所得的结果有一定的差异，使用者可根据情况作出选择。1）主成分方法（principalcomponentmethod）有方差协方差阵S，可对它进行谱分解：S具有特征值、特征向量对(

),则：S=LL’+0=LL’虽然上面的表示是正确的，但这并不有用。这里有p个公共因子，没有给特殊因子留下余地。我们需要地是用较少的公共因子来解释，而不是像上述那样的p个因子。解决的方法也很简单，就是当最后p-m个特征值（λi）很小时，可以将其略去。

而ε这一特殊因子也能得到。它的方差可由S-LL’的对角线元素相减而得：由主成分解，对于一个给定的因子，其因子负荷量当因子增加时不会发生改变。当m=2:对角线元素必等于0，非对角线元素就不一定等于0。选择合适的m，使残差阵（residualmatrix）接近于0。如m=1:残差阵：主对角线元素为0，如果其他元素也很小的话，可以认为是找到了一个合适的m，比较理想的情形是，前面少数几个因子对样本方差的贡献是大的，也就包括了较少的因子数，m较小。由第j个公共因子对第i个变数样本方差的贡献是由第j个公共因子对所有变数样本方差的贡献为：第1个公共因子F1，它对变数方差总的贡献是：第2个公共因子F2，它对变数方差总的贡献是：以此类推。第j个公共因子对于总方差的贡献：（由S阵进行的因子分析）或如同主成分分析，由原始数据得到的因子分析结果与由标准化数据得到的结果不一定相同，且往往很不相同，两者不成比例。应用上是否采用标准化，应视情况而定，掌握的原则跟主成分分析相同。

（由相关阵进行的因子分析）。上述的第j因子对于总方差贡献的比例，常作为一种启发手段，以决定合适的因子数目m。在因子模式中因子数目逐步增加，直至样本方差的一个恰当的成分为此模式所说明，便是较合适的因子数。一般而言，当：将前面m个因子入选。在以R为基础的因子分析时，常把λm>1的前m个因子入选。>=80~85%时2）主因子方法——一种改进的方法我们也可以用另外的方法来描述R的析因过程。当然其过程也适用于方差协方差阵S。如果因子分析模式：则m个公共因子应该能说明R的非对角线元素，以及主对角线元素：如果特异因子从主对角线减去：假如一开始，我们就能估计则：R的对角线元素都经这样的相减：这一相减的相关阵应由m个公共因子所说明：选取m个最大特征值及其对应的特征向量组成：由R*进行特征值、特征向量分解得到（）此时，m个公共因子的贡献率需重新计算：特异因子方差也应重新计算：这一数值很可能与不同。与前面的相近时为止。

应将此过程继续进行，直至经t轮迭代后：以主成分的思想，我们可以以（λ1,λ2,...,λm）的大小来决定m的个数，如取前面m个λi>1的主因子入选，后面的p-m个λi<1的主因子剔除。问题的复杂性还在于此时R*不再是正定阵，它的特征值很可能会出现负值，为因子数的决定带来困难。不过，在大多数情况下，出现负数或绝对值大的负数的情况较少见。比较常见的是前面m个大的特征值的，后面的特征值较小，或者为0，或者小的负数，这种情况对因子分析的影响不算很大。若特征值出现很大的负数，那肯定是不合适的，应该重新选取初值并迭代运算。有时（很多情况下），的选取也有讲究，对结果或迭代次数有较大影响。一般而言：指的是不能被公共因子说明的部分。很多人在1/4~1/3之间选取。也可采用对R阵选求算其逆阵（C），利用其逆阵对角线元素：3）极大似然法假定F与ε是正态分布的（这一般也是成立的），L与E可用极大似然法进行估计。x的概率密度函数：

构建似然函数：将Σ=LL’+E代入上式先给一个L,E的初值，采用迭代估计，得到所需的m,L,E。这一过程比较复杂，可用计算机软件完成。例：Linden(1977)对二战后奥林匹克十项全能(decathlon)比赛的成绩进行了因子分析。7届139人得到了160组数据（有一些人是重复参加比赛的）。十项的成绩进行了标准化处理，因子分析是在相关系数阵的基础上进行的。对相关阵分别进行了主成分和极大似然法分解。

100mljumpshotpthjump400m110hdiscusplvaultjavelin1500m100m10.590.350.340.630.400.280.20.11-0.07ljump

10.420.510.490.520.310.360.210.09shotput

10.380.190.360.730.240.44-0.08hjump

10.290.460.270.390.170.18400m

10.340.170.230.130.39110h

10.320.330.180discus

10.240.34-0.02Plvault

10.240.17javelin

101500m

1相关系数表方法principalcomponentsmaximumlikelihood项目factorlodingseifactorloadingseiF1F2F3F4F1F2F3F4100m.691.217-.520-.206.16-.090.341.830-.169.16Ljump.789.184-.193.092.30.065.433.595.275.38Shotpt.702-.535.047-.175.19-.139.990.000.000.00Hjump.694.134.139.396.35.156.406.336.445.50400m.620.551-.084-.419.13.376.245.671-.137.33110h.687.042-.161.345.38-.021.361.425.388.54Discus.621-.521.109-.234.28-.063.728.030.019.46Plvault.538.087.411.440.34.155.264.229.394.70Javelin.434-.439.372-.235.43-.026.441-.100.098.801500m.147.596.658-.279.11.998.059.000.000.00Cumulative0.380.530.640.73

0.120.370.550.61

因子负荷量表

由主成分因子分析，前面4个特征值分别为3.78，1.52，1.11，0.91，前面4个主成分能说明总变异的73%，可由4个因子所说明。极大似然法说明的成份为61%，比4个主成分因子说明的比例略小。

在这一例中，两种方法得到的结果有较大差别。对于主成分来说，所有项目除了1500m以外，都在第1因子上具有较大的负荷，这可以说是generalathleticability(基本运动能力)。因子2主要说明跑跳能力(腿力)对比于投掷力(臂力)。因子3是速度对比于耐力，其中对撑杆跳也有较高的负荷。因子4在撑杆跳、400m和跳高上有较高负荷，与弹跳力有一定关系。

对于极大似然法而言，1500m是仅有的对第1因子具有大负荷的项目，可能是一个耐力因子。因子2主要在铅球和铁饼上有较大负荷，可视为是力量（臂力）因子。因子3可能主要是速度因子，与速度无关的运动在这一因子上几乎没有负荷量。第4因子也难于直观地看出，基本和腿部力量有关。

特异方差在各个项目上占有一定的比重，如主成分方法仍有27%，本还可以再分解出一些公共因子，但在一些项目上的特异方差已经很小了，公共因子数目的增加已不能说明更多的总变异。在一些项目上，特异方差仍较大，说明这些项目不能为公共因子的线性组合所说明，这也说明这些项目如标枪、撑杆跳等需要较高的技术。

极大似然法的结果也较相似。虽然它在说明方差的成分上稍小，但在一些项目上的特异方差趋于0。这些项目基本能由前面4个因子的线性组合所说明，甚至有些项目只需前面两个因子便能完全说明。而在标枪、撑杆跳等项目上，特异方差较大，单靠臂力、腿力、耐力是不可能有好成绩的，它们需要特有的技术。残差阵：R-LL’-E项目100mljumpshotpthjump400m110hdiscusplvaultjave1500100m0.000.000.012.000-.012.004.000-.018.000Ljump-.0750.000.000.002-.002.006-.025-.009.000Shotput-.030-.0100.000.000-.000-.000-.000-.000.000Hjump-.001-.056.0420-.033.001-.034.006-.045.000400m-.047-.077-.020-.0240.028-.002.008.052.000110h-.096-.092-.032-.122.0220.036-.012-.013.000Discus-.027-.011-.031-.001-.017.0140.043.016.000Plvault.114-.042-.034-.215.067-.129.0090.091.000Javelin.051.042-.158-.022.036.041-.254-.0050.0001500m-.016.017.056.200-.091.076.062-.109-.1120注：主对角线左下为主成分法，主对角线右上为极大似然法。4、因子旋转(factorrotation)上面提到了如果得到一个因子负荷矩阵后，若进行正交变换：对于L的这样一种变换，并不改变因子分析的性质，协方差分解仍保持原有的结构。特异方差ei不变，共因子方差h2不改变，残差阵R-LL-E，S-LL-E也没有改变。那进行变换有什么意义呢？主要在于由主成分或极大似然法得到的因子负荷量矩阵，对因子的解释不够明确，结构不那么简单，不那么清晰易懂。通过正交旋转，则可以获得一个简单、易于解释的结构。这一原理如同调节显微镜的焦距，使细节易于看清一样：Therationaleisverymuchakintosharpeningthefocusofmicroscopeinordertoseethedetailmoreclearly.

较为理想的情形是这样的负荷量矩阵，每一变数只在一个因子上有较大的负荷，而在其它因子上的负荷较小。而实际获得的结果并非如此，因而可以采用旋转的方法帮助这个结果的实现。项目principalcomponents