《非参数统计》课件 易丹辉 第6-10章 二维列联表- 非参数回归

第6章二维列联表

6.1列联分析的基本问题在医学统计中，对变量较多而病例较少的数据进行处理时，可以先使用列联表寻找与被关注变量显著相关的变量。首先看一个例子，对患者按病程分组，521例患者按病程长短分组得到表6-1的结果。病程小于1年1至5年5至10年10至15年15年以上合计人数（人）9087134104106521百分比（%）17.2716.725.7219.9620.35100表6-1从表6-1可以看出，521例患者病程分布相对比较均匀，5至10年的人数略多。研究还希望了解这些患者病情的严重程度，可以按病情严重程度分组，得到表6-2。病情人数（人）百分比（%）较轻529.98一般较重合计16330652131.2958.73100表6-2从表6-2可以看出，521例患者中病情严重的人数最多，占近60%。但是，病情严重是否就是病程长的？它们之间有什么关系？分布规律如何？将上面两张表的数据进行列联分析，可以得到表6-3的结果。从表6-3可以得到比表6-1和表6-2更多的信息。病程10年以上的患者总的来说，病情较重。210名病程10年以上患者中病情较重的有165人，占78.57%；在病情较重的306名患者中，病程在10年以上的有165人，占53.93%。表6-3中有两个变量，分别以X和Y表示病程和病情。病程病情小于1年1至5年5至10年10至15年15年以上合计较轻1618135052一般3733532911163较重3736687095306合计9087134104106521列联表中的X和Y是否独立，是二维列联表独立性检验要解决的问题。如果两者独立就没有必要做过多的分析；如果不独立，那么两者相依程度有多大？X影响Y和Y影响X的程度是否相同？这是二维列联表相关性度量要解决的问题。如果两个变量没有因果关系，或研究中不考察是X影响Y，还是Y影响X，称为对称关系。这种情况下哪个变量为X，哪个为Y均可，两者地位平等。但若两个变量之间存在因果关系，并且要研究这种关系，则需将考察或预测的变量记作Y，称为因变量，常列入表的横行，如表6-3中的病情；另一变量称为自变量，记作X，列入表的纵列，如表6-3中的病程。这种情况是考察病程对病情的影响，称为不对称关系。列联表可以清楚地反映在X变量条件下，Y的次数分布情况。因此，列联表又称作条件次数表。表的最下端是每列的总次数，称为行边缘次数，表的最右列是每行的总次数，称为列边缘次数。表中的次数，称为条件次数。表示在自变量每个条件下，因变量各个值的数目。例如，在表6-1中，病程小于1年的有90人，这是边缘次数，从表6-3可知，其中病情较轻的有16人，这是条件次数。由表6-3中可以看出，病程确实对病情有影响。不同病程的患者,其病情有所不同。病程长的患者，更多的病情较重，而病程较短的患者，相对来说，则更多的倾向于病情较轻。表6-3是一个二维联列表，只有两个变量，变量病情有三个类别，是三行，变量病程有五个类别，是五列，构成二维的3×5列联表，最简单的列联表是2行和2列,称为2×2列联表，也称为四格表。

合计合计表6-4：从条件次数表虽然可以知道在X条件下，Y变量值的次数，但难于比较不同条件下的次数分布，因为作为基数的边缘次数不相同。如表6-3中，病程小于1年的患者有16人的病情较轻，而病程在1至5年的患者有18人，这是否表明，病程较长的患者病情倾向于较轻？观察边缘次数发现，病程1年和1至5年总人数不同，比较的基数不同，因而不宜作出结论。为了能在相同的基础上比较，使列联表的数据提供更多的信息，可以将绝对次数转化成以百分数表示的相对次数，即将条件次数变为百分数。这样的表称为条件百分表，如表6-5。病程病情小于1年1至5年5至10年10至15年15年以上合计较轻30.77%34.62%25%9.61%0100%一般22.7%20.25%32.51%17.79%6.75%100%较重12.09%11.76%22.22%22.88%31.05%100%合计17.27%16.7%25.72%19.96%20.35%100%表6-5在很多时候，研究的是不对称关系，目的是了解自变量X对因变量Y的影响。因此，条件百分数多按自变量X的方向计算，如表6-3。研究病程对病情的影响，沿自变量X的方向计算百分数，如表6-5。这一结果表明在不同病程水平下病情的变化情况。从表6-5可以看出，病程不同，病情的百分数分布也不同，这就是病程长短对病情的影响。有时由于某种原因可能使因变量在样本内的分布不能代表其在总体内的分布，例如为满足资料分析的需要，抽样时扩大了因变量的某个值的数目，使其样本内的分布不同于总体中的分布。这时，以自变量的方向计算百分数会歪曲数据的结果，需要按因变量的方向计算。条件百分表比条件次数表能够提供更多的信息，因此较为有用。但当r×c很大时，百分数的个数会很多，不容易分析两个变量之间的关系。因此，在列联表的相关测量中有许多更实用的方法可供选择研究。

6.2.1独立性检验和齐性检验1.独立性检验（1）基本方法二维列联表的独立性检验是检验行变量（Y）与列变量（X）是否独立。表6-6为两个随机变量的联合分布表，其中每个格子是联合概率，横行合计是X的边缘分布，纵列合计是Y的边缘分布。建立假设组

该假设组用概率的语言描述为

合计合计表6-6

（2）应用【例6.1】检验表6-3数据行变量与列变量之间的独立性。

2.齐性检验（1）基本方法二维列联表的齐性检验要检验的目标是：给定列，条件行分布是否相同；或者给定行，条件列分布是否相同。

