R软件应用多元分析(II)

9.应用多元分析(II)安徽师范大学数学计算机科学学院丁新涛现在是1页\一共有75页\编辑于星期一9.1主成分分析9.1.1总体主成分主成分的定义与导出假定你是一个公司的财务经理，掌握了公司的所有数据，比如固定资产、流动资金、每一笔借贷的数额和期限、各种税费、工资支出、原料消耗、产值、利润、折旧、职工人数、职工的分工和教育程度等等。如果让你向上面介绍公司状况，你能够把这些指标和数字都原封不动地摆出去吗？当然不能。你必须要把各个方面作出高度概括，用一两个指标简单明了地把情况说清楚。本章介绍两种把变量维数降低以便于描述、理解和分析的方法：主成分分析（principalcomponentanalysis）和因子分析（factoranalysis）。实际上主成分分析可以说是因子分析的一个特例。现在是2页\一共有75页\编辑于星期一例子:成绩数据100个学生的数学、物理、化学、语文、历史、英语的成绩如下表（部分）。目前的问题是，能不能把这个数据的6个变量用一两个综合变量来表示呢？现在是3页\一共有75页\编辑于星期一主成分分析例中的的数据点是六维的；也就是说，每个观测值是6维空间中的一个点。我们希望把6维空间用低维空间表示。以二维为例,如果这些数据形成一个椭圆形状的点阵(比如二维正态分布)在短轴Z2方向上，数据变化很少(方差小);进一步,短轴如果退化成一点，则长轴Z1即可解释这些点的变化;这样,由二维到一维的降维就自然完成了.这相当于在平面上做一个坐标变换.Z1,Z2是X1,X2的特殊线性组合.X1X2Z1Z2现在是4页\一共有75页\编辑于星期一推广到p维情况:对于多维变量的情况和二维类似，也有高维的椭球，只不过无法直观地看见罢了。首先把高维椭球的主轴找出来，再用代表大多数数据信息的最长(方差大)的几个轴作为新变量；这样，主成分分析就基本完成了。注意，和二维情况类似，高维椭球的主轴也是互相垂直的。这些互相正交的新变量是原先变量的线性组合，叫做主成分(principalcomponent)。正如二维椭圆有两个主轴，三维椭球有三个主轴一样，有几个变量，就有几个主成分。选择越少的主成分，降维就越好。什么是标准呢？那就是这些被选的主成分所代表的主轴的长度之和占了主轴长度总和的大部分。有些文献建议，所选的主轴总长度占所有主轴长度之和的大约85%即可，其实，这只是一个大体的说法；具体选几个，要看实际情况而定。现在是5页\一共有75页\编辑于星期一定义:设X是p维随机变量,μ=E(X),Σ=var(X)考虑线性变换(p维坐标主轴→p维椭球主轴):我们希望Z1的方差达到最大(Z1是椭球主轴中最长的一个主轴),则a1满足:max=λmax，a1是Σ最大特征值(λ1)的特征向量.Z1=a1TX称为第一主成分.现在是6页\一共有75页\编辑于星期一特征方程，λ是Σ的特征值a是对应的特征向量现在是7页\一共有75页\编辑于星期一类似地,希望Z2的方差达到最大,为保证椭球的主轴也是互相垂直的,a2与a1正交,即:cov(Z1,Z2)=a1T

Σa2=0,类似地,a2是Σ的第二大特征值(设为λ2),Z2=a2TX称为第二主成分.一般地,即为所求线性变换矩阵A.Zi=aiTX即为第i个主成分.,Q是正交阵.现在是8页\一共有75页\编辑于星期一2.主成分的性质主成分的均值和协方差阵.主成分的总方差.主成分分析是把p个原始变量X1,X2,…,Xp的总方差分解成p个不相关变量(cov(Z1,Z2)=0,)Z1,Z2,…,Zp的方差之和.总方差中第i主成分Zi的比例称为主成分Zi的贡献率.总方差中前m个主成分Zi的贡献率之和称Z1,Z2,…,

Zm的累积贡献率.Xj与Zi之间的相关系数.现在是9页\一共有75页\编辑于星期一2.主成分的性质m个主成分对原始变量的贡献率.总方差中前m个主成分Zi的贡献率之和称Z1,Z2,…,

Zm对X1,X2,…,Xp的累积贡献率.Z1,Z2,…,

