数学建模中的统计分析问题(样本比较-置信度评估)

编号专用页

赛区评阅编号（由赛区组委会评阅前进行编号）：

全国统一编号（由赛区组委会送交全国前编号）：全国评阅编号（由全国组委会评阅前进行编号）：

白血病临床治疗的统计分析问题

摘要

一、问题重述

为研究某药物6-MP是否有治疗以缓解病痛的作用，研究者在持续1年的急性白血病治疗的临床试验中，将42位急性白血病患者（进入项目的时间有先后）随机地分成两组（各21人）。对一组病人用药物6-MP治疗以缓解病痛，而另一组病人用安慰剂。安慰剂的外形和颜色与药物完全相同，但不含任何药物，病人自己并不知道实际服用的是药物还是安慰剂。研究者记录下每个病人病痛缓解的持续时间（以周为单位），持续时间越长则疗效越好。数据见表1。

表1数据后面有+者表示，当项目结束时缓解仍在持续。例如，处理组中的20+表示：该病人在项目结束前20小时进入临床治疗，使用6-MP后，缓解持续到项目结束。因此，该病人的实际缓解持续时间至少为20周，很可能大于20周。这种数据在统计学中称为删失数据。

我们需要回答的问题是：

问题1. 6-MP能否显著延长缓解的持续时间？

问题2.如果问题一不能得到肯定的回答，则对该药物没有必要进一步研究；反之，如果结论是肯定的，预测以后的病人在使用6-MP后的缓解持续时间的有关参数，对6-MP的效果给出有足够置信度的量化评估。

二、模型假设

1．假设在项目期间的食物，生存环境，其他药物等外界因素对6-MP药效及病人病痛无影响；

2．假设在项目期间各个阶段病人的6-MP药物服用量充足，治疗方式恰当；3．假设在项目期间无其他病痛误判，粗心等原因引起数据记录失误；

4．假设对同一个病人使用药物的效果始终一样；

5．假设每个病人的身体、精神素质都是相当的，不会因此而使药物的效果变化；6．假设病人自己并不知道实际服用的是药物还是安慰剂；

7．假设病人的年龄、性别对试验无影响；

8．假设两组病人是随机分配的。

三、符号说明及其概念解释

3.1符号说明

3.2概念解释

生存时间：疾病治疗的预后情况，一方面看结局好坏，另一方面还要看出现这种结局所经历的时间长短。所经历的时间称为生存时间。

完全与不完全数据：一部分研究对象可观察到死亡，从而得到准确的生存时间，所提供的信息是完全的，称为完全数据；另一部分病人由于失访、意外事故、或到观察结束时仍存活等原因，无法知道确切的生存时间，它提供了不完全的信息，称为不完全数据（截尾数据、删失数据）。

生存分析：生存时间一般是通过随访收集。不完全数据提供了部分信息。须要用专门的方法进行统计处理，这类统计方法起源于对寿命资料的统计分析，故称为生存分析。

死亡概率：指已活满t时刻的个体，在此后一段时期内（t至t t

+?）死亡的可能性。

生存概率：表示在某单位时段开始时存活的个体到该时段结束时仍存活的可能性大小。

四、 问题分析

本文研究者在持续1年的急性白血病治疗的临床试验中，对42位急性白血病患者随机均等地分成两组，一组病人用药物6-MP 治疗以缓解病痛，而另一组病人用安慰剂。通过对两组病人病痛缓解的持续时间进行对照比较分析，从而研究某药物6-MP 是否有治疗以缓解病痛的作用。

此问题可以转化为生存分析问题，即每个病人的缓解时间可以看每个成个体的寿命，从而可以采用生存分析的相关知识对问题进行分析求解。

由已知可知，每个病人进入项目试验的先后顺序不同，缓解病痛的时效也不同（如图一所示）。

图一：病人的缓解持续时间

通过对表1数据的分析，我们发现有些数据后面有+者表示，当项目结束时缓解仍在持续，这种数据在统计学中称为删失数据，又称截尾数据和不完全数据。对于这类数据的处理，如果我们丢弃删失数据只考虑确切数据，则会损失大量的信息；若将删失数据当作确切数据处理，则会低估了生存时间的平均水平。 用统计学的术语，白血病缓解效果的分析是一个“两样本比较”问题，一般用两正态样本均值比较的t 检验。但现在由于样本分布未知，而且在时间数据的分析中，由于数据分布有很大的偏度，正态分布是一个“坏”的模型。又因为数据是不完全的（有删失数据），常规的、用于完全数据的分析方法不能简单套用。 所以我们引入生存分析这一概念对本文进行分析求解。