如果满足齐性，则应该有对于给定列，条件行分布都相同，即表6-7的每行分别相等。

合计11表6-7建立二维列联表的假设组：

在原假设成立的条件下，对于表6-7的第i行都相等，记第i行的值为

这说明，二维列联表的齐性检验本质上是独立性检验。因此采用的检验统计量和计算p值的方法都用独立性检验相同。（2）应用【例6.2】对表6-3数据进行齐性检验。对表6-3数据进行齐性检验，即检验在不同的病程条件下，病情的条件分布是否相同。建立假设组：

3.独立性检验和齐性检验独立性检验和齐性检验的检验统计量完全相同，但两者的统计意义不同。独立性检验要检验行变量与列变量是否独立，两个变量地位平等，没有考虑因果关系。齐性检验要检验的是条件分布的齐一性，待检验的两个变量地位是不平等的，包含着因果关系，一般作为条件的变量是自变量，例如病程；另一个变量是因变量，例如病情。在应用中，应根据实际问题选择独立性检验还是选择齐性检验。

YX

2.皮尔逊（Pearson）列联相关系数

当列联表的行和列独立时，C系数的值为0；当列联表的行和列不独立时，C系数不会随着r或c的增大而增大，它的值永远小于1。对列联相关C系数可以进行显著性检验。但检验不是利用C系数的抽样分布，而是利用检验统计量Q。因为在计算C系数值的过程中先计算了Q的值，Q值可以作为C系数显著性的一个简单而合适的指标。即只要检验了Q的显著性，就等于检验了由Q计算的C的显著性。检验方法就是前面介绍的独立性检验，这里不再赘述。列联相关C系数的局限性一般而言，一个令人满意的相关系数至少应该满足两个特点：变量完全相关时，系数为1；变量完全不相关时，系数为0。但C不是理想的相关系数。一是变量完全相关时，C值并不等于1。二是两个列联相关C系数不能比较，除非是从同样大小的列联表获得的数据。

三是C系数不能与其它相关系数进行比较，例如斯皮尔曼秩相关系数和肯德尔秩相关系数等。3.克拉默V相关系数（Cramer’sVCoefficientofAssociation）

当行和列完全不相关时，V=0；当行和列完全相关时，V=1。4.应用【例6.3】分析青少年犯罪行为与家庭状况的相关性。数据如下表.青少年行为家庭状况合计离异家庭和睦家庭犯罪178272450未犯罪38502540合计216774990

6.3PRE准则下的二维列联表相关性度量当研究X和Y之间的关系，目的是为了用X去预测Y时，往往希望这一预测的效果是很好的，但事实上难免会出现误差。一般来说，两个变量之间的关系越强，以一个变量预测另一个变量的误差就越小，也就是减少的预测误差越多。换言之，消减的误差有多少，也可以反映变量之间的相关程度强弱。因此，在度量两个变量的相关程度时，可以引进消减误差比例的概念。相关性的测量方法有很多种，凡是其统计值具有消减误差比例（ProportionateReductioninError）的意义，均称为PRE测量法。

其公式定义为

PRE的值可以作为度量相关性大小的值。

6.3.1Lambda相关测量法

1.非对称形式的Lambda相关测量当研究的两个变量间存在某种因果关系，自变量X影响因变量Y的变化，而Y不会影响X，这种情况称为非对称关系。例如，研究出生时间和身体发育的关系，一般认为出生时间对身体发育会有影响，而身体发育不会影响出生时间。

根据PRE的定义可得下式

【例6.4】出生季节与身体发育相关程度的分析出生时间对身体发育不够正常有多大的影响，为此进行了调査。对两组身体发育不够正常的人的调査结果如表6-10。身体发育（Y）出生季节（X）合计春季夏季秋季冬季较低1829181277很低2013162069合计38423432146表6-10

2.对称形式的Lambda测量

【例6.5】某市随机抽取100位居民调查得到表6-11的数据，考察性别和是否愿意看中医的态度之间相关情况。态度（Y）性别（X）总数男女很愿意103040不大愿意401050不愿意10010总数6040100表6-11

如果仅研究性别影响态度的程度，可以算得

这个计算结果表明，以性别X预测态度Y，可以消减40%的误差。从表6-11可以看出，性别对看中医的态度是有影响的，男性的态度基本倾向不大愿意，而女性则较多的认为很愿意。3.Lambda相关测量法的特点

第三，由于Lambda相关测量法是以众数作为预测的准则，没有考虑其它的条件次数，因此，当众数集中在条件次数表的某行或某列时，Lambda系数会等于0。但这并不一定真是

X与Y完全无关。

故Lambda相关测量法的敏感性有问题。4.显著性检验

6.3.2Goodman-KruskalTau相关测量法Goodman-KruskalTau相关测量法是由古德曼和克鲁斯卡尔提出的，采用Tau系数测定两个定类变量间的相关程度。Tau系数是对Lambda系数的改进。它不再用众数对Y进行预测，而是利用边缘次数提供的比例进行预测。1.非对称形式的Tau相关测量

根据PRE准则，可得Tau相关测量系数为（6.14）式。

【例6.6】随机从城镇、乡村两个地区抽取10800户家庭调査，结果如表6-12。讨论城乡地区与健康知识传播渠道的相关。传播渠道（Y）地区（X）城镇乡村卫生机构宣传319111204311网络宣传237011833553商品推介12865281814其他媒介93518711227782301810800表6-12分析：表中提供的数据是两个定类变量的值，测定城乡地区与健康知识传播渠道的相关程度，是为考察城乡不同地区信息传播渠道是否不同，对两个不同地区是否应采用不同的传播形式。自变量是不同地区，因变量是不同的信息渠道。这是非对称的关系，可以采用Tau相关测量法。