Zm对Xj的累积贡献率:原始变量对主成分的影响.qji称为第i个主成分在第j个原始变量Xj上的载荷,它度量了Xj对Zi的重要程度.现在是10页\一共有75页\编辑于星期一3.从相关矩阵出发求主成分的方差矩阵是X的相关矩阵R.p个主成分为:相关矩阵R的主成分性质:E(Z*)=0;var(Z*)=Λ*,其中Λ*=diag(λ1*,λ2*,…,λp*).变量Xj*与主成分Zi*之间的相关系数为:主成分Z1*,Z2*,…,Zm*对Xj*的贡献率为:.设是相关矩阵R的p个特征值，是相应的单位特征向量。现在是11页\一共有75页\编辑于星期一9.1.2样本主成分总体的参数在实际问题中，通常是未知的，则通过样本来估计。样本变量样本方差矩阵:样本相关矩阵R:相关记号:现在是12页\一共有75页\编辑于星期一1.从S出发求主成分S的特征值:λ1≧λ2≧…≧λp对应的特征向量(标准化):a1,a2,…,ap,ai称为主轴向量(简称主轴);则第i个主成分zi=aiTX=aiT(x1,x2,…,xp)T令:z=(z1,z2,…,zp)T=(a1,a2,…,ap)Tx=QTx则样本主成分:z(k)=QTX(k)每一个样本的主成分样本投影到新坐标系下的坐标X的第2个主成分实际问题中,经常将样本进行中心化.协方差矩阵不变现在是13页\一共有75页\编辑于星期一样本主成分性质:样本主成分的总方差等于原变量样本的总方差Xj与Zi的样本相关系数.现在是14页\一共有75页\编辑于星期一2.从R出发求主成分样本相关矩阵R的特征值:λ1*≧λ2*≧…≧λp*对应的特征向量(标准化):a1*,a2*,…,ap*,z(k)*=QTX(k)*性质:略现在是15页\一共有75页\编辑于星期一9.1.3相关R函数以及实例princomp{stats}Description:princompperformsaprincipalcomponentsanalysisonthegivennumericdatamatrixandreturnstheresultsasanobjectofclassprincomp.Usageprincomp(x,...)##S3methodforclass'formula':princomp(formula,data=NULL,subset,na.action,...)Formula:aformulawithnoresponsevariable,referringonlytonumericvariables.##DefaultS3method:princomp(x,cor=FALSE,scores=TRUE,covmat=NULL,subset=rep(TRUE,nrow(as.matrix(x))),...)X:anumericmatrixordataframewhichprovidesthedatafortheprincipalcomponentsanalysis.Cor:alogicalvalueindicatingwhetherthecalculationshouldusethecorrelationmatrixorthecovariancematrix.(Thecorrelationmatrixcanonlybeusediftherearenoconstantvariables.)现在是16页\一共有75页\编辑于星期一3-63.Loadings():Extractorprintloadingsinfactoranalysis(orprincipalcomponentsanalysis).Usage:loadings(x)x:anobjectofclass"factanal"or"princomp"ortheloadingscomponentofsuchanobject.predict():predictisagenericfunctionforpredictionsfromtheresultsofvariousmodelfittingfunctions.Thefunctioninvokesparticularmethodswhichdependontheclassofthefirstargument.Usage:predict(object,...)screeplot():screeplot.defaultplotsthevariancesagainstthenumberoftheprincipalcomponent.Thisisalsotheplotmethodforclasses"princomp"and"prcomp".Usage:screeplot(x,npcs=min(10,length(x$sdev)),type=c("barplot","lines"),main=deparse(substitute(x)),...)npcs:thenumberofcomponentstobeplotted.biplot():画出数据关于主成分的散点图和原坐标在主成分下的方向.Biplot(x,choices=1:2,scale=1,pc.biplot=FALSE,…)Choices:是选择的主成分,默认是第1，2主成分.现在是17页\一共有75页\编辑于星期一7.实例序号X1X2X3X4114841727821393471763160497786414936677951594580866142316676715343768381504377799151427780101393168741114029647412161477884131584978831414033677715137316673161523573791714947827918145357077191604774872015644788521151427382221473873782315739688024147306575251574880882615136748027144366876281413067762913932687330148387078现在是18页\一共有75页\编辑于星期一R实现:student=read.table('dataexample901.txt')student.pr=princomp(student,cor=TRUE)summary(student.pr,loadings=TRUE)Importanceofcomponents:Comp.1Comp.2Comp.3Comp.4Standarddeviation1.88178050.559806360.281795940.25711844ProportionofVariance0.88527450.078345790.019852240.01652747CumulativeProportion0.88527450.963620290.983472531.00000000Loadings:Comp.1Comp.2Comp.3Comp.4x1-0.4970.543-0.4500.506x2-0.515-0.210-0.462-0.691x3-0.481-0.7250.1750.461x4-0.5070.3680.744-0.232#采用从R出发求主成分主成分贡献率:Q由于前2个主轴的贡献率达到96.4%,这个例子的主成分可以认为是comp.1(Z1*),Comp.2(Z2*)现在是19页\一共有75页\编辑于星期一预测:求Z*predict(student.pr)

