多元統(tǒng)計分析-基于R（第3版） 課件 費宇 第12、13章-典型相關(guān)分析、多維標度分析

2024/9/5主編：費宇2第12章典型相關(guān)分析相關(guān)系數(shù)可以衡量兩個變量間的相關(guān)關(guān)系，但兩組變量之間的相關(guān)關(guān)系如何來度量呢？本章討論的典型相關(guān)分析(canonicalcorrelationanalysis)就是研究兩組變量之間相關(guān)關(guān)系的一種多元統(tǒng)計分析方法，它利用主成分的思想來討論兩組隨機變量的相關(guān)性問題，分別對兩組變量提取主成分，通過它們的相關(guān)性來度量兩組變量整體的線性相關(guān)關(guān)系.典型相關(guān)分析的思想首先由Hotelling于1936年提出，現(xiàn)在其已經(jīng)成為一種常用的分析兩組變量相關(guān)性的多元分析方法，在實際中應用廣泛.2024/9/5主編：費宇3第12章典型相關(guān)分析12.1

典型相關(guān)分析基本理論12.2

總體典型相關(guān)變量的概念及其解法12.3典型相關(guān)變量的性質(zhì)12.4原始變量與典型相關(guān)變量的相關(guān)系數(shù)12.5簡單相關(guān)、復相關(guān)和典型相關(guān)之間的關(guān)系12.6分量的標準化處理2024/9/5主編：費宇4第12章典型相關(guān)分析12.7樣本典型相關(guān)系數(shù)及其對應典型相關(guān)變量的計算12.8典型相關(guān)系數(shù)的顯著性檢驗12.9被解釋樣本方差的比例12.10案例分析及R操作2024/9/5主編：費宇512.1

典型相關(guān)分析基本理論典型相關(guān)分析是研究兩組變量之間整體的相關(guān)關(guān)系,它將每一組變量作為一個整體來進行研究,所研究的兩組變量可以是一組變量是自變量,另一組變量是因變量;當然,也可以兩組變量處于同等地位.2024/9/5主編：費宇612.1

典型相關(guān)分析基本理論典型相關(guān)分析的基本原理借助主成分分析的思想,在每組變量中找出變量的線性組合即新的綜合變量,使生成的綜合變量能代表原始變量的主要信息,同時,與由另一組變量生成的新的綜合變量的相關(guān)程度最大,這樣得到的一組新變量稱為第一對典型相關(guān)變量;同樣的方法可以找到第二對典型相關(guān)變量,第三對典型相關(guān)變量…，要求各對典型相關(guān)變量之間互不相關(guān).典型相關(guān)變量間的相關(guān)系數(shù)稱為典型相關(guān)系數(shù),它度量了這兩組變量之間關(guān)系的強度.此項最大化技術(shù)是努力將兩組變量間的一個高維關(guān)系濃縮到用少數(shù)幾個典型變量來表現(xiàn)．2024/9/5主編：費宇712.2總體典型相關(guān)變量的概念及其解法1.總體典型相關(guān)變量假設(shè)有兩組變量,一組變量為x=(x1,x2,…,xp

)T

,另一組變量為y=(y1,y2,…,yq

)T,且p≤q

,變量x與變量y的協(xié)方差陣為:2024/9/5主編：費宇812.2總體典型相關(guān)變量的概念及其解法為研究變量x與變量y之間的線性相關(guān)關(guān)系,我們考慮它們之間的線性組合u和v的方差和協(xié)方差分別為2024/9/5主編：費宇912.2總體典型相關(guān)變量的概念及其解法兩個新變量u和v之間的相關(guān)系數(shù)（即典型相關(guān)系數(shù)）為由于變量u和v乘以不為零常數(shù)不改變它們之間的相關(guān)性,即對任意常數(shù)c≠0,d≠0,有Corr(cu,cv)=Corr(u,v),所以通常需對a和b附加約束條件,使變量u和v不必要的重復,最好的約束條件是2024/9/5主編：費宇1012.2總體典型相關(guān)變量的概念及其解法我們的問題就變成在上述約束條件下求a和b,使得達到最大，于是有以下定義．定義12.1設(shè)，，維隨機向量的均值向量為0，協(xié)方差矩陣

（不妨設(shè)）.如果存在2024/9/5主編：費宇1112.2總體典型相關(guān)變量的概念及其解法和，令，，使得這樣得出的和稱為x，y的第一對（組）典型相關(guān)變量，稱為第一個典型相關(guān)系數(shù)；如果存在和使得1），和前面i-1對典型相關(guān)變量不相關(guān)；2）,;3）與的相關(guān)系數(shù)最大；2024/9/5主編：費宇1212.2總體典型相關(guān)變量的概念及其解法則稱，是x，y的第i對（組）典型相關(guān)變量，它們之間的相關(guān)系數(shù)稱為第i個典型相關(guān)系數(shù)（）．由拉格朗日乘數(shù)法,這個問題等價于求a和b使達到最大,其中μ1和μ2是拉格朗日乘數(shù).2024/9/5主編：費宇1312.2總體典型相關(guān)變量的概念及其解法將(12.7)兩邊分別對向量a和b求導,并令其為0,得方程組以aT和bT分別左乘(12.8)兩式得但,所以μ1=μ2=ρ,即μ1恰好就是u和v的相關(guān)系數(shù).兩邊左乘以

得,同理可得,記,則得2024/9/5主編：費宇1412.2總體典型相關(guān)變量的概念及其解法由方程組(12.8)的第二式得,將其代入(12.8)的第一式得

M1和M2的非零特征值皆為正數(shù)，，既是M1的特征根又是M2的特征根,a和b分別是M1和M2相相應的特征向量.于是求

和a,b的問題就轉(zhuǎn)化為求矩陣M1和M2的特征根和特征向量的問題.設(shè)ai是M1的屬于的特征向量，令2024/9/5主編：費宇1512.2總體典型相關(guān)變量的概念及其解法有則是M2的屬于的特征向量.2024/9/5主編：費宇1612.2總體典型相關(guān)變量的概念及其解法設(shè)M1的m個正特征根為,相對應的特征向量分別為a1,a2,…,am由式(12.10)得出且正交化，b1,b2,…,bm由式(12.11)得出且正交化,

