劉應興的部落格

如果我們有原始資料 Y(ij), i = 1, ..., n, 代表個案， j = 1, ..., k, 代表變項，即資料項，則有兩種熟知的方式來縮減資料維度，即變項數，並使資料結構較清晰。這兩種方式就是主成分分析與因子分析，前者是變項的正交變換

W = P Y

之後忽略不重要，也就是變異數偏低的成分，可以 Y 的共變異矩陣為基礎；或以相關矩陣為基礎，相當於考慮 Y 的標準化變量。後一種是假設原資料變項 Y 符合因子模型：

Y = μ + L F  + ε

在此模型中 F 是各 (Y) 變項的共同因子，ε 各成分則是 Y 各成分變數所特有的，稱特殊因子。由於此模型的共同因子理論上不可能唯一決定，因此在「抽取因子」之後還會透過轉軸，讓因子軸 F 與觀測變項 Y 之間形成緊密而單純的結構，即每一觀測變項只「歸屬於」一個共同因子，使因子（軸）的意義較清晰。

如果我們只關心觀測變項與主成分變項間的關係，也就是 P 矩陣，則主成分分析只需要觀測變項的共變異矩陣，或相關矩陣；對因子模型而言，只要有 Y 的相關矩陣，就可以計算標準化型態的因子負荷量矩陣 L。

現在假如有 k 個物件，或稱「刺激（物）」, 例如商業上某種商品 k 個品牌，我們有這 k 物件間的偏好度、相似度或距離資料，假設形成這種結果是因為極少變數或因素在其中影響所致，例如第 j 個物件有隱藏變數值 (Z(j1),...,Z(jr)), j = 1, ..., k, 故物件 j, j' 之間有距離，例如歐氏距離 (Euclidean distance)

d(j, j') = √[Σ_t (Z(jt) - Z(j't))^2]

以向量形式來表示，

d(i, j)^2 = (Zi - Zj)'(Zi - Zj)

式中 v' 代表行向量 v 的轉置。將 d(i, j)^2 雙重置中 (doubly centered):

b(i,j) = d(i,j)^2 - Σ_i d(i,j)^2/k - Σ_j d(i,j)^2/k + ΣΣ d(i,j)^2/k^2

則 - b(i,j)/2 = (Zi - Z.)'(Zj - Z.)，其中 Z. = Σ Zj/k，各 Zi 是各物件在空間的位置，Z. 即所有物件的中心。此結果也就是說：令 B = [ -b(i,j)/2 ]， 則 B 是一個正半定矩陣 (positive semidefinite matrix)，因為 B 的各行或各列加總都是 0，因此不是正確定矩陣。現在，給予一個距離矩陣 D = [d(i,j)]，當然，d(j,j) = 0，令 B = [b(i, j)]，其中

b(i,j) =  -(d(i,j)^2 - Σ_i d(i,j)^2/k - Σ_j d(i,j)^2/k + ΣΣ d(i,j)^2/k^2)/2

這裡直接把前面的 b(i,j) 變號並折半以符合所需。如果 B 是一個正半定矩陣，則存在 r×k 矩陣 Z 使 B = Z'Z = [b(i,j)]，則 d(i,j) 為 Z 之第 i, j 兩行所代表的歐氏空間之點的歐氏距離；而 Z 的各行所代表的物件所在，圍繞著座標原點為其中心。

根據上面的結果而進行的分析稱為古典多元尺度法，也稱「主座標分析 (principal coordinates analysis, PCoA)」。其程序是：

1. 令 D" = [ d(i,j)^2 ],  d(i,j) 為物件 i 與 j 之距離，或稱臨近度。

2. 雙重置中：令 B = -(1/2) (I - J/k) D" (I - J/k), 其中 J 為元素皆為 1 的 k×k 矩陣。

3. 做 B 之正交對角線化，B = E Λ E', 其中 Λ 為特徵值構成的對角線矩陣，E 之各行為對應之單位特徵向量，各行相互正交。

4. 取最大特徵根 r 個，構成對角線矩陣 Λ*；並取其對應之特徵向量，構成 k×r 矩陣 E*，令

Z = ( E* Λ*^{1/2} )'