Comp.1Comp.2Comp.3Comp.410.06990950-0.23813701-0.35509248-0.26612013921.59526340-0.718473990.32813232-0.1180566463-2.847931510.38956679-0.09731731-0.27948248740.759969880.80604335-0.04945722-0.1629492985-2.739667770.017180870.360126150.35865304462.105831680.322843930.18600422-0.0364560847-1.42105591-0.060531650.21093321-0.0442230928-0.82583977-0.78102576-0.275577980.0572885729-0.93464402-0.58469242-0.088141360.181037746102.36463820-0.365321990.088404760.045520127112.837419160.348758410.03310423-0.03114693012-2.608512240.21278728-0.333980370.21015757413-2.44253342-0.16769496-0.46918095-0.162987830141.866306690.050213840.37720280-0.358821916152.81347421-0.31790107-0.03291329-0.222035112160.063929830.207184480.043343400.70353362417-1.55561022-1.70439674-0.331264060.007551879181.07392251-0.067634180.022836480.04860668019-2.521742120.972743010.12164633-0.39066799120-2.140723770.022178810.374109720.12954896021-0.796244220.163078870.12781270-0.294140762220.28708321-0.35744666-0.039621160.08099198923-0.251510751.25555188-0.556173250.109068939242.057060320.78894494-0.265521090.38808864325-3.08596855-0.057753180.62110421-0.21893961226-0.163675550.043179320.244818500.560248997271.372650530.02220972-0.23378320-0.257399715282.160977780.137332330.355897390.093123683292.40434827-0.48613137-0.16154441-0.007914021300.502874680.14734317-0.20590831-0.122078819Comp1(Z1*)对应的系数符号都是负号,反映了学生身材的魁梧程度,身材高大的学生,他的4个部分值都比较大,则comp1的值就比较小,反之,身材’矮小的学生comp1的值比较大;Comp2(Z2*)是X1,X4与X2,X3的差,即:纵向高度与横向围度的差,所以,”细高”的同学comp2值越大,反之,该值越小说明样本越”矮胖”.蓝色样本身材魁梧;粉色样本身材瘦小;红色样本身材细高;绿色样本身材矮胖;现在是20页\一共有75页\编辑于星期一结果图示化:screeplot(student.pr,type='lines')碎石图:纵坐标:(1.88178050.559806360.281795940.25711844)2通过碎石图可以直观的观察出样本在变换主轴方向的波动程度biplot(student.pr,choice=1:2)横坐标是Comp1纵坐标是comp2矮胖细长魁梧瘦小现在是21页\一共有75页\编辑于星期一9.1.4主成分分析的应用1.主成分分类X1X2X3X4X5X6X7X8X9X10X11X12X13X14X15X16X11X20.791X30.360.311X40.960.740.381X50.890.580.310.91X60.790.580.30.780.791X70.760.550.350.750.740.731X80.260.190.580.250.250.180.241X90.210.070.280.20.180.180.29-0.041X100.260.160.330.220.230.230.250.49-0.341X110.070.210.380.08-0.0200.10.44-0.160.231X120.520.410.350.530.480.380.440.3-0.050.50.241X130.770.470.410.790.790.690.670.320.230.310.10.621X140.250.170.640.270.270.140.160.510.210.150.310.170.261X150.510.350.580.570.510.260.380.510.150.290.280.410.50.631X160.210.160.510.260.2300.120.380.180.140.310.180.240.50.651现在是22页\一共有75页\编辑于星期一R实现:x=scan('dataexample902.txt')names=c('x1','x2','x3','x4','x5','x6','x7','x8','x9','x10','x11','x12','x13','x14','x15','x16')r=matrix(0,nrow=16,ncol=16,dimnames=list(names,names))for(iin1:16){for(jin1:i){r[i,j]=x[(i-1)*i/2+j]r[j,i]=r[i,j]}}pr=princomp(covmat=r)load=loadings(pr)plot(load[,1:2])text(load[,1],load[,2],adj=c(-0.4,0.3))#rij=x[1+2+…+(i-1)+j]长类围类体型指标>summary(pr)Importanceofcomponents:Comp.1Comp.2Comp.3Standarddeviation2.6521.61679711.2775386ProportionofVariance0.43977980.16337700.1020066CumulativeProportion0.43977980.60315690.7051634前7个主成分占贡献率87%现在是23页\一共有75页\编辑于星期一主成分聚类:Comp.1Comp.2x1-0.341770990.20040027x2-0.264991630.14320222x3-0.23415224-0.32862480x4-0.344232670.18112449x5-0.326117710.19965028x6-0.285913570.26980664x7-0.295261390.19214958x8-0.18927312-0.37026699x9-0.084792910.06747164x10-0.15429508-0.17424610x11-0.09835526-0.34784952x12-0.24254582-0.01766472x13-0.317158240.11191444x14-0.18011330-0.37135294x15-0.26635929-0.27122509x16-0.15833266-0.3628239316个样本,每个样本是2维向量;可以用样本的距离对其实施聚类;write(load[,1:2],'dataexample902load12.txt')x=scan('dataexample902load12.txt')x=as.matrix(x);dim(x)=c(16,2)d=dist(scale(x));hc=hclust(d)plot(hc,hang=-1);re3=rect.hclust(hc,k=3)现在是24页\一共有75页\编辑于星期一用相关矩阵直接聚类d=as.dist(1-r);hc=hclust(d)plot(hc,hang=-1);re3=rect.hclust(hc,k=3)长类围类体型指标长类主成分聚类图现在是25页\一共有75页\编辑于星期一2.主成分回归例9.3考虑进口总额Y与3个自变量:国内总产值X1,存储量X2和总消费量X3(单位为10亿法郎)之间的关系.现收集了1949-1959年共11年数据,如表,试对此数据做经典回归分析和主成分回归分析.序号X1X2X3Y1149.34.2108.115.92161.24.1114.816.43171.53.1123.2194175.53.1126.919.15180.81.1132.118.86190.72.2137.720.47202.12.114622.78212.45.6154.126.59226.15162.328.110231.95.1164.327.6112390.7167.626.3现在是26页\一共有75页\编辑于星期一R实现:conomy=read.table('dataexample903.txt')lm.sol=lm(Y~X1+X2+X3,data=conomy)#普通线性回归summary(lm.sol)Coefficients:EstimatePr(>|t|)(Intercept)-10.127996.9e-05***X1-0.051400.488344X20.586950.000444***X30.286850.026277*Residualstandarderror:0.4889MultipleR-squared:0.9919AdjustedR-squared:0.9884F-statistic:285.6p-value:1.112e-07Y=-10.12799-0.0514X1+0.58695X2+0.28685X3lm.sol=lm(Y~X1+X2,data=conomy)Coefficients:EstimatePr(>|t|)(Intercept)-8.440140.000370***X10.145313.14e-08***X20.622480.001243**Residualstandarderror:0.6667MultipleR-squared:0.9828,AdjustedR-squared:0.9785F-statistic:228.3p-value:8.796e-08Y=-8.44014+0.14531X1+0.62248X2现在是27页\一共有75页\编辑于星期一残差图:Y=-8.44014+0.14531X1+0.62248X2H0:同方差,不是小概率事件接受HO:同方差;异方差检验:H0:残差不相关,不是小概率事件接受HO:残差不相关;不相关检验:现在是28页\一共有75页\编辑于星期一分析:λ3≈0X1,X2共线性Y=-8.44014+0.14531X1+0.62248X2后果很严重:回归系数不稳定总之:回归结果仍不可信(多重共线性)Y=-10.12799-0.0514X1+0.58695X2+0.28685X3回归系数不能通过显著性t检验,所以该回归结果不可信.Y=-8.44014+0.14531X1+0.62248X2回归诊断：误差项是否满足独立性、等方差性、正态性；选择线性模型是否合适；是否存在异常样本；回归分析的结果是否对某些样本的依赖过重，即回归模型是否具备稳定性；变量之间是否存在高度相关，即是否有多重共线性问题存在；(不通过)回归残差通过白噪声检验怎么办?现在是29页\一共有75页\编辑于星期一多重共线性的处理方法:删去模型中次要的或可替代的解释变量Y~x2不能通过系数显著性t检验;Y~x1:y=-6.558101+0.146199X1,