從而可得m對線性組合每一對變量稱為一對典型變量,其中u1和v1稱為第一對典型變量,它們之間的相關(guān)系數(shù)ρ1即為第一典型相關(guān)系數(shù).

ui和vi稱為第i對典型變量,它們之間的相關(guān)系數(shù)ρi即為第i典型相關(guān)系數(shù).2024/9/5主編：費宇1712.3典型相關(guān)變量的性質(zhì)我們給出典型變量以下四個性質(zhì)（證明見本章附錄12）:(1)每一對典型變量ui及vi(i=1,2,…,m)的標準差為1.(2)同一組的任意兩個典型變量ui(i=1,2,…,m)彼此不相關(guān),典型相關(guān)變量vi(i=1,2,…,m)彼此不相關(guān),,，

.(3)不同組的任意兩個典型變量ui,vj(i=1,2,…,m;2024/9/5主編：費宇1812.3典型相關(guān)變量的性質(zhì)j=1,2,…,m)的關(guān)系為：(4)典型變量ui及vi的相關(guān)系數(shù)為ρi(i=1,2,…,m)，典型相關(guān)系數(shù)滿足關(guān)系式．2024/9/5主編：費宇1912.3典型相關(guān)變量的性質(zhì)在理論上,典型變量的對數(shù)和相對應的典型相關(guān)系數(shù)的個數(shù)可以等于兩組變量中數(shù)目較少的那一組變量的個數(shù),其中,u1及v1的相關(guān)系數(shù)ρ1反映的相關(guān)成分最多,所以稱為u1,v1第一對典型變量;u2及v2的相關(guān)系數(shù)ρ2反映的相關(guān)成分次之,所以稱u2,v2為第二對典型變量;以此類推.2024/9/5主編：費宇2012.4原始變量與典型相關(guān)變量的相關(guān)系數(shù)記2024/9/5主編：費宇2112.4原始變量與典型相關(guān)變量的相關(guān)系數(shù)2024/9/5主編：費宇2212.4原始變量與典型相關(guān)變量的相關(guān)系數(shù)則上面四個等式可以表示為其中.

2024/9/5主編：費宇2312.4原始變量與典型相關(guān)變量的相關(guān)系數(shù)其中.

其中.

2024/9/5主編：費宇2412.4原始變量與典型相關(guān)變量的相關(guān)系數(shù)其中.

其中.

2024/9/5主編：費宇2512.4原始變量與典型相關(guān)變量的相關(guān)系數(shù)其中.

其中.

2024/9/5主編：費宇2612.5簡單相關(guān)、復相關(guān)和典型相關(guān)之間的關(guān)系當p=q=1時，x與y之間的（唯一）典型相關(guān)就是它們之間的簡單相關(guān)；當p=1或q=1時，x與y之間的（唯一）典型相關(guān)就是它們之間的復相關(guān)；可見，復相關(guān)是典型相關(guān)的一個特例，而簡單相關(guān)是復相關(guān)的一個特例.第一個典型相關(guān)系數(shù)至少同x（或y）的任一分量與y（或x）的復相關(guān)系數(shù)一樣大，即使所有這些2024/9/5主編：費宇2712.5簡單相關(guān)、復相關(guān)和典型相關(guān)之間的關(guān)系復相關(guān)系數(shù)都較小，第一個典型相關(guān)系數(shù)仍可能很大；同樣，從復相關(guān)的定義也可以看出，當p=1（或q=1）時，x（或y）與y（或x）之間的復相關(guān)系數(shù)也不會小于x（或y）與y（或x）的任一分量之間的相關(guān)系數(shù)，即使所有這些相關(guān)系數(shù)都較小，復相關(guān)系數(shù)仍可能很大.2024/9/5主編：費宇2812.6分量的標準化處理一般來說，典型變量是人為定義的，也就是說它沒有實質(zhì)意義.如果使用原始變量，那么典型系數(shù)a，b的單位與x和y的單位成比例．而x和y的各分量的單位往往不全相同.我們希望在對各分量作標準化變換之后再作典型相關(guān)分析，這樣原始變量就有零均值和單位方差，典型變量就沒有測量值單位．2024/9/5主編：費宇2912.6分量的標準化處理記，，

，.為的相關(guān)矩陣.對x和y的各分量作標準化變換，即令.現(xiàn)在來求和的典型相關(guān)變量：2024/9/5主編：費宇3012.6分量的標準化處理于是2024/9/5主編：費宇3112.6分量的標準化處理因為，

，所以其中，.同理2024/9/5主編：費宇3212.6分量的標準化處理其中，.由此可見，為和的第i對典型系數(shù)，其第i個典型相關(guān)系數(shù)仍為ρi，在標準化變換下具有不變性，這一點與主成分分析有所不同．和的第i對典型變量具有零均值，且與x和y的第i對典型變量只相差一個常數(shù)．2024/9/5主編：費宇3312.7樣本典型相關(guān)系數(shù)及其

對應典型相關(guān)變量的計算前面我們是從變量x與變量y的協(xié)方差陣Σ出發(fā)考慮x與y的典型相關(guān)變量,這稱為總體典型相關(guān)變量,但在實際例子中一般并不知道Σ,因此通常采用樣本協(xié)方差陣S代替Σ.由12.6節(jié)的分析可知，在大多數(shù)情況下，我們在進行典型相關(guān)分析時，需將數(shù)據(jù)標準化，這時樣本協(xié)方差矩陣S即為樣本相關(guān)陣.根據(jù)樣本相關(guān)陣計算得到的典型相關(guān)變量，稱為樣本典型相關(guān)變量，具體計算過程如下．設(shè)容量為n的樣本來自正態(tài)總體,兩組變量的觀測值分別記為x=(x1,x2,…,xp

)T和y=(y1,y2,…,yq

)T,不妨設(shè)p≤q,則樣本數(shù)據(jù)矩陣為2024/9/5主編：費宇3412.7樣本典型相關(guān)系數(shù)及其

對應典型相關(guān)變量的計算2024/9/5主編：費宇35（1）計算樣本相關(guān)系數(shù)（1）計算樣本相關(guān)系數(shù)陣

