R语言各种假设检验实例整理

一、正态分布参数检验

   例1. 某种原件的寿命X(以小时计)服从正态分布N(μ, σ)其中μ, σ2均未知。现测得16只元件的寿命如下:

                 159 280 101 212 224 379 179 264
                  222 362 168 250 149 260 485 170

    问是否有理由认为元件的平均寿命大于255小时?

    解:按题意,需检验

                     H0: μ ≤ 225H1: μ >225

      此问题属于单边检验问题

      可以使用R语言t.test

      t.test(x,y=NULL,

          alternative=c(“two.sided”,”less”,”greater”),

          mu=0,paired=FALSE,var.equal=FALSE,

           conf.level=0.95)

      其中x,y是又数据构成e向量,(如果只提供x,则作单个正态总体的均值检验,如果提供x,y则作两个总体的均值检验),alternative表示被则假设,two.sided(缺省),双边检验(H1:μ≠H0),less表示单边检验(H1:μ<μ0),greater表示单边检验(H1:μ>μ0),mu表示原假设μ0,conf.level置信水平,即1-α,通常是0.95,var.equal是逻辑变量,var.equal=TRUE表示两样品方差相同,var.equal=FALSE(缺省)表示两样本方差不同。

      R代码:

       X<-c(159, 280, 101, 212, 224, 379, 179, 264,
              222, 362, 168, 250, 149, 260, 485, 170)
       t.test(X,alternative = “greater”,mu=225)

      结果:

         可见P值为0.257 > 0.05 ,不能拒绝原假设,接受H0,即平均寿命不大于225小时。

    例2.在平炉上进行的一项试验以确定改变操作方法的建议是否会增加刚的得率,试验时在同一个平炉上进行的,每炼一炉刚时除操作方法外,其它条件都尽可能做到相同,先用标准方法炼一炉,然后用新方法炼一炉,以后交替进行,各炼了10炉,其得率分别为

                           标准方法    78.1 72.4 76.2 74.3 77.4 78.4 76.0 75.5 76.7 77.3
                           新方法        79.1 81.0 77.3 79.1 80.0 79.1 79.1 77.3 80.2 82.1

设这两个样本相互独立,且分别来自正态总体N(μ1, σ2)和N(μ2, σ2),其中μ1,μ2和σ2未知。问新的操作能否提高得率?(取α=0.05)

    解1:根据题意,需要假设

              H0: μ1 ≥  μ2     H1: μ1 <  μ2  

    这里假定σ12=σ22=σ2,因此选择t.test,var.equal=TRUE

    R代码:

    X<-c(78.1,72.4,76.2,74.3,77.4,78.4,76.0,75.5,76.7,77.3)
    Y<-c(79.1,81.0,77.3,79.1,80.0,79.1,79.1,77.3,80.2,82.1)
    t.test(X,Y,var.equal = TRUE,alternative = “less”)

   结果:

    可见P值<0.05,接受备择假设,即新的操作能够提高得率。

    解2:

    因为数据是成对出现的,所以采用成对数据t检验比上述的双样本均值检验更准确。所谓成对t检验就是Zi=Xi-Yi,再对Z进行单样本均值检验

    R代码:

    X<-c(78.1,72.4,76.2,74.3,77.4,78.4,76.0,75.5,76.7,77.3)
    Y<-c(79.1,81.0,77.3,79.1,80.0,79.1,79.1,77.3,80.2,82.1)
    t.test(X-Y, alternative = “less”)

    结果:

