劉應興的部落格

前幾天談的「重複量測與相關」談的是一個 X(i) 觀測值對應數個 Y(ij) 觀測值的問題，基本模型是

Y(ij) = α + β X(i) + ε(i,j),   i = 1, ..., k,  j = 1, ..., n(i)

其中假設 誤差項 ε(ij) 是 i.i.d. (0, σ^2)。本文要談另一種惰形：有 m 個個體，每個個體各有 n(i) 個重複測量的數據對 (X(ij), Y(ij))，那麼如何測量變數 X 與 Y 之間的相關？

這問題其實不是表面上那麼簡單：如果我們只布 (X(i), Y(i)) 觀測值，也就是說每一個體就是一組 (X, Y) 成對觀測值，毫無疑問地，我們看的就是不同個體有不同 X, Y 值，問的是不同個體間，不同 X, Y 觀測值之間的關聯；每個個體一個 X(i) 對應 n(i) 個觀測值 Y(ij)，而諸 Y(ij) 只是對應 X(i) 之單純重複觀測時，也一樣，表示每個個體有一個 Y 的理論值 μ(i)，我們的目標是看 μ(i) 與 X(i) 間的直線關聯。但同一個體有多個 (X(ij), Y(ij)) 配對資料時卻不一樣，我們可能要看的是：

(1) 排除個體差異，X 和 Y 純粹的關聯；
(2) 個體間 X 理論值 ξ(i) 與 Y 理論值 η(i) 之間的關聯；
(3) 不理會個體差異問題，X 與 Y 的總和關聯。

如果我們假設模型：

X(ij) = ξ(i) + u(ij), Y(ij) = η(i) + v(ij)

或簡單地表示為 X = ξ + u, Y = η + v，第一項的相關是 corr(u, v), 第二項是 corr(ξ,η), 第三項才是 corr(X, Y), 但第三項或許不是研究者真正關心的。

前段假設的模型可以說是測量誤差模型；不過就本文要談的東西來說，或可稱之為「個別差異模型」。不管如何稱呼，重要的是模型假設與想做的事。基本上 u 與 ξ, v 與 η 之間分別是零相關，甚至可以假設是獨立的；u 與 v 或諸 u(ij) 與諸 v(ij) 分別可能假設是 i.i.d. (0, τ^2) 與 (0, σ^2) 的。如果我們的目的是 corr(ξ, η), 也就是把 (X(ij), V(ij)) 看成是測量誤差模型的觀測值，則

Y(ij) = η(i) + v(ij)
      = α + β ξ(i) + v(ij)
      = α + β X(ij) - β u(ij) + v(ij)

關鍵是 v(ij) 與 X(ij) 可能透過 v(ij) 與 u(ij) 之關聯而不符合線性模型要求的：誤差項與自變項獨立，至少是零相關。也就是說：若我們的目標是 corr(ξ, η), 那麼計算 corr(X, Y) 可能是有偏的，除非 v(ij) 與 u(ij) 間零相關甚或相互獨立。但如果計算 x(i) = Σ_j X(ij)/n(i), y(i) = Σ_j Y(ij)/n(i), 由於它們分別是 ξ(i) 與 η(i) 之合理估計，似乎‵以 (x(i), y(i)) 來估計 ξ(i) 與 η(i) 間的相關是合理的。

在模型 X = ξ + u, Y = η + v 中, 假設 u 和 ξ 獨立，v 和 η 獨立，資料 X(ij) 是 ξ 的隨機觀測值 ξ(i) 與 u 的隨機觀測值 u(ij) 組合而成的，但我們看到的只有 X(ij) 而無法覲測 ξ(i) 和 u(ij)；變數 Y  之觀測值 Y(ij) 與 η(i), v(ij) 類同  X(ij) 與 ξ(i), u(ij)。如果我們要查探的是 corr(u, v), 令 X(i.) = Σ_j X(ij)/n(i), Y(i.) = Σ_j Y(ij)/n(i), 則