,并將

剖分為其中,

是第一組變量x的關(guān)系數(shù)陣,

是第二組變量y的相關(guān)系數(shù)陣,而

、

(=)為變量x與變量y的相關(guān)系數(shù)陣.（2）計算典型相關(guān)系數(shù)及典型變量（2）計算典型相關(guān)系數(shù)及典型變量設(shè)首先求

的特征根，，…，（），并求對應的特征向量

,它是a1,a2,…,am的估計值;再求的特征根對應的特征向量,它是b1,b2,…,bm的估計值.這里

稱為樣本典型相關(guān)系數(shù),而

稱為樣本典型相關(guān)變量.2024/9/5主編：費宇3612.7樣本典型相關(guān)系數(shù)及其

對應典型相關(guān)變量的計算（3）記，由式(12.12第一等式)和(12.13第二等式)得2024/9/5主編：費宇3712.8典型相關(guān)系數(shù)的顯著性檢驗典型相關(guān)系數(shù)是否顯著的不為零,可以通過Bartlett大樣本卡方檢驗來完成.設(shè)的m個特征根為,則典型相關(guān)系數(shù)λ1的顯著性檢驗等價于以下檢驗H0:λ1=0,H1:λ1≠0.檢驗統(tǒng)計量為:式中，2024/9/5主編：費宇38式中，12.8典型相關(guān)系數(shù)的顯著性檢驗在檢驗水平α下,如果，則拒絕原假設(shè),認為第一對典型變量顯著相關(guān).一般,若前j-1個典型相關(guān)系數(shù)在水平α下是顯著的,則當檢驗第j個典型相關(guān)系數(shù)的顯著性時,檢驗統(tǒng)計量為

2024/9/539主編：費宇12.8典型相關(guān)系數(shù)的顯著性檢驗需要指出的是,在實際應用上,通常通過典型相關(guān)系數(shù)的顯著性檢驗以及典型變量和典型相關(guān)系數(shù)的實際解釋,來確定究竟保留幾對典型變量.所求得的典型變量的對數(shù)愈少愈容易解釋,最好是第一對典型變量就能反映足夠多的相關(guān)成分,這樣只保留一對典型變量便比較理想.2024/9/540主編：費宇12.9被解釋樣本方差的比例在進行樣本典型相關(guān)分析時，我們也想了解每組變量提取出的典型變量所能解釋的該組樣本總方差的比例，由此定量出典型變量所包含的原始信息量的大小．對于經(jīng)標準化變換后的樣本數(shù)據(jù)，第一組變量的樣本總方差為，第二組變量的樣本總方差為．稱為樣本典型相關(guān)變量，其中分別是原始變量x,y的標準化結(jié)果．2024/9/541主編：費宇12.9被解釋樣本方差的比例前r對典型相關(guān)變量對樣本總方差的貢獻為：2024/9/542主編：費宇其中，可依據(jù)式（12.19）計算．則第一組變量樣本方差由前r個典型變量解釋的比例為：12.9被解釋樣本方差的比例同理，第二組變量樣本方差由前r個典型變量解釋的比例為：2024/9/543主編：費宇其中，可依據(jù)式(12.22)計算．例12.1數(shù)據(jù)文件為exam12.1例12.1(數(shù)據(jù)文件為exam12.1)康復俱樂部對20名中年人測量了體重(x1)、腰圍(x2)、脈搏(x3)三個生理指標和引體向上次數(shù)(y1)、仰臥起坐次數(shù)(y2)、跳高(y3)三個訓練指標，數(shù)據(jù)詳見表12-1，分析生理指標與訓練指標的相關(guān)性．解:先讀取數(shù)據(jù)，求樣本相關(guān)系數(shù)矩陣．R程序和運行結(jié)果如下：9/5/2024主編：費宇44例12.1數(shù)據(jù)文件為exam12.1表12-1康復俱樂部數(shù)據(jù)9/5/2024主編：費宇45例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇46#exam12.1康復俱樂部數(shù)據(jù)的典型相關(guān)分析#打開數(shù)據(jù)文件exam12.1.xls，選取B1:G21區(qū)域，然后復制data12.1<-read.table("clipboard",header=T)#將exam12.1.xls數(shù)據(jù)讀入到data12.1中R=round(cor(data12.1),3);R#求樣本相關(guān)系數(shù)矩陣，保留三位小數(shù)x1x2x3y1y2y3x11.0000.870-0.366-0.390-0.493-0.226x20.8701.000-0.353-0.552-0.646-0.191x3-0.366-0.3531.0000.1510.2250.035y1-0.390-0.5520.1511.0000.6960.496y2-0.493-0.6460.2250.6961.0000.669y3-0.226-0.1910.0350.4960.6691.000例12.1數(shù)據(jù)文件為exam12.1生理指標和訓練指標之間的相關(guān)性強度中等，其中腰圍和仰臥起坐次數(shù)的相關(guān)系數(shù)最大為-0.646；組內(nèi)較大的是體重和腰圍的相關(guān)系數(shù)為0.87；引體向上次數(shù)和起坐次數(shù)的相關(guān)系數(shù)為0.696，仰臥起坐次數(shù)和跳高的相關(guān)系數(shù)為0.669．作典型相關(guān)分析，求典型相關(guān)系數(shù)和對應的典型變量的系數(shù)，R程序和運行結(jié)果如下：9/5/2024主編：費宇47例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇48X=scale(data12.1)#對數(shù)據(jù)進行標準化處理x=X[,1:3]#指定一組變量數(shù)據(jù)y=X[,4:6]#指定另一組變量數(shù)據(jù)library(CCA)#載入典型相關(guān)分析所用CCA包CCA=cc(x,y)#進行典型相關(guān)分析CCA$cor#輸出典型相關(guān)系數(shù)[1]0.79560.20060.0726CCA$xcoef#輸出x的典型載荷

[,1][,2][,3]x10.77541.8844-0.1910x2-1.5793-1.18060.5060x30.05910.23111.0508例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇49CCA$ycoef#輸出y的典型載荷

