第六章 计数资料的统计描述

日期: 2020-11-24 作者：wxhyihuan

随机变量有连续型和离散型之分，相应的概率分布就可分为连续型分布(如 u分布、t分布和F分布)和离散型分布，常用的离散型分布有 二项分布、泊松分布，负二项分布，多项分布。

6.1 二项分布

二项分布(Binomial distribution)是\(n\)个独立的 成功/失败(或阴/阳，或是/非) 试验中“成功”的次数\(X\)(即取值[0~n]之间)的离散概率分布，其中每次试验的“成功”概率为\(p\)，记为B(n,p)。这样的单次成功/失败试验又称为伯努利试验，即当\(n = 1\)时，二项分布就是伯努利分布，0-1分布。二项分布是显著性差异的二项试验的基础。一般地，如果随机变量\(X\)(即取值[0~n]之间)服从参数为\(n\)和\(p\)的二项分布时，那么\(n\)次试验中正好得到\(x\)次成功的概率质量函数即为： \[\begin{aligned} f(x,n,p)=P(x,n,p) & = \binom{x}{n}p^x(1-p)^{n-x} \\ & = \frac{n!}{x!(n-x)!}p^x(1-p)^{n-x}, x\in[0,n] \end{aligned}\] P(x)就是二项式\([p+(1-p)]^n\)展开式的通项，\(n!/x!(n-x)!\)为二项系数，且\(\sum_{x=0}^nP(x)=1\)。

注意

对二项分布的资料进行抽样时，如果是有返回(放回)随机抽去n个个体，则依旧是刚才的\(P(x)\)分布函数；如果时无返回(不放回)的随机抽样，自由当n远小于总体个数N(如n < N/10)时，也可做近似\(P(x)\)分布函数处理。

适用的条件

  1. 每次实验结果只有发生和不发生两种情形(如成功/失败)，且每次实验中，成功/失败事件发生的概率都相同
  2. 重复实验是相互独立

6.1.1 二项分布的性质

6.1.1.1 均数与标准差

1.“成功”次数的性质

X的总体均数标准差

在n次独立重复实验中，出现“成功”次数X的总体均数，方差和标准差分别是， X的总体均数： \[μ=np\] X的总方差： \[σ^2=np(1-p)\] X的总标准差： \[σ=\sqrt{np(1-p)}\]

率p的总体均数标准差 在n次独立重复实验中，出现“成功”次数的率p(p=0/n,1/n,2/n…,n/n])的总体均数，方差和标准差分别是， “成功”次数的率p的总体均数： \[μ_p=p\]

“成功”次数的率p的总方差： \[σ^2=\frac{p(1-p)}{n}=Var_p\]

“成功”次数的率p的总标准差： \[σ=\sqrt{\frac{p(1-p)}{n}}=S_p\]

一般的，总体的p往往不知道，因此用样本资料计算样本的“成功”率\(p=X/n\)，来计算\(Var_p和S_p\)。

6.1.1.2 二项分布的图像

二项分布的图像偏态主要由p和n决定：

  1. 当p=0.5时，分布呈对称分布
  2. 当p≠0.5时，图形呈偏态，随着n增大，分布呈逐渐趋向于对称
  3. 当p不太靠近0或1，n->∞时，分布呈近似于正态分布

6.1.1.3 与其他分布的关系

  1. 伯努利分布是二项分布在n=1时的特殊情况。X~ B(1,p)与X~ Bern(p)的意思是相同的。相反，任何二项分布B(n,p)都是n次独立伯努利试验的和，每次试验成功的概率为p。
  2. 泊松近似 当试验的次数n->∞，而乘积np固定时，二项分布收敛于泊松分布。因此参数为λ=np的泊松分布可以作为二项分布B(n,p)的近似，近似成立的前提要求n足够大，而p足够小，np不是很小。
  3. 正态近似 如果n足够大，那么分布的偏度就比较小。在这种情况下，如果使用适当的连续性校正，那么B(n,p)的一个很好的近似是正态分布:N(np,np(1-p))。当n越大（至少20）且p不接近0或1时近似效果更好。不同的经验法则可以用来决定n是否足够大,以及p是否距离0或1足够远,其中一个常用的规则是np和n(1 −p)都必须大于5。