生存时间经常服从的基线分布有指数分布、Weibull 分布、对数正态分布、对数logistic 分布和Gamma 分布。由于缓解持续时间不长，因此年龄、体质等可能影响缓解持续时间的因素作用不大，可以认为在任何时刻缓解持续的结束是随机的。又指数分布具有恒定危险率的特点，所以可假设生存时间服从指数分布,由次进行检验。

对于问题二，要预测持续时间参数，则先要给出其相关的参数，在对其置信

区间进行预测。

五、模型建立与求解

5.1问题一的模型

参数回归模型

本案例中样本容量不大，我们事先根据其生存函数曲线将分布假定为指数分布，所以可采用参数回归模型的分析方法，首先对指数分布进行检验，然后使用点估计的方法分别对两组数据的参数进行估计，接着对其有花都进行检验。用参数假设检验来判断处理组与控制组的缓解时间分布是否有显著差别从而判断药物6-MP能否显著延长缓解的持续时间。

首先分别对处理组和控制组的数据进行分析处理，由此拟合一个满意的参数分布，再用参数假设检验来判断处理组和控制组的缓解时间分布是否有显著差别，从而回答问题一：6-MP能否显著延长缓解的持续时间。

5.1.1模型的准备

我们假设每个病人病痛缓解的持续时间为生存时间。在这批数据中，其中控制组（使用安慰剂）的数据是完全的，没有删失数据。因为完全数据的分析比较简单，所以我们先对控制组数据进行处理。在刻画时间分布模型的特征方面，“生存函数”和“危险率函数”是两个重要的函数。

对控制组数据进行处理

通过计算机吃力可得，控制组生存函数图像为：

生存函数又称可靠性函数。是个体寿命超过某个时刻的概率。X记为个体的生存期，生存函数定义为：

()() S t P T t

=>=t时刻仍存活的例数观察病人总数

累积生存函数为：12()()i i k S t P T t p p p =>=…

危险率函数也称为风险函数、瞬时死亡率、年龄别死亡率、条件死亡率，常用h （t ）表示，它表示已存活到t 时刻的一个体，死于（t ，t t +?）小区间内的概率的极限。

0(|)

()lim

t P t T t t t T h t t

?→<≤+?<=? 0()()lim

(ln ())()t S t S t t d

S t t S t dx

?→-+?==-?? 累积危险函数为：0

()()ln ()t

H t h u du S t t λ==-=?

λ为指数分布的危险率，或称为尺度参数，其大小决定了生存时间的长短，

危险率越大，生存率下降越快；危险率越小，生存时间越长。在指数分布模型中，

λ是常数，与时间t 无关。

因为本文所给的生存资料分布具有不规则、不确定或未知分布的特点，所以采用非参数法估计生存率。 根据本文的样本含量为小样本，所以选择乘积极限估计法（Kaplan-meier ）来出来数据。

乘积极限估计法简称积限法或PL 法,是直接用概率乘法原理估计生存率，它是由统计学家Kaplan 和Meier 于1958年首先提出的,因此又称为Kaplan-Meier 法。此法计算生存率时，先将每个个体的生存时间按照由小到大的顺序排列，排序时若截尾值与非截尾值的观察时间相同，则规定非截尾值小于截尾值，排在截尾值之前。然后依次计算出各时段的死亡概率、生存概率，进而计算出从观察开始至各时刻的生存率(如图一、二所示)。

1.对处理组和控制组的累计生存率取自然对数，即取其累计危险函数：y ln ()S t =；

2.以个体寿命t 为横坐标，y ln ()S t =为纵坐标画折线图。

由上图可明显看出其图形走势近似直线。从而可以证明我们构建指数模型的思路是可行的。 因此设

()-ln ()H t S t t λ==

即：

()t S t e λ-=

由此可以直观的认为指数分布对数据的拟合是比较满意。

5.1.2 指数回归模型的构建： 指数分布准备知识：

若随机变量T 具有概率密度函数为：()t f t e λλ-=，则称 T 服从参数为λ（为尺度参数）的指数分布，简记为~()X E λ。 T 的总体分布函数为：0()1t

u t F t e du e λλλ--==-?

由此可导出以下公式：

()1()t S t F t e λ-=-= ()()/()h t f t S t λ==