2.Tau相关测量法的特点

3.显著性检验

【例6.7】从某地区随机抽取620位居民进行调查，调查结果如表6-13所示，欲进行是否定期体检与小病处理方式的相关分析。是否定期体检（Y）小病处理方式（X）去医院

自己买药

忍耐

合计检验结果定期体检非定期体检71182513042624743316合计9742994620表6-13分析：先进行独立性检验，建立假设组：

由R语言可以算得p值很小，所以拒绝源假设，认为定期体检与小病处理方式有关。具体的关系如何？接下来利用Tau相关测量法进一步分析。

6.3.3Gamma相关测量法

1.同序对和异序对同序对是X变量中的数值与Y变量中的数值变化方向一致的序对。表6-14是三个职工受教育程度与经济收入的列联表，由表6-14可以看出，甲、乙、丙三人在受教育程度和经济收入方面的位次有如下的关系：受教育程度：乙>甲；乙>丙；甲>丙，经济收入：甲>乙；乙>丙；甲>丙。经济收入受教育程度高中低高

甲

中乙

低

丙表6-14

2.对称关系的Gamma相关测量Gamma相关测量法具有消减误差比例的意义，它利用同序对和异序对定义系数G。

由PRE准则，可以定义G系数为下式：

G系数没有考虑Y对X的影响，还是X对Y的影响，因而是对称关系的相关测量。

YX高中低高中低表6-15

异序对可以采用类似的方法得到，即将同行同列的同分对舍去，再舍去某一频数右下方的同序对，在某一频数左下方的都可以构成异序对。表6-15中的异序对总数为

【例6.8】从医院某病患者随机抽取55名通过评分测得头晕和头重的不同等级状况如表6-16。是否能够测量头晕和头重的相关性。头重头晕轻中重轻17140中6151重002表6-16

表明两个变量间存在正相关，相关程度不算太高。由于G系数具有消減误差比例的性质，因此G=0.6456意味着，以头晕的相对等级解释头重的相对等级可以消减64.56%的误差。由于G系数是对称关系的相关测量，因而也可以说，以头重的相对等级解释头晕的相对等级可以消减64.56%的误差。头晕与头重呈正相关，表明一个变量等级愈高，另一变量等级也愈高。4.Gamma相关测量法的特点

G系数测量对称关系的相关程度。从定义可以看出，G系数只考虑同序和异序的关系，因此，无论用X预测Y，还是Y预测X，计算的G值都是一样的。

YX表6-16

【例6.9】新药是否能预防肝炎？采用配对样本进行研究，抽取60对，每对中随机指定一人服用新药，另一人服用原有药，经过一段时间的观察，结果如表6-17.YX服用新药服用原有药未患肝炎6048患肝炎012表6-17

5.Gamma系数的检验利用随机样本数据计算的G系数，是否能用以推断总体，必须进行统计检验。建立的假设组为

为判定假设，需要采用随机抽样获得数据，数据至少是定序尺度测量。定义的检验统计量为

【例6.10】表6-16中从医院某病患者随机抽取55名通过评分测得头晕和头重的不同等级状况得到G=0.6456，能否将这一结论推断到所有该病的患者，也就是该系数是否能够说明总体具有的性质，需要进行检验。分析：如果研究的是患者头越晕就感觉头越重，建立单侧备择，建立的假设为

结果表明p=0.007。在0.05的显著性水平下，拒绝原假设，说明可以用55名患者头晕与头重的数据推断得出患者头越晕则头越重的结论。6.3.4Somer'sd相关测量法Somer'sd相关测量法亦称d相关测量法，是通过计算d系数测定变量相关程度的方法。适用于两个定序变量间相关的测量。1.非对称关系的d系数G系数测定具有对称关系的两个定序变量间的相关程度，而当两个变量X与Y为非对称关系，即X对Y的影响和Y对X的影响不一样时，采用G系数不够严谨，而宜采用d系数。d的计算公式为

于是d系数的计算公式为

2.列联表中同分对计算如表6-15将两个变量按等级顺序分别排列，在Y变量上的同分对数目是表中同一横行每两格次数乘积的和，即

在X变量上的同分对数目是表中同一纵列每两格次数乘积的和，即

头晕（Y）心悸（X）合计轻中重轻0202中811524重1281434合计20211960表6-18

3.d系数的显著性检验

对假设作出判定，需要计算检验统计量Z。Z的计算公式为

式中，S'是S的修正值

转化为正态性检验，可得p=0.2444，表明：患者的心悸与头晕成正比，但是几乎没有什么相关。第7章高维列联表与对数线性模型本章内容：7.1三维列联表及其命名法7.2三维表的独立性检验7.3对数线性模型7.1三维列联表及其命名法

年龄（i）吸烟状况（j）呼吸情况（k）合计正常尚可异常<40从不吸烟1615536吸烟73434440-59从不吸烟1315吸烟18312合计25601297表7-1呼吸情况与年龄、吸烟状况调查表单位：人

对于表7-1的数据，可以求得

7.2三维表的独立性检验在二维列联表中，由于只涉及两个变量，因而独立性检验的假定也仅涉及两个变量。高维列联表中涉及的变量多于两个，因而可能研究的是所有变量间是否独立，也可能希望检验某些变量与其他一些变量相互独立，或是某一特定变量与其余变量无关。7.2.1三个变量相互独立性的检验若三个变量分别记为X、Y、Z，则检验三个变量相互独立的假设组为

原假设的意思很明确，但是备择假设的情况就比二维情形复杂。