R-squared 0.931761模型的检验都能通过.是一个”能用的”模型.实际上,一个更好的模型是y~x2+x3:y=-9.74+0.596X1+0.212X3R-squared 0.991277R2比y~x1的要高一些.是一个”比较好的”模型.现在是30页\一共有75页\编辑于星期一主成分法:conomy.pr=princomp(~X1+X2+X3,data=conomy,cor=T)summary(conomy.pr,loadings=TRUE)Importanceofcomponents:Comp.1Comp.2Comp.3Standarddeviation1.4139150.99907670.0518737839ProportionofVariance0.6663850.33271810.0008969632CumulativeProportion0.6663850.99910301.0000000000Loadings:Comp.1Comp.2Comp.3X10.7060.707X2-0.999

X30.707-0.707λ3=0.05187378392≈0X1,X2共线性#Z1,Z2的贡献率已达99.9%小于0.05的数没有显示出来可以通过conomy.pr$loadings[][[2]]等查看0.706330.0356890.7069820.043501-0.999030.0069710.7065440.02583-0.7072现在是31页\一共有75页\编辑于星期一X1starx2starx3starZ1star-1.509720.545705-1.53319-2.12589-1.113050.485071-1.20848-1.61893-0.76971-0.12127-0.8014-1.11517-0.63637-0.12127-0.62209-0.8943-0.4597-1.33395-0.37008-0.64421-0.1297-0.66697-0.09869-0.190350.250307-0.727610.3035530.3596220.5936461.394580.6961010.9718021.050321.0307761.0934961.5593161.2436571.091411.1904221.7669951.480327-1.576481.3503491.931103pre=predict(conomy.pr)conomy$Z1=pre[,1]conomy$Z2=pre[,2]X*=-0.70720.025830.7065440.006971-0.999030.0435010.7069820.0356890.70633Q=(q1,q2,q3)=-2.2296493-1.6979452-1.1695976-0.9379462-0.6756511-0.19964230.37717461.01923441.63542431.85324012.0253583Z1*=X*q1lm.sol=lm(Y~Z1+Z2,data=conomy)summary(lm.sol)Coefficients:EstimatePr(>|t|)(Intercept)21.89091.21e-14***Z12.98926.02e-09***Z2-0.82880.00106**

