基于R语言的假设检验（二）

上一期文章《基于R语言的假设检验（一）》中通过参数方法对假设检验做了一下梳理，包括单样本和两样本的均值、方差假设检验。然而实际应用中并不知道总体服从何种分布，对于这类数据就不可以使用参数方法，而要使用下文中即将讨论的几种非参数方法。

皮尔逊拟合优度卡方检验

拟合优度检验是检验来自总体中的一类数据其分布是否与某种理论分布相一致的统计方法。其核心思想包含两步，其一是将总体中的样本数据离散化（原本为离散数据就不需要离散化了）；其二比较各离散组中实际观测数与理论观测数之间的差异。该卡方统计量的数学表达式为：

其中ni为第i组中实际观测到的数量，n为样本总量，pi为第i组的理论概率。

使用皮尔逊拟合优度卡方检验可以做如下几件事：

1）判别某个总体的不同水平下是否存在显著差异（类似于方差分析）

2）检验某个总体是否服从某种理论的分布

3）检验两个离散总体（因子变量）间是否独立。

用于皮尔逊拟合优度卡方检验的R函数是chisq.test()，该函数的语法和参数如下：

chisq.test(x, y = NULL, correct = TRUE,

p = rep(1/length(x), length(x)),

rescale.p = FALSE,simulate.p.value = FALSE,

B = 2000)

x为分析的向量或矩阵，x和y也可以是因子向量，当x为2维矩阵时（二维列联表），该函数可以用来实现离散变量的独立性检验，但前提是矩阵内没有空单元，且各单元的频数至少为5；

correct为逻辑参数，指明是否用于连续修正，默认为修正；

p为各个离散组的理论概率，默认各组的理论概率相等；

simulate.p.value指明是否需要仿真计算P值，默认值为FALSE；

B表示仿真计算P值的次数，默认为2000次。

关于皮尔逊拟合优度卡方检验使用的几种场景做一个讲解和应用。

一、判别某个总体的不同水平下是否存在显著差异

1、提出原假设和备择假设

H0：各水平之间不存在显著差异，即可以理解为各水平下的理观测数（或理论概率）相同；

H1：各水平之间存在显著差异。

2、应用

为了了解市场中的5种品牌手机是否存在差异，现随机抽取1500人选择自己喜欢的品牌，收集回来的数据如下图所示，根据这5个数据是否能得到手机品牌之间是否存在差异？

按照皮尔逊拟合优度卡方检验的思想，做如下程序设计（自定义函数）：

my_chisq.test <- function(x, alpha = 0.05, pi = Inf){

df <- length(x)-1

ni <- x

if (pi == Inf) pi <- rep(1/length(x), length(x))

else pi = pi

npi <- sum(x)*pi

K <- sum((ni-npi)^2/npi)

P_Value <- 1-pchisq(q = K, df = df)

return(data.frame(K_Stat = K, df = df,

P_Value = P_Value))

}

该自定义函数就是针对离散情况下的显著性检验。

x <- c(231,300,450,108,411)

my_chisq.test(x)

从自定义返回的结果看，5种品牌间的手机还是存在显著差异的。当然可以直接使用R中的chisq.test()函数实现以上自定义函数功能。

chisq.test(x)

得到与自定义函数一致的结果。

二、检验某个总体是否服从某种理论的分布

这里用以下几种分布为例，检验样本是否服从均匀分布、正态分布、泊松分布和指数分布。

set.seed(1234)

x <- rnorm(1000)

y <- runif(1000)

z <- rpois(1000, lambda = 3)

e <- rexp(1000, rate = 2.4)

查看最小值和最大值，用于下面的分组

x.max <- max(x); x.min <- min(x);x.min;x.max

y.max <- max(y); y.min <- min(y);y.min;y.max

e.max <- max(e); e.min <- min(e);e.min;e.max

使用cut函数对数据分组

x.cut <- cut(x, breaks = c(-4,-3,-2,-1,0,1,2,3,4))

y.cut <- cut(y, breaks = c(0,0.2,0.4,0.6,0.8,1))

e.cut <- cut(e,breaks = c(0,1,2,3))

查看分组情况

x.table <- table(x.cut)

y.table <- table(y.cut)

e.table <- table(e.cut)

计算x和y各组的理论概率

均匀分布时

x.p <- punif(c(-3,-2,-1,0,1,2,3,4),

min = x.min, max = x.max)

p.x <- c(x.p[1], x.p[2]-x.p[1], x.p[3]-x.p[2],

x.p[4]-x.p[3],x.p[5]-x.p[4],

x.p[6]-x.p[5],x.p[7]-x.p[6],1-x.p[7])

chisq.test(x =x.table, p = p.x)

也可以使用ks.test()函数来验证连续变量是否服从某种分布，前提是需要知道总体理论分布。

ks.test(x,'punif',y.min,y.max)

确实验证了x不服从均匀分布

y.p <- punif(c(0.2,0.4,0.6,0.8,1),

min = y.min, max = y.max)

p.y <- c(y.p[1], y.p[2]-y.p[1], y.p[3]-y.p[2],

y.p[4]-y.p[3],1-y.p[4])

chisq.test(x =y.table, p = p.y)

ks.test(y,'punif',y.min,y.max)

验证y确实服从均匀分布，这里P值不一致是因为chisq.test()函数中涉及到了认为的分组，不同的分组会得到不一致的P值。