上列程序第四步，也是最後一步，E* Λ*^{1/2} 是 k×r 矩陣，每一行是 B 的一個特徵值的平方根乘以對應的單位特徵向量；而 Z 是 r×k 矩陣，每一行代表一個物件在 r 維空間的位置，這就是取轉置的原因。這程序是假設原始距離矩陣 [d(i,j)] 是存在誤差的，因此我們只期望 Z' Z 近似 B 矩陣，而不希求完全一致。或者，可能要完全複製出距離矩陣，需要較大的 r，也就是需要較高維度空間的點來表示物件位置，但我們希望在較低維度，如 r ≦ 3，使我們能圖示 k 個物件在空間中的相對位置。為衡量用 Z 之各行描述物件相對位置適當與否，一個指標是

Strain(D, Z) = [ ΣΣ (b(i,j) - Zi'Zj)^2/ΣΣ b(i,j)^2 ]^{1/2}

式中 D = [ d(i,j) ] 是原來的距離矩陣，Zi 是 Z 矩陣第 i 行，代表物件 i 的位置。古典多元尺度法的目標是使上列 Strain（張力？）最小，相當於由諸點 Zj 的歐氏距離計算出來的 b(i,j) 近似值，與與原距離 d(i,j) 得到的 b(i,j) 值差異的平方和最小。按 B 可以表示為

B = Σ_t λ_t u_t u'_t = [b(i,j)] = [ Σ_t λ_t u(it) u(jt) ]

式中 u(it) 為特徵值 λ_t 對應的特徵向量 u_t 的第 i 成分。前面的程序取前 r 個最大特徹向量構成 Z' 矩陣，故

ΣΣ (b(i,j) - Zi'Zj)^2 = Σ_i Σ_j (Σ_{t>r} λ_t u(it) u(jt))^2
       = Σ_i Σ_j Σ_{s>r} Σ_{t>r} λ_s λ_t u(it) u(jt) u(is) u(js)
       = Σ_{s>r} Σ_{t>r} λ_s λ_t Σ_i u(is) u(it) Σ_j u(js) u(jt)
       = Σ_{s>r} Σ_{t>r} λ_s λ_t δ(s,t)^2
       = Σ_{t>r} λ_t^2

如果不採用上述古典 MDS 算法，考慮 r = 1，我們希望找一個向量 x = [x_1, ..., x_k]' 使

Q(x) = Σ_i Σ_j (b(i,j) - x_i x_j)^2

最小化。由於有任一 x_i 趨於 ±∞ 時 Q(x) 都趨於 ∞，因此極小值都發生在臨界點，即對 x 微分為 0 的地方，經計算得

Bx = x x'x

由於 x'x 是純量，是 x 的模 (norm) 的平方，所以 x 是 B 的特徹向量，且 x'x 為特徵值，換言之，x 的模是對應之特徵值的（正）平方根。而我們的目標是極小化 Q(x)，所以取 B 的最大特徵值及其對應的特徵向量。接下來要取第二個向量 x*，在已有 x 的情況下，等於找 B 的第二大特徵值及對應的特徵向量。以此類推，要取 r 個向量定位 k 個物件，也就是說尋找 Z 使 Strain(D,Z) 極小化，就是前述古典 MDS 的算法。至於 r 的決定，類似主成分分析那棣用特徵值平方和比例，或 scree plot 等方法即可。

古典多元尺度法是所謂計量（或譯：度量）(metric) MDS 的一種。所謂計量多元尺度法，其根據就是數值型距離矩陣；另有非計量多元尺度法 (non-metric MDS)，其「距離」或「差異性 (dissimilarity)」或差異度矩陣只是順序性的：

D = [ δ(i,j) ],  δ(i,i) = 0, δ(i,j) > 0 if i ≠ j

當 k 個物件被定位在 r 維度空間時，其在此空間的距離順序與 δ(i,j) 的順序一致，也就是說兩者需維持單調 (monotonic) 關係。

首先，就計量多元尺度法而言，以前述  Strain 為評估配適狀況並求其極小化，即是古典多元尺度法；但可用其他準則，例如壓力優化 (stress majorization)

Stress(D, Z) = √[ΣΣ(d(i,j) - ||Zi - Zj||)^2]

之極小化。前面 Strain 是距離 d(ij) 經過平方，雙重置中，以及變號並折半後原始值與「配適值」的差異平方和；此處 Stress 則是 d(i,j) 原始值與配適值的差異平方和。（也有定義 Stress 是根號內再除以 d(i,j) 之平方和的，但如同在古典 MDS，並不影響 Z 矩陣之計算。另外，也有對各平方項加權的，例如用 w(i,j) = d(i,j) 加權。可惜，不論加權與否，壓力優化並無公式解，只能用迭代方式來極小化 Stress 指標。

如果我們的資料不是物件間的差異度或距離，而是相似度 (similarity) s(i,j), 在歐氏空間，假設「相似度」其實就是一種內積，例如

s(i,j) = X'i Xj ;  或  s(i,j) = X'i A Xj,  A 正確定

則 Xi 與 Xj 間之距離可以計算如下：

d(i,j)^2 = s(i,i) + s(j,j) - 2 s(i,j)

此轉換式適用於實向量之任一種內積。實際上如果有 [s(i,j)] 矩陣，而 d(i,j)^2 如上定義，則將距離矩陣 [d(i,j)] 用前面的雙置中程序轉換成 B = [b(i,j)]，則

b(i,j) = s(i,j) - s(i.) - s(.j) + s(..)

式中 s(i.) 如前，代表的是對對應註標的平均，所以等於把相似矩陣的元素做了雙重置中，S* = (I-J/k)S(I-J/k), S=[s(i,j)], 如果 S 是非負確定，則雙重置中後也是非負確定；但因 M = (I-J/k) 不可逆，並不能推論說 S* 非負確定則 S 也是。注意 M 是一個垂直（正交）投影矩陣，M^2 = M。並且，

設 u = M v 是 S 的特徵向量，則它也是 S* 對應同一特徵值的特徵向量。

(證) 設 S u = λ u，則 S* u = M S M u = M S u = M λ u = λ u.

另一方面，如果 S 的特徵向量 u 不在 M 的行空間，也就是 u 的元素和不為 0（u 在 c(M) 中其充要條件）, 則 u 不是 S* 的特徵向量。所以一般而言，尋找 Z（r×k 各列正交的矩陣）使 Z' Z ＝S 和使 Z' Z = S* 或只是「近似」, 並不是等價的問題。但 S = X'X，而 S* = (XM)'(XM) 也就是將各評價變數 Xi 減去所有物件之平均 X. 然後以交叉乘積的方式計算相似度。但 Xi - X. 等於使諸代表物件 i 的點 Xi 圍繞原點為其中心，或說是將各物件點 Xi 中心點設為座標原點，只是座標的平移。以此來看，內積做為相似度的指標並不算合理；再回顧古典多元尺度法將 d(i,j)^2 做雙重置中，雖是解決了用一組座標系統表現各物件，但似乎「雙重置中」只是一種技巧，卻缺乏底層邏輯。或者，依古典 MDS 將距離矩陣雙重置中的結果，也應對相似矩陣雙重置中而後進行 Z 矩陣的計算，無需先轉換成 d(i,j)^2。

如果將相似度矩陣視為非計量資料，彼時只考慮順序關係，相似度與差異度正好是相反順序，是以基本上沒有問題。

非計量多元尺度法的差異度矩陣 Δ = [ δ(i,j) ] 其中 δ(i,j) 無數值意義，僅有順序意義，而「真實的」距離 d(i,j) 可以表示成 δ(i,j) 的單調函數

d(i,j) = f(δ(i,j)),  f 為單調遞增實數值函數

非計量 MDS 通常用的評估準則是

Stress-1(Δ, Z) = √[ΣΣ(d(i,j) - d*(i,j))^2/ΣΣ d*(i,j)^2]

式中 d(i,j) = f(δ(i,j)), d*(i,j) = ||Zi - Zj||，實際計算程序如下：

1. 給予 Z = [ Zi ] 的起始值，計算 d*(i,j)；

2. 用單調迴歸法利用 d*(i,j) 估計 d(i,j) 也就是 f(δ(i,j))；

3. 極小化 Stress-1 得到新的 Z 矩陣，即諸 Zi 向量。

4. 如果結果不滿意，回到 2.；如果不收斂，設法修改 Z 再重複 2-3；如果收斂良好，結束計算。

上面談的都是給予一個矩陣表示物件之間的相似度、差異度或距離，但這些評價在實務上經常是多人提供的綜合結果，目的只是在低維度，最好是二維空間上定點物件相對位置，以此來分析這些物件的特性，例如哪些物件構成一個集群，哪些物件與另一些物件的差別何在等。但有時候我們想把提供評價上述資料的個體差異納入模型：

d(i,j,s)^2 ≒ Σ_t w(s,t)(Z(it) - Z(jt))^2 = d*(i,j,s)^2

此式意謂每個提供物件間距離評價的差異，是在於對各評估維度有不同權重。軟體 ALSCAL 取諸點 Zi 和權量 w(s,t) 最小化目標

Σ_s Σ_i Σ_j (d(i,j,s)^2 - d*(i,j,s)^2)^2

其中 d(i,j,s) 是 s 來源（個案 s）評價的物件 i, j 之間的距離，d*(i,j,s) 是根據物件座標矩陣 Zi 和權量（對角線矩陣）Ws 計算的距離。另一軟體 INDSCAL 採用如同前面古典 MDS 所做的，將平方距離做轉換：

b(i,j,s) = - (d(i,j,s)^2 - d(i,.,s)^2 - d(.,j,s)^2 + d(.,.,s)^2)/2

其中 d(.,j,s)^2 就是諸 d(i,j,k)^2 對 i 做平均，餘類推。而欲極小化的目標函數是

Σ_s || Bs - Z' Ws Z ||^2

式中 Bs 是 k×k 矩陣，元素為 b(i,j,s)；Ws = diag(W(s,1), ..., W(s,r))；Z 是 r×k 矩陣，其第 i 行是物件 i 的座標位置向量。矩陣的模 ||．|| 如果假設是其所有元素平方和的平方根，

||A||^2 = ΣΣ a(i,j)^2 = tr(A A') = tr(A' A)

除了因有多來源（個體）的評估，以及各維度權量不同外，形式上類似古典 ＭＤS。當

 d(i,j,s)^2 = Σ_t w(s,t)(Z(it) - Z(jt))^2 = (Zi - Zj)' Ws (Zi - Zj) 

時，b(i,j,s) = (Zi - Z.)' Ws (Zj - Z.)，如同古典 MDS，因為權重 W(s,t) 都是正值，諸 Bs 都是正半定矩陣。如果只有一個評估來源，假設 Bs = E Λ E'，在權量 Ws 已知的情況下，取最大的 r 個特徵值及其對應的特徵向量，則

Z' = E Λ^{1/2} W^{-1/2}

得解，而任意正權重分配不影響 || Bs - Z' Ws Z ||^2。在評價來源有 n > 1 個時，不能每個評價來源找一組 Zs，而是共同的一組座標 Z，另方面在物件座標矩陣 Z 已知的情況下，可以選擇 Ws 使各殘差平方和

|| Bs - Z' Ws Z ||^2 = Σ_i Σ_j (b(i,j,s) -  Zi' Ws Zj)^2

或即 tr((Bs - Z' Ws Z)^2) 最小，把 Z' 表示為 [ u_1, ..., u_r ], 則目標函數成為

tr((Bs - Z' Ws Z)^2) = tr((Bs - Σ_t W(s,t) u_t u'_t)^2)

如果 u_t 間相互正交，則目標函數對 W(s,t) 的偏導數是

-2 tr(Bs u_t u'_t - W(s,t) u_t (u'_t u_t) u'_t)

令偏導數為 0，解得 W(s,t) = u'_t Bs u_t/(u'_t u_t)^2。若 u_t 間並非正交，則需解聯立方程式：

Σ_t' (u_t u_t')^2 W(s,t') = u'_t Bs u_t,   t = 1, ..., r.

在諸 W(s,t) 已知的條件下，要取 Z 使

Σ_s Σ_i Σ_j (b(i,j,s) -  Zi' Ws Zj)^2

最小化，設 Z(ti) 代表 Zi 的 t 元素，即物件 i 在 t 座標的值，令 e_t 代表第 t 元素為 1 其餘 r-1 個元素都是 0 的標準單位向量，則上列目標函數對 Z(ti) 做偏微並令其為 0，得方程式

Σ_s Σ_j b(i,j,s) Z'j Ws e_t = Σ_s Σ_j Z'i Ws Zj Zj' Ws e_t

變動 i = 1, ..., k, 聯合寫成矩陣式

Σ_s Bs Z' Ws e_t = Σ_s Z' Ws Z Z' Ws e_t

又變動 t, 則 e_1, ..., e_r 姐成單位矩陣，故得

Σ_s Bs Z' Ws = Σ_s Z' Ws Z Z' Ws

其解為何，筆者尚不知曉，或許需迭代計算。

以上只是有個別評價（相似度或相異度）資料的一種模型兩種配適方法，當然還有其他模型或配適準則，例如 PROXSCAL 考慮

Σ_s ΣΣ_{i<j} w(ijs) ( φ(δ(ijs) - d(ij)(Xs))^2

之極小化。此式意謂每個個體所做原相似度或差異度評價 δ(ijs) 是非計量的；並且每個個體最後配適一個知覺空間 Xs, 是 k×r 矩陣，代表個體 s 所認為的，k 個物件在 r 維度空間的位置；最後殘差平方和所用權量更是允許隨所考慮之兩物件，及做出評價的個體而變。不過該模型卻又對 Xk 做了一些限制。

前面 ALSCAL 和 INDSCAL 所假設的模型是 k 個物件在所有評價者（個案）心中有一共同定位 Z 或其圖示，如同前面只有單一距離、差異或相似矩陣時一般，Z 或其圖示稱為「知覺圖」，由於輸入只有相互間的「距離」，因此並沒有方向，理論上知覺圖可以做任意平移及正交變換（旋轉及鏡射）；不過在多評估者的情形，因為有 W(s,t) 賦予知覺圖座標軸不同權重，其座標軸等於被定住了，否則在庢標軸旋轉時，各軸權重可能無意義或至少應配合旋轉方式重算。在多人評價的情形，如果知覺圖各軸被賦予有意義的維度名稱，把各評價者依其權重值 W(s,t) 標示在圖上，可以清楚看出各評價者在不同維度的偏重與否，也可以看出評價者整體（如果是有代表性的樣本）是否集中偏重某個維度。

一些實務上多元尺度分析的原始資料是由評價者按具體標準給予各物件評分，而非直接給予距離或差異評價。由評分資料 X(s,i,q)，原變數 Xq, q = 1, ..., m，可得出眾物件兩兩之間的距離 d(i,j,s)，而最後的知覺圖數據 Z(it),  新變數標示為 Ut, t = 1, ..., r，其在 k 個物件的值構成 u_t，矩陣 Z' 的第 t 行。將變數維度（個數）從 m 縮減為 r，因此多元尺度法一如主成分分析及因子模型，被認為是資料維度縮減的方法。由於 Xq 之間單位不同，相互非為正交，因此由 X(s,i,q) 資料得出 s 個體對物件 i, j 間距離的評價，歐氏距離並不恰當，一般認為 Mahalanobis distance

d(i,j,s) = √[ (X(s,i) - X(s,j))' V^{-1} (X(s,i) - X(s,j)) ]

更為適當。其中 X(s,i) 指諸 Xq 總合成一個行向量變數 X 在個體 s 物件 i 的值，即 [X(s,i,1), ..., X(s,i,m)]'；而 V 是 X 在不同物件的共同的共變異矩陣，可用前項樣本資料，即 X(s,i), s = 1, ..., n; i = 1, ..., k，做聯合估計；當不同物件間 X 的共變異矩陣實際上並不相等時，這樣得到的 V 等於所有不同 Vi, i = 1, ..., k, 的平均。另外，如果諸變數 Xq 是數值的，可以對諸 Ut 做迴歸

Xq = α(q) + β(1,q) U1 + . . . + β(r,q) Ur + error

迴歸函數與諸 Ut 相關的部分 β(1,q) U1 + . . . + β(r,q) Ur 代表 U-空間上的變化對 Xq 的影響，也代表 Xq 變數視為一個向量投影到 U-空間的結果，是從 U-座標原點到 (β(1,q), ..., β(r,q)) 的向量，可以用來表示變數 Xq 在知覺圖座標系統上的方向，例如 0×U1 + 1×U2 表示 X 變數與 U2 軸一致，則 X 變數在知覺座標系上就應指著與 U2 軸相同的方向。（這是筆者個人的理解，實際上是怎麼做的筆者並不清楚。）或者，可以將迴歸係數向量化為大小是複判定係數的正平方根的向量。兩種圖示都表示了 Xq 變數在知覺圖上的方向，但前者向量的大小表示了知覺圖座標對 Xq 的絕對作用力大小，但這很受 Xq 計量單位的影響（可將 Xq 標準化解決），也不能呈現 Ut 等對解釋 Xq 的力量強弱；而後者其向量大小則顯示了知覺圖變數 Ut 等對 Xq 的相對影響力，但無法顯現 Ut 等對 Xq 作用的絕對效果。知覺圖的座標軸要命名可以根據 k 個物件所在位置及物件本身特性做主觀命名；而若多元尺度分析是根據 n 個評價者的評價，特別是以諸變數 Xq 為基礎建立差異評價資料時，Xq 的意義可以很客觀地描述知覺圖座標軸的意義，所以將 Xq 在知覺圖的方向標示出來，讓我們可以較客觀地命名各軸，也更容易分析各物件在知覺圖上定位的含義。

偏好資料的多元尺度分析結果和距離或差異度的分析結果是不同的，也許依差異度建立的知覺圖物件 i 和 j 很接近，偏好卻相差不小；或差異度知覺圖相差甚遠的物件，其偏好知覺卻相近。從偏好知覺圖反向來看差異度知覺圖也一樣。其實這也很合理，如果差異度的評價是多維度的，有任一個維度兩物品有差距，則其差異度不能被其他維度所泯滅；而偏好可以因某一維度的表現而決定，可以既喜歡某個維度某種表現，又偏愛另一個維度另一種表現；至於差異度知覺圖看起來相差不遠的兩物件，卻可能因某個臨界點而在偏好上有了截然相反的結果。距離或差異度可能標準較客觀，而偏好則是較主觀的選擇，不一定都是客觀指標的平滑函數。

偏好資料通常不會只有單一評價，而是由 n 個個體對 k 個物件的評價所構成。評價方式可能是直接對 k 個物件排出一個喜好順序，或給予偏好評分；也可能兩兩或三三比較定出兩物或三物的順序。如何評價偏好及如何彙整，這是複雜的問題，這裡只想談談如何由偏好資料轉變成 k 個物件的偏好構圖，即偏好資料的多元尺度分析。

偏好資料的多元尺度分析也有許多方法，基本上可以由 n 個評價者對 k 物件的偏好評價，轉化為 k 物件間的距離。這其中就有多種假設，也就是多種轉換方式。那麼，k 物件的知覺圖，即 Zi, i = 1, ..., k, 可以古典、計量，或非計量方式建立。偏好資料的 MDS，最重要的是對每一評價者，可以定出一個「理想點」，這有外部法及內部法兩類。外部法是物件的知覺圖建立後，再根據蒐集到的每一個評價者對 k 個物件的評價資料決定其理想點。內部法是定義評價者理想點 Xs 與物件 i 之間的距離，例如

d(s,i) = d(Xs, Zi) = √[ Σ_t (X(s,t) - Z(i,t))^2 ]

在適當優化準則，例如極小化 stress 函數的平方：

Σ_i [ Σ_s (d(s,i) - d*(s,i) / Σ_s (d(s,i) - d(s.)^2 ]

之下得出理想點 Xs 及物件位置 Zi，其中 d(s,i) 必須符合評價者 s 對眾物件的偏好評價：愈偏好之物件 i 則 d(s,i) 愈小。上式中 d(s,i) 是根據 Xs, Zi 算出之距離，d*(s,i) 是 d(s,i) 對偏好順序做單調迴歸的結果，d(s.) 是 Xs 和諸 Zi 之間的平均距離，d(s.) = Σ_i d(s,i)/k。

當然計算從來都不是任一種統計分析的重點，特別是在這個各種套裝軟體能提供各種各樣的計算的時代，即使完全不知道如何由原始資料計算出許多漂亮結果，甚至繪出完美圖形，也一樣可進行各種統計分析。在多元尺度分析，上述結果之後對知覺圖座標軸的命名（解釋）；物件的知覺構圖及分群解釋；評價者在差異度資料分析之側重，在偏好分析之理想點，及消貴者分群之結果闡釋等，有些需要關於這統計方法的概念理解，更重要的是對所分析標的的瞭解。