Residualstandarderror:0.55MultipleR-squared:0.9883AdjustedR-squared:0.9853F-statistic:337.2p-value:

1.888e-08

#回归结果比较理想.Y=21.89+2.99Z1*-0.83Z2*现在是32页\一共有75页\编辑于星期一将Z*表达成X:beta=coef(lm.sol)a=loadings(conomy.pr)x.bar=conomy.pr$centerx.sd=conomy.pr$scalecoef=(beta[2]*a[,1]+beta[3]*a[,2])/x.sdbeta0=beta[1]-sum(x.bar*coef)c(beta0,coef)(Intercept)Z1Z221.89090912.9891518-0.8287678

(Intercept)X1X2X3-9.130107820.072779810.609220120.10625939主成分法回归方程为:Y=-9.13010782+0.07277981X1+0.60922012X2+0.10625939X3现在是33页\一共有75页\编辑于星期一9.2.2因子分析模型一般地，设X=(x1,x2,…,xp)T为可观测的随机变量，且有f=(f1,f2,…,fm)’为公共（共性）因子(commonfactor)，简称因子(factor);e=(e1,e2,…,ep)’为特殊因子（specificfactor）f和e均为不可直接观测的随机变量;μ=(μ1,μ2,…,μp)’为随机变量x的总体均值;A=(aij)p*m为因子负荷（载荷）（factorloading）矩阵因子分析的模型为:假设通常先对x作标准化处理，使标准化得到的新变量均值为零，方差为１．现在是34页\一共有75页\编辑于星期一如果再满足（４）fi与fj相互独立（i≠j），则称该因子模型为正交因子模型。正交因子模型具有如下特性：x的方差可表示为设(1)hi2是m个公共因子对第i个变量的贡献，称为第i个共同度（communality）或共性方差，公因子方差（commonvariance）(2)δi称为特殊方差（specificvariance），是不能由公共因子解释的部分(3)因子载荷（负荷）aij是随机变量xi与公共因子fj的相关系数，系数aij是用来度量Xi可由f1,f2,…,fm线性组合表示的程度。现在是35页\一共有75页\编辑于星期一设称gj2为公共因子fj对x的“贡献”，是衡量公共因子fj重要性的一个指标。现在是36页\一共有75页\编辑于星期一三、因子分析的步骤输入原始数据xn*p，计算样本均值和方差，进行标准化计算（处理）；求样本相关系数矩阵R=(rij)p*p；求相关系数矩阵的特征根λi

(λ1,λ2,…,λp>0)和相应的标准正交的特征向量li；确定公共因子数；计算公共因子的共性方差hi2;对载荷矩阵进行旋转，以求能更好地解释公共因子；对公共因子作出专业性的解释。现在是37页\一共有75页\编辑于星期一9.2.3参数估计(提取因子的方法)1.主成分法（principalcomponentfactor）