u'(ij) ≡ X(ij) - X(i.) = u(ij) - u(i.),   v'(ij) ≡ Y(ij) - Y(i.) = v(ij) - v(i.)

則由離差配對 (u'(ij), v'(ij)) 計算出來的相關係數，有如諸 (u(ij), v(ij)) 資料之樣本相關係數。是 corr(u, v) 的合理估計，只不過在計算中不是以諸 u(ij), v(ij) 之總樣本平均做置中，而是以分組樣本平均 u(i.) 與 v(i.) 做置中而已。對於這項相關係數，可以稱之為「消除個別差異後 X 與 Y 之間的相關係數，ξ(i), η(i) 即是 X, Y 變數中的個別差異成分。

如前面所述，用諸 (X(i.), Y(i.)) 的相關係數來估計 corr(ξ, η) 是合理的方案，也符合 Y 有重複觀測值 Y(ij) 而 X 僅有單一觀測值 X(i), 以及每一個體都只有單一 (X(i), Y(i)) 成對觀測值時的算法；不過，在每個個體只有單一 X(i) 時，我們只能認為 X(i) 就是 ξ(i) 而無測量誤差，因此 corr(ξ, η) = corr(X, η)。同時，在 X(i) 單一而 Y(ij) 有重複時，我們是要估計 corr(Y, X) 或 corr(ξ, η) 也決定了計算方式。在「重複量測與相關」一文中只談及計算得到的相關係數差異及對檢定的影響，其實蛛涉及模型的假設：如果 corr(v, X), 即 Y 變數的重複測量誤差項與 X 變數值無關，不‵管把 X(i) 複製成 Z(ij) 與 Y(ij) 配對計算相關係數，或把 Y(ij) 取平均值 Y(i.) 與 X(i) 配對計算相關係數，它們都是對 Corr(ξ, η) 的適當估計，只是用於估計的資料自由度不同，實際計算結果當然有差異。而基於統計經驗，資料自由度或誤差自由度愈大，統計效率愈高，或說平均誤差愈小，因此該文的總結意見認為：先對 Y(ij) 做平均，不如把 X(i) 複製擴充為 Z(ij)。而當重複資料成為 (X(ij), Y(ij)) 時，先計算平均，縮減資料為 (X(i.), Y(i.)) ，即使假設 Cov(u, η) = 0 = Cov(v, ξ)，則

Cov(X(i.), Y(i.)) = Cov(ξ, η) + Cov(u, v)/n(i) 
Var(X(i.)) = Var(ξ) + τ^2/n(i)
Var(Y(i.)) = Var(η) + σ^2/n(i)

因此，用平均值資料 (X(i.), Y(i.)) 來估計 corr(ξ, η) 仍不是很好的估計方式。或許更適當的做法，是先估計 Cov(u, v), Var(u), Var(v)，而後由相關係數計算式的三個成分：XY 共變異數，X 變異數，及 Y 變異數，其其估計式的期望值，扣去 Cov(u, v) 及 τ^2, σ^2 成分後再計算 corr(ξ, η) 的估計值。

在諸 n(i) 不等的情況，不考慮上述扣除 Cov(u, v), τ^2 及 σ^2 成分，先平均再計算相關有兩種算法：一是加權法：

r = Σ n(i) (X(i.)-X(..))(Y(i.)-Y(..))/√[Σ n(i) (X(i.)-X(..))^2 Σ n(i) (Y(i.)-Y(..))^2]

其中 X(..) = Σ n(i) X(i.)/N,  Y(..) = Σ n(i) Y(i.)/N,  N = Σ n(i)。另一是不加權法