【例7.2】例7.1中的呼吸情况与年龄、吸烟状况是否有关，数据见表7-1。

7.2.2局部独立性检验三个变量间相互独立性的假设被拒绝，并不意味着所有变量之间都存在显著的联系。可能是两个变量间相关，而第三个变量完全独立，即局部独立；也可能是两个变量在第三个变量的每一水平上是独立的，但两个变量的每一个都与第三个变量相关，即当给定第三个变量的水平时，前两个变量是条件独立的。零假设记号假设组X与（Y，Z）独立Y与（X，Z）独立Z与（X，Y）独立三维表的局部独立性及其假设组为了检验假设，可以按照与之前相同的方式进行。这里以检验X与（Y，Z）独立为例说明检验统计量的推导与检验过程，其余情形完全类似。假设组为

检验的p值计算与前面完全类似。【例7.3】例7.1中的年龄与吸烟状况和呼吸情况是否无关，数据见表7-1。分析：由例7.2的结论可知，年龄、吸烟状况、呼吸情况三个变量不是相互独立的，本例检验年龄与吸烟状况和呼吸情况是否无关，建立假设组

计算局部独立性检验统计量值和相应p值的R语言函数也是loglin，故可与例7.2完全类似的方法算得统计量Q相应的p值为p=0.288。在0.05的显著性水平下，不拒绝原假设，表明年龄独立于其他两个变量。7.2.3条件独立性检验

假设组为

原假设记号假设组给定X，Y和Z条件独立给定Y，X和Z条件独立给定Z，X和Y条件独立

三维表独立性检验的期望频数与自由度【例7.4】例7.1中给定年龄时，吸烟状况和呼吸情况是否条件独立？

7.3对数线性模型前面关于列联表分析的假设检验，能够有助于认识变量间较为复杂的关系，但是没有将这些关系量化。对数线性模型恰好能解决这个问题。7.3.1对数线性模型的类型和参数估计

基于上述方程组结合（7.10）（总共4个方程）就可以得到

这类似于方差分析中所用的模型，称之为二维列联表两个变量独立时的对数线性模型。上述模型的统计学意义：

2.更复杂的对数线性模型的类型（一）饱和模型当变量间相互不独立时，对数线性模型称为饱和模型。此时，模型中，不仅有各个变量的主效应，还有变量间的交互效应。

多个变量的饱和模型。与高维列联表分析一样，多变量的对数线性模型通常也以三个变量为例。其饱和模型为（7.14）式。

从（7.13）、（7.14）式可看出，饱和模型包含了变量的全部可能的主效应与交互效应。这类模型能够完全拟合数据，即列联表提供的资料的数目恰恰够确定模型的参数，因而误差为零。（但往往没用）（二）非饱和模型非饱和模型即独立模型，模型中没有变量间的交互效应。两个变量的非饱和模型。两个变量X和Y相互独立时，列联表上每个单元的频数只受两个变量各自的影响，因而其模型为（7.15）式。

多个变量的非饱和模型。这里仍以三个变量为例，三个变量的非饱和模型是三个变量间不存在任何交互效应的对数线性模型，表达式为（7.16）式。

非饱和模型是只包含总平均效应以及每个变量主效应的对数线性模型。由于模型的参数少于列联表中单元的数目，模型对数据不能完全拟合，因此需要对模型进行检验。（三）谱系模型（HierarchicalModels）

与前面三维表独立性检验的情况相同，将三个变量的所有关系分为三种类型。类型I为三个变量相互不独立的情况，类型II是三个变量局部独立的情况，类型III是三个变量条件独立的情况，类型IV是三个变量相互独立的情况。三个变量间的谱系模型

3.模型的参数估计对于三变量的饱和模型

各效应间有约束条件：

【例7.5】对例7.1的数据建立对数线性模型进行分析。分析：由例7.3和例7.4的结论可知，年龄独立于吸烟状况和呼吸情况，但是吸烟情况和呼吸情况是不独立的。因此可建立模型如（7.17）式：

采用与饱和模型相同的约束，即假定各个效应参数满足下列关系

所有效应的估计结果如表：

4.模型的检验

二是似然比统计量如（7.22）式。

【例7.6】对例7.1的数据建立的所有对数线性模型进行拟合优度检验。

5.模型的选择

【例7.7】在例7.6的基础上，选择最优模型。

可以看出，模型（1）是模型（2）和模型（3）的子模型，可以比较（1）和（2）、（1）和（3）有没有显著差异。在R语言中，可以通过anova函数完成模型选择。

由表中的结果可知，没有充分理由认为模型（1）与模型（2）、模型（3）有显著差异，最终选择模型（1），因为该模型最简单。（2）自后淘汰法

【例7.8】以例7.1的数据建立对数线性模型，利用向后淘汰法选择模型。

第8章基于经验分布函数的方法本章内容：8.1非参数Bootstrap8.2Kolmogorov-Smirnov检验8.3Kaplan-Meier估计与对数秩检验8.1非参数Bootstrap在实际数据分析中，总体分布往往并不知道。但经常需要估计分布函数F的某些特征值，例如均值、方差、分位数等。此时，非参数Bootstrap是一种可供选择的方法。Bootstrap通常被翻译为自助法。所谓自助，就是指非参数Bootstrap是通过不断对已有样本重复抽样实现的。1.基本原理

普通样本Bootstrap样本样本服从的分布分布特征数统计量估计量普通样本与Bootstrap样本的关系例如，要估计的是总体分布r阶矩，在总体分布F的框架下有（8.1）式，

在经验分布函数的框架下，相应的r阶矩为（8.2）式，