每一个公共因子的载荷系数之平方和等于对应的特征根，即该公共因子的方差。估计因子荷载矩阵A=(aij)p*n和特殊方差矩阵D特殊方差矩阵D：当总体的方差矩阵未知时，用样本协方差矩阵S代替，sii是S对角线上元现在是38页\一共有75页\编辑于星期一factor.analy1factor.analy1<-function(S,m){

p<-nrow(S);diag_S<-diag(S);sum_rank<-sum(diag_S)

rowname<-paste("X",1:p,sep="")colname<-paste("Factor",1:m,sep="")A<-matrix(0,nrow=p,ncol=m,dimnames=list(rowname,colname))eig<-eigen(S)for(iin1:m)A[,i]<-sqrt(eig$values[i])*eig$vectors[,i]h<-diag(A%*%t(A))rowname<-c("SSloadings","ProportionVar","CumulativeVar")B<-matrix(0,nrow=3,ncol=m,dimnames=list(rowname,colname))for(iin1:m){B[1,i]<-sum(A[,i]^2)B[2,i]<-B[1,i]/sum_rankB[3,i]<-sum(B[1,1:i])/sum_rank}method<-c("PrincipalComponentMethod")list(method=method,loadings=A,var=cbind(common=h,spcific=diag_S-h),B=B)}为方便，设定中间变量A=(Factor1,…,FactormX1………0………

………0………Xp)p*m现在是39页\一共有75页\编辑于星期一R实例：例9.7对55个国家和地区的男子径赛记录作统计，每位运动员记录8项指标：100m跑(X1)、200m跑(X2)、400m跑(X3)、800m跑(X4)、1500m跑(X5)、5000m跑(X6)、10000m跑(X7)、马拉松(X8)。8项指标的相关矩阵R如表9.4所示。取m=2，用主成分法估计因子载荷和共性方差等指标。现在是40页\一共有75页\编辑于星期一R实现：x=scan('dataexample907.txt')names=c('x1','x2','x3','x4','x5','x6','x7','x8')r=matrix(0,nrow=8,ncol=8,dimnames=list(names,names))for(iin1:8){for(jin1:i){r[i,j]=x[(i-1)*i/2+j]r[j,i]=r[i,j]}}source('factor.analy1.r')fa=factor.analy1(r,m=2)Fa$method[1]"PrincipalComponentMethod"$loadingsFactor1Factor2X1-0.8171562-0.53123478X2-0.8674064-0.43231147X3-0.9151503-0.23258703X4-0.9487239-0.01184340X5-0.95935870.13153096X6-0.93764640.29276177X7-0.94395700.28715151X8-0.87992530.41074922$varcommonspcificX10.94995460.05004538X20.93928700.06071301X30.89159680.10840319X40.90021730.09978270X50.93766960.06233044X60.96489030.03510972X70.97351080.02648920X80.94298340.05701655$BFactor1Factor2SSloadings6.62258840.8775213ProportionVar0.82782360.1096902CumulativeVar0.82782360.9375137公共因子fj对X1,…,Xp的总方差贡献现在是41页\一共有75页\编辑于星期一9.2.4因子旋转目的：使因子负荷两极分化，要么接近于0，要么接近于1。常用的旋转方法：（1）方差最大正交旋转（varimaxorthogonalrotation）基本思想：使公共因子的相对负荷（lij/hi2）的方差之和最大，且保持原公共因子的正交性和公共方差总和不变。可使每个因子上的具有最大载荷的变量数最小，因此可以简化对因子的解释。（2）斜交旋转（obliquerotation）因子斜交旋转后，各因子负荷发生了较大变化，出现了两极分化。各因子间不再相互独立，而彼此相关。各因子对各变量的贡献的总和也发生了改变。适用于大数据集的因子分析。现在是42页\一共有75页\编辑于星期一1.理论依据设因子模型：令：(是任意m阶正交矩阵).即：对任一正交矩阵Γ，Z=ΓTF也是公因子向量。相应的AΓ是公因子Z的因子载荷矩阵。为此，在因子分析的实际计算中，当求得初始因子载荷矩阵A后，反复右乘正交矩阵Γ，使得AΓ具有更明显的意义，这种变换载荷矩阵的方法，称为因子轴的正交旋转。现在是43页\一共有75页\编辑于星期一2.因子载荷方差设因子模型A=(aij)p*m为因子载荷矩阵，为变量Xi的共同度。A的每一列(因子载荷向量)数值越分散，相应的因子载荷向量的方差越大。令：第j列p个数据的方差因子载荷矩阵A的方差为：Vj越大，A的第j列值越分散，Vj趋于1或0，称公因子Fj具有简单化结构。V越大越好。现在是44页\一共有75页\编辑于星期一3.方差最大的正交旋转所谓最大方差旋转法就是选择正交矩阵，使得矩阵所有个列元素平方的相对方差之和达到最大。当m=2时，设已求出的因子载荷矩阵为