6.1.2 二项分布的应用

利用二项分布及其正态近似性，可进行总体率的区间估计和差异推断，也可用于研究非遗传性疾病的家族集聚性和做群检验。对于二分类事物的构成问题，对其中一类的构成比进行统计推断，其方法则类似总体率。

6.1.2.1 总体率的估计

6.1.2.1.1 精确法/直接法

对于n≤50的小样本资料，用查表法可以获取百分率的可信区间，即了得到其总体的1-α(α一般取0.05，或0.01)的可信区间。

R语言中的二项分布相关的函数与之前介绍的连续性分布(u分布，t分布和f分布)类似，主要有 dbinom(),pbinom(),qbinom()和rbinom()。 用于进行参数估计和假设检验的函数主要有binom.test()。

使用《医学统计学》中的 案例6-2 的数据在R中进行二项分布的精确法/查表法分析测试。

n <- 13
x <- 6
ratio <- x/n
binom.test(x, n, p = 0.5)
binom.test(x, n, p = 0.5, alternative = "two.sided", conf.level = 0.95)
##  Exact binomial test
## 
## data:  x and n
## number of successes = 6, number of trials = 13, p-value = 1
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.1922324 0.7486545
## sample estimates:
## probability of success 
##              0.4615385 

#使用Hmisc包中的binconf()函数
#install.packages("Hmisc")
library(Hmisc)
binconf(x, n, alpha=0.05, method="all")
##             PointEst     Lower     Upper
## Exact      0.4615385 0.1922324 0.7486545
## Wilson     0.4615385 0.2320607 0.7085620
## Asymptotic 0.4615385 0.1905457 0.7325312
6.1.2.1.2 正态近似法

如果n足够大，那么二项分布的偏度就比较小。在这种情况下，如果使用适当的连续性校正， 那么B(n,p)的一个很好的近似是正态分布:N(np,np(1-p))，对应的样本率p也近似正态分布: N(p,p(1-p)/n)。当n越大（至少20）且p不接近0或1时近似效果更好。 不同的经验法则可以用来决定n是否足够大，以及p是否距离0或1足够远， 其中一个常用的规则是np和n(1 −p)都必须大于5时，可利用样本率p的分布近似正态分布来估计总体的1-α(α=0.05或0.01)的可信区间： \[(p-u_{a/2}S_p,p+u_{a/2}S_p) ,\] \(S_p\)为样本资料计算样本的“成功”率\(S_p=\sqrt{{p(1-p)}/{n}}，p=X/N\)，α=0.05时，\(u_{a/2}=1.96\)；α=0.01时，\(u_{a/2}=2.58\)

使用《医学统计学》中的 案例6-3 的数据在R中进行二项分布的正态近似法分析测试。

##用正态分布近似计算
n1 <- 100
x1 <- 55
p1 <- x1/n1
S_p1=sqrt(p1*(1-p1)/n1)
p1 + c(-1, 1)*qnorm(0.975)*S_p1
## [1] 0.452493 0.647507

##用精确法计算
binom.test(x1, n1, 0.5, alternative="two.sided", conf.level=0.95)
##  Exact binomial test
## 
## data:  x1 and n1
## number of successes = 55, number of trials = 100, p-value = 0.3682
## alternative hypothesis: true probability of success is not equal to 0.5
## 95 percent confidence interval:
##  0.4472802 0.6496798
## sample estimates:
## probability of success 
##                   0.55 

#使用Hmisc包中的binconf()函数
binconf(x1, n1, alpha=0.05, method="all")
##            PointEst     Lower     Upper
## Exact          0.55 0.4472802 0.6496798
## Wilson         0.55 0.4524460 0.6438546
## Asymptotic     0.55 0.4524930 0.6475070

6.1.2.2 样本率与总体率的比较

6.1.2.2.1 直接法

使用直接法时，需要注意是单侧检验还是双侧检验。一般的有：

  1. 单侧检验的情况：
    • 若是回答“差/低”的问题计算的“成功”至多是k次的概率，即是\(P(X≤k)\)；
    • 若是回答“优/高”的问题计算的的“成功”至少是k次的概率，即\(P(X≥k)\)的率。
  2. 双侧检验的情况： 主要是回答“有无差别”，及备择假设是\(p≠p_0\)是否成立，因此计算的是实际样本“成功”k次出现的概率，与更背离原假设的事件出现的概率和。