[,1][,2][,3]y10.34950.3755-1.2966y21.0540-0.12351.2368y3-0.7164-1.0621-0.4188因六個變量沒有用相同單位測量，這里用標準化后的系數(shù)進行分析．第一典型相關(guān)系數(shù)為0.796，它比生理指標和訓練指標兩組間的任一其他對的典型相關(guān)系數(shù)都大．調(diào)用相關(guān)系數(shù)檢驗腳本進行典型相關(guān)系數(shù)檢驗，確定典型變量對數(shù)，R程序和運行結(jié)果如下：例12.1數(shù)據(jù)文件為eg12.19/5/2024主編：費宇50source('corcoef_test.R')#調(diào)用典型相關(guān)系數(shù)檢驗腳本,若該腳本不在R的當前工作路徑下,則要將路徑設(shè)置清晰,如source('C:/ProgramFiles/corcoef_test.R')corcoef_test(r=CCA$cor,n=nrow(x),p=ncol(x),q=ncol(y))#進行典型相關(guān)系數(shù)檢驗

rQP[1,]0.795616.25500.0617[2,]0.20060.74500.9457[3,]0.07260.21090.6461檢驗總體中所有典型相關(guān)系數(shù)均為0的零假設(shè)時概率水平為0.062，故在（或）的顯著性水平下，拒絕所有典型相關(guān)系數(shù)均為0的假設(shè).也就是至少有一對典型相關(guān)是顯著的.從后面的檢驗結(jié)果例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇51可知，只有一對典型相關(guān)是顯著的．結(jié)合前面輸出的典型相關(guān)載荷結(jié)果可知，生理指標的第一典型變量為：它近似地是腰圍和體重的加權(quán)和，在腰圍上的權(quán)數(shù)更大些，在脈搏上的權(quán)數(shù)近似為0.來自訓練指標的第一典型變量為：其在仰臥起坐次數(shù)上的權(quán)數(shù)最大．這對典型變量主要是反映腰圍和仰臥起坐的負相關(guān)關(guān)系．例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇52輸出原始變量和典型變量的相關(guān)系數(shù).R程序和運行結(jié)果如下：CCA$scores$corr.X.xscores#輸出第一組典型變量與X組原始變量之間的相關(guān)系數(shù)

[,1][,2][,3]x1-0.62060.7724-0.1350x2-0.92540.3777-0.0310x30.3328-0.04150.9421CCA$scores$corr.Y.xscores#輸出第一組典型變量與Y組原始變量之間的相關(guān)系數(shù)

[,1][,2][,3]y10.5789-0.0475-0.0467y20.6506-0.11490.0040例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇53y30.1290-0.1923-0.0170CCA$scores$corr.X.yscores#輸出第二組典型變量與X組原始變量之間的相關(guān)系數(shù)

[,1][,2][,3]x1-0.49380.15498-0.0098x2-0.73630.07578-0.0022x30.2648-0.00830.0684CCA$scores$corr.Y.yscores#輸出第二組典型變量與Y組原始變量之間的相關(guān)系數(shù)