现选取正交变换矩阵Γ进行因子旋转，Γ可以表示为:这里θ是坐标平面上因子轴按顺时针方向旋转的角度，只要求出θ，也就求出了Γ.现在是45页\一共有75页\编辑于星期一现在是46页\一共有75页\编辑于星期一m>2:当m>2时，我们可以逐次对每两个公共因子和进行上述旋转。对公因子Fl和Fk进行旋转，就是对A的第l和k两列进行正交变换，使这两列元素平方的相对方差之和达到最大，而其余各列不变，其正交变换矩阵为现在是47页\一共有75页\编辑于星期一其中θ是因子轴Fi和Fj的旋转角度，矩阵中其余位置上的元素全为0。m个公共因子两两配对旋转共需要进行m(m-1)/2次，称其为完成了第一次旋转，并记第一轮旋转后的因子载荷矩阵为A(1)。然后再重新开始，进行第二轮的Cm2次配对旋转，新的因子载荷矩阵记为A(2)。这样可以得到一系列的因子载荷矩阵为现在是48页\一共有75页\编辑于星期一9.2.5因子分析的计算函数Factanal：factanal(x,factors,data=NULL,covmat=NULL,n.obs=NA,subset,na.action,start=NULL,scores=c("none","regression","Bartlett"),rotation="varimax",control=NULL,...)ArgumentsX:Aformulaoranumericmatrixoranobjectthatcanbecoercedtoanumericmatrix.Factors:Thenumberoffactorstobefitted.Data:Anoptionaldataframe(orsimilar:seemodel.frame),usedonlyifxisaformula.Bydefaultthevariablesaretakenfromenvironment(formula).Covmat:Acovariancematrix,oracovariancelistasreturnedbycov.wt.Ofcourse,correlationmatricesarecovariancematrices.Scores：Typeofscorestoproduce,ifany.Thedefaultisnone,“regression”givesThompson‘sscores,“Bartlett”givenBartlett’sweightedleast-squaresscores.Partialmatchingallowsthesenamestobeabbreviated.Rotation：character."none"orthenameofafunctiontobeusedtorotatethefactors:itwillbecalledwithfirstargumenttheloadingsmatrix,andshouldreturnalistwithcomponentloadingsgivingtherotatedloadings,orjusttherotatedloadings.现在是49页\一共有75页\编辑于星期一例9.11取m=2，用factanal()函数估计例9.7因子载荷和共性方差等指标，参数选择方差最大。x=scan('dataexample907.txt')names=c('x1','x2','x3','x4','x5','x6','x7','x8')r=matrix(0,nrow=8,ncol=8,dimnames=list(names,names))for(iin1:8){for(jin1:i){r[i,j]=x[(i-1)*i/2+j]r[j,i]=r[i,j]}}fa=factanal(factors=2,covmat=r);faUniquenesses:

x1x2x3x4x5x6x7x8

0.0810.0750.1520.1350.0820.0330.0180.087现在是50页\一共有75页\编辑于星期一Loadings:

Factor1Factor2x10.2910.914

x20.3820.882

x30.5430.744

x40.6910.622

x50.7990.529

x60.9010.393

x70.9070.399

x80.9140.278

$loadingsFactor1Factor2X1-0.8171562-0.53123478X2-0.8674064-0.43231147X3-0.9151503-0.23258703X4-0.9487239-0.01184340X5-0.95935870.13153096X6-0.93764640.29276177X7-0.94395700.28715151X8-0.87992530.41074922Factor1Factor2SSloadings4.1133.224ProportionVar0.5140.403CumulativeVar0.5140.917Thedegreesoffreedomforthemodelis13andthefitwas0.3329现在是51页\一共有75页\编辑于星期一例9.12现有48名应聘者应聘某公司的某职位，公司为这些应聘者的15项指标打分，某指标与得分情况见例3.17，试用因子分析的方法对15项指标作因子分析，在因子分析中选取5个因子。rt=read.table("applicant.txt")factanal(~.,factors=5,data=rt)Loadings:

Factor1Factor2Factor3Factor4Factor5FL0.1270.7220.102-0.117APP0.4510.1340.2700.2060.258

AA0.1290.686

LA0.2220.2460.827

SC0.9170.167LC0.8510.1250.279-0.420

HON0.228-0.2200.777

SMS0.8800.2660.111EXP0.7730.171DRV0.7540.3930.1990.114AMB0.9090.1870.1120.165GSP0.7830.2950.3540.148-0.181POT0.7170.3620.4460.267KJ0.4180.3990.563-0.585

SUIT0.3510.7640.148f1:外露能力f2:经验f3:讨人喜欢f4,f5:专业能力、外貌现在是52页\一共有75页\编辑于星期一实例：表7.5是研究消费者对购买牙膏偏好的调查数据。通过市场的拦截访问，用7级量表询问受访者对以下陈述的认同程度（1表示非常不同意，7表示非常同意）。

V1：购买预防蛀牙的牙膏是重要的；

V2：我喜欢使牙齿亮泽的牙膏；

V3：牙膏应当保护牙龈；

V4：我喜欢使口气清新的牙膏；

V5：预防坏牙不是牙膏提供的一项重要利益；