使用《医学统计学》中的 案例6-4的数据在R中进行二项分布的样本率与总体率的比较的直接法分析测试。

#比较的是至少9次，或者零一对象最多8次的概率，所以：
p_8<-pbinom(8,10,0.55)
1 - p_8
#[1] 0.0232571

根据计算结果，p(x>8)=p(x≥9)=0.0232571＜0.05，因此拒接原假设\(H_0:p=0.55\)，接受\(H_1:p＞0.55\)。

使用《医学统计学》中的 案例6-5的数据在R中进行二项分布的样本率与总体率的比较的直接法分析测试。

n2<-10
x2<-9
p2<-0.6
#比较的是至少9次，或者零一对象最多8次的概率，所以：
p_3<-pbinom(x2,n2,p2)
1 - p_3
# [1] 0.0232571

根据计算结果，p(x>8)=p(x≥9)=0.0232571＜0.05，因此拒接原假设\(H_0:p=0.55\)，接受\(H_1:p＞0.55\)。

n3<-10
x3<-9
p3<-0.6
#计算的是出现9次的概率：
p4<-dbinom(x3,n3,p3)
# [1] 0.04031078
#计算出现更背离原假设的事件及概率，即：
p4_1<-sum(dbinom(0:n3,n3,p3)[which(dbinom(0:n3,n3,p3)<dbinom(x3,n3,p3))])
# [1] 0.01834117
#双侧检验的p值为：
p4+p4_1
# [1] 0.05865196

#或者直接用which()函数一并计算
sum(dbinom(0:n3,n3,p3)[which(dbinom(0:n3,n3,p3)<=dbinom(x3,n3,p3))])
# [1] 0.05865196

根据计算结果，p(x)=0.05865196＞0.05，因此不拒接原假设\(H_0:p=0.6\)，尚不能认为两种处理的效果不同。

6.1.2.2.2 正态近似法

当n较大、p和1-p均不太小，如np和n(1-p)均大于5时，利用样本率的分布近似正态分布的原理， 可作样本所在的总体率p与已知总体率\(p_0\)的比较。检验统计量u值的计算公式为： \[u=\frac{p-p_0}{\sqrt{p_0(1-p_0)/n}}\]

使用《医学统计学》中的 案例6-6的数据在R中进行二项分布的样本率与总体率的比较的正态近似法分析测试。

n4<-180
x4<-117
p4<-x4/n4
p4_0<-0.45

u<-(p4-p4_0)/sqrt( p4_0*(1-p4_0)/n4)
# [1] 5.393599
#根据正态分布的该v累计分布函数计算出p值
pnorm(u,lower.tail = F)*2
# [1] 6.906029e-08

#整理成函数的形式
binom_comp_uAsym<-function(n1,x1,p0){
  p1<-x1/n1
  p_val<-pnorm((p1-p0)/sqrt( p0*(1-p0)/n1),lower.tail = F)*2
  return( p_val)
}

binom_comp_uAsym(n4,x4,p4_0)
# [1] 6.906029e-08

根据计算结果，p(x)=6.906029e-08＜0.01＜0.05，因此拒接原假设\(H_0:p=0.45\)，接受\(H_1\)，两种处理的效果不同，新的比旧的好。

6.1.2.3 两样本率的比较(正态近似)

当n1与n2均较大，p1、p2、1-p1和1-p2均不太小，如n1p1、n1(1-p1)、n2p2、n2(1-p2)均大于5时，可利用样本率的分布 近似正态分布且独立两正态变量之差也服从正态分布的性质，采用近似正态法对两总体率进行统计检验，u的计算公式为： \[\begin{aligned} u &= \frac{p_1-p_2}{S_{p_1-p_2}}\\ S_{p_1-p_2} &= \sqrt{\frac{X_1+X_2}{n_1+n_2}(1-\frac{X_1+X_2}{n_1+n_2})(\frac{1}{n_1}+\frac{1}{n_2})} \end{aligned}\]

