第十七讲 列联表与对应分析离散多元分析
第十七讲 列联表与对应分析 离散多元分析
多项分布假设每次随机试验有m种可能的互斥事件，概率各为p1,,Pm多项满足p1++pm=1,pi≥0,记n次独立随机试验中各个事件发分布生的个数分别为x1.…xm，则其联合概率函数p(x1.., m) = m2xII1 pr, Z1 xi = n, xi ≥ 0我们称x=（x1.,xm）T服从多项分布Mm(n,(p1,...Pm）)。容易证明：命题1.若(x1,..,xm)~Mm(n,(p1,..,Pm)),P1 + ...+ pm=1,则(a)合并一些格子计数之后(概率也同样合并)仍然服从多项分布，(b)对任何下标集S = (i,,ik} {1,.,m}，xs =(xij.,Xik)Xsxi.+xix=N~Mk(N, (pin,.,Pi)/(pin +..+ pin),且xs与x-s条件独立.(c) (x1,..,xm)= (y1..,ym)ly1 +.. + ym = n, 其中y1,..,ym独立，Yi~Poisson(api).(a)的说明：比如5项分布（x1..,xs)~Ms(n,(p1,.,ps)）合并3成项分布(x1 + X5,X2,X3 + x4)~M3(n, (p1 + P5,P2,P3 +p4))
2 多项分布 假设每次随机试验有𝑚种可能的互斥事件，概率各为𝑝1, . , 𝑝𝑚, 满足𝑝1 + ⋯ + 𝑝𝑚 = 1, 𝑝𝑖 ≥ 0, 记𝑛 次独立随机试验中各个事件发 生的个数分别为 𝑥1, . , 𝑥𝑚，则其联合概率函数 𝑝 𝑥1, . , 𝑥𝑚 = 𝑛! ς𝑖=1 𝑚 𝑥𝑖 ! ς𝑖=1 𝑚 𝑝𝑖 𝑥𝑖 , σ𝑖=1 𝑚 𝑥𝑖 = 𝑛, 𝑥𝑖 ≥ 0 我们称𝐱 = 𝑥1, . , 𝑥𝑚 ⊤服从多项分布𝑀𝑚(𝑛, 𝑝1, . , 𝑝𝑚 )。 命题1 . 若 𝑥1, . , 𝑥𝑚 ~𝑀𝑚(𝑛, 𝑝1, . , 𝑝𝑚 ), 𝑝1 + ⋯ + 𝑝𝑚 = 1, 则 (a) 合并一些格子计数之后(概率也同样合并)仍然服从多项分布. (b) 对任何下标集𝑆 = 𝑖1, . , 𝑖𝑘 ⊂ 1, . , 𝑚 ，𝑥𝑆 = 𝑥𝑖1 , . , 𝑥𝑖𝑘 𝑥𝑆 |𝑥𝑖1 +⋯+𝑥𝑖𝑘 =𝑁~𝑀𝑘(𝑁, 𝑝𝑖1 , . , 𝑝𝑖𝑘 /(𝑝𝑖1 + ⋯ + 𝑝𝑖𝑘 )), 且𝑥𝑆与𝑥−𝑆条件独立. (c) 𝑥1, . , 𝑥𝑚 = 𝑦1, . , 𝑦𝑚 |𝑦1 + ⋯ + 𝑦𝑚 = 𝑛, 其中𝑦1, . , 𝑦𝑚独 立，𝑦𝑖~𝑃𝑜𝑖𝑠𝑠𝑜𝑛(𝜆𝑝𝑖). 容易证明: (a)的说明： 比如5项分布 𝑥1, . , 𝑥5 ~𝑀5(𝑛, 𝑝1, . , 𝑝5 ) 合并3成项分布: 𝑥1 + 𝑥5, 𝑥2, 𝑥3 + 𝑥4 ~𝑀3(𝑛, 𝑝1 + 𝑝5, 𝑝2, 𝑝3 + 𝑝4 )。 多项 分布
命题2.假设x= (x1,.,xm)T~Mm(n,p),p= (p1,.,Pm)T,则(a) E(x) = np, 即E(xi) = npi;(p:(1-p)-PiP2-PiPmP2(1- P2)-P2Pm(b) var(xi) = npi(1 - pi),-P2PlN=::Bcov(xi,x) = -npipj, i ± j;pm(1-pm)-pmPi- PmP2..即var(x) = n[diag(p) - ppT] ≤ nZ,(C)一(x- np) Nm(0,Z),n → 80. (退化多元正态,rank() = m-1)(a) (x- p)(m2)(αx-np)=Z m x-1 npi证明：(c)中心极限定理；1是的广义逆：(d)下面验证A=diagDdiadiag(p) = Im, diag=1mP0Pm) [diag(p) -ppT] = Im - 1mpT→ AZ = diag→ ZAZ= [diag(p) -ppT](Im-1mpT)diag(p)1m = p= diag(p) -ppT - diag(p)1mpT + ppT1mpTpT1m = 1= diag(p) - ppT = Z.3
3 命题2. 假设𝐱 = 𝑥1, . , 𝑥𝑚 ⊤~𝑀𝑚 𝑛, 𝐩 , 𝐩 = 𝑝1, . , 𝑝𝑚 ⊤ ,则 (a) 𝐸 𝐱 = 𝑛𝐩, 即𝐸 𝑥𝑖 = 𝑛𝑝𝑖 ; (b) var 𝑥𝑖 = 𝑛𝑝𝑖(1 − 𝑝𝑖), cov 𝑥𝑖 , 𝑥𝑗 = −𝑛𝑝𝑖𝑝𝑗 , 𝑖 ≠ 𝑗; 即𝑣𝑎𝑟 𝐱 = 𝑛 𝑑𝑖𝑎𝑔 𝐩 − 𝐩𝐩 ⊤ ≜ 𝑛Σ, (c) 1 𝑛 𝐱 − 𝑛𝐩 → 𝑁𝑚 0, Σ , 𝑛 → ∞. d (𝐱 − 𝑛𝐩) ⊤ 𝑛Σ − 𝐱 − 𝑛𝐩 = σ𝑖=1 𝑚 (𝑥𝑖−𝑛𝑝𝑖 ) 2 𝑛𝑝𝑖 → 𝜒𝑚−1 2 . 证明: (c)中心极限定理; (d) 下面验证𝐴 = diag 1 𝑝1 , . , 1 𝑝𝑚 是Σ的广义逆： diag 1 𝑝1 , . , 1 𝑝𝑚 𝑑𝑖𝑎𝑔 𝐩 = 𝐼𝑚, diag 1 𝑝1 , . , 1 𝑝𝑚 𝐩 = 𝟏𝑚 ⇒ 𝐴Σ = diag 1 𝑝1 , . , 1 𝑝𝑚 𝑑𝑖𝑎𝑔 𝐩 − 𝐩𝐩 ⊤ = 𝐼𝑚 − 𝟏𝑚𝐩 ⊤ ⇒ Σ𝐴Σ = 𝑑𝑖𝑎𝑔 𝐩 − 𝐩𝐩 ⊤ (𝐼𝑚−𝟏𝑚𝐩 ⊤) = 𝑑𝑖𝑎𝑔 𝐩 − 𝐩𝐩 ⊤ − 𝑑𝑖𝑎𝑔 𝐩 𝟏𝑚𝐩 ⊤ + 𝐩𝐩 ⊤𝟏𝑚𝐩 ⊤ = 𝑑𝑖𝑎𝑔 𝐩 − 𝐩𝐩 ⊤ = Σ. （退化多元正态 , 𝑟𝑎𝑛𝑘(Σ) = 𝑚 − 1） (1 ) (1 ) (1 ) 1 2 2 1 2 2 2 1 1 1 2 1 m m m m m m p p p p p p p p p p p p p p p p p p 𝑑 𝑑𝑖𝑎𝑔 𝐩 𝟏𝑚 = 𝐩 𝐩 ⊤𝟏𝑚 = 1 𝑑
拟合优度检Ho:Pi=pio，Pio已知．Ho成立时，验：Pearsond(Oi-Ei)2X2 = Zn, moXm-1,Zm卡方Einpio称为Pearson拟合优度卡方检验（Pearson,1900）。Pearson（19o0）提出了多项分布的拟合优度Pearson卡方检验，标志着现代统计的开端。Pearson(1904)提出列联表两个属性变量独立性的Pearson卡方检验K.Pearson (1goo).OntheCriterion that a Given System ofDeviationsfromtheProbableintheCaseofaCorrelatedSystemofVariablesis SuchthatitCanbeReasonablySupposedtohaveArisenfromRandomSampling.拟合优度K.Pearson(19o4).OnContingencyanditsRelationtoAssociationandNormalCorrelation.列联表1
4 𝐻0: 𝑝𝑖 = 𝑝𝑖0, 𝑝𝑖0已知. 𝐻0成立时， 𝑋 2 = σ𝑖=1 𝑚 (𝑥𝑖−𝑛𝑝𝑖0) 2 𝑛𝑝𝑖0 ≜ σ𝑖=1 𝑚 (𝑂𝑖−𝐸𝑖 ) 2 𝐸𝑖 → 𝜒𝑚−1 2 , 称为Pearson拟合优度卡方检验 (Pearson,1900) 。 d Pearson (1900) 提出了多项分布的拟合优度Pearson卡方检验，标 志着现代统计的开端。Pearson(1904)提出列联表两个属性变量独 立性的Pearson卡方检验. 拟合优度检 验: Pearson 卡方 • K. Pearson (1900). On the Criterion that a Given System of Deviations from the Probable in the Case of a Correlated System of Variables is Such that it Can be Reasonably Supposed to have Arisen from Random Sampling. 拟合优度 • K. Pearson (1904). On Contingency and its Relation to Association and Normal Correlation. 列联表
注意 Z (mp加项分母是均值E=npi而不是方差Vi=npi(1一pi)，npi这是Pearson卡方令人疑惑的地方。X2为什么不是标准化后的平方和：2xi-Exi-npiZmZm1npi(1-a1．简单理由：由命题1(c),x1,,xm可看作是独立的xi~Pois(npi),则E(xi) =var(xi) = npix-E(确是Poisson变量的标准化。故xi-npiVnpiVvar(xi)2.复杂理由：约束Zm=1xi=n导致分母上是均值。多项分布中x;的标准化确实应该是“,但约束Z1xi =n导致不同xi/Vi之间负相关，需要将这些相关性考虑在内（即命题2（d)中的-）。5
5 注意 σ𝑖=1 𝑚 (𝑥𝑖−𝑛𝑝𝑖 ) 2 𝑛𝑝𝑖 加项分母是均值 𝐸𝑖 = 𝑛𝑝𝑖 而不是方差𝑉𝑖 = 𝑛𝑝𝑖 1 − 𝑝𝑖 ， 这是Pearson卡方令人疑惑的地方。 𝑋 2为什么不是标准化后的平方和： σ𝑖=1 𝑚 𝑥𝑖−𝐸𝑖 𝑉𝑖 2 = σ𝑖=1 𝑚 𝑥𝑖−𝑛𝑝𝑖 𝑛𝑝𝑖 (1−𝑝𝑖 ) 2 ? 2. 复杂理由：约束 σ𝑖=1 𝑚 𝑥𝑖 = 𝑛 导致分母上是均值。 多项分布中𝑥𝑖 的标准化确实应该是 𝑥𝑖−𝐸𝑖 𝑉𝑖 , 但约束 σ𝑖=1 𝑚 𝑥𝑖 = 𝑛导致不同 𝑥𝑖 之间负相关， 需要将这些相关性考虑在内（即命题2 (d)中的Σ −）。 1. 简单理由：由命题1(c), 𝑥1, . , 𝑥𝑚可看作是独立的 𝑥𝑖~𝑃𝑜𝑖𝑠(𝑛𝑝𝑖), 则 𝐸 𝑥𝑖 = 𝑣𝑎𝑟 𝑥𝑖 = 𝑛𝑝𝑖 , 故 𝑥𝑖−𝑛𝑝𝑖 𝑛𝑝𝑖 = 𝑥𝑖−𝐸(𝑥𝑖 ) 𝑣𝑎𝑟(𝑥𝑖 ) 确是Poisson变量的标准化
例：二项分布（m=2情形）假设（x1,x2)~M2(n,(p1,p2)),通常记作x1~B(n,p1),则Pearson卡方(x1-np1)2 + (x2-np2)2 _ (x1-np1)2 + (n-x1-n(1-p1)2np1np2n(1-p1)npiX2=n-X1= (x1 - np1)2p2 = 1 - p1n(1-p1(np1(x1-np1)2np1(1-p1)xi-npi≥ z2/npi(1-pi)z是x1的通常的中心标准化，而左端-mp+2-mp2"中的npinp2X1 -np1 X2 -np2VnpiVnp2是约束x1+x2=n条件下x1，x2的同时中心标准化6
6 例：二项分布（𝑚 = 2情形） 假设 𝑥1, 𝑥2 ~𝑀2 𝑛, 𝑝1, 𝑝2 , 通常记作 𝑥1~𝐵(𝑛, 𝑝1),则Pearson卡方 (𝑥1−𝑛𝑝1) 2 𝑛𝑝1 + (𝑥2−𝑛𝑝2) 2 𝑛𝑝2 = (𝑥1−𝑛𝑝1) 2 𝑛𝑝1 + (𝑛−𝑥1−𝑛(1−𝑝1)) 2 𝑛(1−𝑝1) = (𝑥1 − 𝑛𝑝1) 2 1 𝑛𝑝1 + 1 𝑛(1−𝑝1) = (𝑥1−𝑛𝑝1) 2 𝑛𝑝1(1−𝑝1) = 𝑥1−𝑛𝑝1 𝑛𝑝1(1−𝑝1) 2 ≜ 𝑧 2 𝑧是𝑥1的通常的中心标准化，而左端(𝑥1−𝑛𝑝1) 2 𝑛𝑝1 + (𝑥2−𝑛𝑝2) 2 𝑛𝑝2 中的 𝑥1 − 𝑛𝑝1 𝑛𝑝1 , 𝑥2 − 𝑛𝑝2 𝑛𝑝2 是约束 𝑥1 + 𝑥2 = 𝑛 条件下𝑥1，𝑥2的同时中心标准化. 𝑥2 = 𝑛 − 𝑥1 𝑝2 = 1 − 𝑝1
列联表与Pearson卡方如果多项分布的格子由两个属性变量划分而成，则是列联表问题列联表列联表分析主要考察两个属性变量的关系。为了与连续变量情形的相关系数概念区分，Pearson将属性变量之间的关联性度量称为contingency，将属性变量交又分类得到的计数表格称为列联表contingencytable.Contingency：一个事件可能引发的事件（依赖性）、应急、临时、附加条款。假设属性变量x,y各有p,q个类别或水平。记wii为x取水平i，y取水平j的样本频数，Wpxg=（wii)称为p×q列联表(contingencytable）或交叉分类表(crosstabulation)。关心的问题：(1)x,y是否独立；(2)如果不独立，x与y哪个水平之间关联性最强?y=jWijrix=iCj
7 列联表与Pearson卡方 列联表 假设属性变量𝑥, 𝑦各有𝑝, 𝑞个类别或水平。记 𝑤𝑖𝑗为 𝑥 取水平𝑖， 𝑦 取 水平 𝑗 的样本频数， 𝑊𝑝×𝑞 = 𝑤𝑖𝑗 称为𝑝 × 𝑞 列联表(contingency table) 或交叉分类表(cross tabulation)。关心的问题： (1) 𝑥, 𝑦是否独立; (2) 如果不独立，𝑥 与 𝑦 哪个水平之间关联性最强? 如果多项分布的格子由两个属性变量划分而成，则是列联表问题. 列联表分析主要考察两个属性变量的关系。为了与连续变量情形 的相关系数概念区分，Pearson将属性变量之间的关联性度量称为 contingency，将属性变量交叉分类得到的计数表格称为列联表 contingency table 。 Contingency：一个事件可能引发的事件(依 赖性)、应急、临时、附加条款。 𝑥 = 𝑖 𝑤𝑖𝑗 𝑦 = 𝑗 𝑟𝑖 𝑐𝑗
例1.x=i与y=i共同发生的次数wii代表了水平i与i的关联度，但考察计数wi的同时应参考其所在行、列的其它计数，特别地需要考虑各行各列的边际总数。y=3x=i：x取值为第i个水平，6这里的代表水平类别，不9是通常的实数1511x=214520表中W23=5，X=2与y=3的关联大还是小？x23=5在第2行最大，似乎x=2与y=3关联程度较高。但观察第3列，5在第3列中并不大，因此x=2与y=3关联程度可能并不高。另一方面，x23=5在第3列中中等大小，这似乎说明y=3与x=2关联程度不高，但它在第2行是最大的计数。如何综合地考察wi相对于所在行总和、列总和的大小并不是简单的问题，也没有唯一答案。8
8 表中𝑤23 = 5， x = 2 与𝑦 = 3 的关联大还是小？ 𝑥23 = 5在第2行最大, 似乎x = 2 与𝑦 = 3关联程度较高。但观察 第 3 列， 5在第 3 列中并不大，因此 x = 2 与𝑦 = 3 关联程度可 能并不高。 另一方面，𝑥23 = 5在第3列中中等大小，这似乎说明𝑦 = 3 与𝑥 = 2关联程度不高，但它在第2行是最大的计数。 如何综合地考察𝑤𝑖𝑗相对于所在行总和、列总和的大小并不是简 单的问题，也没有唯一答案。 𝑥 = 2 𝑦 = 3 6 1 1 5 1 1 4 5 例1. 𝑥 = 𝑖 与 𝑦 = 𝑗 共同发生的次数 𝑤𝑖𝑗 代表了水平 𝑖 与 𝑗 的关 联度，但考察计数𝑤𝑖𝑗的同时应参考其所在行、列的其它计数， 特别地需要考虑各行各列的边际总数。 9 20 𝑥 = 𝑖 ：𝑥取值为第𝑖个水平， 这里的代表水平/类别，不 是通常的实数
例2.n=781件出土陶器按考古地点(0-6)和类型(a-d)两个属性变量交又分类得到下述列联表。（i,j)格子元素xi为地点i出土的第j类陶器的个数/计数。bd不同考古地点代表不同的时期总计ac或文化，因此01010398930534162751陶器类型分布相似的考古地4127311161点，其年代可能接近；32061431通过研究陶器类型分布判断446363713132考古地点的年代次序关系和54565910120交流情况。616281695218总计2839133374781列联表矩阵与一般的数据矩阵略有不同，其元素非负，代表了关联大小，与此类似的是强度矩阵（intensity）、丰度矩阵（abundance）网络的邻接矩阵等矩阵，矩阵(i，i)元素xii代表行标i与列标的联系紧密程度，恰当的标准化对于此类矩阵分析尤其重要
9 例2. 𝑛 = 781件出土陶器按考古地点(0-6)和类型(a-d)两个属性变量交 叉分类得到下述列联表。 (𝑖,𝑗)格子元素𝑥𝑖𝑗为地点𝑖出土的第𝑗类陶器的 个数/计数。 a b c d 总计 0 30 10 10 39 89 1 53 4 16 2 75 2 73 1 41 1 116 3 20 6 1 4 31 4 46 36 37 13 132 5 45 6 59 10 120 6 16 28 169 5 218 总计 283 91 333 74 781 不同考古地点代表不同的时期 或文化，因此 • 陶器类型分布相似的考古地 点，其年代可能接近； • 通过研究陶器类型分布判断 考古地点的年代次序关系和 交流情况。 列联表矩阵与一般的数据矩阵略有不同，其元素非负，代表了关联大 小，与此类似的是强度矩阵（intensity）、丰度矩阵(abundance）、 网络的邻接矩阵等矩阵，矩阵(𝑖,𝑗)元素𝑥𝑖𝑗代表行标𝑖与列标𝑗的联系紧 密程度，恰当的标准化对于此类矩阵分析尤其重要
如上所述，考察wi大小的时候需要综合考虑其所在行、列的其它计算，需要考虑wii相对于行计数总和ri与列计数总和c;的大小，并作某种归一化或标准化。记号假设p×q列联表X=（wi)，定义→行和:r = W1g = (r,.,rp)T列和:c = WT1p = (Ci,,Cq)T Dr = diag(r); Dc= diag(c)行归一化：P=（wii/r）=D-1W，其每一行的和都是1列归一化：Q=（wij/c）=WD=1，其每一列的和都是1W行归一化后，P的每一行的总和为1，各行之间具有可比性我们认为x各水平的特征刻画，比如第i行Pi = (Wi, ., Wiq)/ri是属性x的的水平i的特征刻画。同样O的第i列qj = (W1j..,Wpj)/cj是属性的的水平i的特征刻画10
10 如上所述，考察𝑤𝑖𝑗大小的时候需要综合考虑其所在行、列的其它 计算，需要考虑 𝑤𝑖𝑗 相对于行计数总和 𝑟𝑖 与列计数总和 𝑐𝑗 的大小， 并作某种归一化或标准化。 假设 𝑝 × 𝑞 列联表 𝑋 = 𝑤𝑖𝑗 ，定义 行和: 𝐫 = 𝑊𝟏𝑞 = (𝑟1, . , 𝑟𝑝) ⊤ 列和: 𝐜 = 𝑊⊤𝟏𝑝 = (𝑐1, . , 𝑐𝑞) ⊤ 𝐷𝑟 = 𝑑𝑖𝑎𝑔 𝐫 ; 𝐷𝑐= 𝑑𝑖𝑎𝑔 𝐜 行归一化: 𝑃 = 𝑤𝑖𝑗 /𝑟𝑖 = 𝐷𝑟 −1𝑊，其每一行的和都是1 列归一化: 𝑄 = 𝑤𝑖𝑗 /𝑐𝑗 = 𝑊𝐷𝑐 −1 ，其每一列的和都是1 记号 𝑊行归一化后， 𝑃的每一行的总和为1，各行之间具有可比性， 我们认为 𝑥 各水平的特征刻画，比如第 𝑖 行 𝐩𝑖 = 𝑤𝑖1, . , 𝑤𝑖𝑞 /𝑟𝑖 是属性 𝑥 的的水平 𝑖 的特征刻画。同样𝑄的第 𝑗 列 𝐪𝑗 = 𝑤1𝑗, . , 𝑤𝑝𝑗 /𝑐𝑗 是属性 𝑦 的的水平 𝑗 的特征刻画