可见P值<0.05,接受备择假设,即新的操作能够提高得率。并且P值更小可见比双样本均值检验更准确

  例3.对例2进行方差检验,方差是否相同

    解:根据题意,需检验

     H0: σ12 =  σ22H1: σ12 ≠  σ22

    方差检验可以用var.test

    var.test(x, y, ratio = 1,
            alternative = c(“two.sided”, “less”, “greater”),
            conf.level = 0.95, …)

    x,y是来自两样本数据构成的向量,ratio是方差比的原假设,缺省值为1.alternative是备择假设,two.sided表示双边检验(H1:σ1222<ratio),greater表示单边检验(H1:σ12

    R代码:

    X<-c(78.1,72.4,76.2,74.3,77.4,78.4,76.0,75.5,76.7,77.3)
    Y<-c(79.1,81.0,77.3,79.1,80.0,79.1,79.1,77.3,80.2,82.1)  

    var.test(X,Y)

    结果:

   可见P值为0.559>0.05,接受原假设,认为两者方差相同

二、二项分布参数检验

    例4.有一批蔬菜种子的平均发芽率p0=0.85,现随即抽取500粒,用种衣剂进行浸种处理,结果有445粒发芽。试检验种衣剂对种子发芽率有无效果。

    解:根据题意,所检验的问题为

                H0:p=p0=0.85, H1:p≠p0

            可以用R语言的binom.test

           binom.test(x, n, p = 0.5,
                    alternative = c(“two.sided”, “less”, “greater”),
                     conf.level = 0.95)

     其中x是成功的次数;或是一个由成功数和失败数组成的二维向量。n是试验总数,当x是二维向量时,此值无效。P是原假设的概率。

    R语言代码:

   binom.test(445,500,p=0.85)

   结果:

    可知P值0.01207<0.05,拒绝原假设,说明种衣剂对种子的发芽率有显著效果。

三、其它重要的非参数检验法

  3.1.理论分布完全已知的情况下

     3.1.1.皮尔森拟合优度检验

    例5.某消费者协会为了确定市场上消费者对5种品牌啤酒的喜好情况,随即抽取了1000名啤酒爱好者作为样品进行试验:每个人得到5种品牌的啤酒各一瓶,但未标明牌子。这5种啤酒分别按着A、B、C、D、E字母的5张纸片随即的顺序送给每一个人。下表是根据样本资料整理的各种品牌啤酒爱好者的频数分布。试根据这些数据判断消费者对这5种品牌啤酒的爱好有无明显差异?

       最喜欢的牌子    A          B          C          D          E

        人数X               210      312       170      85        223

    解:如果消费者对5种品牌的啤酒无显著差异,那么,就可以认为喜好这5种拍品啤酒的人呈均匀分布,即5种品牌啤酒爱好者人数各占20%。据此假设

                  H0:喜好5种啤酒的人数分布均匀

     可以使用Pearson χ2拟合优度检验,R语言中调用chisq.test(X)

     chisq.test(x, y = NULL, correct = TRUE,
                   p = rep(1/length(x), length(x)), rescale.p = FALSE,
                   simulate.p.value = FALSE, B = 2000)

      其中x是由观测数据构成的向量或者矩阵,y是数据向量(当x为矩阵时,y无效)。correct是逻辑变量,标明是否用于连续修正,TRUE(缺省值)表示修正,FALSE表示不修正。p是原假设落在小区间的理论概率,缺省值表示均匀分布,rescale.p是逻辑变量,选择FALSE(缺省值)时,要求输入的p满足和等于1;选择TRUE时,并不要求这一点,程序将重新计算p值。simulate.p.value逻辑变量(缺省值为FALSE),当为TRUE,将用仿真的方法计算p值,此时,B表示仿真的此值。

    R语言代码:

   X<-c(210, 312, 170, 85, 223)

   chisq.test(X)

    结果:

    例6.为研究电话总机在某段时间内接到的呼叫次数是否服从Poisson分布,现收集了42个数据,如下表所示,通过对数据的分析,问能否确认在某段时间内接到的呼叫次数服从Poisson分布(α = 0.1)?

    解:R语言代码:

    #输入数据
    X<-0:6; Y<-c(7, 10, 12, 8, 3, 2, 0)
    #计算理论分布,其中mean(rep(X,Y))为样本均值
    q<-ppois(X, mean(rep(X,Y))); n<-length(Y)
    p=rep(0,n)
    p[1]<-q[1]; p[n]<-1-q[n-1]
    for (i in 2:(n-1))
       p[i]<-q[i]-q[i-1]
    #作检验
    chisq.test(Y,p=p)

提示结果可能不准确,因为皮尔森卡方拟合由度检验要求分组后每组的频数至少要大于等于5,而后三组中出现的频率分别为3,2,0,均小于5,解决问题的方法是将后三组合成一组,此时的频数为5,满足要求,重写R语言代码

    R语言代码:

   #输入数据
    X<-0:6; Y<-c(7, 10, 12, 8, 3, 2, 0)
   #计算理论分布,其中mean(rep(X,Y))为样本均值
    q<-ppois(X, mean(rep(X,Y))); n<-length(Y)
    p<-rep(0,n)
    p[1]<-q[1]; p[n]<-1-q[n-1]
    for (i in 2:(n-1))
      p[i]<-q[i]-q[i-1]
    #重新分组
    Z<-c(7, 10, 12, 8, 5)
    #重新计算理论分布
    n<-length(Z); p<-p[1:n-1]; p[n]<-1-q[n-1]
    #作检验
    chisq.test(Z,p=p)

可见P值>>0.1,可以确认在某段时间之内接到的电话次数服从Poisson 分布

 3.1.2.正态W检验

    例7.已知15名学生体重如下,问是否服从正态分布

    解:

    R语言代码:

    w <- c(75.0, 64.0, 47.4, 66.9, 62.2, 62.2, 58.7, 63.5,
       66.6, 64.0, 57.0, 69.0, 56.9, 50.0, 72.0)
    shapiro.test(w)

   P值>0.05,接受原假设,认为来自正态分布总体。

3.2.理论分布依赖于若干个未知参数的情况

    3.2.1Kolmogorov-Smirnov 检验

    例8.对一台设备进行寿命检验,记录10次无故障工作时间,并按从小到大的次序排列如下:(单位)

       420  500 920  1380  1510  1650  1760 2100 2300 2350

    试用Kolmogorov-Smirnov K 检验方法检验此设备无故障工作时间分布是否服从λ = 1/1500的指数分布?

   解: R语言进行Kolmogorov-Smirnov K 检验使用ks.test( )

    ks.test(x, y, …,
    alternative = c(“two.sided”, “less”, “greater”),
    exact = NULL)    # x是待检测的样品构成的向量,y是原假设的数据向量或是原假设的字符串。

    R语言代码:

        X<-c(420, 500, 920, 1380, 1510, 1650, 1760, 2100, 2300, 2350)
        ks.test(X, “pexp”, 1/1500)

    P值大于0.05,无法拒绝原假设,因此认为此设备无故障工作时间的分布服从λ = 1/1500的指数分布。

    例9.假定从分布函数未知的F(x)和G(x)的总体中分别抽出25个和20个观察值的随即样品,其数据由下表所示。现检验F(x)和G(x)是否相同。

     R语言代码:

        X<-scan( )
        0.61 0.29 0.06 0.59 -1.73 -0.74 0.51 -0.56
        1.64 0.05 -0.06 0.64 -0.82 0.37 1.77
        2.36 1.31 1.05 -0.32 -0.40 1.06 -2.47
        0.39 1.09 -1.28
        Y<-scan( )
        2.20 1.66 1.38 0.20 0.36 0.00
        0.96 1.56 0.44  1.50 -0.30 0.66       

        2.31 3.29 -0.27 -0.37 0.38 0.70
        0.52 -0.71
        ks.test(X,Y)

P值>0.05,无法拒绝原假设,说明F(x)和G(x)分布函数相同。

    3.2.2.列联表数据的检验

    例10.为了研究吸烟是否与患肺癌相关,对63位肺癌患者及43名非肺癌患者(对照组)调查了其中的吸烟人数,得到2×2列联表,如下表所示

   解:

        进行Pearson卡方检验

        R语言代码:

        x<-c(60, 3, 32, 11)
        dim(x)<- c(2,2)
        chisq.test(x,correct = F)

P值<0.05,拒绝原假设,认为吸烟与患肺癌相关。

    例11.某医师为研究乙肝免疫球蛋白预防胎儿宫内感染HBV的结果,将33例HBsAg阳性孕妇随即分为预防注射组和对照组,结果由下表所示,问两组新生儿的HBV总体感染率有无差别?

  解: 最小期望值T=(11*9)/33=3 < 5 ,因此不能使用卡方检验,使用fisher检验,在R语言中使用fisher.test( )

    fisher.test(x, y = NULL, workspace = 200000, hybrid = FALSE,
        control = list(), or = 1, alternative = “two.sided”,
         conf.int = TRUE, conf.level = 0.95)

    #其中x是具有二维列联表形式的矩阵或是由因子构成的对象,y是由因子构成的对象,当x是矩阵时,此值无效。workspace的输入值时一整数,其整数表示用于网络算法空间的大小。hybrid为逻辑变量,FALSE(缺省值)表示精确计算概率,TRUE表示用混合算法计算概率。alternative为备择,有”two.sided”(缺失值)双边,”less”单边小于,”greater”单边大于,conf.int逻辑变量,当conf.int=TRUE(缺省值),给出 区间估计。conf.level为置信水平,缺省值为0.95,其余参数见在线说明。

    R语言代码:

    x<-c(4,5,18,6); dim(x)<-c(2,2)
    fisher.test(x)

  可见P值>0.05,接受原假设,认为两变量是独立的,即两组新生儿的HBV总体感染率无差别   

    例12.某胸科医院同时用甲乙两种方法测定202份痰样本中的抗酸杆菌,结果如下表所示,问甲、乙两种方法检出率有无差异。

解:因为是在相同个体上进行的两次检验,因此使用McNemar检验,

    H0:对相同痰样本测定中,甲乙两种方法检出率没有差异。

    在R语言中进行McNemar检验用到函数mcnemar.test( )

    mcnemar.test(x, y = NULL, correct = TRUE)  #其中x是具有二维列联表形式的矩阵或是由因子构成的对象。y是由因子构成的对象,当x是矩阵时,此值无效。correct是逻辑变量,TRUE(缺省值)表示在计算检验统计量时用连续修正,FALSE是不用修正。

    R语言代码:

    X<-c(49, 21, 25, 107); dim(X)<-c(2,2)
    mcnemar.test(X,correct=FALSE)

P值> 0.05,不能认为两种检测方法有差异。

3.2.3.符号检验

    例13.联合国人员在世界上66个大城市的生活花费指数(以纽约市1996年12月为100)按自小至大的次序排列如下(这里北京的指数为99):

       66 75 78 80 81 81 82 83 83 83 83
       84 85 85 86 86 86 86 87 87 88 88
       88 88 88 89 89 89 89 90 90 91 91
       91 91 92 93 93 96 96 96 97 99 100
       101 102 103 103 104 104 104 105 106 109 109
       110 110 110 111 113 115 116 117 118 155 192

    假设这个样品是从世界许多大城市中随即抽样得到的。试用符号检验分析,北京是在中位数之上,还是在中位数之下。

    解:样本的中位数(M)作为城市生活水平的中间值,因此需要检验:

                     H 0 : M ≥ 99,  H 1 : M < 99.

     输入数据,作二项检验。

     R语言代码:

     X <- c(66,75, 78 ,80 ,81 ,81 ,82, 83, 83, 83, 83,
           84 , 85, 85, 86, 86, 86, 86, 87 ,87, 88, 88,
           88, 88, 88, 89 ,89, 89, 89, 90 ,90 ,91 ,91,
           91 ,91, 92, 93, 93, 96, 96, 96, 97, 99, 100,
          101, 102, 103, 103 ,104, 104, 104 ,105, 106, 109, 109,
           110 ,110, 110, 111, 113, 115, 116, 117 ,118, 155 ,192)
     binom.test(sum(X>99), length(X), al=”l”)

    在程序中,sum(x>99)表示样本中大于99的个数。al是alternative的缩写,”l”是”less”的缩写。计算出的P值小于0.05,拒绝原假设,也就是说,北京的生活水平高于世界的中位水平。

例14.用两种不同的饲料养猪,其增重情况如下表所示。试分析两种饲料养猪有无显著差异。

     R语言代码:

     x<-scan()
     25 30 28 23 27 35 30 28 32 29 30 30 31 16

     y<-scan()
     19 32 21 19 25 31 31 26 30 25 28 31 25 25

     binom.test(sum(x<y), length(x))

sum(x < y)表示样品X小于样品Y的个数。计算出P值大于0.05,无法拒绝原假设,可以认为两种饲料养猪无显著差异。

例15.某饮料店为了解顾客对饮料的爱好情感,进一步改进他们的工作,对顾客喜欢咖啡还是喜欢奶茶,或者两者同样爱好进行了调查,顾客喜欢咖啡超过奶茶用正号表示,喜欢奶茶超过咖啡用负号表示,两者同样爱好用0表示。现将调查的结果列在下表中。试分析顾客是喜欢咖啡还是喜欢奶茶。

解:根据题意可检验如下假设:

    H0:顾客喜欢咖啡等于喜欢奶茶;      H1:顾客喜欢咖啡超过奶茶。

    以上资料中有以人(即6号顾客)表示对咖啡和奶茶有同样爱好,用0表示,因此在样本容量中不加计算,所以实际上N=12.如果H0假设为真,那么符合p为1/2的二项分布,如果H1为真,那么顾客喜欢奶茶的人数小于理论值,al=”l”,因此用R软件进行计算,显著性水平取α = 0.10,

    R语言代码:

    binom.test(3,12,p=1/2, al=”l”, conf.level = 0.90)

    可见P值 < 0.1 ,置信区间也不包括0.5,因此拒绝原假设人口喜欢咖啡的人超过喜欢奶茶的人。

    在符号检验法中,只计算符号的个数,而不考虑每个符号差所包含的绝对值的大小,因此常常使用弥补了这个缺点的wilcoxon符号秩检验。

3.3.4.符号秩检验

    例16.假定某电池厂宣称该厂生产的某种型号电池寿命的中位数为140安培小时。为了检验改厂生产的电池是否符合其规定的标准,现从新近生产的一批电池中抽取了随即样本,并对这20个电池的寿命进行了测试,其结果如下(单位:安培小时):

                 137.0 140.0 138.3 139.0 144.3 139.1 141.7 137.3 133.5 138.2

                 141.1 139.2 136.5 136.5 135.6 138.0 140.9 140.6 136.3 134.1

试用Wilcoxon符号秩检验分析该厂生产的电池是否符合其标准。

    解:根据题意假设:

                 H0:电池中位数M≥ 140安培小时;

                 H1:电池中位数<140安培小时。

          在R语言中进行符号秩检验可以使用wilcox.test( )

          wilcox.test(x, y = NULL,
              alternative = c(“two.sided”, “less”, “greater”),
              mu = 0, paired = FALSE, exact = NULL, correct = TRUE,
              conf.int = FALSE, conf.level = 0.95, …)

           其中x,y是观察数据构成的数据向量。alternative是备择假设,有单侧检验和双侧检验,mu待检参数,如中位数M0.paired是逻辑变量,说明变量x,y是否为成对数据。exact是逻辑变量,说明是否精确计算P值,当样本量较小时,此参数起作用,当样本两较大时,软件采用正态分布近似计算P值。correct是逻辑变量,说明是否对P值的计算采用连续性修正,相同秩次较多时,统计量要校正。conf.int是逻辑变量,说明是否给出相应的置信区间。

         R语言代码:

         X<-scan()
        137.0 140.0 138.3 139.0 144.3 139.1 141.7 137.3 133.5 138.2 141.1  139.2 136.5 136.5 135.6 138.0 140.9 140.6 136.3 134.1

         wilcox.test(X, mu=140, alternative=”less”,
               exact=FALSE,correct=FALSE, conf.int=TRUE)

    这里V=34是wicoxon的统计量,P值<0.05,即拒绝原假设,接受备择假设,中位值小于小于140安培小时。

    例17. 为了检验一种新的复合肥和原来使用的肥料相比是否显著提高了小麦的产量,在一个农场中选择了10块田地,每块等分为两部分,其中任指定一部分使用新的复合肥料,另一部分使用原肥料。小麦成熟后称得各部分小麦产量如下表所示。试用Wilcoxon符号检验法检验新复合肥是否会显著提高小麦的产量,并与符号检验作比较(α = 0.05)。

    解:根据题意作如下假设:

                     H0:新复合肥的产量与原肥料的产量相同,

                     H1:新复合肥的产量高于原肥料的产量。

           符号秩检验R语言代码:

           x<-c(459, 367, 303, 392, 310, 342, 421, 446, 430, 412)
           y<-c(414, 306, 321, 443, 281, 301, 353, 391, 405, 390)
           wilcox.test(x, y, alternative = “greater”, paired = TRUE)

         可见P值<0.05拒绝原假设,即新复合肥能显著提高小麦产量。

         符号检验R语言代码:

         x<-c(459, 367, 303, 392, 310, 342, 421, 446, 430, 412)
         y<-c(414, 306, 321, 443, 281, 301, 353, 391, 405, 390)
         binom.test(sum(x>y), length(x), alternative = “greater”)        

         用符号检验P值>0.05,因此在α = 0.05的水平下,就所给数据而言,符号检验还不足以区分两种化肥对提高小麦的产量产生差异。

例18.今测得10名非铅作业工人和7名铅作业工人的血铅值,如下表所示。试用Wilcoxon秩和检验分析两组工人血铅值有无差异。

    解:进行Wilcoxon秩和检验R语言同样可以使用wilcox.test( )

    R语言代码:

    x<-c(24, 26, 29, 34, 43, 58, 63, 72, 87, 101)

    y<-c(82, 87, 97, 121, 164, 208, 213)

    wilcox.test(x,y,alternative=”less”,exact=FALSE,correct=FALSE)

   P值小于0.05,拒绝原假设,即铅作业工人血铅值高于非作业工人。

例19.某医院用某种药物治疗两型慢性支气管炎患者共216例,疗效由下表所示,试分析该药物对两型慢性支气管炎的治疗是否相同。

    解:我们想象各病人的疗效用4个不同的值表示(1表示最好,4表示最差),这样就可以位这216名排序,因此,可用Wilcoxon秩和检验来分析问题。

    R语言代码:

    x<-rep(1:4, c(62, 41, 14,11)); y<-rep(1:4, c(20, 37, 16, 15))
    wilcox.test(x, y, exact=FALSE)

    P值<0.05,拒绝原假设,即认为该药物对两型慢性支气管炎的治疗是不相同的。因为数据有结点存在,故无法精确计算P值,其参数为exact=FALSE。 

 3.3.5.二元数据相关检验

    例20.某种矿石中两种有用成分A,B,取10个样品,每个样品中成分A的含量百分数x(%),及B的含量百分数y(%)的数据下表所示,对两组数据进行相关性检验。

    解:进行相关性检验,在R语言中可以使用cor.test( )

    cor.test(x, y,
         alternative = c(“two.sided”, “less”, “greater”),
         method = c(“pearson”, “kendall”, “spearman”),
         exact = NULL, conf.level = 0.95, …)

    #其中x,y是数据长度相同的向量,alternative是备择假设,缺省值为”two.sided”,method是检验方法,缺省值是Pearson检验,conf.level是置信区间水平,缺省值为0.95

    cor.test( )还有另一种使用格式

        cor.test(formula, data, subset, na.action, …) #其中formula是公式,形如’~u+v’ , ‘u’, ‘v’ 必须是具有相同长度的数值向量,data是数据框,subset是可选择向量,表示观察值的子集。   

   假设此例中两组数据均来自正态分布,使用pearson相关性检验,

   R语言代码:

    ore<-data.frame(
         x=c(67, 54, 72, 64, 39, 22, 58, 43, 46, 34),
         y=c(24, 15, 23, 19, 16, 11, 20, 16, 17, 13)
    )
    cor.test(ore$x,ore$y)

    可见P值<0.05,拒绝原假设,认为X与Y相关。

   例21.一项有六个人参加表演的竞赛,有两人进行评定,评定结果用下表所示,试用Spearman秩相关检验方法检验这两个评定员对等级评定有无相关关系。

    解:

    R语言代码:

    x<-c(1,2,3,4,5,6); y<-c(6,5,4,3,2,1)
    cor.test(x, y, method = “spearman”)

可见P值<0.05,拒绝原假设,认为x与y相关,rs=-1,表示这两个量是完全负相关,即两人的结论有关系,但完全相反。

参考资料

统计建模与R软件(上册)