使用《医学统计学》中的 案例6-7的数据在R中进行二项分布的样本率与总体率的比较的正态近似法分析测试。

n5<-120
n6<-110
x5<-36
x6<-22

p5<-x5/n5
p6<-x6/n6

S_5_6<-sqrt((x5+x6)/(n5+n6)*(1-(x5+x6)/(n5+n6))*(1/n5+1/n6))
u<-(p5-p6)/S_5_6
pnorm(u,lower.tail = F)*2
# [1] 0.08107076

#整理成函数的形式
binom_twosamplecomp_uAsym<-function(n1,x1,n2,x2){
  p1<-x1/n1
  p2<-x2/n2
  S_val<-sqrt((x1+x2)/(n1+n2)*(1-(x1+x2)/(n1+n2))*(1/n1+1/n2))
  p_val<-pnorm((p1-p2)/S_val,lower.tail = F)*2
  return( p_val)
}
binom_twosamplecomp_uAsym(n5,x5,n6,x6)
# [1] 0.08107076

根据计算结果，p值=0.08107076＞0.05＞0.01，因此不拒接原假设\(H_0\)，尚不能认为两种处理的效果不同。

6.1.2.4 非遗传疾病的家族聚集性

非遗传性疾病的家族集聚性(Clustering in familes)，系指该疾病的发生在家族成员间是否有传染性，如果没有传染性，则家族成员间的患病是独立的，该疾病无家族集聚性。 否则，家族成员间的患病是非独立的，该疾病存在家族集聚性。

当这种遗传性疾病无家族集聚性时，以家族为样本，在n个成员中，出现X个成员患病的概率分布可认为服从二项分布; 否则，使不服从二项分布。 实际资料与二项分布进行拟合优度(Goodness of fit)的卡方检验得出p值。

6.1.2.5 群检验

为了解决检验大批标本的阳性率问题，把N个标本分为n个群，每个群m个标本，即N=n*m。检验每个群是否为阳性群（一旦检测到阳性就停止检测当前群），只有阴性群才需要检测所有标本，可以大大减少检测数目。 记每个标本阳性率为p，则\((1-p)^m\)为m个标本均为阴性的概率，\(1-(1-p)^m\)为一个群为阳性群的概率。 通过阳性群率计算阳性率：假设受检的n个群中，X个群为阳性群，则X/n可作为一个群为阳性群概率的估计值，则有：

\[(1-p)^m=\frac{X}{n}\]

从而估计出：

\[p=1-\sqrt[m]{1-X/n}\]

使用《医学统计学》中的 案例6-9的数据在R中进行二项分布的群检验分析测试。

N<-6000
n<-300
m<-20
X<-270
1-(1-X/n)^(1/m)
# [1] 0.1087491

#整理成函数的形式
binom_commtest<-function(N,n,m,X){
  return( 1-(1-X/n)^(1/m))
}
binom_commtest(N,n,m,X)
# [1] 0.1087491

根据计算结果，该地区的钉螺血吸虫的感染率的估计值为10.87%。

6.2 泊松分布

泊松分布，也叫Poisson分布，是二项分布的一种极端(即n很大,p很凶啊)情况，已发展为描述小概率事件发生规律的一种重要分布，适合于描述单位时间内随机事件发生的次数的概率分布。如某一服务设施在一定时间内受到的服务请求的次数，电话交换机接到呼叫的次数、汽车站台的候客人数、机器出现的故障数、DNA序列的变异数、放射性原子核的衰变数、发病率低的非传染性疾病发病或人数分布等、单位时间或面积、空间某罕见事物的分布等。如果随机变量X服从总体均数为\(\lambda\)的泊松分布，X取值为0，1，2，···的对应概率质量函数为： \[P(X)=\frac{e^{-\lambda}\lambda^X}{X!}, X=0,1,2,\dots\] 其中，\(\lambda\)为总体均数，e=2.71828为自然常数，有\(\sum{P(X)}=1\)。