2.估计量偏差和标准误差的Bootstrap估计

具体做法如下：

对于偏差的估计，重复抽样的次数B需要大一些，一般设为2000。对标准误差的估计，重复抽样的次数B不用太大，一般50次就够了。这些都是经验法则，仅供参考。【例8.1】

有82个大学法学院的LSAT成绩（法学院入学测试成绩的平均分）和GPA成绩，从中随机抽取容量为15的样本，LSAT576635558578666580555661GPA3.393.32.813.033.443.0733.43LSAT651605653575545572594

GPA3.363.133.122.742.762.882.96

利用这个15个数据估计两种成绩之间的相关系数，并计算估计的偏差和标准误差。分析：在这个问题中，可以用皮尔逊相关系数、肯德尔秩相关系数、斯皮尔曼秩相关计算两类成绩的相关系数，由于不知道总体相关系数是多少，标准误差无法直接计算。因此尝试用非参数Bootstrap方法估计标准误差。R代码见教材。15个样本数据的皮尔逊相关系数为0.776，全部数据的皮尔逊相关系数为0.760。

3.置信区间的非参数Bootstrap估计置信区间的Bootstrap估计内容十分丰富，这里介绍百分位数Bootstrap置信区间，此方法体现了非参数Bootstrap估计的独特思想。

【例8.2】有一批人进行打靶，分别记录了他们培训前和培训后的打靶成绩，数据如下：

打中次数总次数培训前9811034培训后11911037根据数据回答培训是否有效果。

4.假设检验的非参数Bootstrap方法

当原假设成立时的p值为

【例8.3】利用非参数Bootstrap重新做例3.3。幼儿园的生活对孩子的社会知识是否有影响？有人认为儿童上幼儿园有助于其认识社会，有人则认为儿童在家一样可以获得社会知识。为了解两种方式对孩子社会认知是否存在显著差异，对8个同性孪生儿童进行实验，随机指定8对中一个上幼儿园，另一个则在家。经过一个时期后，通过对他们的询问。给他们分别作出相应的评价，数据如下：配对编号12345678上幼儿园儿童7870678176728583在家儿童6258637780738278分析：建立假设组

采用Wilcoxon符号秩检验完成这个问题检验。将上幼儿园儿童的得分减去在家儿童的得分，得到差值，再进行有放回等概率抽样得到Bootstrap样本，计算Bootstrap样本下的统计量的值以及全部样本下的统计量值，得到p值。结果表明p=0.0624。这个结果与例3.3的结果略有不同。此处是不依赖正态近似得到的p值，但最终结论与例3.3是一致的。8.2Kolmogorov-Smirnov检验Kolmogorov-Smirnov检验简写为K-S检验，常译成柯尔莫哥洛夫-斯米尔诺夫检验，是一种通过经验分布函数与理论分布的接近程度进行的分布检验。检验的目的是考察所获取的样本数据是否来自具有某一理论分布的总体。1．单样本K-S检验

建立假设组

如果是单侧检验，检验统计量需要进行修正，对于右侧检验

检验统计量为

对于左侧检验

检验统计量为

当原假设成立时，这些统计量都有很复杂的精确分布和渐近分布，一般而言，当样本量少于100且没有打结时，用精确分布。在R语言中有函数ks.test可以对前面的各种情况输出精确检验和渐近检验的p值。【例8.4】公共交通设施适合性的研究——公共汽车到达时间是否服从正态分布。公共汽车按计划每15分钟通过一个商店旁。然而，由于交通条件，乘客数目等影响，汽车实际到达的时间可能会有很大不同。通过一天随机的观察，获得的数据如下表。比计划提前到达的为负值，取大的整数，如提前1分10秒到达，记作-1；比计划晚到的为正值，也取大的整数，如迟到1分10秒，记作+2。公共汽车到达时间是否服从标准差为3的正态分布？到达时间（x）-5-3-1012478观测频数（f）112155311

根据统计量的定义，从表中可知，D的值为0.1970。在某些书上有专门的统计表查找相应的临界值。在R语言中，利用ks.test函数可以很方便得到p值，结果表明，精确检验的p值为0.3707，渐近检验的p值为0.4197。都不能拒绝原假设。2.两样本K-S检验

单侧检验假设组为

单侧检验也需要进行相应的修正，感兴趣的读者可以自己推导。【例8.5】城郊县是否比边远县有较低的人口增长速度？根据城镇和农村的划分，将城市的近郊县与边远的远郊县分为两类地区。研究这两类地区的人口相对增长速度的分布是否存在显著差异，是否远郊县的相对增长速度高于城近郊县。随机地从两类地区分别抽取9个远郊县和7个近郊县调查，依据人口普查的资料得知，近10年人口的年均增长速度（单位：％）分别如下表。远郊县5.33.75.83.54.85.66.89.310.3近郊县3.12.18.22.73.47.93.2

结果表明精确p值为0.0180，拒绝原假设，认为存在远郊县比近郊县有较低的人口增长速度的情况。8.3Kaplan-Meier估计与对数秩检验