[,1][,2][,3]y10.7276-0.2370-0.6438y20.8177-0.57300.0544y30.1622-0.9586-0.2339例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇54整理后得表12-2．表12-2原始變量與第一對典型變量的相關(guān)系數(shù)由表12-2可知來自生理指標的第一典型變量u1與腰圍的相關(guān)系數(shù)為-0.925，與體重的相關(guān)系數(shù)為-0.621，它們都是負的．但在典型變量中體重的載荷為正（0.775），即體重在中的載荷和它與的相關(guān)系數(shù)反號．來自例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇55訓練指標的第一典型變量與三個訓練指標的相關(guān)系數(shù)都是正數(shù)，其中跳高在中的載荷（-0.716）和它與的相關(guān)系數(shù)（0.1622）反號；因此，體重和跳高在這組變量中分別是一個校正（或抑制）變量．一個變量具有同典型變量的相關(guān)系數(shù)相反符號的載荷似乎是矛盾的.為了理解這是怎樣發(fā)生的，考慮簡單的情況：用多元回歸方法由腰圍和體重來預測仰臥起坐次數(shù).一般來說，胖的人比瘦的人仰臥起坐次數(shù)少，這似乎是有道理的.假定這組樣本中沒有非常高的人，于是腰圍和體重之間的相關(guān)系數(shù)（0.87）是很大的．檢驗肥胖同自變量之間的相關(guān)性：例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇56腰圍大的人傾向于比腰圍小的人胖，因此腰圍與仰臥起坐次數(shù)之間的相關(guān)為負相關(guān)．體重大的人傾向于比體重小的人胖，于是體重與起坐起坐次數(shù)之間的相關(guān)為負相關(guān)．固定體重的值，腰圍大的人傾向于較強壯和較胖，于是腰圍的多元回歸系數(shù)應是負的．固定腰圍的值，體重大的人傾向于比較高和比較瘦，因此體重的多元回歸系數(shù)應為正的，它與體重和仰臥起坐次數(shù)間的相關(guān)反號．因此，第一典型相關(guān)一般解釋為以體重和跳高作為例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇57抑制變量來提高腰圍和起坐次數(shù)之間的相關(guān)性，但樣本的大小對于得出確定的結(jié)論還不夠大．計算典型變量解釋原變量方差的比例，R程序和運行結(jié)果如下：apply(CCA$scores$corr.X.xscores,2,function(x){mean(x^2)})#第一組典型變量解釋原第一組變量方差的比例[1]0.45080.24700.3022apply(CCA$scores$corr.Y.xscores,2,function(x){mean(x^2)})#第一組典型變量解釋原第二組變量方差的比例[1]0.25840.01750.0008例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇58第一對典型變量中u1解釋生理指標的標準方差的比例為0.451，第一對典型變量中v1解釋訓練指標的標準方差的比例為0.408，但兩者都不能很好地全面預測對應的那組變量.因為來自生理指標的標準方差被對方第一個典型變量v1解釋的方差比例為0.285，而來自訓練指標的標準方差被對方第一典型變量u1解apply(CCA$scores$corr.X.yscores,2,function(x){mean(x^2)})#第二組典型變量解釋原第一組變量方差的比例[1]0.28540.00990.0016apply(CCA$scores$corr.Y.yscores,2,function(x){mean(x^2)})#第二組典型變量解釋原第二組變量方差的比例[1]0.40810.43450.1574例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇59釋的方差比例為0.258．計算得分，并繪制得分等值平面圖.R程序如下：u<-as.matrix(x)%*%CCA$xcoef#計算得分v<-as.matrix(y)%*%CCA$ycoef#計算得分plot(u[,1],v[,1],xlab="u1",ylab="v1")#繪制第一對典型變量得分的散點圖,x軸名稱為u1,y軸名稱為v1，見圖12-1abline(0,1)#在散點圖上添加一條y等于x的線,以查看散點分布情況圖12-1康復俱樂部數(shù)據(jù)第一對典型相關(guān)得分等值平面圖例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇60例12.1數(shù)據(jù)文件為exam12.19/5/2024主編：費宇61通過作第一對典型相關(guān)變量得分等值平面圖可以看出，散點在一條近似直線上，雖然有偏離情況發(fā)生，但總體還是呈現(xiàn)出了線性相關(guān)關(guān)系.綜合來看，生理指標與訓練指標之間的關(guān)系雖有波動，但從整體來看較為明顯．12.10案例分析及R操作案例12.1(數(shù)據(jù)文件為case12.1)表12-3給出了2008-2016年我國科技活動和經(jīng)濟發(fā)展的部分代表指標.其中，科技活動指標：x1為R&D人員全時當量（單位：萬人年），x2為R&D經(jīng)費支出（單位：億元），x3為R&D項目(課題)數(shù)（單位：項），x4為發(fā)表科技論文數(shù)（單位：篇），x5為專利申請授權(quán)數(shù)（單位：件）；經(jīng)濟發(fā)展指標：y1為國內(nèi)生產(chǎn)總值（單位：億元），y2為城鎮(zhèn)居民家庭人均可支配收入（單位：元），y3為農(nóng)村居民家2024/9/562主編：費宇案例12.1數(shù)據(jù)文件為case12.1庭人均純收入（單位：元）；利用這些數(shù)據(jù)進行典型相關(guān)分析來分析我國科技活動和經(jīng)濟發(fā)展的關(guān)系．2024/9/563主編：費宇案例12.1數(shù)據(jù)文件為case12.1解:先讀取數(shù)據(jù)，求樣本相關(guān)系數(shù)矩陣．R程序和運行結(jié)果如下：2024/9/564主編：費宇#case12.1我國科技活動和經(jīng)濟發(fā)展的典型相關(guān)分析#打開數(shù)據(jù)文件case12.1.xls,選取B1:I10區(qū)域,然后復制case12.1<-read.table("clipboard",header=T)#將case12.1.xls數(shù)據(jù)讀入到case12.1中R=round(cor(case12.1),3);R#求樣本相關(guān)系數(shù)矩陣,保留三位小數(shù)案例12.1數(shù)據(jù)文件為case12.12024/9/565主編：費宇x1x2x3x4x5y1y2y3x11.0000.9880.9870.9950.9690.9900.9840.979x20.9881.0000.9990.9840.9920.9950.9970.995x30.9870.9991.0000.9810.9920.9930.9950.993x40.9950.9840.9811.0000.9690.9850.9820.979x50.9690.9920.9920.9691.0000.9870.9950.997y10.9900.9950.9930.9850.9871.0000.9980.993y20.9840.9970.9950.9820.9950.9981.0000.997y30.9790.9950.9930.9790.9970.9930.9971.000科技活動指標和經(jīng)濟發(fā)展指標之間的相關(guān)性很強，組內(nèi)相關(guān)性也很強.作典型相關(guān)分析，求典型相關(guān)系數(shù)和對應的典型變量的系數(shù)，R程序和運行結(jié)果如下：案例12.1數(shù)據(jù)文件為case12.12024/9/566主編：費宇X=scale(case12.1)#對數(shù)據(jù)進行標準化處理x=X[,1:5]#指定一組變量數(shù)據(jù)y=X[,6:8]#指定另一組變量數(shù)據(jù)library(CCA)#載入作典型相關(guān)分析所用CCA包CCAc12.2=cc(x,y)#進行典型相關(guān)分析CCAc12.2$cor#輸出典型相關(guān)系數(shù)[1]0.999980.867060.29489CCAc12.2$xcoef#輸出x的典型載荷

[,1][,2][,3]x10.1188-7.5900.1671x2-0.7634-1.44013.9738x30.65382.655-19.2193x4-0.16523.6411.6708x5-0.84572.6453.3987案例12.1數(shù)據(jù)文件為case12.12024/9/567主編：費宇CCAc12.2$ycoef#輸出y的典型載荷

[,1][,2][,3]y10.6016-13.1187.092y2-0.89369.993-19.181y3-0.70583.03612.106因六個變量沒有用相同單位測量，這里用標準化后的系數(shù)進行分析.第一典型相關(guān)系數(shù)為0.99998，它比科技活動指標和經(jīng)濟發(fā)展指標間的任一相關(guān)系數(shù)都大．調(diào)用相關(guān)系數(shù)檢驗腳本進行典型相關(guān)系數(shù)檢驗，確定典型變量對數(shù)，R程序和運行結(jié)果如下：案例12.1數(shù)據(jù)文件為case12.12024/9/568主編：費宇source('corcoef_test.R')#調(diào)用典型相關(guān)系數(shù)檢驗腳本,若該腳本不在R的當前工作路徑下,則要將路徑設(shè)置清晰,如source('C:/ProgramFiles/corcoef_test.R')corcoef_test(r=CCAc12.2$cor,n=nrow(x),p=ncol(x),q=ncol(y))#進行典型相關(guān)系數(shù)檢驗