适用的条件 1. 普通性 在充分小的观测单位上X的取值最多为1。简单而言，就是在试验次数n足够大时，每次试验可看作是一个“充分小的观测单位”，且每次试验只会发生两种互斥的可能结果之一(阳性或阴性)，这样阳性数X的取值最多为1。 2. 独立增量性 在某个观测单位上的取值与前面各观测单位上的取值无关。简而言之，就是前 面的试验结果不影响下一次的试验结果，各次试验具有独立性。 3. 平稳性 X的取值只与观测单位的大小有关，而与观测单位的位置无关。简单而言，就是每一次 试验阳性事件发生的概率都应相同，为x=/n，这样阳性数 的取值只与重复试验的次数有关，为合计的阳性数，可看作是大量独立试验的总结果。

6.2.1 泊松分布的性质

6.2.1.1 均数与标准差

  1. 总体均数\(\lambda\)与总体方差\(\sigma^2\)相等
  2. n很大p很小，且\(np=\lambda\)为常数时，二项分布近似Poisson分布

6.2.1.2 泊松分布的图像

  1. Poisson分布的图形是右偏分布，均值越大时，其对称性越好
  2. 当\(\lambda≤1\)时，随X取值越大，P(X)取值反而越小；当\(\lambda＞1\)时，随X取值越大，P(X)取值先正达后变小；
  3. 当\(\lambda\)为整数时，则在\(X=\lambda\)和\(X=\lambda-1\)处有最大概率值
  4. \(\lambda\)越大，Poisson分布渐近正态分布，\(\lambda>=20\)时可近似为正态分布
  5. Possion分布具有可加性(和正态分布类似)，但不具有可乘性(可由X取值和均数，方差看出)

6.2.2 泊松分布的应用

利用泊松分布及其正态近似性，可进行总体率的区间估计和差异推断。

6.2.2.1 总体均数的区间估计

利用服从泊松分布的样本资料，可进行总体均属的1-α可信区间，α取值为0.05或0.01。

6.2.2.1.1 直接法

对于n≤50的服从泊松分布的小样本资料，用查泊松分布表可以获取百分率的可信区间，即了得到其总体的1-α的可信区间。

R语言中的泊松分布相关的函数与之前介绍的分布(如二项分布，t分布和f分布)类似，主要有 dpois(),ppois(),qpois()和rpois()。 用于进行参数估计和假设检验的函数主要有poisson.test()。

使用《医学统计学》中的 案例6-10 的数据在R中进行二项分布的精确法/查表法分析测试。

n <- 1
x <- 21
poisson.test(x, n,conf.level=0.95)
##  Exact Poisson test
## 
## data:  21 time base: 1
## number of events = 21, time base = 1, p-value = 1
## alternative hypothesis: true event rate is not equal to 21
## 95 percent confidence interval:
##  12.99933 32.10073
## sample estimates:
## event rate 
##         21 
poisson.test(x, n,conf.level=0.99)
##  Exact Poisson test
## 
## data:  21 time base: 1
## number of events = 21, time base = 1, p-value = 1
## alternative hypothesis: true event rate is not equal to 21
## 99 percent confidence interval:
##  11.06923 35.94628
## sample estimates:
## event rate 
##         21 

#使用exactci包中的poisson.exact()方法，可以使用plot参数绘制出区间图
#install.packages("exactci")
library("exactci")
poisson.exact(21, plot=T, conf.level=0.95)
poisson.exact(21, plot=T, conf.level=0.99)

根据计算结果，95%和99%的可信区间分别是(12.99933,32.10073)和(11.06923,35.94628)。

6.2.2.1.2 正态近似法

对于n＞50的服从泊松分布的样本资料，可以采用正态近似法获取百分率的可信区间，即了得到其总体的1-α的可信区间。 \[(X-u_{\alpha/2}\sqrt{X},X+u_{\alpha/2}\sqrt{X})\]

使用《医学统计学》中的 案例6-11 的数据在R中进行泊松分布的正态近似法分析测试。

x<-68
alpha1<-0.05
alpha2<-0.01

posi_uAsym<-function(X,alpha) {
  Low<-X-qnorm(alpha/2)*sqrt(X)
  Up<-X+qnorm(alpha/2)*sqrt(X)
  return(list(Low,Up))
}

posi_uAsym(x,alpha1)
# 84.16228
# 51.83772
posi_uAsym(x,alpha2)
# 89.24083
# 46.75917

