第十二章 简介定性资料的统计分析
本章不是全面的介绍这方面的理论、方法和应用,而是初步反映一下这方面的主要内容,目的是展示进一步可学的知识,以便更好地解决实际问题。
§12.1 定性变量数量化
前面几章所介绍的各种统计方法,主要是研究与定量变量(或称间隔尺度变量)有关的问题,但在实际应用中,往往不可避免地要涉及到定性变量(或称名义尺度变量),例如人的性别、职业、天气状态,经济工作中选择的政策以及地层的构成类型等等,这些变量都只有各种状态的区别,而没有数量之区别。若定性变量不进入数学关系式,则会丢失信息,若要进入,又难于直接参加运算,于是从20世纪五十年代起开始发展了数量化理论,首先应用于“计量社会学”,六十年代后,逐步应用于各种学科,随着电子计算机的普及和发展,数量化理论将会在自然科学和社会科学的许多方面发挥出更大的作用。
如何对定性变量给以相应的数值描述,从而进行有关的统计分析,这就是数量化理论所研究的主要内容。
数量化理论已有专著出版,本节为了应用上的需要,仅介绍常用的0-1赋值法。 例如定性变量是性别,记为X,如此赋值:
当性别为女⎧1, 当性别为男⎧1,
X=⎨ 或X=⎨
⎩0, 当性别为女⎩0, 当性别为男
如此赋值的理由是简单,并没有任何数量大小的意义,它仅仅用来说明观察单位的特征
或属性,因此不同特性或属性的观察单位应取不同的值。
例如:天气可取晴、阴、雨三类,则用两个变量(X1,X2)表示天气,如此赋值:
当天气晴⎧(0,0),
⎪
(X1,X2)=⎨(1,0), 当天气阴
⎪(0,1), 当天气雨⎩
例如:有多种有害物污染了大气,由于有害物的结构不同,将污染物分为五类地区;甲、
乙、丙、丁、成戊将地区用4个变量(X1, X2, X3, X4)来表示,如此赋值:
甲类地区⎧(0,0,0,0),
⎪
乙类地区⎪(1,0,0,0),
⎪
(X1,X2,X3,X4)=⎨(0,1,0,0), 丙类地区
⎪(0,0,1,0), 丁类地区⎪⎪戊类地区⎩(0,0,0,1),
综上所述,推广为一般的赋值法如下:若某定性变量可取K类,则用K-1个变量表示,
如此赋值:
⎧(0,0,0, ,0), 第一类⎪
第二类⎪(1,0,0, ,0),
⎪第三类⎪(0,1,0, ,0),
(X1,X2, ,Xk-1)=⎨
第四类⎪(0,0,1, ,0),
⎪ ⎪⎪第K类⎩(0,0,0, ,1),
以上K个类的次序可以交换。
对于取K个类的定性变量,为什么用K-1个变量而不用K个变量表现?例如某定性变量可取甲、乙、丙、西四个类,可否如下赋值:
⎧(0,0,0,0), 取甲类⎪
取乙类⎪(1,0,0,0),
(X1,X2,X3,X4)=⎨
取丙类⎪(0,1,0,0),
⎪(0,0,1,0), 取丁类⎩
易知,如此赋值将使X1+X2+X3+X4 =1,不论是第几次观测,也不论定性变量取哪一类,皆使上式成立,即4个变量之和有稳定的线性关系式,知道其中任意三个就可推知另一个。
定性变量数量化后,就可以全部作为定量变量来统一处理进行预测或分类等研究。
§12.2 列联表
主要介绍二维列联表,对于三维以上的列联表只要在形式上稍加改变就能适用于高维表,原则上是一样的,只不过高维列联表符号更复杂一些,也增加些分析的难度。
1 列联表的概念
列联表讨论的主要是定性资料,此处介绍二维列联表的目的,不是将其数量化,而是直接进行分析并给出两个定性变量之间是否独立性检验。
先看一个简单例子:研讨吸烟与患肺癌的关系,这里用A表示一个人是否患肺癌,用B表示一个人是否吸烟,从一批被调查的对象中得到的统计表如下:
研讨患肺癌是否与吸烟有关?
这张统计表称为2×2列联表,表中考察两个定性变量A和B,每个变量有两类,即A分为患肺癌与未患肺癌两类,B分为吸烟与不吸烟两类,表中间的数值是频数,每一个被抽到的人,都可确定他的(AiBj)取值,比如表中数值60,表示被抽人群中吸烟又患肺癌的人数,数值32表示吸烟示患肺癌的人数。
一般2×2列联表形式如下:
其中nij(i,j=1,2)表示第i行Ai和第j列Bj的样品出现的频数,一般nij可取任意非负整数。
这是一个最简单的列联表,如果两个定性变量分别考察r和c类,则相应的列联表为r⨯c表(r和c可以不等)有如下形式:
如果一个问题涉及到很多的定性变量,相应的频数表就是一个高维列联表。
在概率统计中描述两个随机变量的相关程度是用线性相关系数,为了避免术语上的混淆,描述两个一性随机变量之间的相关性是指广义的相关性,称为关联性,两个定性随机变量之间的关联程度在某种意义上就是指的“不独立性”,它与独立的情形差距越大,就表明彼此的关系越密切,这种关系不一定是线性关系,然而在实际问题中,重要的是判断变量之间是否独立,因为不独立就意味着是关联的。如何判断是否独立有很多方法,这里仅介绍一种常用的皮尔逊拟合优度x2检验。
2×2列联表,对应一个多项分布,检验A与B是否独立,等价于检验:
H0:pij=pi.p.j
其中pij表示A为i、B为j的样品概率,pi.和p.j是相应的边缘概率,当独立性成立时, 理论频数为:npij=npi.p.j 其中n=
2
2
∑∑n
i=1j=1
ij
实际频数为:nij
运用x2检验作判定,需要知道列联表中实际频数与相应的理论频数。用估计量
nn
ˆ.j=.j代替pi.和p.j。基实际频数与理论频数有差异,这时可用其差值的大小ˆi.=i., pp
nn
来度量两个变量相关程度。相差愈大,表明H0为真的可能性愈小,即A与B无关的可能性愈小。相反差值愈小,即二愈接近,H0为真的可能性愈大,A与B之间相关的可能性愈小。为避免实际频数与理论频数的差值出现正负抵消,可采用差值的加权平方和来检验,于是给
2
出皮尔逊的拟合优度x统计量为:
2
x=
2
∑∑
i=1
=
∑∑
i=1
2
⎛nn⎫ nij-ni.j.⎪2 nn⎪⎝⎭
nn.jj=1
ni.
nn
2nn-nn2
iji..j
nnni..jj=1
2
()
它的极限分布是自由度为1的x2分布,根据给定的显著性水平a,查x2分布表得到临界值λa。若x2≥λa则拒绝H0,表示A与B之间不独立,存在相关,若x2
H0,表明A与B之间独立,不存在相关。
将前面的例子作x2检验:
计算
(106⨯60-63⨯92)2(106⨯32-43⨯92)2
x=+
106⨯63⨯92106⨯43⨯92(106⨯3-63⨯14)2(106⨯11-43⨯14)2++
106⨯63⨯14106⨯43⨯14(6360-5796)2(3392-3956)2=+
[1**********]6(318-882)2(1166-602)2++
9349263812
=0.75857+0.5775+4.98489+3.40239=9.66360
2
取显著性水平a=0.05,自由度为1,查x2分布表,临界值λa=3.84。显然x2=9.66360>3.84,表明在5%的显著性水平上,拒绝H0即说明吸烟与肺癌不独立,而是存在相关的。
如果列联表中变量间存在相关,那么如何度量变量间的相关程度?又如何从一个变量去预测另一变量呢?解决这类总是还有很多方法,已超出本书范围,不再详述,有兴趣的读者可查阅这方面的参考书。
§12.3 对数线性模型
如前所述,列联表能够反映定性变量之间的关系,但能否像定量变量那样建立起数学模型如方差分析模型、回归分析模型等以便进一步描述定性变量之间的复杂关系呢?对数线性模型和Logistic回归模型就是解决这一问题的极为有效的方法,它们从不同角度出发导出不同的处理方法。
对数线性模型,近十年来是国外实际工作者常用的方法,它的主要优点是可以把方差分析和线性模型的一些方法系数地移植过来,在概念和理解上均可进行对比,对数线性模型能够估计模型中各个参数,而这些参数值使各个变量的效应和变量间的交互作用效应得以数量化。下面即将看到这些结论。
1 模型
对数线性模型又分为很多种类型,常用的模型有:饱和模型(当变量间相互不独立时),非饱和型(变量间相互独立),谱系模型(包含高阶效应)等。
下面从2×2的频数表与概率表出发,推导对数线性模型: (频数表)(概率表)
将概率取对数后进行分解处理,使处理后的变量有较好的数学、统计的性质。
⎛pij⎫ ⎪ μij=lnpij=lnpi.p.j pi.p.j⎪⎝⎭
pij
=lnpi.+lnp.j+ln
pi.p.j
记Ai=lnpi., Bj=lnp.j, (AB)ij=ln由上式可写成
pijpi.p.j
μij=Ai+Bj+(AB)ij i,j=1.2
显然上式的结构类似于两因子有交互作用,各因子均为二水平的方差分析模型,于是令
μi.=∑μij, μ.j=∑μij, μ=∑∑μij
j=1
i=1
i=1j=1
2222
然后再进行平均,对i, j =1,2
i.=μi., .j=μ.j ..=μ..
记
1
21212
αi=μi.-μ..
βj=j.-..
λij=μij-i.-.j+..
则有关系式:
⎧⎪
⎪μij=..+αi+βj+λij⎪22⎪
i,j=1,2 ⎨且αi=0, βj=0,
j=1⎪i=1
⎪22⎪λ=λ=0 ⎪i=1ijj=1ij⎩
可见通过上边分解处理,可以完全化成与方差分析模型有同样的结构,因此借助于方差分析的术语,上式中μ..表示“总平均效益”,αi表示A属性的“主效应”,βj表示B属性
∑∑
∑∑
的“主效应”,λij表示,A, B的 “交互作用效应”,直观可以理解当交互作用效应为0,即等价于A、B独立。上式模型称为对数线性模型的饱和模型,当λij=0时,称为非饱和模型。
n.jni.
ˆp=,.j,这时就可以看
nnn
到对数线性模型是将列联表上每个单元的频数作为因变量,表上所有变量作为自变量,建立各个自变量的效应与每个单元频数的对应之间的函数关系。因而可以用它分析列联表上的各个变量的关系。主效应αi或βj若大于0,表明效应为正;若小于0,表明效应为负。αi是
ˆij=在实际应用时概率可用其估计量代替,即p
ˆi.=,p
第一个变量的第i个水平对总平均效应μ..的增减量;βj是第二个变量的第j个水平对总平均效应μ..的增减量,λij代表变量1和变量2在各自的第i个水平和第j个水平之间交互作用效应,是其交互作用对总平均效应的增减量。若λij
nij
1
i.=
2
∑
1μij=
2j=1
n
∑
j=1
2
(ln
nijn
) „ 第i行频数对数的平均
1μ.j=
2
∑
1μij=
2i=1
n
n
∑
i=1
2
(ln
nijn
) „ 第j列频数对数的平均
11
μ..=μ..=
44
∑∑
i=1
1
μij=
4j=1
2
∑∑(ln
i
j
22
nijn
) „ 各个观测值对数的总平均即总平均效
应对本章前面的例子,按上述模型估计各效应参数。
各单元的频数对数表:
计算:
α1=1.-..=2.5964-2.5141=0.0823
α2=2.-..=2.4318-2.5141=-0.0823
β1=.1-..=3.2800-2.5141=0.7659
β2=.2-..=1.7482-2.5141=-0.7659
主效应估计值:
计算:
λ12=μ12-1.-.2+..=1.0986-2.5964-1.7482
+2.5141=3.6127-4.3426=-0.7319
λ21=μ21-μ2.-μ.1+μ..=2.4657-2.4318-3.2800
+2.5141=4.9798-5.7188=-0.7320
λ11=μ11-1.-.1+..=4.0943-2.5964-3.2800
+2.5141=6.6083-0.8319=0.8319
λ22=μ22-2.-.2+..=2.3979-2.4318-1.7482
+2.5141=4.9120.1800=0.7320
变量间交互作用效应估计值:
λ11=0.8319 λ12=-0.7399 λ22=0.7320
λ21=-0.7317
主效应大于0,表明效应为正,如α1=0.0832>0是因为患肺癌比未患肺癌的人多;主效应小于0,表明效应为负,如β2=-0.76599
§12.4 Logistic回归
对数线性模型是将列联表中每格的概率(或理论频数)取对数后分解参数获得的,Logistic回归模型是将概率比取对数后,再时行参数化而获得的。研究概率比这样的量在不少问题中是常常遇到的,当列联表中因变量是一个多级分类的变量时,就需要考虑两两比较的情况。
Logistic回归要解决的问题,类似于普通回归所要解决的很大一类问题。比如在医药行业中,因变量y取0, 1, „, g这g+1个不同的值,y=0表示正常情况类型,y=1….,g表示不同用药后的反应,显然它与药的剂量x1 ,性别x2,年龄x3,体重x4,血压x5,„等等自然变量有关,这里因变量是定性的,自变量有定性的也有定量的,问这些自变量对一个定性变量的关系是否独立?不独立又会具有什么形式的联系?是线性的还是非线性的等等。
1 Logit变换
为了给出Logistic回归模型,先介绍Logit变换。在现实问题中,人们常常要研究某一事件A发生的概率p以及p值的大小与某些因素的关系,但由于p对x的变化在p = 0或p =1的的附近是很敏感的,或说是缓慢的比如像可靠系统,可靠度p已经是0.988了,即使再改善条件和系统结构,它的可靠度增长只能在小数点后面的第三位或第四位,于是自然希望寻找一个p的函数θ(p),使它在p = 0或p =1附近变化幅度较大,而且函数的形式也不要太复杂,根据数学上导数的意义,提出用时希望p = 0或p =1,
dθ(p)
来反映θ(p)在p附近的变化是很合适的,同dp
dθ(p)
有较大的值,因此取 dp
dθ(p)111
==+ dpp(1-p)p1-p
即
θ(p)=ln
p
1-p
称上式为Logit变换。
p
由于θ(p)=ln,因此p也可用θ表示:
1-p
eθ
p= θ
1+e
如果θ是某些自变量x1, „, xk,的线性函数
∑a
i=1
k
k
xi,则p就x1, „, xk的下述函数:
p=
ei=1
∑aixi∑aixi
k
k
1+i=1
显然,如果如果p对x不是线性关系,则θ对x就可以是线性关系了,这就是Logit变换带来的方便。
2 Logistic回归模型
设因变量y为二值定性变量,用0,1表示取两个不同的状态,y=1的概率p=P(y=1)是我们要研究的对象。如果有很多因素影响y的取值,这些因素就是自变量,记为x1,„, xk,这其中既有定性变量也有定量变量,最重要的一个条件是
pln=a0+a1x1+ +akxk 1-p适合于上式线性函数称为Logistic线性回归。如果有已知函数g(x1, ,xk),其中含有若干待定的参数,且
pln=g(x1, ,xk) 1-p称上式为非线性Logistic回归。
Logistic线性模型相当于广义线性模型,因此可以系统地应用线性模型的方法,在处理时就比较方便。
如果某一事件A发生的概率p依赖于一些自变量x1, ,xk(定性定量均可),对x=(x1, ,xk)'观测了m组结果,在第a组中,共试验了na次,A发生了ra次,于是A发
ˆa=生的概率pa可用p
ra
来估计。假定pa适合于Logistic线性回归关系式,即 napaln=β0+β1xa1+ +βkxak a=1, ,m 1-pa
ˆa代入上式中pa就有关系式: 其中xai表示xi在第a组所取的值。用p
ln
ˆap
=β0+β1xa1+ +βkxak+εa a=1, ,m ˆa1-p
其中εa是随机误差项。记
⎡1x11 x1k⎤
⎢1x⎥ xˆap212k⎥, X=⎢ ya=ln
m⨯(k+1)⎢ ⎥ˆ 1-pa
⎢⎥⎣1xml xmk⎦
y=(y1, ,ym)',假定E(εa)=0,V(εa)=va,ε1, ,εm相互独立,于是就有
⎡β0⎤
0⎤⎡v1⎢⎥
β⎥
Ey=X⎢1⎥∆Xβ Var(y)=V=⎢ ⎢⎥⎢ ⎥
⎢vm⎥⎢⎥⎣0⎦β⎣k⎦
ˆ 这就是典型的线性模型,因此β的最小二乘估计(因V不是单位阵,应该相应地加权)β
有公式:
ˆ=(X'βV-1X)-1X'V-1y
这样就求得了β的估计值。要讨论某些xi是否对A发生的概率有无影响,也即要检验xi相应的回归系数βi=0这一假设是否成立,这时要搬用线性模型中已知的结论时,必须知道y
是否服从正态分布以及V的估计,数学上可以证明:
lnˆa=v
ˆa⎛⎫ppa1
的渐近分布为正态的,即N ln, 1-pnp(1-p)⎪⎪V是的va的估计值为ˆa1-paaaa⎭⎝
11
ˆˆnapa(1-pa)
如果在m组试验结果中,有的ra=0或ra=na,此时
ln
ˆapra
=ln ˆa1-pna-ra
会取-∞或∞的值,ya就不是一个有限的值,上述方法就会行不通,于是要进行修正,修正
rapa
的目的是使ln尽可能接近ln,可以证明下面的修正是合理的。
na-ra1-pa
ra+0.5
a=1, ,m
na-ra+0.5(na+1)(na+2)
ˆa=v a=1, ,m
na(ra+1)(na-ra+1)Za=ln
上式Za称为经验的Logistic变换,相应的线性模型也称为经验的Logistic回归模型:
⎧⎡EZ1⎤⎪⎢⎥EZ= ⎪⎢⎥=Xβ⎪⎢⎪⎣EZm⎥⎦
⎨
ˆ0⎤⎡v1⎪
⎥⎪Var(Z)=⎢ ⎢⎥⎪
⎢ˆm⎥v⎪⎣0⎦⎩
ˆ的估计值如下: ˆa为修正后的表达式,β其中Za,v
ˆ=(X'Vβ
-11-2V2
X)-1X'V
-
11
-2V2
~。 用加权最小二乘估计回归系数时,第a组的权系数是va
引入经验Logistic变换后,Logistic回归就和普通的回归一样了,而且用的是加权最小
二乘,在普通回归中遇到的一些检验的问题,如回归系数的显著性检验、共线性问题、„„等等,在这里都有同样的意义,而且处理方法也相似,经验的Logistic回归是Logistic回归中最常用的方法。
3 实例
对本章所例举的实例,讨论吸烟与肺癌的关系,作经验的Logistic变换并再次检验吸烟与肺癌是否相关。
将吸烟人作为一类,不吸烟人作为一类,此时m =2, n1=92, n2=14, r1=60, r2=3。 作经验Logistic变换:
r+0.560.5
Z1=ln1=ln=0.6214
n1-r1+0.532.5
-
1
2
r2+0.53.5
=ln=-1.1896
n2-r2+0.511.5(n1+1)(n1+2)(92+1)(92+2)8742ˆ1=v===0.047204
n1(r1+1)(n1-r1+1)92⨯61⨯33185196Z2=ln
ˆ2=v
(n2+1)(n2+2)(14+1)(14+2)240
===0.3571428
n2(r2+1)(n2-r2+1)14⨯4⨯12672
为了作假设检验,设吸烟人患肺癌的概率是p1,未患肺癌的概率就是1-p1,不吸烟的人
患肺癌的概率是p2,未患肺癌的概率就是1-p2,作Logit变换:
pp
θ1=ln1, θ2=ln2
1-p11-p2设θ2=θ,则
θ1=θ2+(θ1-θ2)=θ+∆
因此患肺癌是否与吸烟有关,就归结为检验
H0:∆=0
是否成立。
ˆ1+vˆ2),由于Z1,Z2可以看成是相互独立的正态变量。因此,H0成立时,Z1-Z2~N(0,v
即
Z1-Z21+v2v
计算
~N(0,1)
Z1-Z21+v2v
=
1.81100.0472+0.357
=
1.8110.4043
=
1.8110
=2.8483
0.6358
取a=0.05,查N(0,1)表得临界值λa=0.834。
2.8483>0.834,应拒绝H0,即认为吸烟与不吸烟对于患肺癌是有影响的。
综上所述,对于列联表,既可以用对数线性模型也可以用Logistic模型,要根据实际问题的需要来选择模型。
定性资料统计分析的内容是很丰富的,这里不再详述,有兴趣的读者可查有关资料。
附录:矩阵代数
矩阵和行列式是研究多元统计分析的重要工具,这里针对本书的需要,对有关矩阵的知识作复习,同时再简单的补充矩阵的分块和矩阵的微商。(因为在部分章节中要用到)。
§1 矩阵及基本运算
1 矩阵的定义
将n×p个实数,a11, a12, „, a1p, a21, a22, „, a2p, „, an1, an2, „, anp排成如下形式的矩形表,记为A
⎡a11a12 a1p⎤⎢a⎥a a21222p⎥ A=⎢⎢ ⎥⎢⎥aa a⎢⎥n1n2np⎣⎦
则称A为n⨯p阶矩阵,一般记为A=(aij)n⨯p,其中aij是A中元素,它也可以是复数,但本书的aij均为实数。当n = p时,称A为n阶方阵;若p =1,A只有一列,称A为列向量,
记为
⎡a11⎤⎢a⎥a=⎢21⎥
⎢ ⎥⎢⎥⎣an1⎦
当n=1时,A只有一行,称A为行向量,记为
a'=(a11,a12, ,a1p)
当A为n阶方阵,称a11,a22, ,ann为A的对角线元素,其它元素称为非对角元素,若方阵A的非对角元素全为0,称A为对角阵,记为
O⎤⎡a11
⎢⎥a22⎥=diag(a,a, ,a) A=⎢1122nn⎢⎥ ⎢⎥Oann⎦⎣
进一步若a11 = a22 =„= ann =1,称
O⎤⎡1
⎢⎥1⎥ A=⎢⎢⎥ ⎢⎥O1⎣⎦
为n阶单位阵,记为In或A=I(I的阶数从上下文可明确)。
若A为n×p阶矩阵,它的转置A′是阶p×n阶矩阵,即
⎡a11a21 an1⎤⎢a⎥a a1222n2⎥ A'=⎢⎢ ⎥⎢⎥aa a⎢2pnp⎥⎣1p⎦
若A是方阵,且A′=A,则称A为对称阵;若方阵A=(aij)的元素aij = 0,对一切i<j成立,则称
⎡a11⎢aA=⎢21
⎢ ⎢⎣an1
a22 an2
O⎤⎥⎥ ⎥
⎥
ann⎦
为下三角阵;若A′为下三角阵,则称A为上三角阵。
2 矩阵的运算
若A与B是n×m阶阵,则A与B的和定义为:
A+B=(aij+bij)n⨯m
若a为一常数,它与A的积定义为
aA=(aaij)n⨯m
若A=(aik)p⨯q,B=(bkj)q⨯r,则A与B的积定义为: AB=(
q
∑a
k=1
ikkj
b)p⨯r
在一般情况下AB≠BA。
容易验证上述运算符合下面运算规律:
A+(-1)A=0
(AB)′=B′A′ (A′)′=A
(A+B)′=A′+B′ A(BC)=(AB)C A(B+C)=AB+AC (A+B)C=AC+BC AI=IA=A
a(A+B)=aA+aB a(AB)=(aA)B=A(aB)
若A为方阵,满足A′A=AA′=I,则称A为正交阵。
§2 行列式、逆矩阵和矩阵的秩
1 行列式
一个p阶方阵A=(aij)p⨯p对应一个数,记为|A|
⎡a11⎢a21
A=⎢
⎢ ⎢⎢⎣ap1
=
j1j2 jp
a12 a1p⎤
a22 a2p⎥⎥ ⎥
⎥
ap2 app⎥⎦
τ(j1 jp)
∑(-1)
j1 jp
a1j1,a2j2, ,apjp
称|A|为A的行列式或记为detA。这里∑表示对所有p元排列求和。τ(j1 jp)表示,若jt>js,则称这两个数组成一个j1 jp的逆序数(在一个排列j1j2 ,jt, ,js jp)逆序。一个排列中逆序的总数称为此排列的逆序数)。
不难证明:
A=
∑(-1)
j=1
p
i+j
,aijMij
其中Mij是在A中去掉第i行、第j列而形成的p-1阶方阵。称(-1)i+jMij为aij的代数余子式,通常记为Aij。
直接由行列式的定义来计算行列式是很麻烦的,通常利用行列式的如下一些性质可以简化行列式的计算:
(1)若A的某行(或列)为零,则|A|=0; (2)|A|=|A′|;
(3)将A的某行(或列)乘以数a所得矩阵的行列式等于a|A|; (4)若A的两行(或列)相同,则|A|=0;
(5)若将A的两行(两列)互换所得矩阵的行列式等于-|A|;
(6)若将A的某一行(或列)乘上一个常数后加到另一行相应的元素上,所得矩阵的行列式不变,仍等于|A|。
2 逆矩阵
设A为p阶方阵,若|A|≠0,则称A是非退化阵或称非奇异阵,若|A|=0,则称A是退化阵或称奇异阵。
若A是p阶非退化阵,则存在唯一的矩阵B,使得AB=Ip,B称为A的逆矩阵,记为B=A-1,不难证明:
A-1
⎡A11⎢⎢A⎢A12=⎢A⎢⎢ ⎢A1p⎢⎢⎣A
A21
AA22A A2pA
Ap1⎤
⎥A⎥Ap2⎥
⎥A⎥ ⎥App⎥
⎥A⎥⎦
其中Aij是aij的代数余子式。
一般情况,这个求逆公式只有理论的价值,在多元分析中求逆矩阵是通过消去变换来实现的,并且同时可求得该矩阵的行列式。消去变换在后面介绍。
逆矩阵的基本性质如下:
(1)AA
-1
=A-1A=I
(2)(A')-1=(A-1)'
(3)若A和C均为p阶非退化阵,则
(AC)-1=C-1A-1
(4)设A为p阶非退化阵,b和a为p维列向量,则方程:
Ab=a
的解为
-1
b=Aa
(5)A-1=A
-1
(6)若A是正交阵,由于AA'=I推知A-1=A';若A是对角阵,A=diag(a11,a22,„,app)
-1-11
且aij≠0,i=1, ,p,则A-1=diag(a11,a22, ,a-pp)。
3 矩阵的秩
设A为p×q阶阵,若存在它的一个r阶子方阵的行列式不为零,而A的一切(r+1)阶子方阵的行列式均为零,则称A的秩为r,记作rk(A)=r。它有如下基本性质:
(1)rk(A)=0,当且仅当A=0;
(2)若A为p×q阶阵,则0≤rk(A)≤min(p,q); (3)rk(A)=rk(A');
(4)rk(AB)≤min(rk(A),rk(B)); (5)rk(A+B)≤rk(A)+rk(B);
(6)若A和C为非退化阵,则rk(ABC)=rk(B)。
§3 特征根、特征向量和矩阵的迹
1 特征根和特征向量
设A为p阶方阵,则方程A-λIp=0是λ的p次多项式,由多项式理论知道必有p个根(可以有重根),记为λ1,λ2, ,λp,称为A的特征根或称特征值。
若存在一个p维向量li,使得(A-λiIp)li=0,则称li为对应于λi的A的特征向量。今后总假设li' li=1。
特征根有如下性质:
(1)若A为实对称阵,则A的特征根全为实数,故可按大小次序排列成λ1≥λ2≥ ≥λp,若λi≠λj,则相应的特征向量li与lj必正交。
(2)A和A′有相同的特征根。
(3)若A与B分别是p×q与q×p阶阵,则AB与BA有相同的非零特征根。 (4)若A为三角阵(上三角或下三角),则A的特征根为其对角元素。
-11-1(5)若λ1,λ2, ,λp是A的特征根,A可逆,则A-1的特征根为λ1,λ-2, ,λp。
2 矩阵的迹
若A是p阶方阵,它的对角元素之和称为A的迹,记为tr(A)=若A是p阶方阵,它的特征根为λ1,λ2, λp,则tr(A)=特征根和迹有如下关系: (1)tr(AB)=tr(BA) (2)tr(A)=tr(A′)
(3)tr(A+B)=tr(A)+tr(B) (4)tr(aA)=atr(A)
∑a
i=1
p
ii
;
∑λ
i=1
p
i
。
§4 二次型与正定阵
称表达式
Q=
∑∑a
i=1j=1
pp
ijxixj
为二次型,其中aij=aji是实常数;x1,x2, ,xp是p个实变量。
若A=(aij)p⨯p为对称阵,X=(x1, ,xp)' 则
p
p
Q=
∑∑a
i=1j=1
ijxixj
=X'AX
若方阵A对一切X≠0,都有X'AX>0,则称A与其相应的二次型是正定的,记为A>0;若对一切X≠0,都有X'AX≥0,则称A与二次型是非负定的,记为A≥0。
记A≥B,表示A-B≥0;记A>B,表示A-B>0。 正定阵和非负定阵有如下性质:
(1)一个对称阵是正(非负)定的当且仅当它的特征根为正(非负);
(2)若A>0,则A-1>0;
(3)若A>0,则CA>0,其中C为正数; (4)若A≥0,因它是对称阵,则必存在一个正交阵Γ,使Γ'AΓ=agdi其中λ1, ,λp为A的特征根,Γ的列向量为相应的特征向量,于是
A=ΓΛΓ'
(5)由性质1,λ1, ,λp均非负,即Λ≥0。记f(Λ)=diag(f(λ1), ,f(λp)),f(A)=
(λ1,λ2, ,λp)=Λ
Γf(Λ)Γ'
1
特别Λ2
12
=diag(λ
12
1, ,
λ
12) p
1A2
=
1ΓΛ2Γ'
称A为A的平方根;
(6)若A≥0 (>0),则存在
1A2
≥0 (>0),使得A=
1A21A2
。
§5 消去变换
在多元分析中经常要解线性方程组,求矩阵的逆和行列式,或进行某种逆推运算,通过消去变换可以达到目的。
*
设A=(aij)是n×p阶阵,若aij≠0,将A变换为A*=(aij),使得
*
aαβ
当a≠i, β≠j⎧aαβ-aiβaαjaij
⎪
当a≠i,β=j⎪-aajaij
=⎨
当β≠j,α=i⎪aiβaij
⎪a 当a=i,β=j⎩ij
即A*阵为如下形式:
-a1jaij⎡⎤
⎢*⎥* **⎢⎥⎢⎥-ai-1jaij⎢⎥
aijaij+1aij aimaij⎥ A*=⎢ai1aij aij-1aij
⎢⎥-ai+1jaij⎢⎥
* **⎥⎢*
⎢⎥-anjaij⎣⎦
其中**部分第(α,β)位置的元素是aαβ-aαjaiβaij。这个变换称为以(i, j)为枢轴
的消去变换,记作A*=Tij(A)。
消去交换具有如下性质:
(1)Tij(TijA)=A,即对A连续施行两次(i,j)消去交换,其结果仍是A不变。 (2)若i≠k,j≠l,则Tij(TklA)=Tkl(TijA),它表明Tij在某种意义下的可交换性。
§6 矩阵的分块和矩阵的微商
1 矩阵的分块
对于任意一个p×q阶矩阵A,可以用纵线和横线按某种需要将它们划分成若干块低阶的矩阵,也可以看作是以所分成的子块为元素的矩阵,称为分块矩阵,即:
⎡a11a12 a1q⎤⎢a⎥a a21222q⎥ A=⎢⎢ ⎥⎢⎥aa a⎢⎥p1p2pq⎣⎦
A12⎤⎡A
写成 A=⎢11⎥
AA22⎦⎣21
其中
A11=[aij] i=1, ,m;j=1, ,n
A12=[aij] i=1, ,m;j=n+1, ,q
A21=[aij] i=m+1, ,p;j=1, ,n A22=[aij] i=m+1, ,p;j=n+1, ,q
矩阵的分块是相当任意的,同一个矩阵可以根据不同的需要,划分成不同的子块,构成
不同的分块矩阵。
值得注意的是,在分块矩阵中,每行(列)上各元素(指原矩阵的子块)具有相同的行(列)数(指在原矩阵中的行(列)数)。
分块矩阵也满足平常矩阵的加法、乘法等运算规律。 不难证明:
若A=⎢
⎡A11⎣A21'A21⎤⎡A11
', 则A=⎢
'A22⎥⎦⎣A12'⎤A21
。 ⎥'⎦A22
若A11, A22是 方阵且是非退化阵,则
-1-1
A=A11A22-A21A11A12=A22A11-A12A22A21
2 矩阵的微商
设x=(x1, ,xp)'为实向量,y=f(x)为x的实函数。 则f(x)关于x的微商定义为:
⎡∂f⎤⎢⎥∂x1⎥
∂f(x)⎢
∆⎢ ⎥ ∂x⎢⎥
∂f⎢⎥⎢∂xp⎥⎣⎦
若
则定义
⎡x11 x1p⎤⎢⎥
X=⎢ ⎥
⎢xn1 xnp⎥⎣⎦
∂f⎤⎡∂f
⎢∂x∂x1p⎥11⎥∂f(X)⎢
∆⎢ ⎥ ∂X⎢∂f∂f⎥
⎢∂x⎥∂xnp⎦⎣n1
由上述定义不难推出以下公式:
(1)若x=(x1, ,xp)', A=(a1, ,ap)' 则
∂(x'A)
=A ∂X
∂(x'x)
=2x ∂x
(3)若x=(x1, ,xp)', B=(bij)p⨯p对称阵
(2)若x=(x1, ,xp)', 则则
∂(x'Bx)
=2Bx ∂x
(4)若y=tr(X'AX),式中X为n×p阶阵,A为n×n阶阵,则 则
∂tr(X'AX)
=(A+A')X
∂X
若A为对称阵,则
∂tr(X'AX)
=2AX
∂X
参 考 文 献
[1] 张尧庭,方开泰(1982):《多元统计分析引论》,科学出版社。 [2] 王学仁(1982):《地质数据的多变量统计分析》,科学出版社。 [3] 方开泰,潘恩沛(1982):《聚类分析》,地质出版社。 [4] [英]M. 肯德尔(1983):《多元分析》,科学出版社。
[5] Morris L. Eaton(1983), Multivariate Statistical A Vactor Space Approach.
[6] Anderson, T. W. (1984), An Introduction to Multivariate Statistical Analysis (2nd Edition), Wiley.
[7] 田中丰·垂水共之·胁本和昌编(1984):《 卜2统计解 II多变量解析编》,共立出版株式会社。
[8] 余金生、李裕伟(1985):《地质因子分析》,地质出版社。 [9] 陈希孺,王松桂(1987):《近代回归分析—原理方法及应用》,安徽教育出版社。 [10] 罗积玉、邢瑛(1987):《经济统计分析方法及预测—附实用计算机程序》,清华大学出版社。
[11] 周光亚等(1988):《多元统计分析》,吉林大学出版社。 [12] 方开泰(1989):《实用多元统计分析》,华东师范大学出版社。 [13] 孙尚拱(1990):《实用多变量统计方法与计算程序》,北京医科大学、中国协和医科大学联合出版社。
[14] 孙尚拱、潘恩沛(1990):《实用判别分析》,科学出版社。 [15] 王学仁、王松桂(1990):《实用多元统计分析》,上海科学出版社。 [16] 胡国定、张润楚(1990):《多元数据分析方法》,南开大学出版社。 [17] 张明立、于秀林(1991):《多元统计分析方法及程序——在体育科学中的应用》,北京体育学院出版社。
[18] 张尧庭等(1991):《定性资料的统计分析》,广西师范大学出版社。 [19] 于秀林(1993):《多元统计分析及程序》,中国统计出版社。 [20] 王国梁、何晓群(1993):《多变量经济数据统计分析》,陕西科学技术出版社。 [21] Alvin C.Rencher (1995), Methods of Multivariate Analysis. [22] [美]David Freedman等著,魏宗舒等译(1997):《统计学》,中国统计出版社。
第十二章 简介定性资料的统计分析
本章不是全面的介绍这方面的理论、方法和应用,而是初步反映一下这方面的主要内容,目的是展示进一步可学的知识,以便更好地解决实际问题。
§12.1 定性变量数量化
前面几章所介绍的各种统计方法,主要是研究与定量变量(或称间隔尺度变量)有关的问题,但在实际应用中,往往不可避免地要涉及到定性变量(或称名义尺度变量),例如人的性别、职业、天气状态,经济工作中选择的政策以及地层的构成类型等等,这些变量都只有各种状态的区别,而没有数量之区别。若定性变量不进入数学关系式,则会丢失信息,若要进入,又难于直接参加运算,于是从20世纪五十年代起开始发展了数量化理论,首先应用于“计量社会学”,六十年代后,逐步应用于各种学科,随着电子计算机的普及和发展,数量化理论将会在自然科学和社会科学的许多方面发挥出更大的作用。
如何对定性变量给以相应的数值描述,从而进行有关的统计分析,这就是数量化理论所研究的主要内容。
数量化理论已有专著出版,本节为了应用上的需要,仅介绍常用的0-1赋值法。 例如定性变量是性别,记为X,如此赋值:
当性别为女⎧1, 当性别为男⎧1,
X=⎨ 或X=⎨
⎩0, 当性别为女⎩0, 当性别为男
如此赋值的理由是简单,并没有任何数量大小的意义,它仅仅用来说明观察单位的特征
或属性,因此不同特性或属性的观察单位应取不同的值。
例如:天气可取晴、阴、雨三类,则用两个变量(X1,X2)表示天气,如此赋值:
当天气晴⎧(0,0),
⎪
(X1,X2)=⎨(1,0), 当天气阴
⎪(0,1), 当天气雨⎩
例如:有多种有害物污染了大气,由于有害物的结构不同,将污染物分为五类地区;甲、
乙、丙、丁、成戊将地区用4个变量(X1, X2, X3, X4)来表示,如此赋值:
甲类地区⎧(0,0,0,0),
⎪
乙类地区⎪(1,0,0,0),
⎪
(X1,X2,X3,X4)=⎨(0,1,0,0), 丙类地区
⎪(0,0,1,0), 丁类地区⎪⎪戊类地区⎩(0,0,0,1),
综上所述,推广为一般的赋值法如下:若某定性变量可取K类,则用K-1个变量表示,
如此赋值:
⎧(0,0,0, ,0), 第一类⎪
第二类⎪(1,0,0, ,0),
⎪第三类⎪(0,1,0, ,0),
(X1,X2, ,Xk-1)=⎨
第四类⎪(0,0,1, ,0),
⎪ ⎪⎪第K类⎩(0,0,0, ,1),
以上K个类的次序可以交换。
对于取K个类的定性变量,为什么用K-1个变量而不用K个变量表现?例如某定性变量可取甲、乙、丙、西四个类,可否如下赋值:
⎧(0,0,0,0), 取甲类⎪
取乙类⎪(1,0,0,0),
(X1,X2,X3,X4)=⎨
取丙类⎪(0,1,0,0),
⎪(0,0,1,0), 取丁类⎩
易知,如此赋值将使X1+X2+X3+X4 =1,不论是第几次观测,也不论定性变量取哪一类,皆使上式成立,即4个变量之和有稳定的线性关系式,知道其中任意三个就可推知另一个。
定性变量数量化后,就可以全部作为定量变量来统一处理进行预测或分类等研究。
§12.2 列联表
主要介绍二维列联表,对于三维以上的列联表只要在形式上稍加改变就能适用于高维表,原则上是一样的,只不过高维列联表符号更复杂一些,也增加些分析的难度。
1 列联表的概念
列联表讨论的主要是定性资料,此处介绍二维列联表的目的,不是将其数量化,而是直接进行分析并给出两个定性变量之间是否独立性检验。
先看一个简单例子:研讨吸烟与患肺癌的关系,这里用A表示一个人是否患肺癌,用B表示一个人是否吸烟,从一批被调查的对象中得到的统计表如下:
研讨患肺癌是否与吸烟有关?
这张统计表称为2×2列联表,表中考察两个定性变量A和B,每个变量有两类,即A分为患肺癌与未患肺癌两类,B分为吸烟与不吸烟两类,表中间的数值是频数,每一个被抽到的人,都可确定他的(AiBj)取值,比如表中数值60,表示被抽人群中吸烟又患肺癌的人数,数值32表示吸烟示患肺癌的人数。
一般2×2列联表形式如下:
其中nij(i,j=1,2)表示第i行Ai和第j列Bj的样品出现的频数,一般nij可取任意非负整数。
这是一个最简单的列联表,如果两个定性变量分别考察r和c类,则相应的列联表为r⨯c表(r和c可以不等)有如下形式:
如果一个问题涉及到很多的定性变量,相应的频数表就是一个高维列联表。
在概率统计中描述两个随机变量的相关程度是用线性相关系数,为了避免术语上的混淆,描述两个一性随机变量之间的相关性是指广义的相关性,称为关联性,两个定性随机变量之间的关联程度在某种意义上就是指的“不独立性”,它与独立的情形差距越大,就表明彼此的关系越密切,这种关系不一定是线性关系,然而在实际问题中,重要的是判断变量之间是否独立,因为不独立就意味着是关联的。如何判断是否独立有很多方法,这里仅介绍一种常用的皮尔逊拟合优度x2检验。
2×2列联表,对应一个多项分布,检验A与B是否独立,等价于检验:
H0:pij=pi.p.j
其中pij表示A为i、B为j的样品概率,pi.和p.j是相应的边缘概率,当独立性成立时, 理论频数为:npij=npi.p.j 其中n=
2
2
∑∑n
i=1j=1
ij
实际频数为:nij
运用x2检验作判定,需要知道列联表中实际频数与相应的理论频数。用估计量
nn
ˆ.j=.j代替pi.和p.j。基实际频数与理论频数有差异,这时可用其差值的大小ˆi.=i., pp
nn
来度量两个变量相关程度。相差愈大,表明H0为真的可能性愈小,即A与B无关的可能性愈小。相反差值愈小,即二愈接近,H0为真的可能性愈大,A与B之间相关的可能性愈小。为避免实际频数与理论频数的差值出现正负抵消,可采用差值的加权平方和来检验,于是给
2
出皮尔逊的拟合优度x统计量为:
2
x=
2
∑∑
i=1
=
∑∑
i=1
2
⎛nn⎫ nij-ni.j.⎪2 nn⎪⎝⎭
nn.jj=1
ni.
nn
2nn-nn2
iji..j
nnni..jj=1
2
()
它的极限分布是自由度为1的x2分布,根据给定的显著性水平a,查x2分布表得到临界值λa。若x2≥λa则拒绝H0,表示A与B之间不独立,存在相关,若x2
H0,表明A与B之间独立,不存在相关。
将前面的例子作x2检验:
计算
(106⨯60-63⨯92)2(106⨯32-43⨯92)2
x=+
106⨯63⨯92106⨯43⨯92(106⨯3-63⨯14)2(106⨯11-43⨯14)2++
106⨯63⨯14106⨯43⨯14(6360-5796)2(3392-3956)2=+
[1**********]6(318-882)2(1166-602)2++
9349263812
=0.75857+0.5775+4.98489+3.40239=9.66360
2
取显著性水平a=0.05,自由度为1,查x2分布表,临界值λa=3.84。显然x2=9.66360>3.84,表明在5%的显著性水平上,拒绝H0即说明吸烟与肺癌不独立,而是存在相关的。
如果列联表中变量间存在相关,那么如何度量变量间的相关程度?又如何从一个变量去预测另一变量呢?解决这类总是还有很多方法,已超出本书范围,不再详述,有兴趣的读者可查阅这方面的参考书。
§12.3 对数线性模型
如前所述,列联表能够反映定性变量之间的关系,但能否像定量变量那样建立起数学模型如方差分析模型、回归分析模型等以便进一步描述定性变量之间的复杂关系呢?对数线性模型和Logistic回归模型就是解决这一问题的极为有效的方法,它们从不同角度出发导出不同的处理方法。
对数线性模型,近十年来是国外实际工作者常用的方法,它的主要优点是可以把方差分析和线性模型的一些方法系数地移植过来,在概念和理解上均可进行对比,对数线性模型能够估计模型中各个参数,而这些参数值使各个变量的效应和变量间的交互作用效应得以数量化。下面即将看到这些结论。
1 模型
对数线性模型又分为很多种类型,常用的模型有:饱和模型(当变量间相互不独立时),非饱和型(变量间相互独立),谱系模型(包含高阶效应)等。
下面从2×2的频数表与概率表出发,推导对数线性模型: (频数表)(概率表)
将概率取对数后进行分解处理,使处理后的变量有较好的数学、统计的性质。
⎛pij⎫ ⎪ μij=lnpij=lnpi.p.j pi.p.j⎪⎝⎭
pij
=lnpi.+lnp.j+ln
pi.p.j
记Ai=lnpi., Bj=lnp.j, (AB)ij=ln由上式可写成
pijpi.p.j
μij=Ai+Bj+(AB)ij i,j=1.2
显然上式的结构类似于两因子有交互作用,各因子均为二水平的方差分析模型,于是令
μi.=∑μij, μ.j=∑μij, μ=∑∑μij
j=1
i=1
i=1j=1
2222
然后再进行平均,对i, j =1,2
i.=μi., .j=μ.j ..=μ..
记
1
21212
αi=μi.-μ..
βj=j.-..
λij=μij-i.-.j+..
则有关系式:
⎧⎪
⎪μij=..+αi+βj+λij⎪22⎪
i,j=1,2 ⎨且αi=0, βj=0,
j=1⎪i=1
⎪22⎪λ=λ=0 ⎪i=1ijj=1ij⎩
可见通过上边分解处理,可以完全化成与方差分析模型有同样的结构,因此借助于方差分析的术语,上式中μ..表示“总平均效益”,αi表示A属性的“主效应”,βj表示B属性
∑∑
∑∑
的“主效应”,λij表示,A, B的 “交互作用效应”,直观可以理解当交互作用效应为0,即等价于A、B独立。上式模型称为对数线性模型的饱和模型,当λij=0时,称为非饱和模型。
n.jni.
ˆp=,.j,这时就可以看
nnn
到对数线性模型是将列联表上每个单元的频数作为因变量,表上所有变量作为自变量,建立各个自变量的效应与每个单元频数的对应之间的函数关系。因而可以用它分析列联表上的各个变量的关系。主效应αi或βj若大于0,表明效应为正;若小于0,表明效应为负。αi是
ˆij=在实际应用时概率可用其估计量代替,即p
ˆi.=,p
第一个变量的第i个水平对总平均效应μ..的增减量;βj是第二个变量的第j个水平对总平均效应μ..的增减量,λij代表变量1和变量2在各自的第i个水平和第j个水平之间交互作用效应,是其交互作用对总平均效应的增减量。若λij
nij
1
i.=
2
∑
1μij=
2j=1
n
∑
j=1
2
(ln
nijn
) „ 第i行频数对数的平均
1μ.j=
2
∑
1μij=
2i=1
n
n
∑
i=1
2
(ln
nijn
) „ 第j列频数对数的平均
11
μ..=μ..=
44
∑∑
i=1
1
μij=
4j=1
2
∑∑(ln
i
j
22
nijn
) „ 各个观测值对数的总平均即总平均效
应对本章前面的例子,按上述模型估计各效应参数。
各单元的频数对数表:
计算:
α1=1.-..=2.5964-2.5141=0.0823
α2=2.-..=2.4318-2.5141=-0.0823
β1=.1-..=3.2800-2.5141=0.7659
β2=.2-..=1.7482-2.5141=-0.7659
主效应估计值:
计算:
λ12=μ12-1.-.2+..=1.0986-2.5964-1.7482
+2.5141=3.6127-4.3426=-0.7319
λ21=μ21-μ2.-μ.1+μ..=2.4657-2.4318-3.2800
+2.5141=4.9798-5.7188=-0.7320
λ11=μ11-1.-.1+..=4.0943-2.5964-3.2800
+2.5141=6.6083-0.8319=0.8319
λ22=μ22-2.-.2+..=2.3979-2.4318-1.7482
+2.5141=4.9120.1800=0.7320
变量间交互作用效应估计值:
λ11=0.8319 λ12=-0.7399 λ22=0.7320
λ21=-0.7317
主效应大于0,表明效应为正,如α1=0.0832>0是因为患肺癌比未患肺癌的人多;主效应小于0,表明效应为负,如β2=-0.76599
§12.4 Logistic回归
对数线性模型是将列联表中每格的概率(或理论频数)取对数后分解参数获得的,Logistic回归模型是将概率比取对数后,再时行参数化而获得的。研究概率比这样的量在不少问题中是常常遇到的,当列联表中因变量是一个多级分类的变量时,就需要考虑两两比较的情况。
Logistic回归要解决的问题,类似于普通回归所要解决的很大一类问题。比如在医药行业中,因变量y取0, 1, „, g这g+1个不同的值,y=0表示正常情况类型,y=1….,g表示不同用药后的反应,显然它与药的剂量x1 ,性别x2,年龄x3,体重x4,血压x5,„等等自然变量有关,这里因变量是定性的,自变量有定性的也有定量的,问这些自变量对一个定性变量的关系是否独立?不独立又会具有什么形式的联系?是线性的还是非线性的等等。
1 Logit变换
为了给出Logistic回归模型,先介绍Logit变换。在现实问题中,人们常常要研究某一事件A发生的概率p以及p值的大小与某些因素的关系,但由于p对x的变化在p = 0或p =1的的附近是很敏感的,或说是缓慢的比如像可靠系统,可靠度p已经是0.988了,即使再改善条件和系统结构,它的可靠度增长只能在小数点后面的第三位或第四位,于是自然希望寻找一个p的函数θ(p),使它在p = 0或p =1附近变化幅度较大,而且函数的形式也不要太复杂,根据数学上导数的意义,提出用时希望p = 0或p =1,
dθ(p)
来反映θ(p)在p附近的变化是很合适的,同dp
dθ(p)
有较大的值,因此取 dp
dθ(p)111
==+ dpp(1-p)p1-p
即
θ(p)=ln
p
1-p
称上式为Logit变换。
p
由于θ(p)=ln,因此p也可用θ表示:
1-p
eθ
p= θ
1+e
如果θ是某些自变量x1, „, xk,的线性函数
∑a
i=1
k
k
xi,则p就x1, „, xk的下述函数:
p=
ei=1
∑aixi∑aixi
k
k
1+i=1
显然,如果如果p对x不是线性关系,则θ对x就可以是线性关系了,这就是Logit变换带来的方便。
2 Logistic回归模型
设因变量y为二值定性变量,用0,1表示取两个不同的状态,y=1的概率p=P(y=1)是我们要研究的对象。如果有很多因素影响y的取值,这些因素就是自变量,记为x1,„, xk,这其中既有定性变量也有定量变量,最重要的一个条件是
pln=a0+a1x1+ +akxk 1-p适合于上式线性函数称为Logistic线性回归。如果有已知函数g(x1, ,xk),其中含有若干待定的参数,且
pln=g(x1, ,xk) 1-p称上式为非线性Logistic回归。
Logistic线性模型相当于广义线性模型,因此可以系统地应用线性模型的方法,在处理时就比较方便。
如果某一事件A发生的概率p依赖于一些自变量x1, ,xk(定性定量均可),对x=(x1, ,xk)'观测了m组结果,在第a组中,共试验了na次,A发生了ra次,于是A发
ˆa=生的概率pa可用p
ra
来估计。假定pa适合于Logistic线性回归关系式,即 napaln=β0+β1xa1+ +βkxak a=1, ,m 1-pa
ˆa代入上式中pa就有关系式: 其中xai表示xi在第a组所取的值。用p
ln
ˆap
=β0+β1xa1+ +βkxak+εa a=1, ,m ˆa1-p
其中εa是随机误差项。记
⎡1x11 x1k⎤
⎢1x⎥ xˆap212k⎥, X=⎢ ya=ln
m⨯(k+1)⎢ ⎥ˆ 1-pa
⎢⎥⎣1xml xmk⎦
y=(y1, ,ym)',假定E(εa)=0,V(εa)=va,ε1, ,εm相互独立,于是就有
⎡β0⎤
0⎤⎡v1⎢⎥
β⎥
Ey=X⎢1⎥∆Xβ Var(y)=V=⎢ ⎢⎥⎢ ⎥
⎢vm⎥⎢⎥⎣0⎦β⎣k⎦
ˆ 这就是典型的线性模型,因此β的最小二乘估计(因V不是单位阵,应该相应地加权)β
有公式:
ˆ=(X'βV-1X)-1X'V-1y
这样就求得了β的估计值。要讨论某些xi是否对A发生的概率有无影响,也即要检验xi相应的回归系数βi=0这一假设是否成立,这时要搬用线性模型中已知的结论时,必须知道y
是否服从正态分布以及V的估计,数学上可以证明:
lnˆa=v
ˆa⎛⎫ppa1
的渐近分布为正态的,即N ln, 1-pnp(1-p)⎪⎪V是的va的估计值为ˆa1-paaaa⎭⎝
11
ˆˆnapa(1-pa)
如果在m组试验结果中,有的ra=0或ra=na,此时
ln
ˆapra
=ln ˆa1-pna-ra
会取-∞或∞的值,ya就不是一个有限的值,上述方法就会行不通,于是要进行修正,修正
rapa
的目的是使ln尽可能接近ln,可以证明下面的修正是合理的。
na-ra1-pa
ra+0.5
a=1, ,m
na-ra+0.5(na+1)(na+2)
ˆa=v a=1, ,m
na(ra+1)(na-ra+1)Za=ln
上式Za称为经验的Logistic变换,相应的线性模型也称为经验的Logistic回归模型:
⎧⎡EZ1⎤⎪⎢⎥EZ= ⎪⎢⎥=Xβ⎪⎢⎪⎣EZm⎥⎦
⎨
ˆ0⎤⎡v1⎪
⎥⎪Var(Z)=⎢ ⎢⎥⎪
⎢ˆm⎥v⎪⎣0⎦⎩
ˆ的估计值如下: ˆa为修正后的表达式,β其中Za,v
ˆ=(X'Vβ
-11-2V2
X)-1X'V
-
11
-2V2
~。 用加权最小二乘估计回归系数时,第a组的权系数是va
引入经验Logistic变换后,Logistic回归就和普通的回归一样了,而且用的是加权最小
二乘,在普通回归中遇到的一些检验的问题,如回归系数的显著性检验、共线性问题、„„等等,在这里都有同样的意义,而且处理方法也相似,经验的Logistic回归是Logistic回归中最常用的方法。
3 实例
对本章所例举的实例,讨论吸烟与肺癌的关系,作经验的Logistic变换并再次检验吸烟与肺癌是否相关。
将吸烟人作为一类,不吸烟人作为一类,此时m =2, n1=92, n2=14, r1=60, r2=3。 作经验Logistic变换:
r+0.560.5
Z1=ln1=ln=0.6214
n1-r1+0.532.5
-
1
2
r2+0.53.5
=ln=-1.1896
n2-r2+0.511.5(n1+1)(n1+2)(92+1)(92+2)8742ˆ1=v===0.047204
n1(r1+1)(n1-r1+1)92⨯61⨯33185196Z2=ln
ˆ2=v
(n2+1)(n2+2)(14+1)(14+2)240
===0.3571428
n2(r2+1)(n2-r2+1)14⨯4⨯12672
为了作假设检验,设吸烟人患肺癌的概率是p1,未患肺癌的概率就是1-p1,不吸烟的人
患肺癌的概率是p2,未患肺癌的概率就是1-p2,作Logit变换:
pp
θ1=ln1, θ2=ln2
1-p11-p2设θ2=θ,则
θ1=θ2+(θ1-θ2)=θ+∆
因此患肺癌是否与吸烟有关,就归结为检验
H0:∆=0
是否成立。
ˆ1+vˆ2),由于Z1,Z2可以看成是相互独立的正态变量。因此,H0成立时,Z1-Z2~N(0,v
即
Z1-Z21+v2v
计算
~N(0,1)
Z1-Z21+v2v
=
1.81100.0472+0.357
=
1.8110.4043
=
1.8110
=2.8483
0.6358
取a=0.05,查N(0,1)表得临界值λa=0.834。
2.8483>0.834,应拒绝H0,即认为吸烟与不吸烟对于患肺癌是有影响的。
综上所述,对于列联表,既可以用对数线性模型也可以用Logistic模型,要根据实际问题的需要来选择模型。
定性资料统计分析的内容是很丰富的,这里不再详述,有兴趣的读者可查有关资料。
附录:矩阵代数
矩阵和行列式是研究多元统计分析的重要工具,这里针对本书的需要,对有关矩阵的知识作复习,同时再简单的补充矩阵的分块和矩阵的微商。(因为在部分章节中要用到)。
§1 矩阵及基本运算
1 矩阵的定义
将n×p个实数,a11, a12, „, a1p, a21, a22, „, a2p, „, an1, an2, „, anp排成如下形式的矩形表,记为A
⎡a11a12 a1p⎤⎢a⎥a a21222p⎥ A=⎢⎢ ⎥⎢⎥aa a⎢⎥n1n2np⎣⎦
则称A为n⨯p阶矩阵,一般记为A=(aij)n⨯p,其中aij是A中元素,它也可以是复数,但本书的aij均为实数。当n = p时,称A为n阶方阵;若p =1,A只有一列,称A为列向量,
记为
⎡a11⎤⎢a⎥a=⎢21⎥
⎢ ⎥⎢⎥⎣an1⎦
当n=1时,A只有一行,称A为行向量,记为
a'=(a11,a12, ,a1p)
当A为n阶方阵,称a11,a22, ,ann为A的对角线元素,其它元素称为非对角元素,若方阵A的非对角元素全为0,称A为对角阵,记为
O⎤⎡a11
⎢⎥a22⎥=diag(a,a, ,a) A=⎢1122nn⎢⎥ ⎢⎥Oann⎦⎣
进一步若a11 = a22 =„= ann =1,称
O⎤⎡1
⎢⎥1⎥ A=⎢⎢⎥ ⎢⎥O1⎣⎦
为n阶单位阵,记为In或A=I(I的阶数从上下文可明确)。
若A为n×p阶矩阵,它的转置A′是阶p×n阶矩阵,即
⎡a11a21 an1⎤⎢a⎥a a1222n2⎥ A'=⎢⎢ ⎥⎢⎥aa a⎢2pnp⎥⎣1p⎦
若A是方阵,且A′=A,则称A为对称阵;若方阵A=(aij)的元素aij = 0,对一切i<j成立,则称
⎡a11⎢aA=⎢21
⎢ ⎢⎣an1
a22 an2
O⎤⎥⎥ ⎥
⎥
ann⎦
为下三角阵;若A′为下三角阵,则称A为上三角阵。
2 矩阵的运算
若A与B是n×m阶阵,则A与B的和定义为:
A+B=(aij+bij)n⨯m
若a为一常数,它与A的积定义为
aA=(aaij)n⨯m
若A=(aik)p⨯q,B=(bkj)q⨯r,则A与B的积定义为: AB=(
q
∑a
k=1
ikkj
b)p⨯r
在一般情况下AB≠BA。
容易验证上述运算符合下面运算规律:
A+(-1)A=0
(AB)′=B′A′ (A′)′=A
(A+B)′=A′+B′ A(BC)=(AB)C A(B+C)=AB+AC (A+B)C=AC+BC AI=IA=A
a(A+B)=aA+aB a(AB)=(aA)B=A(aB)
若A为方阵,满足A′A=AA′=I,则称A为正交阵。
§2 行列式、逆矩阵和矩阵的秩
1 行列式
一个p阶方阵A=(aij)p⨯p对应一个数,记为|A|
⎡a11⎢a21
A=⎢
⎢ ⎢⎢⎣ap1
=
j1j2 jp
a12 a1p⎤
a22 a2p⎥⎥ ⎥
⎥
ap2 app⎥⎦
τ(j1 jp)
∑(-1)
j1 jp
a1j1,a2j2, ,apjp
称|A|为A的行列式或记为detA。这里∑表示对所有p元排列求和。τ(j1 jp)表示,若jt>js,则称这两个数组成一个j1 jp的逆序数(在一个排列j1j2 ,jt, ,js jp)逆序。一个排列中逆序的总数称为此排列的逆序数)。
不难证明:
A=
∑(-1)
j=1
p
i+j
,aijMij
其中Mij是在A中去掉第i行、第j列而形成的p-1阶方阵。称(-1)i+jMij为aij的代数余子式,通常记为Aij。
直接由行列式的定义来计算行列式是很麻烦的,通常利用行列式的如下一些性质可以简化行列式的计算:
(1)若A的某行(或列)为零,则|A|=0; (2)|A|=|A′|;
(3)将A的某行(或列)乘以数a所得矩阵的行列式等于a|A|; (4)若A的两行(或列)相同,则|A|=0;
(5)若将A的两行(两列)互换所得矩阵的行列式等于-|A|;
(6)若将A的某一行(或列)乘上一个常数后加到另一行相应的元素上,所得矩阵的行列式不变,仍等于|A|。
2 逆矩阵
设A为p阶方阵,若|A|≠0,则称A是非退化阵或称非奇异阵,若|A|=0,则称A是退化阵或称奇异阵。
若A是p阶非退化阵,则存在唯一的矩阵B,使得AB=Ip,B称为A的逆矩阵,记为B=A-1,不难证明:
A-1
⎡A11⎢⎢A⎢A12=⎢A⎢⎢ ⎢A1p⎢⎢⎣A
A21
AA22A A2pA
Ap1⎤
⎥A⎥Ap2⎥
⎥A⎥ ⎥App⎥
⎥A⎥⎦
其中Aij是aij的代数余子式。
一般情况,这个求逆公式只有理论的价值,在多元分析中求逆矩阵是通过消去变换来实现的,并且同时可求得该矩阵的行列式。消去变换在后面介绍。
逆矩阵的基本性质如下:
(1)AA
-1
=A-1A=I
(2)(A')-1=(A-1)'
(3)若A和C均为p阶非退化阵,则
(AC)-1=C-1A-1
(4)设A为p阶非退化阵,b和a为p维列向量,则方程:
Ab=a
的解为
-1
b=Aa
(5)A-1=A
-1
(6)若A是正交阵,由于AA'=I推知A-1=A';若A是对角阵,A=diag(a11,a22,„,app)
-1-11
且aij≠0,i=1, ,p,则A-1=diag(a11,a22, ,a-pp)。
3 矩阵的秩
设A为p×q阶阵,若存在它的一个r阶子方阵的行列式不为零,而A的一切(r+1)阶子方阵的行列式均为零,则称A的秩为r,记作rk(A)=r。它有如下基本性质:
(1)rk(A)=0,当且仅当A=0;
(2)若A为p×q阶阵,则0≤rk(A)≤min(p,q); (3)rk(A)=rk(A');
(4)rk(AB)≤min(rk(A),rk(B)); (5)rk(A+B)≤rk(A)+rk(B);
(6)若A和C为非退化阵,则rk(ABC)=rk(B)。
§3 特征根、特征向量和矩阵的迹
1 特征根和特征向量
设A为p阶方阵,则方程A-λIp=0是λ的p次多项式,由多项式理论知道必有p个根(可以有重根),记为λ1,λ2, ,λp,称为A的特征根或称特征值。
若存在一个p维向量li,使得(A-λiIp)li=0,则称li为对应于λi的A的特征向量。今后总假设li' li=1。
特征根有如下性质:
(1)若A为实对称阵,则A的特征根全为实数,故可按大小次序排列成λ1≥λ2≥ ≥λp,若λi≠λj,则相应的特征向量li与lj必正交。
(2)A和A′有相同的特征根。
(3)若A与B分别是p×q与q×p阶阵,则AB与BA有相同的非零特征根。 (4)若A为三角阵(上三角或下三角),则A的特征根为其对角元素。
-11-1(5)若λ1,λ2, ,λp是A的特征根,A可逆,则A-1的特征根为λ1,λ-2, ,λp。
2 矩阵的迹
若A是p阶方阵,它的对角元素之和称为A的迹,记为tr(A)=若A是p阶方阵,它的特征根为λ1,λ2, λp,则tr(A)=特征根和迹有如下关系: (1)tr(AB)=tr(BA) (2)tr(A)=tr(A′)
(3)tr(A+B)=tr(A)+tr(B) (4)tr(aA)=atr(A)
∑a
i=1
p
ii
;
∑λ
i=1
p
i
。
§4 二次型与正定阵
称表达式
Q=
∑∑a
i=1j=1
pp
ijxixj
为二次型,其中aij=aji是实常数;x1,x2, ,xp是p个实变量。
若A=(aij)p⨯p为对称阵,X=(x1, ,xp)' 则
p
p
Q=
∑∑a
i=1j=1
ijxixj
=X'AX
若方阵A对一切X≠0,都有X'AX>0,则称A与其相应的二次型是正定的,记为A>0;若对一切X≠0,都有X'AX≥0,则称A与二次型是非负定的,记为A≥0。
记A≥B,表示A-B≥0;记A>B,表示A-B>0。 正定阵和非负定阵有如下性质:
(1)一个对称阵是正(非负)定的当且仅当它的特征根为正(非负);
(2)若A>0,则A-1>0;
(3)若A>0,则CA>0,其中C为正数; (4)若A≥0,因它是对称阵,则必存在一个正交阵Γ,使Γ'AΓ=agdi其中λ1, ,λp为A的特征根,Γ的列向量为相应的特征向量,于是
A=ΓΛΓ'
(5)由性质1,λ1, ,λp均非负,即Λ≥0。记f(Λ)=diag(f(λ1), ,f(λp)),f(A)=
(λ1,λ2, ,λp)=Λ
Γf(Λ)Γ'
1
特别Λ2
12
=diag(λ
12
1, ,
λ
12) p
1A2
=
1ΓΛ2Γ'
称A为A的平方根;
(6)若A≥0 (>0),则存在
1A2
≥0 (>0),使得A=
1A21A2
。
§5 消去变换
在多元分析中经常要解线性方程组,求矩阵的逆和行列式,或进行某种逆推运算,通过消去变换可以达到目的。
*
设A=(aij)是n×p阶阵,若aij≠0,将A变换为A*=(aij),使得
*
aαβ
当a≠i, β≠j⎧aαβ-aiβaαjaij
⎪
当a≠i,β=j⎪-aajaij
=⎨
当β≠j,α=i⎪aiβaij
⎪a 当a=i,β=j⎩ij
即A*阵为如下形式:
-a1jaij⎡⎤
⎢*⎥* **⎢⎥⎢⎥-ai-1jaij⎢⎥
aijaij+1aij aimaij⎥ A*=⎢ai1aij aij-1aij
⎢⎥-ai+1jaij⎢⎥
* **⎥⎢*
⎢⎥-anjaij⎣⎦
其中**部分第(α,β)位置的元素是aαβ-aαjaiβaij。这个变换称为以(i, j)为枢轴
的消去变换,记作A*=Tij(A)。
消去交换具有如下性质:
(1)Tij(TijA)=A,即对A连续施行两次(i,j)消去交换,其结果仍是A不变。 (2)若i≠k,j≠l,则Tij(TklA)=Tkl(TijA),它表明Tij在某种意义下的可交换性。
§6 矩阵的分块和矩阵的微商
1 矩阵的分块
对于任意一个p×q阶矩阵A,可以用纵线和横线按某种需要将它们划分成若干块低阶的矩阵,也可以看作是以所分成的子块为元素的矩阵,称为分块矩阵,即:
⎡a11a12 a1q⎤⎢a⎥a a21222q⎥ A=⎢⎢ ⎥⎢⎥aa a⎢⎥p1p2pq⎣⎦
A12⎤⎡A
写成 A=⎢11⎥
AA22⎦⎣21
其中
A11=[aij] i=1, ,m;j=1, ,n
A12=[aij] i=1, ,m;j=n+1, ,q
A21=[aij] i=m+1, ,p;j=1, ,n A22=[aij] i=m+1, ,p;j=n+1, ,q
矩阵的分块是相当任意的,同一个矩阵可以根据不同的需要,划分成不同的子块,构成
不同的分块矩阵。
值得注意的是,在分块矩阵中,每行(列)上各元素(指原矩阵的子块)具有相同的行(列)数(指在原矩阵中的行(列)数)。
分块矩阵也满足平常矩阵的加法、乘法等运算规律。 不难证明:
若A=⎢
⎡A11⎣A21'A21⎤⎡A11
', 则A=⎢
'A22⎥⎦⎣A12'⎤A21
。 ⎥'⎦A22
若A11, A22是 方阵且是非退化阵,则
-1-1
A=A11A22-A21A11A12=A22A11-A12A22A21
2 矩阵的微商
设x=(x1, ,xp)'为实向量,y=f(x)为x的实函数。 则f(x)关于x的微商定义为:
⎡∂f⎤⎢⎥∂x1⎥
∂f(x)⎢
∆⎢ ⎥ ∂x⎢⎥
∂f⎢⎥⎢∂xp⎥⎣⎦
若
则定义
⎡x11 x1p⎤⎢⎥
X=⎢ ⎥
⎢xn1 xnp⎥⎣⎦
∂f⎤⎡∂f
⎢∂x∂x1p⎥11⎥∂f(X)⎢
∆⎢ ⎥ ∂X⎢∂f∂f⎥
⎢∂x⎥∂xnp⎦⎣n1
由上述定义不难推出以下公式:
(1)若x=(x1, ,xp)', A=(a1, ,ap)' 则
∂(x'A)
=A ∂X
∂(x'x)
=2x ∂x
(3)若x=(x1, ,xp)', B=(bij)p⨯p对称阵
(2)若x=(x1, ,xp)', 则则
∂(x'Bx)
=2Bx ∂x
(4)若y=tr(X'AX),式中X为n×p阶阵,A为n×n阶阵,则 则
∂tr(X'AX)
=(A+A')X
∂X
若A为对称阵,则
∂tr(X'AX)
=2AX
∂X
参 考 文 献
[1] 张尧庭,方开泰(1982):《多元统计分析引论》,科学出版社。 [2] 王学仁(1982):《地质数据的多变量统计分析》,科学出版社。 [3] 方开泰,潘恩沛(1982):《聚类分析》,地质出版社。 [4] [英]M. 肯德尔(1983):《多元分析》,科学出版社。
[5] Morris L. Eaton(1983), Multivariate Statistical A Vactor Space Approach.
[6] Anderson, T. W. (1984), An Introduction to Multivariate Statistical Analysis (2nd Edition), Wiley.
[7] 田中丰·垂水共之·胁本和昌编(1984):《 卜2统计解 II多变量解析编》,共立出版株式会社。
[8] 余金生、李裕伟(1985):《地质因子分析》,地质出版社。 [9] 陈希孺,王松桂(1987):《近代回归分析—原理方法及应用》,安徽教育出版社。 [10] 罗积玉、邢瑛(1987):《经济统计分析方法及预测—附实用计算机程序》,清华大学出版社。
[11] 周光亚等(1988):《多元统计分析》,吉林大学出版社。 [12] 方开泰(1989):《实用多元统计分析》,华东师范大学出版社。 [13] 孙尚拱(1990):《实用多变量统计方法与计算程序》,北京医科大学、中国协和医科大学联合出版社。
[14] 孙尚拱、潘恩沛(1990):《实用判别分析》,科学出版社。 [15] 王学仁、王松桂(1990):《实用多元统计分析》,上海科学出版社。 [16] 胡国定、张润楚(1990):《多元数据分析方法》,南开大学出版社。 [17] 张明立、于秀林(1991):《多元统计分析方法及程序——在体育科学中的应用》,北京体育学院出版社。
[18] 张尧庭等(1991):《定性资料的统计分析》,广西师范大学出版社。 [19] 于秀林(1993):《多元统计分析及程序》,中国统计出版社。 [20] 王国梁、何晓群(1993):《多变量经济数据统计分析》,陕西科学技术出版社。 [21] Alvin C.Rencher (1995), Methods of Multivariate Analysis. [22] [美]David Freedman等著,魏宗舒等译(1997):《统计学》,中国统计出版社。