人们用不变的危险率λ来刻划指数分布的特征，λ为指数分布模型中唯一的参数，其极大似然估计为：

1

=n i

i m m T t λ∧==∑

其中，n 为样本含量；i t 为每个观察对象的生存时间，i =1,..,n 包括完全数

据和截尾数据；m 为数据中完全数据的个数。 指数分布只有中只有一个参数，令控制组参数为1λ

因为控制组的数据是完全的，我们使用如下方法对参数1λ进行估计。 1、没有删失数据的极大似然估计：

1

1

2

n

()n

i

i t t t t n

L e

e

e

e

λ

λλλλλλλλ=----∑==…

两边取对数得：1

ln ()ln n

i i L n t λλλ==-∑

再对其进行求导得：

1

ln ()n =0n

i i d L t d λλλ=-=∑ 最后得：1

=n i

i n

t λ∧=∑

而控制组是完全数据不含删失数据，所以代入该极大似然估计模型得：

11

21==0.115182n i

i n t λ∧==∑

因此我们可以初步得出用0.115()t S t e -=来拟合控制组数据，为进一步确定这一假设，我们将对其拟合度进行优化检验。 用拟合优度对拟合效果进行评估

构造一个度量所假设的分布拟合数据优度的经验统计量D (0)D ≥，D 的值越小则表示拟合的越好。为了更好的确定其拟合的效果，我们令0.5ξ=，近似的认为当D ξ≤时拟合效果较好。

由于该分布为但一分布，则采用kolmogorov 检验法进行优度检验。

首先作出控制组总体分布的经验分布函数()n F t ，总体分布函数为0.115(t)F t e -= 假设控制组21人的病痛缓解的持续时间是总体的一样本，则可以得到控制组经验分布函数13()F t 的观察值为：

130

,12/21,124/21,235/21,347/21,459/21,58()13/21,81115/21,111217/21,121518/21,151719/21,172220/21,22231

,23x x x x x x F t x x x x x x x ≤??<≤?

?<≤?

<≤??<≤?

<≤??

=<≤??<≤?

<≤??

<≤??<≤?

<≤??>?

因为对于任一实数t ，当n →∞时，()n F t 以概率1一致收敛于分布函数()F t ，即：

{}

lim sup ()()01n n x P F t F t →∞-∞<<∞

-==

换句话说，对于任一实数t ，当n 充分大时，经验分布函数的任一个观察值()n F t 与总体分布函数()F t 只有微小的差别，从而在实际上可当作()F t 来使用。 检验统计量sup ()()i n i i t t

D F t F t ∈=-

为了便于求解，我们可简化 'max ()()i n i i t t

D D F t F t ∈≈=-

通过matlab 求解得：'0.1729D =，则'0.5D D ε≈<=故此得到0.115()t S t e -=对控制数据的拟合较好。

对处理组数据的处理：

通过计算机处理得到处理组的生存函数曲线图

对于含有删失数据的分布，我们采用含有删失数据的乘积限估计。 因为处理组中含有删失数据，所以我们对n 个数据合为k 个不同的死亡时期，从小到大排列为：1230k t t t t <<<<…<。在时间i t 上重复的个数为i g ，满足大于或等于i t （死亡、刪失、重复）的个数记为i n 。 假设

①在区间[]1,i i t t +上的刪失数据，其时间不早于1i t +。 ②当1i i t t t +<<时，条件概率(T|)i P t t T ≤≤的经验估计为：

i i i i n g p n ∧-= , i=0,1,…，k

由此得到其经验生存函数()n S t ：

当0i t t t ≤≤，00

()n g S t n -==1

当1i i t t t +<<，由21(T)()(|)i P t P t T P t T t T ≤=≤≤≤… 得：

累积生存函数121

()i

j j

n i j j

n g S t p p p n ∧

∧

∧

=-==∏

…

使用excel 软件对附录一、二的数据做出如下处理：

1、对处理组和控制组的累计生存率取自然对数，即取其累计危险函数：

y ln ()S t =；

2、以个体寿命t 为横坐标，y ln ()S t =为纵坐标画折线图。

由上图可明显看出其图形走势近似直线。从而可以证明我们构建指数模型的思路是可行的。

因此设

()-ln ()H t S t t λ==

即：

()t S t e λ-=

由此可以直观的认为指数分布对数据的拟合是比较满意。