6.2.2.2 样本均数与总体均数的比较

6.2.2.2.1 直接法

当总体均数\(\lambda\)<20\(时，可采用直接计算概率的方式对样本均数与已知总体均数间的差别进行有无 统计学意义的比较。如用于发病率很低的非传染性疾病，则是对以样本计数X为代表的总体率p与已知的总 体率\)p_0$，是否有差别进行推断。

使用《医学统计学》中的 案例6-12 的数据在R中进行样本均数与总体均数的泊松分布-直接法分析测试。

p0<-0.008
n<-120
X<-4
lambda<-n*p0

#计算不小于4的泊松分布概率和
ppois(X-1, lambda=lambda,lower.tail = F)
#或者
1-sum(dpois(0:(X-1), lambda=lambda))
# [1] 0.01663305

根据计算结果，p值为0.0166<0.05，按照α=0.05的检验水准，拒绝原假设\(H_0\)，接受\(H_1\)，即认为母亲吸烟会增大小孩的先天心脏病危险。

6.2.2.2.2 正态近似法

当总体均数\(\lambda\)≥20$时，可采用正态近似法的方式对样本均数与已知总体均数间的差别进行有无 统计学意义的比较，检验统计量u值的计算公式为： \[u=\frac{X-\lambda}{\sqrt{\lambda}}\]

使用《医学统计学》中的 案例6-13的数据在R中进行样本率与总体率的比较的泊松分布-正态近似法分析测试。

p0<-0.003
n=25000
X<-123

posi_comp_uAsym<-function(n,X,p0) {
  u_val<-(X-n*p0)/sqrt(n*p0)
  #print (u_val)
  p_val<-pnorm(u_val,lower.tail = F)*2
  return(p_val)
}

posi_comp_uAsym(n,X,p0)
# [1] 2.980768e-08

根据计算结果，p值为2.980768e-08<0.05，按照α=0.05的检验水准，拒绝原假设\(H_0\)，接受\(H_1\)。

6.2.2.3 两样本均数的比较(正态近似)

对服从 Poisson 分布的样本，其样本计数可看作是样本均数。两个样本均数的比较，目的在于推断两样 本所代表的两总体均数是否有差别。设两个样本计数分别为\(X_1\)和\(X_2\)，可利用正态近似法进行比较，有下面两种情况：

两个样本的观察单位数相等，即\(n_1=n_2\) 当\(X_1+X_2≥20\)时： \[u=\frac{X_1-X_2}{\sqrt{X_1+X_2}}\] 当\(5＜X_1+X_2＜20\)时： \[u=\frac{|X_1-X_2|-1}{\sqrt{X_1+X_2}}\]

两个样本的观察单位数不相等，即\(n_1≠n_2\)： 当\(X_1+X_2≥20\)时： \[u=\frac{\bar{X_1}-\bar{X_2}}{\sqrt{X_1/n_1^2+X_2/n_2^2}}\] 当\(5＜X_1+X_2＜20\)时： \[u=\frac{|\bar{X_1}-\bar{X_2}|-1}{\sqrt{X_1/n_1^2+X_2/n_2^2}}\] 其中的\(\bar{X_1}\)和{X_2}分别是\(X_1/n_1\)和\(X_2/n_2\)。

使用《医学统计学》中的 案例6-14，6-15的数据在R中进行两样本均数的比较的泊松分布-正态近似法分析测试。

n1<-1
n2<-1
x1<-4
x2<-7

posi_twosamplecom_uAsym <-function(n1,n2,x1,x2){
  u_val<-NULL
  if(n1==n2){
    if((x1+x2)>=20)
      u_val<-(x1-x2)/sqrt(x1+x2)
    else if (5<(x1+x2) && (x1+x2)<20) 
      u_val<-(abs(x1-x2)-1)/sqrt(x1+x2)
  }else{
    x1m<-x1/n1
    x2m<-x2/n2
     if((x1+x2)>=20)
      u_val<-(x1m-x2m)/sqrt(x1/n1^2+x2/n2^2)
    else if (5<(x1+x2) && (x1+x2)<20) 
      u_val<-(abs(x1m-x2m)-1)/sqrt(x1/n1^2+x2/n2^2)
  }
  #print(u_val)
  return(pnorm(u_val,lower.tail = F)*2)
}
#案例6-14
posi_twosamplecom_uAsym(n1,n2,x1,x2)
# [1] 0.5464936