8.3.1生存分析的基本概念1.终点事件（outcomeevent）和生存时间（survivaltime）在生存分析中，终点事件或称结局事件。可以是产品下架、政策失败、肿瘤复发、病人死亡等研究者所定义的一个被关注的事件。起点时间，是对所关注对象开始观察的时间。为了使研究对象具有同质性，一般需要定义对象的某个事件发生或者某种干预处理的时间作为起始时间。在生存分析中，如果研究对象的终点事件发生，则终止该对象的观测并且定义终点事件发生的时刻为该对象的终止观测时间。在实际中，有些对象在研究结束时，终点事件还没有发生，比如某些破坏性试验中的对象始终正常工作；有些患者到观测结束仍旧存活；有些研究对象观察期间发生失访，例如病人转院。这些没有观测到终点事件发生，通常以最后一次联系到这个对象的时间为终止观测时间。从起点时间到终止观测时间的时间间隔定义为生存时间。也称为失效时间（failuretime）。对于终点事件发生的生存时间称为完全数据（completedata）。没有观测到终点事件发生的生存时间称为删失数据（censor），也称为截尾数据。删失的生存数据常用“+”表示。【例8.6】脑出血病人的生存曲线。Dupont（2002）讨论了脑出血病人的生存时间受不同基因的影响问题，对应的数据在作者的个人网站中可以找到。这个数据就是生存数据，其中有完全数据和删失数据，总共对71个病人进行了观测，结果如表所示其中，数据带“+”的是删失数据。2.生存概率、生存率、生存函数、生存率曲线（1）生存概率生存概率表示某段观测时间开始时存活的个体到该时段结束时仍存活的可能性大小，用p表示，计算公式为

其中，q表示死亡概率。如果该时段内有删失，则分母须用校正的个体数，校正的个体数为期初观测个体数与1/2删失数之差。（2）生存率和生存函数

（3）生存率曲线生存率曲线也简称生存曲线，是以时间为横轴，生存率为纵轴绘制的曲线。8.3.2生存率曲线的Kaplan-Meier估计当个体中存在删失，Kaplan-Meier估计的基本思想如下。

严格来讲，还需要注意到两个量。

【例8.7】某医院治疗某种疾病后，随访记录患者死前存活月数，结果为：2+，13，8+，11+，6，1，9，3，17，7。估计该疾病的生存率。第一步，将样本数据从小到大排序；第二步，列出各时段初期的存活个体数，以及各时段内死亡个体数；第三步，计算各时段的死亡概率、生存概率；第四步，利用第三步所得生存概率的估计值，连乘得到生存率的估计。具体结果如下表秩观测时间（月）期初病例数死亡人数死亡概率生存概率生存率111010.10000.90000.900022+900.00001.00000.900033810.12500.87500.787546710.14290.85710.675057610.16670.83330.562568+500.00001.00000.562579410.25000.75000.4219811+300.00001.00000.4219913210.50000.50000.21091017111.00000.00000.0000上表的最后一列是生存率，可以以生存率为纵轴，时间为横轴画出曲线，该曲线就是生存率曲线，如下图，其中中间的实线为生存曲线，两侧的虚线是置信带。可以看出，该曲线是阶梯状的，能大致反映生存率随时间的变化趋势。样本量越大，该曲线越光滑。8.3.3生存率曲线的假设检验在实际数据分析中，常常需要比较两条生存率曲线是否显著不同。例如在生物医学统计中，常要比较试验组与对照组生存率曲线的差异。对数秩检验是简单而常用的方法。其基本思想是比较实际失效（死亡）率与期望失效率之间有无显著差别。基本步骤如下：

第二步，将两组待检验的样本混合后统一从小到大排序，利用Kaplan-Meier估计的方法计算合并样本的死亡概率。

【例8.8】本例将估计例8.6数据的生存率曲线，并检验不同基因下的生存率是否显著不同。按照例8.6数据的排列顺序（按行排列），对应的基因分组如下表。0111111101111101NA111101111010000111011000110001110010000111000000100001

其中，0和1分别代表两种基因，分别记为组1和组2。基因数据缺失的患者不参与分析。分析：将例8.6和例8.8的数据结合，分析两种不同基因组（0或1）的生存率，可以画出两组的生存率曲线，用R语言程序可以实现。这个程序有两部分：一是绘制生存曲线，二是进行对数秩检验。这里只显示结果，代码见教材。实线为组1的生存曲线，虚线为组2的生存曲线。从两组的生率曲线可以看出，不同基因组病人的生存曲线显得很不一样，从直观上似乎可以断定他们有不同的生存率。再进行对数秩检验：

第9章

概率密度函数的估计9.1直方图与概率密度函数9.2核密度估计9.3概率密度的样条估计9.1直方图与概率密度函数

第三步，对数据进行分组，统计落入每组的样本点个数，即频数，就能得到频数或者频率分布表，将表画成图就是直方图。【例9.1】在一项研究中，为了研究学生的数学能力和他们解决问题偏好之间的关系，研究人员通过四个可视化测试收集了68个学生空间能力的数据。对于每个学生，他们的数据被合并成一个单一的分数，给出空间能力的总体评价。数据见教材表9-1.1.概率密度的直方图估计

当样本量n足够大时，根据格林文科定理有

即直方图其实是总体概率密度函数的某种近似。2.直方图估计的缺陷与改进

Rosenblatt估计仍然存在如下问题：

9.2核密度估计9.2.1核函数（KernelFunction）与核密度估计

定义9.1的几点注释：

【例9.2】例9.1数据的Rosenblatt估计。

Rosenblatt估计得到概率密度曲线不连续是它的缺点。在实际应用中，可以选择别的核函数，得到连续的概率密度曲线估计。常用的有以下几种：

在R程序包kedd中，提供上面的全部7种核函数。在实际应用中需要选择合适的核函数；同时，对于选定的核函数，需要考虑适当的窗宽。用高斯核对例9.1的数据在不同窗宽下进行概率密度函数的估计如下图：9.2.2最优窗宽的选择1.拇指准则（ThumbRules）法拇指准则法依赖渐近理论，首先不加证明地介绍几个基本的理论结果。