正态分布时

x.p <- pnorm(c(-3,-2,-1,0,1,2,3,4), mean(x), sd(x))

p.x <- c(x.p[1], x.p[2]-x.p[1], x.p[3]-x.p[2],

x.p[4]-x.p[3],x.p[5]-x.p[4],x.p[6]-x.p[5],

x.p[7]-x.p[6],1-x.p[7])

chisq.test(x =x.table, p = p.x)

ks.test(x,'pnorm',mean(x),sd(x))

这里警告信息的产生是因为分组过程中有的组样本量低于5个。

指数分布时(指数分布的参数rate为均值的倒数)

e.p <- pexp(c(1,2,3),rate = 1/mean(e))

p.e <- c(e.p[1],e.p[2]-e.p[1],1-e.p[2])

chisq.test(x =e.table, p = p.e)

ks.test(e,'pexp',1/mean(e))

泊松分布时(泊松分布的参数lambda等于均值)

z.table <- table(z)

z.mean <- sum(as.numeric(names(z.table))*z.table)/length(z)

z.p <- ppois(c(0,1,2,3,4,5,6,10), lambda = z.mean)

p.z <- c(z.p[1],z.p[2]-z.p[1],z.p[3]-z.p[2],

z.p[4]-z.p[3],z.p[5]-z.p[4],

z.p[6]-z.p[5],z.p[7]-z.p[6],1-z.p[7])

chisq.test(c(z.table[1:7],z2 = sum(z.table[8:10])),p.z)

警告信息的参数与之前的原因一致。请注意，对于离散变量，ks.test()函数就无法进行理论分布的假设检验了。

三、检验两个离散总体（因子变量）间是否独立

用chisq.test()函数检验两个离散总体（因子变量）间是否独立的话，需要x变量为一个矩阵。

m <- matrix(c(100,20,34,213,145,45), ncol = 3)

chisq.test(m)

由P值的返回值可知，行与列之间并不是独立的。

chisq.test()和ks.test()函数的优缺点：

使用chisq.test()函数对样本服从某种分布的假设检验时首先需要将样本离散化，这需要人工操作，这是该函数不足的地方，但其优点是可以验证离散变量服从某种分布的假设。ks.test()与之正好相反，即不需要人为的将样本进行分组，而且该函数还可以分析两个总体是否具有相同的分布，但也只局限于连续变量的分布假设。

ks.test()函数的语法和参数：

ks.test(x, y, ...,

alternative = c("two.sided", "less", "greater"),

exact = NULL)

x为分析的数值向量，y即可以是数值向量，也可以是字符串表示的某个累积分布函数(如'pnorm'，'pexp'等)，如果x和y都是数值向量时，ks.test()函数用来检验x和y是否具有相同的分布，如果y表示某个累积分布函数的字符串时，该函数就是用来验证变量x是否服从y表示的分布；

...用来指定分布函数y的参数；

alternative指定备择假设的三种形式，默认为双尾检验；

exact指定是否计算精确的p值。

在《基于R语言的假设检验（一）》一文中假设总体服从或近似服从正态分布时对单样本和两样本均值和方差的假设检验，但如果不知道总体服从何种分布时，可以使用Wilcoxon符号秩检验实现单样本或两样本均值检验。关于Wilcoxon符号秩检验的理论知识可以参考《统计建模与R语言》书中第5.3.7节内容。R中可以使用wilcox.test()函数实现Wilcoxon符号秩检验，有关该函数的语法和参数如下：

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为所需分析的数值向量，当x给定或x和y为配对数据时，mu表示x或x-y关于mu是对称的，默认mu值为0；

alternative指定备择假设假设的三种形式；

paired为逻辑参数，如果分析的两样本为配对数据时需要设置该参数为TRUE；

exact是否进行精确P值的计算；

correct是否进行连续修正；

conf.int是否返回置信区间；

conf.level指定置信水平。

应用

一、检验单样本是否关于mu对称

检验样本x1是否关于mu对称

set.seed(1234)

x1 <- rnorm(1000,mean = 4,sd = 3)

wilcox.test(x1, mu = 3.8)

wilcox.test(x1, mu = 3, alternative = 'less')

wilcox.test(x1, mu = 2.3, alternative = 'greater')

返回结果显示x1关于3.8对称，而且3确实小于4，2.3不大于4。

二、检验配对样本是否存在差异

set.seed(1234)

y1 <- rt(1000, 2)

y2 <- rt(1000, 5)

wilcox.test(y1, y2, paired = TRUE)

y1与y2并无显著差异

z1 <- rexp(1000, 2)

z2 <- rt(1000, 5)

wilcox.test(z1, z2, paired = TRUE)

wilcox.test(z1, z2, paired = TRUE, alternative = 'less')

z1-z2显著不为0，说明两者存在差异，进一步显示z1显著大于z2，即z1-z2>=0。

三、检验非配对样本之间的差异

q1 <- c(20,12,22,34,25,3,67,44,23,10)

q2 <- c(111,130,258,330,120,240)

wilcox.test(q1, q2)

wilcox.test(q1, q2, alternative = 'less')

q1与q2存在显著差异，进一步得到q1显著小于q2，即q1-q2 < 0。

总结：文中涉及到的R包和函数

stats包

chisq.test()

ks.test()

wilcox.test()