#案例6-15
posi_twosamplecom_uAsym(4,3,32,12)
# [1] 0.02845974

根据计算结果，案例6-14的p值为0.5464936>0.05，按照α=0.05的检验水准，不拒绝原假设\(H_0\)。

根据计算结果，案例6-15的p值为0.02845974<0.05，按照α=0.05的检验水准，拒绝原假设\(H_0\)，接受\(H_1H_1\)。

6.3 负二项分布

负二项分布是统计学上一种描述在一系列独立同分布的伯努利试验中， “成功”(或“失败”)次数到达指定次数（记为r）时“失败”(或“成功”)的次数(记为k)的离散概率分布。 根据参数r是否为整数，负二项分布有两种，当r为整数时，也称为帕斯卡分布(Pascal distribution)； 当r为延申到非整数时，也称为波利亚分布(Polya distribution)。 这里仅考虑r为整数时，也就是帕斯卡分布的情况，其定义有下面两种形式：

  1. 在一系列伯努利试验中，失败次数到达指定次数时，成功次数的离散概率分布
  2. 在一系列伯努利试验中，成功次数到达指定次数（指定值，记为r）时，失败次数（观察变量，记为k）的离散概率分布

这两种定义只是将“成功”和“失败”对调，其本质上没差别。一般的，为方便与之前的分布概念保持一直，便于理解， 都采用第二种形式。记r是成功的次数，为固定值；k是失败的次数，为自变量；p是伯努利试验成功的概率，失败概率则为1-p；因此下面将以第二种形式为例。 其概率质量函数为： \[\begin{aligned} f(k,r,p)=P(k,r,p)=p(x=k) & = \binom{k+r-1}{r-1}p^r(1-p)^k\\ & = \frac{(k+r-1)!}{k!(r-1)!}p^r(1-p)^k ,k=0,1,2,\dots\\ & = (-1)^k\binom{-r}{k}p^r(1-p)^k \end{aligned}\] 该表达式可以写成带负值参数\((-1)^k\)的二项系数的形式，解释了“负二项”名称的来源。

注意 有些资料对负二项分布的定义可能与这里的资料略有不同。最常见的变化是随机变量X计算不同的东西。 这些变化可以在这里的表格中看到。

6.3.1 负二项分布的性质

6.3.1.1 均数与标准差

负二项分布的概率分布依赖参数p和r，其总体的均值和方差，标准差为：

X的总体均数： \[μ=r\frac{1-p}{p}\] X的方差： \[σ^2=r\frac{1-p}{p^2}=\frac{μ}{p}=μ+μ^2/r\] X的标准差： \[σ=\frac{\sqrt{r(1-p)}}{p}\]

6.3.1.2 负二项分布的图像

  1. 负二项分布的图像偏态主要由p和n决定：
  2. 负二项分布 与 二项分布 的区别在于：“二项分布”是固定试验总次数N的独立试验中，成功次数k的分布；而“负二项分布”是所有到r次成功时即终止的独立试验中，失败次数k的分布。
  3. 正态近似 如果n足够大，那么分布的偏度就比较小

6.3.1.3 与其他分布的关系

  1. 负二项分布 与 二项分布 的区别在于：“二项分布”是固定试验总次数N的独立试验中，成功次数k的分布；而“负二项分布”是所有到r次成功时即终止的独立试验中，失败次数k的分布。
  2. 当r=1时，负二项分布退化为几何分布

6.3.2 负二项分布的参数估计

负二项分布有两个参数即u和r。r值越大，分布的方差与均数的比值就越接近 1； 而r值越小，分布的方差与均数的比值就越大。为此，可以用r值的大小来衡量 分布的离散程度即聚集 趋向的程度，常称为聚集指数(Cluster index)。

负二项分布的参数一般可用样本均数\(\bar{k}\)作为其估计值，即\(u=\bar{k}\)， 但参数r的估计就复杂一些。关于r的估计方法常用的有矩法、零频数法和最大似然法等。