有删失数据的极大似然估计： 设有n 个急性白血病患者（观察对象）进入急性白血病治疗的临床试验中, 其中有删失数据 m 个，即病人治疗到一半项目结束（即12m 0t t t ≤≤≤≤…）, 而另外n-m 个病人在m t 时能够接受到治疗。利用这一样本，我们用最大似然估计法来估计。

我们可以知道一个观察对象在[]t ,t t i i i d +失效的概率近似为

i i ()t t t i i f t d e d λλ-=（i =1,..,m ），而其余得n-m 个观察对象活过m t 得概率为

()

t =m

m

t

t t e d e λλλ∞--?

即：

()

()

n-m

t

=m

m

n m

t t

t e

d e λλλ∞

---?

故上述观察结果出现的概率近似地为：

()()()()

1

2m n-m

12m t t t m

t t t t n e d e d e d e m λλλλλλλ----?? ???

…

其中12m t t t d d d ，，…，为常数。因忽略一个常数因子不影响q 的最大似然估计,故可取似然函数为：

12[()]()m m t t t n m t m L e λλλ-++++-=…

对数似然函数为:

12ln ()m ln [()]m m L t t t n m t λλλ=-++++-… 令

12ln ()m

=[()]0m m d L t t t n m t d λλλ

-++++-=…

于是得到λ的最大似然估计为：

m =()

m s t λ∧

其中，12()()m m m s t t t t n m t =++++-…称为总观察时间，它表示直至时刻m t 为止n 个病人的试验时间的总和；m 表示出现观察终点的人数。

因为处理组中含有删失数据，所以将数据代入以上极大似然估计模型得：

2m 9===0.025()359m s t λ∧

因此我们可以初步得出用0.025()t S t e -=来拟合控制组数据，为进一步确定这一假设我们同样对其拟合度进行优化检验。（同控制组） 经检验得到0.025()t S t e -=能较好的拟合处理组数据。

由此得，控制组拟合函数为0.115()t S t e -=，处理组拟合函数为0.025()t S t e -= 12λλ>，则控制组的危险率要大于处理组，故同一时期控制组的生存率低。如下图所示，我们可以对两组病人的生存状况进行直观的比较。图中显示了两组病人的生存率随时间的延长成指数下降的趋势，处理组和控制组的下降率分别为0.025%/周和0.115%/周。

控制组和处理组的危险率比较

5.1.3 似然比检验模型

为了回答药物6-MP 能否显著延长缓解的持续时间，我们需构建似然比检验模型进行定性的分析：

原假设 0H ： 12=λλ (两组患者的生存率是相同的) 备择假设 1H ： 12λλ≠ (两组患者的生存率不相同) 因1λ，2λ的极大似然估计分别为：

111

m =T λ , 222

m

=T λ

又根据两样本的合并样本得到合并的尺度参数的估计为：

12

12

c m m T T λ+=+

构造似然函数：

[][]1212121211221212()=exp exp ()m m m m L T T m m λλλλλλλλ--=-+，

[][]12121212()=exp ()exp ()m m m m c c c c L T T m m λλλλ++-+=-+

则似然比为：

12

121212

()==()m m c c m m L L λλλλλλ+∧，

相应的对数似然函数为：

1212112212log()=log[()]log log ()L m m m m λλλλλλ=+-+，，

1212log()=log[()]()log()()c c c L m m m m λλλ=+-+ 如果0H 成立，则似然比统计量为：

[]2121,=2log 2log()log()~c G αλλλχ-∧=--，

即： []21211221,2()log()log log ~c G m m m m αλλλχ=-+-- 取显著性水平0.05α=，按自由度为 1的2χ界值作出决策。 同时：因为21,0.05 3.841χ=，

而统计量2χ的实测值216.75 3.841G χ==>。

由于16.75显著大于3.84，其否定区域，即备择假设成立。所以药物6-MP 能否显著延长缓解的持续时间。

综上所述，6-MP 具有显著延长缓解的持续时间，且其平均延长时间为控制组的4.6倍，即12/ 4.6λλ=。

由于控制组平均缓解时间()8.67i E t =控。则(t )8.67 4.640i E =?≈处周

5.2 问题二的模型

由问题一得出的结果证明了6-MP 能够显著延长缓解的持续时间，所以我们可以预测以后的病人在使用6-MP 后的缓解持续时间的有关参数，对6-MP 的效果给出有足够置信度的量化评估。

5.2.1 模型的建立

由研究者对42位急性白血病患者做持续1年的急性白血病治疗的临床试验可知，研究者不仅希望得到一种有效的药能够对急性白血病患者具有显著延长缓解的持续时间的作用，而且也想通过试验得到该药的有效适用范围。 由此，我们确定以下评估效果的参数：