rQP[1,]0.9999840.491410.00038[2,]0.867065.195750.73646[3,]0.294890.348460.95067檢驗總體中所有典型相關(guān)系數(shù)均為0的零假設(shè)時概率水平遠小于,否定所有典型相關(guān)系數(shù)均為0的假設(shè)，也就是至少有一對典型相關(guān)是顯著案例12.1數(shù)據(jù)文件為case12.12024/9/569主編：費宇的；典型相關(guān)系數(shù)檢驗p值的第二個值為0.736、第三個值為0.951，因此在顯著性水平為0.05的情況下只有一對典型相關(guān)是顯著的．結(jié)合前面輸出的典型相關(guān)載荷結(jié)果來看，科技活動指標的第一典型變量為：它近似地是專利申請授權(quán)數(shù)、R&D經(jīng)費支出和R&D項目（課題）數(shù)的加權(quán)和.在專利申請授權(quán)數(shù)上的權(quán)數(shù)最大，其次是R&D經(jīng)費支出，在R&D項目（課題）數(shù)上的權(quán)數(shù)也較大．案例12.1數(shù)據(jù)文件為case12.12024/9/570主編：費宇來自經(jīng)濟發(fā)展指標的第一典型變量為：它在城鎮(zhèn)居民家庭人均可支配收入上的權(quán)數(shù)最大，其次為農(nóng)村居民家庭人均純收入．輸出原始變量和典型變量的相關(guān)系數(shù)，R程序和運行結(jié)果如下：CCAc12.2$scores$corr.X.xscores#輸出第一組典型變量與X組原始變量之間的相關(guān)系數(shù)

[,1][,2][,3]x1-0.97482-0.20944-0.02677x2-0.99413-0.08663-0.04131案例12.1數(shù)據(jù)文件為case12.12024/9/571主編：費宇x3-0.99275-0.07465-0.08732x4-0.97608-0.157720.01888x5-0.998780.02189-0.03766CCAc12.2$scores$corr.Y.xscores#輸出第一組典型變量與y組原始變量之間的相關(guān)系數(shù)

[,1][,2][,3]y1-0.99059-0.11673-0.00708y2-0.99682-0.05889-0.01212y3-0.99907-0.024930.00931CCAc12.2$scores$corr.X.yscores#輸出第二組典型變量與X組原始變量之間的相關(guān)系數(shù)

[,1][,2][,3]x1-0.97480-0.18160-0.00790案例12.1數(shù)據(jù)文件為case12.12024/9/572主編：費宇x2-0.99411-0.07511-0.01218x3-0.99273-0.06472-0.02575x4-0.97606-0.136750.00557x5-0.998750.01898-0.01111CCAc12.2$scores$corr.Y.yscores#輸出第二組典型變量與y組原始變量之間的相關(guān)系數(shù)

[,1][,2][,3]y1-0.99061-0.13463-0.02399y2-0.99684-0.06792-0.04108y3-0.99909-0.028750.03156整理后得到表12-4.案例12.1數(shù)據(jù)文件為case12.1表12-4原始變量與第一對典型變量的相關(guān)系數(shù)2024/9/573主編：費宇來自科技活動指標的第一典型變量u1與R&D經(jīng)費支出、發(fā)表科技論文數(shù)、專利申請授權(quán)數(shù)的相關(guān)系數(shù)分別為-0.97482、-0.97608、-0.99878，u1與R&D人員全時當量、R&D項目（課題）數(shù)的相關(guān)系數(shù)分別為-0.99413、案例12.1數(shù)據(jù)文件為case12.12024/9/574主編：費宇-0.99275，它們都是負的，因此R&D人員全時當量、R&D項目（課題）數(shù)是抑制變量，其含義是它們在中的載荷（0.119，0.654）和它們與的相關(guān)系數(shù)（-0.97482，-0.99275）反號．來自經(jīng)濟發(fā)展指標的第一典型變量v1與三個經(jīng)濟發(fā)展指標的相關(guān)系數(shù)是負值，因國內(nèi)生產(chǎn)總值在中的系數(shù)和它與的相關(guān)系數(shù)是反號，故國內(nèi)生產(chǎn)總值也是一個抑制變量．案例12.1數(shù)據(jù)文件為case12.12024/9/575主編：費宇計算典型變量解釋原變量方差的比例，第一對典型變量能很好地全面預測對應的那組變量，來自科技活動指標的標準方差被第一個典型變量u1解釋的方差比例為0.9749，第一對典型變量v1解釋經(jīng)濟發(fā)展指標的標準方差的比例為0.9911；而來自經(jīng)濟發(fā)展指標的標準方差被第一典型變量v1解釋的方差比例為0.9748，經(jīng)濟發(fā)展指標的標準方差被對方第一典型變量u1解釋的方差比例為0.991．R程序和運行結(jié)果如下：案例12.1數(shù)據(jù)文件為case12.12024/9/576主編：費宇apply(CCAc12.2$scores$corr.X.xscores,2,function(x){mean(x^2)})#第一組典型變量解釋原第一組變量方差的比例[1]0.97490.01650.0024apply(CCAc12.2$scores$corr.Y.xscores,2,function(x){mean(x^2)})#第一組典型變量解釋原第二組變量方差的比例[1]0.99100.00590.000009apply(CCAc12.2$scores$corr.X.yscores,2,function(x){mean(x^2)})#第二組典型變量解釋原第一組變量方差的比例[1]0.97480.01240.0002apply(CCAc12.2$scores$corr.Y.yscores,2,function(x){mean(x^2)})#第二組典型變量解釋原第二組變量方差的比例[1]0.99110.00790.0011案例12.1數(shù)據(jù)文件為case12.12024/9/577主編：費宇計算得分，并繪制得分等值平面圖.R程序如下：u<-as.matrix(x)%*%CCAc12.2$xcoef#計算得分v<-as.matrix(y)%*%CCAc12.2$ycoef#計算得分plot(u[,1],v[,1],xlab="u1",ylab="v1")#繪制第一對典型變量得分的散點圖,x軸名稱為u1,y軸名稱為v1，見圖12-2abline(0,1)#在散點圖上添加一條y等于x的線,以查看散點分布情況圖12-2科技活動和經(jīng)濟發(fā)展數(shù)據(jù)第一對典型相關(guān)得分等值平面圖案例12.1數(shù)據(jù)文件為case12.19/5/2024主編：費宇78案例12.1數(shù)據(jù)文件為case12.12024/9/579主編：費宇由得分等值平面圖圖12-2可以看出，第一對典型相關(guān)變量得分散點在一條直線上分布，兩者之間呈高度線性相關(guān)關(guān)系，散點圖上沒有離開群體的差異點.這說明我國科技活動與經(jīng)濟發(fā)展之間的關(guān)系很穩(wěn)定，整體波動平穩(wěn)．2024/9/5主編：費宇主編：費宇，魯筠中國人民大學出版社，2024年第13章多維標度分析2024/9/5