当真实的概率密度函数与正态分布相近时，这个方法效果很好。而且因为有显示表达式，这个方法的计算速度很快。但是如果真实的概率密度函数与正态分布相去甚远，效果就很差，甚至误导。这是拇指法的缺点。需有其他方法弥补。2.

留一交叉验证法（Leave-One-OutCross-Validation）（1）交叉验证法的基本思想

（2）最大似然交叉验证法

（3）最小二乘交叉验证法

【例9.3】对例9.1数据进行核密度估计，确定最优窗宽。核函数极大似然交叉验证无偏交叉验证gaussian0.10580.0431epanechnikov0.13980.0877uniform0.13650.0773triangular0.14580.0903triweight0.17380.1185tricube0.15610.1048biweight0.15710.1020cosine0.14320.0894本例通过R程序包kedd实现。首先用拇指法确定最优窗宽，kedd包中的h.amise函数能够计算拇指法的窗宽。但是需要选择kernel为gaussian，并且将deriv.order设为0，这是这个函数默认的选项。最大似然交叉验证法的R函数是h.mlcv，无偏交叉验证法的R函数是h.ucv。7种核函数的最优窗宽分别用两种交叉验证法进行计算，代码见教材。结果表明：拇指法最优窗宽的结果为0.1482。对于高斯核，本例数据的两种交叉验证法所得结果都与拇指法有较大不同，这说明拇指法可能不适合本例。此外，最大似然交叉验证法的窗宽比无偏交叉验证法的窗宽在数值上要偏大。但仅从数值上比较没有太大意义。因为两种交叉验证法基于不同的损失函数，有不同的意义。9.2.3自适应窗宽（最近邻法）

狭义上的最近邻法估计概率密度函数，应该选择的核函数是均匀核。但是均匀核有个致命缺点，就是估计的曲线不是平滑曲线，因此也可以将均匀核换成其他光滑的核函数。例如可以换成高斯核，这时的估计是最近邻法的一种推广，可以得到更加平滑的概率密度曲线估计。【例9.4】利用最近邻法估计例9.1数据的概率密度曲线。

9.3概率密度的样条估计9.3.1样条的基本概念

9.3.2对数样条密度估计对数样条密度估计（LogsplineDensityEstimation）是最大似然估计与样条方法的结合，由Kooperberg和Stone在1991年提出。

该方法是基于指数族推导得出的，对于非指数族分布，有众多改进研究。并且对样条函数的选择，也有很多文献进行了讨论。感兴趣的读者可查阅相关文献，例如BakKY,JhongJH,LeeJJ,etal.Penalizedlogsplinedensityestimationusingtotalvariationpenalty.ComputationalStatistics&DataAnalysis,2020:107060.，该研究利用B样条提出了一个改进的半参数方法。节点的选择对于样条方法很重要，节点越多，越能刻画曲线的细节，但也有过拟合的风险，如何恰当选取节点是一个重要问题。关于模型选择，可以采用AIC、BIC等方法实现。在R中概率密度曲线的对数样条估计的程序包是logspline，其中的oldlogspline函数和logspline函数都可以实现概率密度函数的三次样条估计，后者是前者的改进版本，这两个函数都能够自动确定节点并自动选择模型。

分析：分别采用核密度估计和对数三次样条法估计密度曲线。首先产生200个服从题设分布的随机数，再分别利用核密度估计和对数三次样条估计得到概率密度曲线的估计，画出图形，并将真实的概率密度曲线图也一并画出。其中核密度估计采用拇指法确定窗宽，核函数选择高斯核，其余用到的参数均为R中的默认参数，相应程序见教材。运行结果，如图9-6。从图9-6可以看出，两者都得到平滑的曲线估计，对数三次样条法的结果看起来比核密度估计更接近真实的密度曲线。一般来说，核密度估计对曲线尾部概率高估，导致估计的偏差较大；三次样条估计对曲线的尾部概率低估，导致估计的峰值偏大。本例总体是对数正态的混合，尾部概率很小，所以尾部低估带来的误差比高估带来的误差更小。图9-6两种方法估计混合对数正态密度曲线9.3.3惩罚似然估计

【例9.6】利用例9.5的数据，进行惩罚似然估计。分析：在R中有程序包gss可以实现概率密度函数的惩罚似然估计，具体的函数是ssden。首先利用ssden函数得到惩罚似然估计，并用dssden函数计算估计得到的密度函数值，这样就可以画出估计的曲线，相应程序见教材。运行结果如图9-7。估计时的参数全部采用R中的默认参数。从图9-7可以看出，惩罚似然估计能够较好地识别概率密度函数的细节。对比图9-6的结果可知，对数样条估计和惩罚似然估计看起来都比核密度估计要好，但是对数样条估计与惩罚似然估计孰优孰劣，这里不再展开，有兴趣的读者不妨加以思考。图9-7混合对数正态密度曲线的惩罚似然估计第10章非参数回归10.1基本思想10.2一元核回归10.3一元样条回归10.4可加模型（AdditiveModels）10.1基本思想10.1.1经典回归模型及其不足

图10-1模拟摩托车事故数据建模10.1.2局部拟合思想与K近邻回归从例10.1可以看出，对于这样的数据，线性回归、二次回归、三次回归都不能进行恰当的拟合。原因在于，该数据中两个变量的关系较为复杂，并非为常见的关系，用确定形式的模型强行拟合，势必导致拟合效果不好。

K近邻回归的基本思想