V6：购买牙膏时最重要的考虑是富有魅力的牙齿。现在是53页\一共有75页\编辑于星期一表7.5牙膏属性评分得分表现在是54页\一共有75页\编辑于星期一将表7.5中的数据通过SPSS进行因子分析1.特征根和累计贡献率提取两个因子累计方差贡献率就达到82%，第三个特征根相比下降较快，因此我们选取两个公共因子。现在是55页\一共有75页\编辑于星期一2.因子的含义从因子载荷阵可以看出：因子1与V1（预防蛀牙），V3（保护牙龈），V5（预防坏牙）相关性强，其中V5的载荷是负数，是由于这个陈述是反向询问的；因子2与V2（牙齿亮泽），V4（口气清新），V6（富有魅力）的相关系数相对较高。因此，我们命名因子1为“护牙因子”，是人们对牙齿的保健态度；因子2是“美牙因子”，说明人们“‘通过牙膏美化牙齿’影响社交活动”的重视。从这两方面分析，对牙膏生产企业开发新产品都富有启发意义。旋转后因子载荷矩阵现在是56页\一共有75页\编辑于星期一R实现：rt=read.table("dataexample916.txt")factanal(rt,factors=2)Loadings:Factor1Factor2v10.968v20.749v30.898-0.140v40.784v5-0.887v60.830Factor1Factor2SSloadings2.5421.892ProportionVar0.4240.315CumulativeVar0.4240.739现在是57页\一共有75页\编辑于星期一典型相关性分析典型相关分析（canonicalcorrelationanalysis）是研究两组变量之间相关关系的一种统计分析方法，它能够有效地揭示两组变量之间的相互线性依赖关系。我们知道，在一元统计分析中，用相关系数来衡量两个随机变量之间的线性相关关系；用复相关系数研究一个随机变量和多个随机变量的线性相关关系。然而，这些统计方法在研究两组变量之间的相关关系时却无能为力。比如要研究生理指标与训练指标的关系，居民生活环境与健康状况的关系，人口统计变量（户主年龄、家庭年收入、户主受教育程度）与消费变量（每年去餐馆就餐的频率、每年出外看电影的频率）之间是否具有相关关系？阅读能力变量（阅读速度、阅读才能）与数学运算能力变量（数学运算速度、数学运算才能）是否相关？这些多变量间的相关性如何分析？现在是58页\一共有75页\编辑于星期一9.3.1总体典型相关1936年霍特林（Hotelling）最早就“大学表现”和“入学前成绩”的关系、政府政策变量与经济目标变量的关系等问题进行了研究，提出了典型相关分析技术。之后，Cooley和Hohnes(1971)，Tatsuoka(1971)及Mardia，Kent和Bibby(1979)等人对典型相关分析的应用进行了讨论，Kshirsagar(1972)则从理论上给出了最好的分析。典型相关分析的目的是识别并量化两组变量之间的联系，将两组变量相关关系的分析，转化为一组变量的线性组合与另一组变量线性组合之间的相关关系分析。其基本思想和主成分分析非常相似。首先在每组变量中找出变量的线性组合，使得两组的线性组合之间具有最大的相关系数。然后选取和最初挑选的这对线性组合不相关的线性组合，使其配对，并选取相关系数最大的一对，如此继续下去，直到两组变量之间的相关性被提取完毕为此。被选出的线性组合配对称为典型变量，它们的相关系数称为典型相关系数。典型相关系数度量了这两组变量之间联系的强度。现在是59页\一共有75页\编辑于星期一引入：一般情况，设是两个相互关联的随机向量，分别在两组变量中选取有代表性的综合变量U、V，使得每一个综合变量是原变量的线性组合，即典型变量

我们希望寻找使相关系数达到最大的向量a与b，由于随机向量乘以常数时并不改变它们的相关系数，所以，为防止结果的重复出现，令

现在是60页\一共有75页\编辑于星期一根据条件极值的求法引入Lagrange乘数，将问题转化为求的极大值，其中λ，ν是Lagrange乘数。根据求极值的必要条件得

A和B具有相同的特征根，a，b则是其相应的特征向量.，最大特征根λ2对应的特征向量就是所求的典型变量的系数向量.

2.典型变量和典型相关系数的计算现在是61页\一共有75页\编辑于星期一计算过程：求A的最大特征值和相应的特征向量，令：为第1对典型相关系数，为第1对典型变量。求A的第k个最大特征值和相应的特征向量，令：为第1对典型相关系数，为第k对典型变量。现在是62页\一共有75页\编辑于星期一样本典型相关在实际分析应用中，总体的协差阵通常是未知的，往往需要从研究的总体中随机抽取一个样本，根据样本估计出总体的协差阵，并在此基础上进行典型相关分析。设服从正态分布从该总体中抽取样本容量为n的样本，得到下列数据矩阵：样本均值向量,其中，样本协差阵现在是63页\一共有75页\编辑于星期一

计算过程：求A的第k个最大特征值和相应的特征向量，令：为第1对典型相关系数，为第k对典型变量。现在是64页\一共有75页\编辑于星期一9.3.3典型相关分析的计算R软件函数：cancorcancor(x,y,xcenter=TRUE,ycenter=TRUE)Argumentsxnumericmatrix(n*p1),containingthexcoordinates.ynumericmatrix(n*p2),containingtheycoordinates.xcenterlogicalornumericvectoroflengthp1,describinganycenteringtobedoneonthexvaluesbeforetheanalysis.IfTRUE(default),subtractthecolumnmeans.IfFALSE,donotadjustthecolumns.Otherwise,avectorofvaluestobesubtractedfromthecolumns.ycenteranalogoustoxcenter,butfortheyvalues.现在是65页\一共有75页\编辑于星期一表9.6生理指标和训练指标序号X1X2X3Y1Y2Y311913650516