82多維標度法（MultidimensionalScaling，MDS）是一種維數(shù)縮減方法，把高維的數(shù)據(jù)點映射到一個低維空間上，然后在低維空間中以較高的相似度重新展示這些點的數(shù)據(jù)結(jié)構(gòu)，并由此對原始高維數(shù)據(jù)進行統(tǒng)計分析．多維標度法起源于上世紀40年代的心理測度.1958年

Torgerson

正式提出了這一方法．目前應用廣泛，內(nèi)容豐富，方法較多．兩類方法：度量分析法，非度量分析法.第13章多維標度分析13.1

多維標度法的基本思想13.2

古典多維標度法

多維標度法的幾個基本概念；

已知距離矩陣時CMDS解；

已知相似系數(shù)矩陣時CMDS解；13.3

案例分析與R實現(xiàn)2024/9/5

8313.1

多維標度法的基本思想2024/9/5

84多維標度法是用低維空間

Rk

(k<p)中的n個點去重新標度和展示高維空間Rp中相應的n個點，將高維空間中的研究對象(樣本或變量)簡化到低維空間中進行定位、歸類和分析且有效保留研究對象間原始關(guān)系的多元數(shù)據(jù)分析方法．是一種維數(shù)縮減方法，主要思想是在降維的同時讓新得到的n個低維點與原來的n高維個點保持較高的相似度(如位置關(guān)系、距離、類別等)，使降維后重新標度的數(shù)據(jù)結(jié)構(gòu)發(fā)生的“形變”盡量?。?024/9/5

85多維標度分析與主成分分析的異同：

共同點：均先降維，再做簡明有效的分析；

均依據(jù)特定矩陣的特征值的累積貢獻率來

決定降維后的維數(shù)．

不同點：主成分法按包含信息大小選取主成分；

多維標度法按標度前后距離矩陣盡量接近

或盡量相似的原則來構(gòu)造擬合構(gòu)造點．13.2.1

多維標度法的幾個基本概念2024/9/5

8613.2

古典多維標度法古典多維標度(CMDS)解的定義對于距離陣，多維標度法的目的是要尋找較小的

k和Rk中的n個點,記，表示

xi

與xj

在Rk中的歐氏距離，使得與在某種意義下盡量接近,記稱為的一個古典多維標度(CMDS)解,稱xi

為的一個擬合構(gòu)造點，亦稱為的擬合構(gòu)圖，稱為的擬合距離陣．2024/9/587特別，當時，稱xi

為的一個構(gòu)造點，稱為的構(gòu)圖，注意的構(gòu)圖不唯一2024/9/5

88注：歐式距離陣也稱為歐式型距離陣或Euclid距離陣

如何判斷一個距離陣

是否為歐氏型距離陣呢？幾個相關(guān)矩陣的構(gòu)造：2024/9/5

89

距離陣為歐氏型距離陣的充要條件：

從歐氏距離陣D出發(fā)得到構(gòu)圖

X的步驟：見下面例13.1．2024/9/5

90例13.1(數(shù)據(jù)文件exam13.1)給定距離矩陣D如下，問它是否是歐式型矩陣？并求其的古典多維標度解（CMDS）

X及擬合構(gòu)造點．2024/9/5

91由B的特征值進行判斷并求X

解：在Excel中由矩陣D計算出B很方便：將D的上三角部分補齊輸入單元格區(qū)域A1:D4（如圖13-1

）；然后在A6單元格內(nèi)輸入“=-(A1^2)/2”后拖放填充至D9得到矩陣A；接著如圖13-1在A的左方、下方輸入公式并拖放填充計算出行均值、列均值和總均值；最后在單元格A12內(nèi)輸入公式“=A6-$E6-A$10+$E$10”后拖放填充至D15得到矩陣B．圖13-1在Excel中由矩陣D算出矩陣B（A12:D15）的簡單過程2024/9/5

92

最后在R中計算B的所有特征值來判斷D是否是歐氏型距離陣，并求出CMDS解

X及擬合構(gòu)造點．2024/9/5

93>B=read.table("clipboard",header=F)

#讀入矩陣B（先復制單元格區(qū)域A12:D15）>eig<-eigen(B);eig

#求B的特征值特征向量并顯示eigen()decomposition$values[1]70.6013138636.6152425180.795952899-0.002509279

B的四個特征值均大于等于0（第四個特征值很小可視為0），由定理13.1知D為歐氏型距離陣．為求D的CMDS解

X及擬合構(gòu)造點（實際上可近似看作構(gòu)造點），在R中輸入如下命令:2024/9/5

94>

D=read.table("clipboard",header=F)

#讀入矩陣D(先復制對應單元格區(qū)域A1:D4)>D13.1=cmdscale(D,k=2,eig=T)；D13.1#使用stats中的cmdscale函數(shù)，k取為2$points

[,1]

[,2][1,]2.526493-0.1531671[2,]

-7.110732-0.4721962[3,]3.545859-1.4437017[4,]1.038381

2.0690650$eig[1]7.059702e+016.611734e+007.937158e-012.428205e-15>sum(abs(D13.1$eig[1:2]))/sum(abs(D13.1$eig))

#計算a1.2[1]0.9898245>sum((D13.1$eig[1:2])^2)/sum((D13.1$eig)^2)

#計算a2.2[1]0.9998747D的CMDS解：2024/9/595前兩個特征值的累積絕對貢獻率和累積平方貢獻率均超過98%，說明k取為2是適當?shù)模膫€構(gòu)造點在坐標分別為(2.526,-0.153),(-7.111,-0.472),(3.546,-1.444),(1.038,2.069)．擬合距離陣（可與原始距離陣D對比）和擬合構(gòu)圖如下：圖13-2

距離陣D的擬合構(gòu)圖13.2.2