运行结果如图10-2。由K近邻回归的估计过程可知，这种方法与线性回归等参数回归的最大区别是不假定具体的模型形式，而是根据数据自身的特点局部进行估计，这是典型的非参数回归方法。从图10-2可以看出，K近邻回归可以反映出回归曲线变化的细节，并且能够刻画曲线变化的趋势，这说明K近邻回归比图10-1中的诸方法效果更好。图10-2例10.1数据的K近邻估计K值的选择K近邻回归曲线的形状与K的取值有关。从图10-2可以看出，K值越大，曲线越平滑，如k=50时；反之越陡峭，如k=10。对于K的最优选择，可以使用交叉验证实现。先给出K的候选集合，在候选集合中选定一个K，先将数据随机分成m份，每次留出一份作为测试集，其余m-1份作为训练集，这样不停轮换建立K近邻回归模型，得到m个测试误差，取这m个测试误差的平均，就得到给定K时对应的交叉验证误差（CV误差）。依同样步骤，每个K值都可以得到一个CV误差，这些CV误差中最小那个对应的K就是要选择的K。本例采用5折交叉验证，K的候选为数字1到50，得到最优的K是8，CV误差为117.72。R代码见教材。图10-3利用5折交叉验证选择确定的最优K近邻估计10.2一元核回归

10.2.1局部常数核回归1.基本思想

2.窗宽选择

窗宽选择的目的就是确定一个估计，使得方差和偏差都不太大。和第9章类似，可以通过拇指法和交叉验证法实现窗宽选择。（1）拇指法

（2）留一交叉验证法

【例10.4】例10.1数据的局部常数核回归。分析：通过本例实现例10.1数据的局部常数核回归估计，并针对不同的窗宽，作出相应的回归曲线。在R中有很多程序包可以实现，本例采用stats包中的ksmooth函数实现。窗宽选择采用npbr包中的kern_smooth_bw函数实现，该函数提供交叉验证和BIC两种方法进行窗宽选择，这里仅考虑交叉验证的情形，更多的细节请读者自己查阅相关资料。相应程序见教材。运行结果如图10-4。图10-4用到的核函数为高斯核，读者可以选择其他核函数，但是结果相差不大。图中窗宽分别为0.5、2.5、5，利用交叉验证得到的最优窗宽值14.7。从图中可以看出，窗宽越小，回归曲线越粗糙，窗宽越大，回归曲线越平滑。在实际应用中，有研究者指出：“我们的实践经验是，利用自动方法选择窗宽，如交叉验证，通常效果很好，但有时会产生明显与基于背景知识选出的值不一致估计结果”，也就是说，需要结合数据背景选择窗宽，不可迷信交叉验证。见FarawayJJ.ExtendingtheLinearModelwithR:GeneralizedLinear,MixedEffectsandNonparametricRegressionModels,SecondEdition[M].2016。原文为：Ourpracticalexperiencehasbeenthatautomaticmethods,suchasCV,oftenworkwell,butsometimesproduceestimatesthatareclearlyatoddswiththeamountofsmoothingthatcontextualknowledgewouldsuggest.Forthisreason,weareunwillingtotrustautomaticmethodscompletely.图10-4不同窗宽的局部常数核回归估计10.2.2局部多项式和局部线性核回归

1.局部多项式核回归

2.局部线性核回归

【例10.5】利用局部线性和局部多项式估计例10.1的数据。分析：本例实现例10.1数据的局部多项式回归和局部线性回归。在R中有多个程序包可以实现本节方法，本例采用locpol包。这个包中的regCVBwSelC函数可以实现基于交叉验证的窗宽选择，locPolSmootherC函数可以实现局部多项式和局部线性估计。当采用局部多项式核回归时，需要选择多项式的次数，有研究建议选择奇数次。从理论上讲，多项式的次数越高，估计的精度就越高，但是计算会越复杂。

图10-5局部线性与局部多项式回归曲线10.3一元样条方法10.3.1多项式样条多项式样条不采用逐点逼近，而是用多项式逐段逼近回归曲线。1.基本样条（BaseSpline）基本样条（BaseSpline）的最简单的样条，其基本思想就是逐段地用线性函数逼近回归曲线。

2.M样条（M-spline）

3.自然样条（naturalspline）M样条在边缘区域的估计效果不好，这是它的一个缺陷。为了克服这个缺陷，一种方法是用线性函数去估计边缘区域的回归曲线，这就是自然样条方法。

4.节点的选择在实际中，一般选择三次样条函数足以胜任大部分情形。节点的选择是回归样条的重点和难点，没有统一的选择方法。这里介绍两种经验法则。一是用自变量的分位点作为节点，在R语言的splines程序包中就是用这种方法作为默认的节点确定方法。具体要采用多少个分位点作为节点，可以采用k折交叉验证实现。但是这种方法采用等间隔的分位点不完美，因为曲线弯曲程度大和样本点密集的地方一般需要更多的节点，如果强行使用这个方法可能带来较大误差。二是根据样本点的疏密程度和曲线的弯曲程度选择不同的节点。这种方法设计交叉验证算法的难度较大，主观的程度较高，这是一缺点。有学者发展了一些数据驱动的节点选择方法，感兴趣的读者可以参考论文VivienGoepp等（2018）的文章SplineRegressionwithAutomaticKnotSelection。【例10.6】利用多项式样条和自然样条方法估计例10.1的数据。分析：在R中还有很多程序包可以实现样条方法，感兴趣的读者可以参考综述论文ArisPerperoglou等（2019）的文章AreviewofsplinefunctionproceduresinR。本例采用在R语言中的splines程序包实现回归样条，其中M样条可以用该程序包中的bs函数实现，自然样条可以用ns函数实现。首先通过5折交叉验证确定节点。在splines包中，有两种方法可供选择，一种是通过选择不同的自由度（df）确定