矩法 用样本的均数\(\bar{k}\)和方差\(S^2\)分别作为负二项分布的均数μ和方差\(\sigma^2\)的估计值， 可以利用公式： \[\bar{r}=\frac{\bar{k}^2}{S^2-\bar{k}}, \bar{k}=\sum{fk}/N, S^2=(\sum{fk^2}-(\sum{fk})^2/N)/(N-1)\] 式子中的\(f\)为样本观察到的“失败”数k所对应的频数，N为观察但总数。

使用《医学统计学》中的 案例6-16 的数据在R中进行矩法估计聚集指数r(Cluster index)分析测试。

arr1<-array(dim=c(2,7))
arr1[1,]<-0:6
arr1[2,]<-c(30,14,8,4,2,0,2)

negbino_jufa_restim<-function(k_val,freq){
  N<-sum(freq)
  sum_fk<-sum(k_val * freq)
  sum_fkk<-sum(k_val^2 * freq)
  mean_val<-sum_fk/N
  var_val<-(sum_fkk-sum_fk^2/N)/(N-1)
  #print(paste(N, mean_val ,var_val, sep = "   "))
  r_estim<-mean_val^2/(var_val-mean_val)
  return(r_estim)
}

negbino_jufa_restim(arr1[1,],arr1[2,])
# [1] 1.033333

零频数法

零频法(Zero frequency method)是利用样本计数k=0时，所对应的频数\(f_0\)占观察单位总数的比例，即 \(f_0/N\)来估计r值，即求解满足下面的方程的r值： \[r\lg(1+\bar{k}/r)=\lg(N/f_0)\] 使用此方法要注意满足适当的条件。 - \(\frac{f_0}{N}＞\frac{1}{3}\) - 当\(\bar{r}<10\)时，还必须有充分的\(f_0\)项，以满足不等式：\((\bar{r}_0.17)(\frac{f_0}{N}-0.32)＞0.20\) - 往往需要使用线性内插法(Linear interpolation method)才能求得r的估计值。

最大似然法

最大似然法(Maximum likelihood method)是满足下面是指的z=0的r值： \[z=\sum_{k=0}^m\frac{A_k}{r+k}-N\ln(1+\bar{k}/r)\] 公式中\(m=k_max\)，即样本计数k所取得的最大值；\(A_k=\sum_{i=k+1}^mf_i\)，即样本中所有计数发育k的频数之和。

需要使用线性内插法(Linear interpolation method)才能求得r的估计值，该方法的估计更为精确，但是计算很复杂。

6.3.3 负二项分布的应用

R语言中的二项分布相关的函数与之前介绍的连续性分布(u分布，t分布和f分布)类似，主要有 dnbinom(),pnbinom(),qnbinom()和rnbinom()。

6.3.3.1 拟合优度检验

是通过理论频数和实际频数的比较，对样本分布拟合负二项分布的适合情况进行检验。类似二项分布。此方法需要先进行r值估计，然后计算出 不同样本频数的理论概率和理论频数，因此比较麻烦。目前在R语言中在那时没找到相应的包和函数可以比较好作者部分处理，后期看能否自己开发相应的工具代码吧。

具体案例和内容请先参考《医学统计学》P92页的内容。

6.3.3.2 两样本均数比较

对服从负二项分布的两个样本均数的比较，以推断其所代表的两个总体有无差别。先将两个样本的原始观察数据按公式： \[Y_i=\ln(k_i+0.5r_c)\] 进行转换，然后再对所得到两组转换数据做T检验。公式中，\(k_i\)为每一组独立试验的样本计数， \(r_c\)代表两样本的负二项分布的参数\(r_1\)和\(r_2\)的合并至，采用矩法估计的公式为： \[r_c=\frac{\bar{k1}^2(S_1^2-\bar{k1})+\bar{k2}^2(S_2^2-\bar{k2})}{(S_1^2-\bar{k1})^2+(S_2^2-\bar{k2})^2}\] 公式中，\(\bar{k1}，S_1\)分别是第一个样本计数的均数和方差；\(\bar{k2}，S_2\)分别是第二个样本计数的均数和方差。