已知距離矩陣時CMDS解上面計算CMDS解的過程在R軟件中可使用stats包中的cmdscale()函數(shù)來實現(xiàn)，也可以使用MASS包中處理非度量MDS問題的isoMDS()函數(shù)來實現(xiàn)，但cmdscale()函數(shù)的好處是可以同時計算出的特征值和特征向量以及兩個累積貢獻率a1.k和a2.k的值．

2024/9/5

96

例13.2對表13-1給出的我國六個城市間的距離矩陣D，利用R軟件stats包中的cmdscale函數(shù)求

的CMDS解，給出擬合構(gòu)圖

及擬合構(gòu)造點．2024/9/597

北京濟南青島鄭州上海南京北京0

濟南4390

青島6683620

鄭州7144437720

上海12598867769840

南京10656266177103220表13-1我國六個城市間的道路距離（單位:公里）2024/9/5

98解在R中的程序為:>setwd("C:/data")#設(shè)定工作路徑>exam13.2<-read.csv("exam13.2.csv",header=T)#將exam13.2數(shù)據(jù)讀入>d13.2=exam13.2[,-1]#先去掉exam13.2第一列樣本名稱>rownames(d13.2)=exam13.2[,1]#為d13.2的行重新命名>D13.2=cmdscale(d13.2,k=2,eig=T);D13.2#使用基本包stats中的cmdscale函數(shù),k取為2由R計算結(jié)果可見，矩陣的六個特征值分別為：1051894，311141，59859，10289，0，-11999最后一個特征值為負，表明距離矩陣不是歐氏型距離陣．a(chǎn)1.k=94.3%，a2.k=99.7%，故k=2就可以了．由前兩個特征向量可得六個擬合構(gòu)造點分別為：(612.2,119.4)，(218.2,11.8)，(38.0,319.8)(193.7,430.2)，(646.1,57.4)，(416.0,78.3)2024/9/5

99>x=D13.2$points[,1];>y=D13.2$points[,2]>plot(x,y,xlim=c(-700,800),ylim=c(-300,600))#根據(jù)兩個特征向量的分量大小繪散點圖>text(x,y,labels=s(d13.2),adj=c(0,-0.5),cex=0.8)#將擬合點用行名標出再畫六個城市距離矩陣的擬合構(gòu)圖，并用中文標明（注意：擬合構(gòu)圖主要表示六個城市間的相對距離，和各城市在地圖上的實際位置可能不一致）圖13-3

我國六城市距離陣的擬合構(gòu)圖2024/9/5

1002024/9/5101

易計算出六個擬合構(gòu)造點在R2中的歐氏距離陣，將它們與表13-1中城市間的原始距離數(shù)據(jù)進行對比，可見大多數(shù)距離數(shù)據(jù)擬合較好，有6個城市距離相差公里數(shù)為個位數(shù)，只有兩個城市距離相差在50~60公里．>D1<-dist(D13.2$points[,1:2],method="euclidean",diag=T,p=2)>D1

北京

濟南

青島

鄭州

上海

南京北京

0.0000濟南

408.48960.0000青島

608.1445356.81290.0000鄭州

690.8417442.6457765.94660.0000上海1259.8374865.4658732.6946971.06630.0000南京1047.0937640.5733603.8211704.0036267.07110.000013.2.3

已知相似系數(shù)矩陣時的CMDS解2024/9/5102

例13.3對55個國家和地區(qū)的男子徑賽記錄作統(tǒng)計，每位運動員記錄8項指標：100米、200米、400米、800米、1500米、5000米、10000米和馬拉松．這8項指標的相關(guān)矩陣C如表13-2所示．求C的CMDS解，并給出擬合構(gòu)圖及擬合構(gòu)造點

(數(shù)據(jù)文件為exam13.3)．

2024/9/5

1032024/9/5

104

100米200米400米800米1500米5000米10000米馬拉松100米10.9230.8410.7560.7000.6190.6330.520200米0.92310.8510.8070.7750.6950.6970.596400米0.8410.85110.8700.8350.7790.7870.705800米0.7560.8070.87010.9180.8640.8690.8061500米0.7000.7750.8350.91810.9280.9350.8665000米0.6190.6950.7790.8640.92810.9750.93210000米0.6330.6970.7870.8690.9350.97510.943馬拉松0.5200.5960.7050.8060.8660.9320.9431表13－2運動員徑賽8項指標的相關(guān)系數(shù)矩陣2024/9/5

105>setwd("C:/data")#設(shè)定工作路徑>eg13.3<-read.csv("example13.3.csv",header=T)#將數(shù)據(jù)讀入>c13.3=eg13.3[,-1]#exam13.3的第一列為樣本名稱不是數(shù)值，先去掉>d13.3=round(sqrt(2-2*c13.3),3)#相似陣轉(zhuǎn)換成廣義距離陣,取三位小數(shù)>rownames(d13.3)=exam13.3[,1]>d13.3

X100米X200米X400米

X800米X1500米X5000米X10000米馬拉松100米0.0000.3920.5640.6990.7750.8730.8570.980200米0.3920.0000.5460.6210.6710.7810.7780.899400米0.5640.5460.000

0.5100.5740.6650.6530.768800米0.6990.6210.5100.0000.4050.522

0.5120.6231500米0.7750.6710.574

0.4050.0000.3790.3610.5185000米0.8730.7810.665

0.5220.3790.0000.2240.36910000米0.8570.7780.6530.5120.3610.2240.0000.338馬拉松0.9800.8990.7680.6230.5180.3690.3380.0002024/9/5

106>D13.3=cmdscale(d13.3,k=2,eig=T),D13.3

#k取為2，給出特征向量和特征值$points[,1][,2]100米-0.5273079050-0.4241746300.12620103400米-0.219741318-0.21411107800米0.003249949-0.227485801500米0.146859676-0.084591405000米0.3046913650.0667815610000米0.2994955630.06799955馬拉松0.4169272990.12579130$eig[1]8.842051e-011.650194e-011.258240e-018.178272e-026.736816e-02[6]4.726735e-022.266967e-023.045637e-18……>sum(abs(D13.3$eig[1:2]))/sum(abs(D13.3$eig))#計算a1.2[1]