r' = Σ (X(i.) - X'(..))(Yi.) - Y'(..))/√[Σ (X(i.) - X'(..))^2 Σ (Y(i.) - Y'(..))^2]

其中 X'(..) = Σ X(i.)/m,  Y'(..) = Σ Y(i.)/m。加權算法樣本相關係數 r 之組成部分期望值為

E[Σ n(i) (X(i.)-X(..))(Y(i.)-Y(..))]
     = (N - Σ n(i)^2/N) Cov(ξ, η) + (m-1) Cov(u, v)
E[Σ n(i) (X(i.)-X(..))^2]
     = (N - Σ n(i)^2/N) Var(ξ) + (m-1)τ^2
E[Σ n(i) (Y(i.)-Y(..))^2]
     = (N - Σ n(i)^2/N) Var(η) + (m-1)σ^2

而 Cov(u, v) 可用 ΣΣ u'(ij) v'(ij)/(N-m) 估計，τ^2 用  ΣΣ (u'(ij))^2/(N-m) 估計，σ^2 用  ΣΣ (v'(ij))^2/(N-m) 估計。則上列 r 可修正成

corr(ξ, η) = Cov(ξ, η)/√(Var(ξ) Var(η))

的較合理估計。採不加權法，則 r' 之組成三部分的期望值為

E[Σ(X(i.)-X'(..))(Y(i.)-Y'(..))]
   = (m-1) Cov(ξ, η) + (1-1/m) Σ 1/n(i) Cov(u, v)
E[Σ(X(i.)-X'(..))^2] = (m-1) Var(ξ) + (1-1/m) Σ 1/n(i) τ^2
E[Σ(Y(i.)-Y'(..))^2] = (m-1) Var(η) + (1-1/m) Σ 1/n(i) σ^2

同樣，藉由扣除上列三成分中 Cov(u, v), τ^2 及 σ^2 的估計量，可以把 r' 做修正，使其貌似更合理估計 corr(ξ, η)。再者，直接由 (X(ij), Y(ij)) 明細資料估計 corr(X, Y)，應也可以經校正而估得 corr(ξ, η), 因

E[ΣΣ (X(ij) - X(..))(Y(ij) - Y(..))]
    = (N - Σ n(i)^2/N) Cov(ξ, η) + (N-1) Cov(u, v)
E[ΣΣ (X(ij) - X(..))^2]
    = (N - Σ n(i)^2/N) Var(ξ) + (N-1) τ^2
E[ΣΣ (Y(ij) - Y(..))^2]
    = (N - Σ n(i)^2/N) Var(η) + (N-1) σ^2

注意其與 corr(ξ, η) 相干部分與平均值加權計算法相同，而其與 u, v 相關部分，在 corr(X, Y) 計算中乘數是 N-1, 而平均值加權計算之 r 中乘數為 m-1。所以如果沒有扎與 u, v 相干的部分扣除，平均值加權計算比 corr(X, Y) 更接近 corr(ξ, η)。三種笠法的三個組成成分（XY共變，X變異，及 Y變異）其期望值在我們最簡化設定下都可以分成兩個部分，一是個別差異，即 (ξ,η) 部分；另一是同一個體之重複變異，即 (u, v) 相干部分。三種計算方法造成兩部分乘數的不同：

如果把 Σ n(i)^2/N 看成是諸 n(i) 的加權平均（以自身加權），則明細算法兩部分比值是 N - avg(n) 比 N - 1, 平均值加權算法是  N - av(n) 比 m - 1,, 而平均值不加權算法是 m-1 比 (m-1)avg(1/n), 取平均值計算的兩部分比都大約 1 : 1/n，此處 n 是諸 n(i) 的某種平均，或 1/n 是諸 1/n(i) 的某種平均，因此如果不做校正，把三成分中屬於 (u, v) 的期望值部分以估計量取代替除，則當我們實際關心的是個體間 (X, Y) 之理論值 (ξ, η) 的相關時，顯然應先取平均值，相當於以平均值做理論值的估計。不過，如果我們利用重複觀測資料消除同一個體內之重複變異的效果：Cov(u, v), τ^2 = Var(u) 和 σ^2 = Var(v)，那麼直接修正 corr(X, Y) 之明細計算式，結果可能沒有實質差別。