1．以后的病人使用6-MP 后的期望缓解的持续时间；

2．以后的病人使用6-MP 后，其缓解的持续时间超过半年（26周）的概率； 3．以后的病人使用6-MP 后，具有70%的可能性其期望缓解的持续时间不低于某个下限。

由问题一的求解可知，处理组的指数分布能够较好的拟合其缓解持续时间分布，其指数拟合分布函数0.025()t S t e -=。

⑴ 由指数分布的性质知，以后的病人使用6-MP 后的期望缓解的持续时间最大似然估计值为：

11()400.025

u E x λ∧∧====周；

⑵ 以后的病人使用6-MP 后，其缓解的持续时间超过半年（26周）的概率为：

26260.025(26)0.522P x e e λ--?≥===；

⑶ 采用分为点估测法

设估计分位数为0.1x ，则满足0.30.3()1()0.7P x x P x x ≥=-<=

由此可知

0.30.3()0.7x P x x e λ-≥==

则0.3ln 0.7

x λ

-=

，由此得

0.3ln 0.7

14.270.025

x ∧

-=

=

得出理想值之后，我们需要对其效果给出足够置信度的量化估计，即：

1. 病人期望缓解持续时间的置信区间；

2. 持续时间超过半年（26周）的概率置信下限；

3. 病人使用6-MP 后，具有70%的可能性其期望缓解的持续时间不低于某个

下限的置信下限。

设处理组函数分布的似然函数为()L θ，θ为缓解期望时间估计值()E x 的待估

计量，信息函数2

2()ln ()I L θθθ

?=-?

由问题一可知，我们取的显著性水平0.05α=，则其置信水平10.95p α=-=。根据最大似然估计的渐进正态性有：

~(0,1)N ∧

，取u θ=。

另外，由于处理组服从指数分布，则似然函数123()(,,,...,;)n L L x x x x θθ=。 得：

123[...()]()m m x x x x n m x m L e λθθ-+++++-=，m 为处理组数据中完全数据的个数

同时取似然的自然对数得：

123ln ()ln [...()]m m L m x x x x n m x θθθ=-+++++-

对θ取二次偏微分得：123ln ()[...()]m m L m

x x x x n m x θθθ

?=-+++++-?

222

ln ()L m

θθθ?=-?

由此知，229

()()0.005640m I I u u θ∧====

~(1,0)N ∧

不依赖任何未知参数，根据正态分布α的分位定义得：

2

1P Z αα??<=-???

即：22

1P u Z u u Z ααα∧∧??

<<=-????

所以我们得到期望时间的置信区间为：

（2

2

,u Z u Z αα∧

∧

-）

又因为0.05α=，10.95p α=-=查表得0.0252

1.96Z Z α==

由代入上式可得期望延续时间u 为95%的置信区间

[u u ∧

∧

-+

0.75=

故：（1）u 的置信区间[39.853，40.147]。

由此可知：期望延续时间u 的置信下限为μ-

=39.853，此时的参数

1/1/39.853λμ∧

-

==

所以（2）(6)P X ≥的95%的置信下限为

26/39.853

(6)0.521X e

P --

≥== （3）0.3x 的95%的置信下限为

(ln 0.7)/(ln 0.7)39.85314.215λ-=-?=

5.2.2 模型的求解

5.2.3模型的检验

六、模型的评价与推广

本文所建立的模型具有以下特点：

合理性：每个实际问题都有一定得“背景机理”，而我们选择的统计模型考虑了该问题的背景机理；

合适性：我们的模型能够较好的拟合数据；

简单性：我们建立的模型比较简单，同时又可以较好的拟合数据，能够避免由随机因素造成的模型变形，而且容易从背景机理上解释。

参考文献

附录

附录三：相关的程序及代码

clear,clc

t=[6 7 9 10 11 13 16 17 19 20 22 23 25 32 34 35 ];

p=[0.857142857 0.806722689 0.806722689 0.752941176

0.752941176 0.690196078 0.627450909 0.627450909

0.627450909 0.627450909 0.537815065 0.448179221

0.448179221 0.448179221 0.448179221 0.448179221

x=6:0.01:35;

yi_nearest=interp1(t,p,x,'nearest');

plot(t,p,'ko');hold on;

plot(x,yi_nearest,'g','linewidth',1.5);hold on; grid on;

title('nearest method')