统计检验
注意：这里着重于介绍相应函数的使用，可能会存在函数适用场景把握不准确的情形。 如发现有不妥当的地方，欢迎指教。 具体什么场合应该采用什么检验，可以参考相关的书，比如John H. McDonald的Handbook of Biological Statistics
正态性检验
诸如相关性分析、回归分析、t检验、方差分析等很多参数统计方法的前提的是数据程正态分布，因此有必要对数据的正态性进行查看或者检验。如果数据分布偏离正态分布，那么需要尝试对数据进行转换。如果转换之后更符合正态分布，那么使用转换之后的数据进行统计检验就更为准确。 常用的转换方式有取对数、取倒数、取平方根等，具体参照统计学相关书籍。 一般都认为参数检验对数据分布偏离正态分布不是特别敏感，有些偏离依然可以使用参数检验，除非偏离正态分布特别严重。如果两组都严重偏离正态分布且偏离的方向不一样，那可能选用非参数检验就更合适。然而，到底偏离多少才叫严重，这个没有人给个明确的说法。可以参照自己领域同行的常用选择，比如自己领域的同行都用非参数检验，那么即便你的数据看起来并没有那么偏离正态分布，跟大家一样使用非参数检验可能就更合适。
也有人认为正态性检验不是特别有用1)，因为对于对于小样本来说，正态性检验更趋向于得到不显著的p值，而对于大样本来说，正态性检验更倾向于得到显著的p值。当样本比较大的时候，数据不符合正态分布对诸如t检验等统计检验的结果的影响不大，可以不用进行正态性检验了。至于样本数比较大的标准，有的说每组样本数大于30或402)，有的说每组样本大于503)。
Histogram查看数据分布
可以用Histogram检查数据的正态性，对于偏离正态性比较远的数据，在进行统计检验前，需要进行转换，以便更适合使用相应的统计检验。如下例子为114个男生身高数据，通过Histogram可以看出其接近正态分布。
# 读入数据并转为数值型 sg <- readLines("ftp://ftp.biolab.wang/data/mheight.txt") sg <- as.numeric(sg) # 可以使用hist快速查看数据的分布情况 hist(sg) # 可以把具有相同的均值和方差的正态分布一起展示看看 library(ggplot2) ggplot() + geom_histogram(aes(x=sg, y = ..density..), binwidth=2.5) + stat_function(fun = dnorm, args = list(mean = mean(sg), sd = sd(sg)), colour="red", geom="line", size=2 )
比如，对于t检验来说，如果两个组的分布一致（比如都偏向左边），那么t检验的结果对其不符合正态分布这一现实并不敏感。
Histogram查看需要样本比较大，如果样本比较小，那么Histogram就不适合用来看数据是否符合正态分布了。
Shapiro-Wilk Normality Test
正态性检验可使用Shapiro-Wilk Normality Test，R里采用stats包里的shapiro.test函数，把一组数据组成的向量里用作shapiro.test的参数即可：
shapiro.test(sg) ## ## Shapiro-Wilk normality test ## ## data: sg ## W = 0.9884, p-value = 0.4405
如上例子从结果可以看到，p-value = 0.4405 > 0.05, 因此没有足够的证据表明这些数据的分布偏离了正态分布。
方差齐性检验
首先大家可能需要知道，在某个检验前先做方差齐性检验已经不再被统计学家广泛推荐4)。方差齐性检验有多种: Bartlett's Test, Levene's test, Fligner-Killeen's Test等。
# 比如我们这里随机产生三组数据 set.seed(1); (R1 <- rnorm(6, mean = 0, sd = 1)) ## [1] -0.6264538 0.1836433 -0.8356286 1.5952808 0.3295078 -0.8204684 set.seed(2); (R2 <- rnorm(6, mean = 0, sd = 1)) ## [1] -0.89691455 0.18484918 1.58784533 -1.13037567 -0.08025176 0.13242028 set.seed(2); (R3 <- rnorm(6, mean = 10, sd = 11)) ## [1] 1.030855 11.848492 25.878453 -1.303757 9.197482 11.324203
Bartlett's Test
Bartlett's Test 需要数据符合正态分布，R里使用bartlett.test函数，将要检验的几组数据放到一个list里用作参数即可，比如这里要检验R1、R2、R3三组数据之间的。 当然，也可以像car::leveneTest那样，第一个参数给一个向量，第二参数给对应的分组信息的因子。
X <- list(R1=R1, R2=R2, R3=R3) bartlett.test(X) ## ## Bartlett test of homogeneity of variances ## ## data: X ## Bartlett's K-squared = 27.636, df = 2, p-value = 0.0000009976
Levene's test
Levene's test是最为广泛运用的，相对于Bartlett Test来说对数据偏离正态分布更不敏感。R里使用car::leveneTest函数进行Levene's test，默认用法是第一个参数给一个数值型向量，第二参数为对应的分组的因子。
(Y <- c(R1, R2, R3)) ## [1] -0.62645381 0.18364332 -0.83562861 1.59528080 0.32950777 -0.82046838 -0.89691455 0.18484918 1.58784533 -1.13037567 -0.08025176 0.13242028 1.03085453 11.84849185 25.87845331 -1.30375674 9.19748243 11.32420284 (G <- rep(c("R1", "R2", "R3"), each=6)) ## [1] "R1" "R1" "R1" "R1" "R1" "R1" "R2" "R2" "R2" "R2" "R2" "R2" "R3" "R3" "R3" "R3" "R3" "R3" car::leveneTest(Y, factor(G)) ## Levene's Test for Homogeneity of Variance (center = median) ## Df F value Pr(>F) ## group 2 5.3059 0.01809 * ## 15 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
对于这里这个例子，leveneTest得到的p值要比bartlett.test大得多。
bartlett.test(X[2:3])$p.value ## [1] 0.0001246031 car::leveneTest(Y[7:18], factor(G[7:18]))[["Pr(>F)"]][1] ## [1] 0.04249441
比较明确知道R2组和R3两组数据之间的方差差距是比较大的，可以看到两种方法之间的差别还是比较大的。
Fligner-Killeen's Test
Fligner-Killeen's Test是非参数检验，无需数据符合正态分布。 R里使用fligner.test进行，用法与bartlett.test一致。
fligner.test(X) ## ## Fligner-Killeen test of homogeneity of variances ## ## data: X ## Fligner-Killeen:med chi-squared = 7.8684, df = 2, p-value = 0.01956
F-test
F-test用于两组之间的比较，需要数据符合正态分布。R里使用var.test进行F-test，简单地把两组数据的向量作为前两个参数即可。
var.test(R1, R2) ## ## F test to compare two variances ## ## data: R1 and R2 ## F = 0.95516, num df = 5, denom df = 5, p-value = 0.9611 ## alternative hypothesis: true ratio of variances is not equal to 1 ## 95 percent confidence interval: ## 0.1336565 6.8259430 ## sample estimates: ## ratio of variances ## 0.9551607 var.test(R1, R3) ## ## F test to compare two variances ## ## data: R1 and R3 ## F = 0.0095516, num df = 5, denom df = 5, p-value = 0.00009364 ## alternative hypothesis: true ratio of variances is not equal to 1 ## 95 percent confidence interval: ## 0.001336565 0.068259430 ## sample estimates: ## ratio of variances ## 0.009551607 var.test(R2, R3) ## ## F test to compare two variances ## ## data: R2 and R3 ## F = 0.01, num df = 5, denom df = 5, p-value = 0.0001049 ## alternative hypothesis: true ratio of variances is not equal to 1 ## 95 percent confidence interval: ## 0.00139931 0.07146382 ## sample estimates: ## ratio of variances ## 0.01
T检验
R里使用t.test函数进行T检验，设置不同的参数完成不同类型的T检验。使用? t.test查看具体的用法。两个用法：
## Default S3 method: t.test(x, y = NULL, alternative = c("two.sided", "less", "greater"), mu = 0, paired = FALSE, var.equal = FALSE, conf.level = 0.95, ...)
初学者会这种用法即可。大多数情况下，比较两组数据（分别列为向量X和Y）的均值之间是否有差异，那么用t.text(x=X, y=Y)或者简单地t.text(X, Y)即可。而有的时候需要根据特定的情况设置：
- 如果不设定y的值，那么默认
y = NULL，也就是进行单样本t检验，假设H0为x的均值是mu，如果不设置mu的值，则mu默认为0；
- 默认
alternative = “two.sided”，进行双尾检验，一般都需要进行双尾检验，除非有非常明确的理由可以选择单尾。alternative = “greater”则是单尾检验假设x的均值比y的均值要大；
paired = TRUE则进行配对t检验，x和y里的元素按照顺序一一对应的；
var.equal = FALSE为默认设置，进行Welch近似t检验（Welch’s t-test），这时候将分开计算x和y的方差，并且使用Welch法矫正自由度。当方差不等及两组样本大小不同时，Welch’s t-test更可靠，一般总是使用Welch’s t-test即可。
## S3 method for class 'formula' t.test(formula, data, subset, na.action, ...)
这种用法是把x和y放到一个向量（lhs）里，另外有一个向量（rhs）记录数字对应的分组，于是设置参数formula= lhs ~ rhs即可。这对于非配对t检验有时候会比较方便，对于配对t检验，用这种方法时需要lhs是一个“Pair” class的矩阵，相对就复杂一点，用上方的默认用法即可。
单样本T检验（One-sample Student's t-test）
| 例：某医生测量了36名从事铅作业男性工人的血红蛋白含量（如下表，单位：g/L），问从事铅作业男性工人的血红蛋白含量均数是否不等于正常成年男性的均数（140g/L）? | |||||||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 112 | 137 | 129 | 126 | 88 | 90 | 105 | 178 | 130 | 128 | 126 | 103 | 172 | 116 | 125 | 90 | 96 | 162 |
| 157 | 151 | 135 | 113 | 175 | 129 | 165 | 171 | 128 | 128 | 160 | 110 | 140 | 163 | 100 | 129 | 116 | 127 |
案例及数据引自颜艳, 王彤, 刘洪波, 等. 医学统计学[G]. 5. 人民卫生出版社, 2020.第45页例3-5 |
|||||||||||||||||
这种情况适用单样本T检验：
# 首先构建向量录入数据 X <- c(112, 137, 129, 126, 88, 90, 105, 178, 130, 128, 126, 103, 172, 116, 125, 90, 96, 162, 157, 151, 135, 113, 175, 129, 165, 171, 128, 128, 160, 110, 140, 163, 100, 129, 116, 127)
这里样本容量为36，正态性对检验的结果的影响已经比较小，可以不进行正态性检验了。如果用shapiro.test(X)进行检验得到p-value = 0.09881 > 0.05，我们不认为这组数据不符合正态分布，更加确定可以使用t检验。
# 采用单样本t检验，不给y参数，而是只给 t.test(X, mu=140) ## ## One Sample t-test ## ## data: X ## t = -2.1367, df = 35, p-value = 0.03969 ## alternative hypothesis: true mean is not equal to 140 ## 95 percent confidence interval: ## 122.1238 139.5428 ## sample estimates: ## mean of x ## 130.8333
显然这里得到p = 0.03969 < 0.05, 因此我们认为有足够的证据表明从事铅作业的男性工人血红蛋白含量不同于140g/L。
对于t检验的结果，我们可以从中提取p值，我们可以用str函数查看所得对象的结构，从而取出我们需要值，比如p值(p.value)。
t_test_result <- t.test(X, mu=140) str(t_test_result) ## List of 10 ## $ statistic : Named num -2.14 ## ..- attr(*, "names")= chr "t" ## $ parameter : Named num 35 ## ..- attr(*, "names")= chr "df" ## $ p.value : num 0.0397 ## $ conf.int : num [1:2] 122 140 ## ..- attr(*, "conf.level")= num 0.95 ## $ estimate : Named num 131 ## ..- attr(*, "names")= chr "mean of x" ## $ null.value : Named num 140 ## ..- attr(*, "names")= chr "mean" ## $ stderr : num 4.29 ## $ alternative: chr "two.sided" ## $ method : chr "One Sample t-test" ## $ data.name : chr "X" ## - attr(*, "class")= chr "htest" # 从结果中可以看出来，t.test得到的结果是个list，p.value是其一个元素，很容易就能够取出来 t_test_result$p.value ## [1] 0.03969288
独立样本T检验（Student’s t-test for two independent samples）
通常说“双样本t检验”，或者“独立样本t检验”，指的是Student’s t-test，也称为等方差t检验。这里需要留意，由于方差齐性检验可能不具有足够的检验效力，而使得应该被认为是异方差的没有被认为不是等方差，从而不恰当地使用了等方差t检验，从而增加了I类误差，因此，有人认为始终采用异方差t检验更为合适，详见异方差t检验一节。
| 例：为探讨青藤碱抗兔动脉粥样硬化作用，将16只造模成功的家兔随机分为2组，青藤碱组相较于模型组在饲料中添加青藤碱，其余条件不变，实验结束后测定心脏血液中高密度脂蛋白（HDL，mmol/L），结果见如下表。问能否认为模型组与青藤碱组心脏血液中HDL含量的总体均值不同? | ||||||||
|---|---|---|---|---|---|---|---|---|
| 青藤碱组X1 | 0.66 | 0.76 | 0.79 | 0.88 | 0.78 | 0.66 | 0.75 | 0.88 |
| 模型组X2 | 0.58 | 0.69 | 0.59 | 0.70 | 0.69 | 0.68 | 0.58 | 0.60 |
案例及数据引自颜艳, 王彤, 刘洪波, 等. 医学统计学[G]. 5. 人民卫生出版社, 2020.第50页例3-7，书中采用的是总体方差相等的t检验。 |
||||||||
# 首先构建两个向量录入数据 X1 <- c(0.66, 0.76, 0.79, 0.88, 0.78, 0.66, 0.75, 0.88) X2 <- c(0.58, 0.69, 0.59, 0.70, 0.69, 0.68, 0.58, 0.60) # 用histogram看看数据分布情况 hist(X1, xlab=NULL, breaks=5) hist(X2, xlab=NULL, breaks=5)
从直方图来看，模型组的数据看起来偏离了正态分布。这里先根据教材上的检验先忽略掉这种偏离，采用t检验。
# 采用F检验检验方差齐性 var.test(X1, X2) ## ## F test to compare two variances ## ## data: X1 and X2 ## F = 2.299, num df = 7, denom df = 7, p-value = 0.2944 ## alternative hypothesis: true ratio of variances is not equal to 1 ## 95 percent confidence interval: ## 0.4602708 11.4833515 ## sample estimates: ## ratio of variances ## 2.299011 # 上边F检验得到p>0.05，无足够证据表明两组数据方差不同， # 因此可以采用等方差t检验，参数x=X1, y=X2，检验X1和X2均值是否相等，因此保持mu=0，就不用专门写出来了 t.test(X1, X2, var.equal = TRUE) ## ## Two Sample t-test ## ## data: X1 and X2 ## t = 3.689, df = 14, p-value = 0.00243 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## 0.05494115 0.20755885 ## sample estimates: ## mean of x mean of y ## 0.77000 0.63875
显然这里得到p = 0.00243 < 0.05, 因此可以认为青藤碱组心脏血液中的HDL含量的总体均值高于模型组。
这里进行t检验忽略了两组数据的正态性，shapiro.test(X1)得到p-value = 0.2879，shapiro.test(X2)得到p-value = 0.01726。按照这个结果，认为模型组不符合正态分布。如果用wilcox.test(X1, X2)进行Wilcoxon检验得到p-value = 0.01331，比t.test(X1, X2, var.equal = TRUE)得到的p-value = 0.00243更大。对于这里这个例子，由于模型组偏离了正态分布，同时样本也比较小，根据《医学统计学》的说法，选择使Student's t test似乎有些不合适。当然，也有人认为检查正态性是没有必要的，因为这个检验本身可能也是不准确的。比如有的数据，我们知道它是正态分布的，那是不是自然就适合与t检验了呢？诸如这些问题，目前也没见有人说得很清楚。
配对样本T检验（Student's t-test for paired samples）
| 例：为比较两种方法对乳酸饮料中脂肪含量测定结果是否不同,随机抽取了10份乳酸饮料制品,分别用脂肪酸水解法和哥特里-罗紫法测定,结果如下表（%）。问两法测定结果是否不同? | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| 样本编号 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 |
| 脂肪酸水解法（X1） | 0.840 | 0.591 | 0.674 | 0.632 | 0.687 | 0.978 | 0.750 | 0.730 | 1.200 | 0.870 |
| 哥特里-罗紫法（X2） | 0.580 | 0.509 | 0.500 | 0.316 | 0.337 | 0.517 | 0.454 | 0.512 | 0.997 | 0.506 |
例题引自颜艳, 王彤, 刘洪波, 等. 医学统计学[G]. 5. 人民卫生出版社, 2020.第41页例3-6。原文采用配对T检验 |
||||||||||
# 输入数据 X1 <- c(0.840, 0.591, 0.674, 0.632, 0.687, 0.978, 0.750, 0.730, 1.200, 0.870) X2 <- c(0.580, 0.509, 0.500, 0.316, 0.337, 0.517, 0.454, 0.512, 0.997, 0.506) # 查看正态性 hist(X1-X2, xlab=NULL) # shapiro.test(X1-X2)$p.value ## [1] 0.9992947
t.test(X1, X2, paired=TRUE) ## ## Paired t-test ## ## data: X1 and X2 ## t = 7.926, df = 9, p-value = 0.00002384 ## alternative hypothesis: true mean difference is not equal to 0 ## 95 percent confidence interval: ## 0.1946542 0.3501458 ## sample estimates: ## mean difference ## 0.2724
得到p-value = 0.00002384 < 0.05, 因此，我们拒接了H0，认为两种方法的测定结果不同，X1组的均值要高于X2组的。
配对T检验与用差值进行单样本T检验的结果是一样的：
t.test(X1 - X2) ## ## One Sample t-test ## ## data: X1 - X2 ## t = 7.926, df = 9, p-value = 0.00002384 ## alternative hypothesis: true mean is not equal to 0 ## 95 percent confidence interval: ## 0.1946542 0.3501458 ## sample estimates: ## mean of x ## 0.2724
异方差双样本T检验（Welch’s t-test）
当总体方差不等时，常用的是Welch's T-test，使用Welch–Satterthwaite近似校正自由度。R里的t.test默认就是进行Welch's T-test，即默认情况下设置了var.equal=FALSE。
有文献指出，应该默认就直接使用Welch’s t-test5)，因为Welch’s t-test在方差不相等的情况下能够更好地控制I类错误，而在方差相等的情况下相较Student’s t-test并不会降低检验效能（或者仅有很小的降低）。
| 例：某研究调查了某社区II型糖尿病患者,各测量了25名血糖控制较差(甲组)和25名血糖控制良好(乙组)的血清总胆固醇含量(mmolL) ，如下表。问两组的均值是否不同？ | |||||||
|---|---|---|---|---|---|---|---|
| 甲组 | 5.61 | 4.18 | 3.21 | 2.41 | 2.12 | 6.81 | 3.47 |
| 4.46 | 1.98 | 3.03 | 3.96 | 4.91 | 4.42 | 6.92 | |
| 4.89 | 2.99 | 3.51 | 5.91 | 4.90 | 4.59 | 5.04 | |
| 4.16 | 4.91 | 4.85 | 4.70 | ||||
| 乙组 | 4.76 | 6.54 | 5.65 | 5.30 | 4.96 | 6.30 | 6.17 |
| 4.40 | 5.75 | 4.42 | 4.54 | 5.10 | 5.37 | 4.21 | |
| 5.44 | 5.81 | 4.65 | 6.05 | 4.82 | 5.53 | 5.00 | |
| 4.66 | 4.37 | 4.42 | 4.38 | ||||
案例及数据引自颜艳, 王彤, 刘洪波, 等. 医学统计学[G]. 5. 人民卫生出版社, 2020.第51页例3-8。 |
|||||||
# 首先构建两个向量录入数据 X1 <- c(5.61, 4.18, 3.21, 2.41, 2.12, 6.81, 3.47, 4.46, 1.98, 3.03, 3.96, 4.91, 4.42, 6.92, 4.89, 2.99, 3.51, 5.91, 4.90, 4.59, 5.04, 4.16, 4.91, 4.85, 4.70) X2 <- c(4.76, 6.54, 5.65, 5.30, 4.96, 6.30, 6.17, 4.40, 5.75, 4.42, 4.54, 5.10, 5.37, 4.21, 5.44, 5.81, 4.65, 6.05, 4.82, 5.53, 5.00, 4.66, 4.37, 4.42, 4.38) # 采用F检验检验方差齐性 var.test(X1, X2) ## ## F test to compare two variances ## ## data: X1 and X2 ## F = 3.5022, num df = 24, denom df = 24, p-value = 0.0032 ## alternative hypothesis: true ratio of variances is not equal to 1 ## 95 percent confidence interval: ## 1.543310 7.947462 ## sample estimates: ## ratio of variances ## 3.502199 # 上边F检验得到p < 0.05，可以认为这两组数据方差不同， # 因此采用异方差t检验 t.test(X1, X2, var.equal = FALSE) ## ## Welch Two Sample t-test ## ## data: X1 and X2 ## t = -2.8322, df = 36.672, p-value = 0.007465 ## alternative hypothesis: true difference in means is not equal to 0 ## 95 percent confidence interval: ## -1.4177999 -0.2350001 ## sample estimates: ## mean of x mean of y ## 4.3176 5.1440
这里如果强行采用var.equal = TRUE，得到的p值为0.006738，要更小一些。因此需要注意，如果p值刚好在0.05附近错误地采用等方差t检验，那可能引致不准确的结论。
威尔科克森秩和检验与符号秩检验（Wilcoxon test）
T检验是参数检验，其假设数据连续分布且符合正态分布。而当数据不适合使用t检验的时候，要比较两组样本均值，可以采用威尔科克森检验。威尔科克森检验是非参数检验，不需要假定数据符合任何特定的分布。
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, tol.root = 1e-4, digits.rank = Inf, ...)
其中y、alternative、mu、paired等参数与t.test中的一样。默认paired = FALSE，进行威尔科克森秩和检验（Wilcoxon秩和检验，Wilcoxon rank sum test），用于独立样本；设置paired = TRUE，进行Wilcoxon符号秩检验（Wilcoxon signed-rank test），用于配对样本或者单样本。
威尔科克森秩和检验（Wilcoxon rank sum test）
Wilcoxon rank sum test又被称为Mann-Whitney U test（曼-惠特尼 U 检验）或者Mann-Whitney-Wilcoxon test，或者Wilcoxon-Mann-Whitney test。用于样本非正态分布、样本量比较小(N < 30) 的情况，样本量比较大的时候就能直接用T检验了。
在R里边采用wilcox.test函数加上paired = TFLASE参数进行。
威尔科克森符号秩检验（Wilcoxon signed-rank test）
| 例：采用配对设计，用两种饲料饲喂8对大鼠，然后测定肝脏中维生素A的含量(IU/mg 如下表)，问不同饲料饲喂的大鼠肝脏中的维生素A的含量有无差异？ | ||||||||
|---|---|---|---|---|---|---|---|---|
| 编号 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 |
| 正常饲料组(G1) | 3.55 | 2.00 | 3.00 | 3.95 | 3.80 | 3.75 | 3.45 | 3.05 |
| 维生素E缺乏饲料组(G2) | 2.45 | 2.40 | 1.80 | 3.20 | 3.25 | 2.70 | 2.50 | 1.75 |
| 案例及数据引自颜艳, 王彤, 刘洪波, 等. 医学统计学[G]. 5. 人民卫生出版社, 2020.第126页例8-1，书上采用了Wilcoxon signed-rank test | ||||||||
# 输入数据 G1 <- c(3.55, 2.00, 3.00, 3.95, 3.80, 3.75, 3.45, 3.05) G2 <- c(2.45, 2.40, 1.80, 3.20, 3.25, 2.70, 2.50, 1.75) # 先看看分布情况 hist(G1, xlab=NULL) hist(G2, xlab=NULL)
通过histogram来看，两组数据看起来都偏离了正态分布，且偏向趋势不一致。然而，shapiro.test(G1)得到p-value = 0.1525，而shapiro.test(G2)得到p-value = 0.399，都未显示显著偏离正态分布，这或许正如一些人所说的正态性检验对小样本效力不够，偏向于得到更大的p值，对于大样本又过于有效，偏向于给出更小的p值。shapiro.test(G1-G2)得到p-value = 0.03705，他们之间的差偏离了正态分布。这里因为差值偏离了正态分布而选择了非参数检验。
wilcox.test(G1, G2, paired=TRUE) ## ## Wilcoxon signed rank exact test ## ## data: G1 and G2 ## V = 35, p-value = 0.01563 ## alternative hypothesis: true location shift is not equal to 0 wilcox.test(G1 - G2) ## ## Wilcoxon signed rank exact test ## ## data: G1 - G2 ## V = 35, p-value = 0.01563 ## alternative hypothesis: true location is not equal to 0
配对Wilcoxon检验和对两组的差值进行单样本Wilcoxon检验得到的结果是一样的。
使用t.test(G1, G2, paired=TRUE)检验得到p-value = 0.004001，t检验的p值要更小一些。一般来说，Wilcoxon检验得到的P值相较T检验偏大，Wilcoxon检验的检验效力要低于T检验的。如果数据本身适合与T检验，但是却采用了Wilcoxon，那么可能把一些本来显著的错误地认为是不显著的，也就是增加了II类误差。而反过来，如果不适合T检验却使用了T检验，那么p值偏小，就会错误地把不显著的当做显著的了。有的时候Wilcoxon检验得到的p值会更小，因为Wilcoxon检验针对的是中位数，而T检验针对的是均值，所以这两种检验的结果并不都是正相关的。
方差分析（ANOVA）
可以使用使用aov函数进行方差分析，TukeyHSD函数进行Post-hoc检验进行两两比较。分别用? aov和? TukeyHSD查看使用手册。
aov(formula, data = NULL, projections = FALSE, qr = TRUE, contrasts = NULL, ...)
这里的formula和lm函数的一致，格式为 response ~ factor，因变量response和自变量factor之间用波浪线~链接。对于自变量超过一个时，放在右边用+连接、:表示交叉，factor1*factor2则等同于factor1 + factor2 + factor1:factor2。factorA*factorB*factorC就等同于factorA + factorB + factorC + factorA:factorB + factorA:factorC + factorB:factorC + FactorA:FactorB:FactorC。factor1/factor2则等同于factor1 + factor1:factor2。详见aov和lm说明文档。
aov进行的是Type I ANOVA，公式中变量的顺序可能会影响结果，aov先确定第一个变量能够解释多少方差并进行显著性检验，然后才轮到下一个变量解释剩下的方差，以此类推。当自变量之间有相关性的时候，改变顺序会得到不同的结果。因此要做好实验的设计，或者要了解自变量之间的关系。
TukeyHSD(x, which, ordered = FALSE, conf.level = 0.95, ...)
x就是aov函数运行得到的对象，其余参数一般保持默认即可。
单因素方差分析（One Way ANOVA）
单因素方差分析是t检验的推广，同时用于多个组间的比较。其H0为：各个组来自均值相等的群体。并不能使用对任何可能的两两组合进行t检验来替代方差分析，因为检验越多，也可能仅仅因为随机因素而得到显著的p值。
| 例：为了研究3种降血脂药物的临床疗效,按统一纳人标准选择120例高血脂患者采用完全随机设计方法将患者等分为3组,随机安排服用的药物,进行双盲试验。以用药6周前后甘油三酯降低量为评价指标（如下表，单位：mmol/L）。问3个组患者的甘油三酯降低量总体均数有无差别? | ||||||||||
|---|---|---|---|---|---|---|---|---|---|---|
| 甲药组 | 1.77 | 1.56 | 2.68 | 2.34 | 4.32 | 2.95 | 3.36 | 3.11 | 1.81 | 2.42 |
| 2.97 | 2.65 | 2.17 | 2.93 | 2.86 | 2.72 | 2.63 | 2.22 | 2.9 | 1.98 | |
| 2.96 | 3.52 | 3.85 | 3.56 | 3.53 | 4.07 | 4.04 | 3.93 | 4.19 | 3.3 | |
| 2.31 | 2.8 | 2.98 | 2.27 | 2.52 | 3.72 | 2.56 | 3.57 | 4.02 | 2.36 | |
| 乙药组 | 2.68 | 3.21 | 2.48 | 2.28 | 2.39 | 2.28 | 2.28 | 3.23 | 2.32 | 2.86 |
| 4.02 | 2.66 | 2.58 | 3.64 | 2.61 | 2.65 | 2.32 | 2.68 | 1.65 | 2.66 | |
| 3.25 | 2.64 | 4.59 | 1.66 | 3.34 | 3.13 | 3.59 | 2.56 | 3.5 | 3.53 | |
| 2.68 | 3.04 | 4.29 | 2.66 | 2.41 | 3.7 | 2.42 | 1.81 | 2.97 | 3.48 | |
| 丙药组 | 2.59 | 2.96 | 3 | 2.98 | 2.33 | 3.55 | 3.93 | 3.3 | 3.16 | 1.37 |
| 2.52 | 0.94 | 2.97 | 3.37 | 2.16 | 1.69 | 1.74 | 2.11 | 2.81 | 1.98 | |
| 1.72 | 1.19 | 1.63 | 1.27 | 1.08 | 1.89 | 1.06 | 2.17 | 2.28 | 0.89 | |
| 3.71 | 2.47 | 3.19 | 1.41 | 1.88 | 1.92 | 2.51 | 1.02 | 2.1 | 1.31 | |
例题引自颜艳, 王彤, 刘洪波, 等. 医学统计学[G]. 5. 人民卫生出版社, 2020.第58页例4-2。 |
||||||||||
# 构建一个数据框来录入数据 X1 <- data.frame( 甘油三酯降低量=c( 1.77, 1.56, 2.68, 2.34, 4.32, 2.95, 3.36, 3.11, 1.81, 2.42, 2.97, 2.65, 2.17, 2.93, 2.86, 2.72, 2.63, 2.22, 2.90, 1.98, 2.96, 3.52, 3.85, 3.56, 3.53, 4.07, 4.04, 3.93, 4.19, 3.30, 2.31, 2.80, 2.98, 2.27, 2.52, 3.72, 2.56, 3.57, 4.02, 2.36), 分组 = "甲药组") X2 <- data.frame( 甘油三酯降低量=c( 2.68, 3.21, 2.48, 2.28, 2.39, 2.28, 2.28, 3.23, 2.32, 2.86, 4.02, 2.66, 2.58, 3.64, 2.61, 2.65, 2.32, 2.68, 1.65, 2.66, 3.25, 2.64, 4.59, 1.66, 3.34, 3.13, 3.59, 2.56, 3.50, 3.53, 2.68, 3.04, 4.29, 2.66, 2.41, 3.70, 2.42, 1.81, 2.97, 3.48), 分组 = "乙药组") X3 <- data.frame( 甘油三酯降低量=c( 2.59, 2.96, 3.00, 2.98, 2.33, 3.55, 3.93, 3.30, 3.16, 1.37, 2.52, 0.94, 2.97, 3.37, 2.16, 1.69, 1.74, 2.11, 2.81, 1.98, 1.72, 1.19, 1.63, 1.27, 1.08, 1.89, 1.06, 2.17, 2.28, 0.89, 3.71, 2.47, 3.19, 1.41, 1.88, 1.92, 2.51, 1.02, 2.10, 1.31), 分组 = "丙药组") X <- rbind(X1, X2, X3) head(X) ## 甘油三酯降低量 分组 ## 1 1.77 甲药组 ## 2 1.56 甲药组 ## 3 2.68 甲药组 ## 4 2.34 甲药组 ## 5 4.32 甲药组 ## 6 2.95 甲药组
# 方差齐性检验 car::leveneTest(X$甘油三酯降低量, factor(X$分组)) Levene's Test for Homogeneity of Variance (center = median) Df F value Pr(>F) group 2 1.8226 0.1662 117
Levene's test得到p-value = 0.1662 > 0.05，因此我们不认为这三组数据的方差之间存在差异，可以进行ANOVA。
# 进行模型拟合 fit <- aov(甘油三酯降低量 ~ 分组, data=X) # 查看概要信息 summary(fit) ## Df Sum Sq Mean Sq F value Pr(>F) ## 分组 2 13.62 6.811 12.27 0.0000145 *** ## Residuals 117 64.92 0.555 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 str(summary(fit)) ## List of 1 ## $ :Classes ‘anova’ and 'data.frame': 2 obs. of 5 variables: ## ..$ Df : num [1:2] 2 117 ## ..$ Sum Sq : num [1:2] 13.6 64.9 ## ..$ Mean Sq: num [1:2] 6.811 0.555 ## ..$ F value: num [1:2] 12.3 NA ## ..$ Pr(>F) : num [1:2] 0.0000145 NA ## - attr(*, "class")= chr [1:2] "summary.aov" "listof" summary(fit)[[1]][["Pr(>F)"]][1] ## [1] 0.00001448701
# 诊断图共四个，plot前先用par把作图界面设成两行两列 par(mfrow=c(2, 2)) plot(fit)
Residuals vs Fitted图用于查看方差齐性，可以看到Residuals（残差，即观察值和拟合值之差）与拟合的值之间没有显著的关联，我们可以认为三组之间符合方差齐性。这与上边使用Levene's test得到的p-value是相一致的。 如果发现图上的红线不是水平的，比较严重地偏离了水平，那么说明各组之间不是等方差的，这时候可能就需要考虑使用Welch's ANOVA了，Welch's ANOVA不假定各组之间是等方差的。同样也有人建议说应该始终使用Welch's ANOVA6)。
Scale-Location图 与Residuals vs Fitted图差不多，只是纵坐标换成了标准化的Residuals的绝对值的平方根，这样更便于观察误差分布范围。
Q-Q Residuals图用于查看正态性，这里可以看到，数据点都落在参考线上或者附近，说明可以认为符合正态假设。如果正态性不满足，有可能就要考虑使用非参数检验Kruskal-Wallis test。当然，不是特别严重的偏离正态分布对参数检验的结果的影响不大，而偏离多少才算严重确实是一个没有人给过量化指标的问题。
Constant Leverage: Residuals vs Factor Levels图以分组作为横坐标，也是一样显示方差齐性。
ANOVA告诉我们，至少有两组之间存在了均值的差异，那么具体是哪两组之间存在差异呢？就需要进行Post-hoc检验，可以采用TukeyHSD来完成：
TukeyHSD(fit) ## Tukey multiple comparisons of means ## 95% family-wise confidence level ## ## Fit: aov(formula = 甘油三酯降低量 ~ 分组, data = X) ## ## $分组 ## diff lwr upr p adj ## 甲药组-丙药组 0.75625 0.3608367 1.1516633 0.0000407 ## 乙药组-丙药组 0.66425 0.2688367 1.0596633 0.0003409 ## 乙药组-甲药组 -0.09200 -0.4874133 0.3034133 0.8454854 # 作图展示结果 plot(TukeyHSD(fit))
有可能ANOVA检验得到p < 0.05，但是最后两两比较的时候却没有得到任何两组之间的p < 0.05。这也是可以理解的，因为ANOVA检验的是组间的差异，而不是至少有两组存在差异。ANOVA检验发现了差异，但是Post-hoc并没有识别出哪两组之间存在差异，这也是合理的。
Welch's ANOVA
R里进行Welch's ANOVA使用oneway.test函数，Post-hoc检验使用Games-Howell test，可以使用rstatix包里的games_howell_test函数。 使用? oneway.test和? rstatix::games_howell_test查看函数的使用手册。
同样使用上边的数据，可以看到，结果会有些差异：
oneway.test(甘油三酯降低量 ~ 分组, data=X) ## ## One-way analysis of means (not assuming equal variances) ## ## data: 甘油三酯降低量 and 分组 ## F = 10.677, num df = 2.000, denom df = 77.323, p-value = 0.00008043
rstatix::games_howell_test(甘油三酯降低量 ~ 分组, data=X) ## # A tibble: 3 × 8 ## .y. group1 group2 estimate conf.low conf.high p.adj p.adj.signif ## * <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <chr> ## 1 甘油三酯降低量 丙药组 甲药组 0.756 0.338 1.17 0.000136 *** ## 2 甘油三酯降低量 丙药组 乙药组 0.664 0.260 1.07 0.000549 *** ## 3 甘油三酯降低量 甲药组 乙药组 -0.0920 -0.462 0.278 0.824 ns
Kruskal-Wallis Rank Sum Test
Kruskal-Wallis Rank Sum Test是ANOVA对应的非参数检验。R里边使用kruskal.test函数来完成（在stats包）。
? kruskal.test
双因素方差分析（Two Way ANOVA）
例：比如说，现在有一个项研究，研究四种饲料对某种动物体重的影响，每种饲料用于喂养64头猪，专门有16个饲养员来负责。假如采用嵌套设计，就是每个饲养员分到1种饲料，喂养16头猪，最后得到的数据录入如下：
dat <- data.frame( 体重 = c( 98, 97, 104, 102, 105, 103, 102, 105, 99, 102, 97, 101, 103, 101, 96, 101, 104, 95, 98, 102, 104, 104, 105, 104, 101, 95, 105, 104, 96, 102, 97, 103, 104, 104, 100, 101, 99, 97, 100, 97, 99, 95, 103, 99, 103, 105, 95, 105, 99, 98, 105, 98, 101, 100, 103, 97, 104, 101, 97, 102, 104, 101, 102, 97, 108, 108, 112, 115, 110, 115, 106, 106, 110, 113, 115, 109, 109, 105, 114, 106, 105, 105, 106, 111, 107, 112, 109, 114, 107, 114, 114, 112, 113, 114, 108, 113, 106, 106, 113, 107, 105, 105, 112, 105, 112, 115, 108, 105, 114, 114, 113, 108, 115, 115, 115, 108, 114, 105, 111, 111, 114, 107, 113, 110, 105, 113, 109, 105, 127, 134, 137, 134, 144, 143, 129, 134, 137, 132, 136, 142, 145, 130, 133, 126, 134, 135, 144, 137, 131, 134, 144, 128, 125, 142, 135, 128, 142, 135, 130, 125, 140, 130, 141, 128, 127, 135, 138, 127, 139, 138, 137, 128, 133, 145, 129, 144, 129, 141, 132, 137, 131, 139, 138, 143, 139, 127, 134, 127, 138, 130, 141, 143, 151, 154, 152, 152, 152, 155, 148, 149, 153, 151, 148, 153, 146, 145, 151, 154, 149, 151, 147, 146, 146, 152, 153, 155, 149, 152, 153, 154, 150, 153, 147, 147, 152, 153, 148, 146, 146, 145, 146, 146, 150, 151, 153, 155, 154, 151, 148, 155, 149, 145, 149, 151, 153, 155, 149, 148, 153, 147, 147, 148, 146, 152, 148, 149), 饲料 = rep(c("甲", "乙", "丙", "丁"), each=64), 饲者 = rep(LETTERS[1:16], each=16))
首先我们看一下这个数据的分布情况
library(ggplot2) ggplot(dat, aes(x=饲者, y=体重, fill=饲料)) + geom_boxplot() + theme_classic() + theme(legend.position="inside", legend.position.inside=c(0, 1), legend.justification.inside=c(-0.1, 1))
从这个图上明显可以看出，饲料影响比较大，而用同一种饲料的饲者养出来的动物体重都差不多。那么饲料和饲者的交叉是否影响呢。这时候做分析可以如下：
fit <- aov(体重~饲料 + 饲料:饲者, data=dat) summary(fit) ## Df Sum Sq Mean Sq F value Pr(>F) ## 饲料 3 98284 32761 1905.328 <0.0000000000000002 *** ## 饲料:饲者 12 57 5 0.274 0.993 ## Residuals 240 4127 17 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
可以看到，饲料与饲者的互作会有一些影响，但是并不显著。饲料的影响是显著的。
像这样的数据，就需要格外小心，不同formula的构建得到的结论会非常不一样：
fit01 <- aov(体重~饲料*饲者, data=dat) summary(fit01) ## Df Sum Sq Mean Sq F value Pr(>F) ## 饲料 3 98284 32761 1905.328 <0.0000000000000002 *** ## 饲者 12 57 5 0.274 0.993 ## Residuals 240 4127 17 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 fit02 <- aov(体重~饲者*饲料, data=dat) summary(fit02) ## Df Sum Sq Mean Sq F value Pr(>F) ## 饲者 15 98340 6556 381.3 <0.0000000000000002 *** ## Residuals 240 4127 17 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
对于fit01，在分析的时候，把饲料写在前边，就会先用饲料来解释方差，然后用饲者来解释方差，饲料和饲者就解释了100%的方差，因此看不到饲料:饲者互作一行。而如果把饲者写在前边，饲者就足够解释所有的方差了，所以可以看到饲料和饲料:饲者互作就没有了。足以见得，要非常小心，不合适的设计可能得到完全错误的结论。
fit03 <- aov(体重~饲料/饲者, data=dat) summary(fit03) ## Df Sum Sq Mean Sq F value Pr(>F) ## 饲料 3 98284 32761 1905.328 <0.0000000000000002 *** ## 饲料:饲者 12 57 5 0.274 0.993 ## Residuals 240 4127 17 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 fit04 <- aov(体重~饲者/饲料, data=dat) summary(fit04) ## Df Sum Sq Mean Sq F value Pr(>F) ## 饲者 15 98340 6556 381.3 <0.0000000000000002 *** ## Residuals 240 4127 17 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 fit05 <- aov(体重~饲者, data=dat) summary(fit05) ## Df Sum Sq Mean Sq F value Pr(>F) ## 饲者 15 98340 6556 381.3 <0.0000000000000002 *** ## Residuals 240 4127 17 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 fit06 <- aov(体重~饲者, data=dat) summary(fit06) ## Df Sum Sq Mean Sq F value Pr(>F) ## 饲者 15 98340 6556 381.3 <0.0000000000000002 *** ## Residuals 240 4127 17 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 fit07 <- aov(体重~饲者:饲料, data=dat) summary(fit07) ## Df Sum Sq Mean Sq F value Pr(>F) ## 饲者:饲料 15 98340 6556 381.3 <0.0000000000000002 *** ## Residuals 240 4127 17 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 fit08 <- aov(体重~饲料:饲者, data=dat) summary(fit08) ## Df Sum Sq Mean Sq F value Pr(>F) ## 饲料:饲者 15 98340 6556 381.3 <0.0000000000000002 *** ## Residuals 240 4127 17 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
例：比如说前边提到的饲养实验，研究四种饲料对某种动物体重的影响，每种饲料用于喂养64头猪，专门有16个饲养员来负责。假如采取平衡的设计，就是每个饲养员分到4种饲料，每种饲料，喂养4头猪，假如最终猪的体重还是那样的，但是猪和饲养员的对应关系不一样了，最后得到的数据录入如下：
Dat <- data.frame( 体重 = c( 98, 97, 104, 102, 105, 103, 102, 105, 99, 102, 97, 101, 103, 101, 96, 101, 104, 95, 98, 102, 104, 104, 105, 104, 101, 95, 105, 104, 96, 102, 97, 103, 104, 104, 100, 101, 99, 97, 100, 97, 99, 95, 103, 99, 103, 105, 95, 105, 99, 98, 105, 98, 101, 100, 103, 97, 104, 101, 97, 102, 104, 101, 102, 97, 108, 108, 112, 115, 110, 115, 106, 106, 110, 113, 115, 109, 109, 105, 114, 106, 105, 105, 106, 111, 107, 112, 109, 114, 107, 114, 114, 112, 113, 114, 108, 113, 106, 106, 113, 107, 105, 105, 112, 105, 112, 115, 108, 105, 114, 114, 113, 108, 115, 115, 115, 108, 114, 105, 111, 111, 114, 107, 113, 110, 105, 113, 109, 105, 127, 134, 137, 134, 144, 143, 129, 134, 137, 132, 136, 142, 145, 130, 133, 126, 134, 135, 144, 137, 131, 134, 144, 128, 125, 142, 135, 128, 142, 135, 130, 125, 140, 130, 141, 128, 127, 135, 138, 127, 139, 138, 137, 128, 133, 145, 129, 144, 129, 141, 132, 137, 131, 139, 138, 143, 139, 127, 134, 127, 138, 130, 141, 143, 151, 154, 152, 152, 152, 155, 148, 149, 153, 151, 148, 153, 146, 145, 151, 154, 149, 151, 147, 146, 146, 152, 153, 155, 149, 152, 153, 154, 150, 153, 147, 147, 152, 153, 148, 146, 146, 145, 146, 146, 150, 151, 153, 155, 154, 151, 148, 155, 149, 145, 149, 151, 153, 155, 149, 148, 153, 147, 147, 148, 146, 152, 148, 149), 饲料 = rep(c("甲", "乙", "丙", "丁"), each=64), 饲者 = rep(LETTERS[1:16], time=16))
首先我们看一下这个数据的分布情况
library(ggplot2) ggplot(Dat, aes(x=饲者, y=体重, color=饲料)) + geom_boxplot() + theme_classic() + theme(legend.position="none", legend.position.inside=c(0, 1), legend.justification.inside=c(0, 1))
Fit01 <- aov(体重~饲料*饲者, data=Dat) summary(Fit01) ## Df Sum Sq Mean Sq F value Pr(>F) ## 饲料 3 98284 32761 1775.501 <0.0000000000000002 *** ## 饲者 15 163 11 0.590 0.881 ## 饲料:饲者 45 477 11 0.575 0.985 ## Residuals 192 3543 18 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Fit02 <- aov(体重~饲者*饲料, data=Dat) summary(Fit02) ## Df Sum Sq Mean Sq F value Pr(>F) ## 饲者 15 163 11 0.590 0.881 ## 饲料 3 98284 32761 1775.501 <0.0000000000000002 *** ## 饲者:饲料 45 477 11 0.575 0.985 ## Residuals 192 3543 18 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
很显然，对于这样的一个设计，公式中自变量的顺序就不影响最终的结果了。这样的设计显然非常明确了饲养员对最终结果的影响不显著，饲养员与饲料的互作也没有显著的影响。
多因素方差分析（Multi-Factor ANOVA）
假如前文中的养猪实验再加上一个因素，比如是否添加某种益生菌，我们假如益生菌与几种饲料的组合能够促进生长，这样我们模拟一组示例数据，考虑到饲料会影响到体重，益生菌在对体重增加较小的两个组中会增加体重，如下：
# 随机产生一些体重数据 set.seed(1); W1 <- sample(95:105, 128, replace=TRUE) set.seed(2); W2 <- sample(105:115, 128, replace=TRUE) set.seed(3); W3 <- sample(125:145, 128, replace=TRUE) set.seed(4); W4 <- sample(145:155, 128, replace=TRUE) # 随机产生一些益生菌增加体重的值，比如对于W1对应增加多一些，对于W2增加少一些，对W3和W4无影响 set.seed(5); A1 <- sample(5:10, 64, replace=TRUE) set.seed(6); A2 <- sample(0:5, 64, replace=TRUE) DAT <- data.frame( 体重 = c(W1, W2, W3, W4), 饲料 = rep(c("甲", "乙", "丙", "丁"), each=128), 饲者 = rep(LETTERS[1:16], time=32), 加菌 = rep(c(rep("是", 16), rep("否", 16)), time=8)) DAT[DAT$饲料 == "甲" & DAT$加菌 == "是", ]$体重 <- A1 + DAT[DAT$饲料 == "甲" & DAT$加菌 == "是", ]$体重 DAT[DAT$饲料 == "乙" & DAT$加菌 == "是", ]$体重 <- A2 + DAT[DAT$饲料 == "乙" & DAT$加菌 == "是", ]$体重
完全可以忽略掉上边这个产生示例数据的过程，从boxplot看出数据的情况：
library(ggplot2) ggplot(DAT, aes(x=饲者, y=体重, color=饲料)) + geom_boxplot() + facet_wrap(~加菌) + theme_classic() + theme(legend.position="none", legend.position.inside=c(0, 1), legend.justification.inside=c(0, 1))
FIT01 <- aov(体重~饲料*饲者*加菌, data=DAT) summary(FIT01) ## Df Sum Sq Mean Sq F value Pr(>F) ## 饲料 3 172855 57618 3753.404 < 0.0000000000000002 *** ## 饲者 15 357 24 1.548 0.0857 . ## 加菌 1 758 758 49.382 0.00000000000966 *** ## 饲料:饲者 45 829 18 1.200 0.1855 ## 饲料:加菌 3 1536 512 33.353 < 0.0000000000000002 *** ## 饲者:加菌 15 186 12 0.807 0.6696 ## 饲料:饲者:加菌 45 605 13 0.876 0.6997 ## Residuals 384 5895 15 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 FIT02 <- aov(体重~饲料*加菌*饲者, data=DAT) summary(FIT02) ## Df Sum Sq Mean Sq F value Pr(>F) ## 饲料 3 172855 57618 3753.404 < 0.0000000000000002 *** ## 加菌 1 758 758 49.382 0.00000000000966 *** ## 饲者 15 357 24 1.548 0.0857 . ## 饲料:加菌 3 1536 512 33.353 < 0.0000000000000002 *** ## 饲料:饲者 45 829 18 1.200 0.1855 ## 加菌:饲者 15 186 12 0.807 0.6696 ## 饲料:加菌:饲者 45 605 13 0.876 0.6997 ## Residuals 384 5895 15 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 FIT03 <- aov(体重~加菌*饲料*饲者, data=DAT) summary(FIT03) ## Df Sum Sq Mean Sq F value Pr(>F) ## 加菌 1 758 758 49.382 0.00000000000966 *** ## 饲料 3 172855 57618 3753.404 < 0.0000000000000002 *** ## 饲者 15 357 24 1.548 0.0857 . ## 加菌:饲料 3 1536 512 33.353 < 0.0000000000000002 *** ## 加菌:饲者 15 186 12 0.807 0.6696 ## 饲料:饲者 45 829 18 1.200 0.1855 ## 加菌:饲料:饲者 45 605 13 0.876 0.6997 ## Residuals 384 5895 15 ## --- ## Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
可以看到，正如我们所设想的那样，饲料影响显著，加菌影响显著，加菌和饲料互作有影响。可以看出来，因子的顺序并不影响最终的结果。
psh <- TukeyHSD(FIT01) head(psh[[1]]) ## diff lwr upr p adj ## 丁-丙 15.82031 14.556595 17.084030 0 ## 甲-丙 -30.24219 -31.505905 -28.978470 0 ## 乙-丙 -23.24219 -24.505905 -21.978470 0 ## 甲-丁 -46.06250 -47.326217 -44.798783 0 ## 乙-丁 -39.06250 -40.326217 -37.798783 0 ## 乙-甲 7.00000 5.736283 8.263717 0 head(psh[[7]]) ## diff lwr upr p adj ## 丁:A:否-丙:A:否 15.75 3.412191 28.087809 0.0001889579 ## 甲:A:否-丙:A:否 -37.75 -50.087809 -25.412191 0.0000000000 ## 乙:A:否-丙:A:否 -25.25 -37.587809 -12.912191 0.0000000000 ## 丙:B:否-丙:A:否 -9.50 -21.837809 2.837809 0.6967895859 ## 丁:B:否-丙:A:否 8.25 -4.087809 20.587809 0.9615631558 ## 甲:B:否-丙:A:否 -36.00 -48.337809 -23.662191 0.0000000000
卡方检验（χ2检验, Chi-square test）
四格表资料的卡方检验（Pearson's Chi-squared test）
R 里的卡方检验使用chisq.test函数来进行? 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是一个至少有2行和2列的矩阵时，进行皮尔森卡方检验；或者当x和y都是因子时，从x和y计算出一个列联表，再进行皮尔森卡方检验；当x是向量或者只有一行或一列的矩阵，而y又没有给定时，进行拟合优度检验（详见后文）。chisq.test默认对卡方统计量进行Yates连续性校正，即默认correct = TRUE。默认simulate.p.value = FALSE，设置simulate.p.value = TRUE进行蒙特卡洛模拟，通常当运行出现警告信息“ Chi-squared近似算法有可能不准”时加上这一设置。
卡方检验的适用适用条件包括：（1）样本数据是随机抽取的；（2）一个观测值只属于一个类别；（3）观测值相互独立，即一个观测值不会影响其它观测值的结果；（4）样本量大（通常认为n>40），至少80%的理论频数要大于5，对于四格表（2×2），要求所有单元格理论频数都大于5。
在实际工作中，对于四格表资料，通常规定：
（1）当四格总数n≥40，且所有理论频数T≥5时，设置correct = FALSE，计算皮尔逊卡方统计量，但当得到的P≈α (检验水准，通常为0.05)时，改用Fisher精确检验；
（2）当n≥40，所有T≥1，存在T<5时，设置correct = TRUE，或者改用Fisher精确检验；
（3）当n<40，或存在T<1时，改用Fisher精确检验。也有人说，设置simulate.p.value = TRUE即可，但蒙特卡洛模拟是个随机过程，每次运行的结果都会存在一定的差异，既然适用fisher.test也使用很方便，直接使用fisher.test即可。
| 例：某研究者欲比较膳食干预(试验组)和普通健康教育(对照组)对糖尿病患者血糖的控制达标情况,将171例糖尿病患者随机分为两组,干预6个月后结果见下表。问两组控制空腹血糖的总体达标率有无差别? | ||||||
| 组别 | 达标 | 未达标 | ||||
|---|---|---|---|---|---|---|
| 对照组 | 30 | 56 | ||||
| 实验组 | 43 | 42 | ||||
例题引自颜艳, 王彤, 刘洪波, 等. 医学统计学[G]. 5. 人民卫生出版社, 2020.第107页例7-1。 |
||||||
dat01 <- " 组别 达标 未达标 对照组 30 56 实验组 43 42" dat01 <- read.table(textConnection(dat01), row.names=1, header=TRUE) dat01 <- as.matrix(dat01) dat01 ## 达标 未达标 ## 对照组 30 56 ## 实验组 43 42 chisq_result <- chisq.test(dat01, correct = FALSE) chisq_result ## ## Pearson's Chi-squared test ## ## data: dat01 ## X-squared = 4.3094, df = 1, p-value = 0.0379 # 看看理论频数 chisq_result$expected ## 达标 未达标 ## 对照组 36.71345 49.28655 ## 实验组 36.28655 48.71345
《医学统计学》教材上采用的检验是未经过矫正的，即这里设置了correct = FALSE。这里四格总数大于40，且所有格子中的理论频数都大于5，可以不使用校正。
| 例：某医师欲比较胞磷胆碱与神经节苷脂治疗脑血管疾病的疗效,将 78 例脑血管疾病患者随机分为两组,结果如下表。问两种药物治疗脑血管疾病的有效率是否相? | ||||||
| 组别 | 有效 | 无效 | ||||
|---|---|---|---|---|---|---|
| 胞磷胆碱组 | 46 | 6 | ||||
| 神经节苷脂组 | 18 | 8 | ||||
例题引自颜艳, 王彤, 刘洪波, 等. 医学统计学[G]. 5. 人民卫生出版社, 2020.第109页例7-2。 |
||||||
dat02 <- " 组别 有效 无效 胞磷胆碱组 46 6 神经节苷脂组 18 8" dat02 <- read.table(textConnection(dat02), row.names=1, header=TRUE) dat02 <- as.matrix(dat02) dat02 ## 有效 无效 ## 胞磷胆碱组 46 6 ## 神经节苷脂组 18 8 chisq.test(dat02, correct = TRUE) ## ## Pearson's Chi-squared test with Yates' continuity correction ## ## data: dat02 ## X-squared = 3.1448, df = 1, p-value = 0.07617 ## ## 警告信息: ## In chisq.test(dat02, correct = TRUE) : Chi-squared近似算法有可能不准 fisher.test(dat02) ## ## Fisher's Exact Test for Count Data ## ## data: dat02 ## p-value = 0.05844 ## alternative hypothesis: true odds ratio is not equal to 1 ## 95 percent confidence interval: ## 0.879042 13.548216 ## sample estimates: ## odds ratio ## 3.347519 P <- c(); for(i in 1:20) P <- c(P, chisq.test(dat02, simulate.p.value = TRUE)$p.value) P ## [1] 0.06046977 0.06596702 0.06796602 0.05547226 0.06896552 0.06346827 0.05347326 0.05747126 0.05647176 0.05747126 ## [11] 0.05597201 0.06646677 0.06296852 0.06796602 0.05997001 0.05497251 0.06496752 0.06196902 0.04797601 0.05547226 median(P) ## [1] 0.06021989
《医学统计学》中使用了校正的卡方检验。这里出现了警告信息，或许使用fisher.test会更好。可以看到fisher.test得到的p值有些差异，但也比较接近。使用蒙特卡洛模拟得到的p值也接近，但是每次运行得到的结果都不一样，还有可能出现p<0.05这样导致完全相反结论的可能，多进行几次取其中位数就比较合适，但这样显然比较麻烦，而且要实现可重复性，每次运行前都要set.seed，要设置不同的seed才能得到不同的p值。对于这样的一个数据，采用fisher.test可能是最合适的。
但是如果不进行校正，chisq.test(dat02, correct = FALSE)得到p-value = 0.03695，得到的结论可能就相反了，结果就非常不可靠了，因此这样来做不可取，这点大概是没有争议的。
配对四格表资料的卡方检验（McNemar's Chi-squared test）
R里用stats::mcnemar.test函数来做配对配对四格表资料的卡方检验，此检验用于配对四格表资料。比如同样或配对的材料，用两种方法进行了检测，每个个体或对都有两个结果，同时是阳性的频数为ａ，同时是阴性的频数为ｄ，结果不一致的频数分别为ｂ和ｃ，例如下表:
方法二 阳性 阴性 方法一 阳性 a b 阴性 c d
mcnemar.test(x, y = NULL, correct = TRUE)
x可以是二维矩阵或因子。当x是二维矩阵时忽略掉y，当x是因子时，y也要是因子。默认进行连续矫正（correct = TRUE）。当期望频数(b+c)/2 > 5时，检验才有效7)。当b+c<40时，需要进行连续矫正8)。
| 例：某实验室分别用乳胶凝集法和免疫荧光法对58例可疑系统性红斑狼疮患者血清中抗核抗体进行测定,结果如下表。问两种方法的检测结果有无差别? | ||||
| 乳胶凝集法 | 例题引自颜艳, 王彤, 刘洪波, 等. 医学统计学[G]. 5. 人民卫生出版社, 2020.第110页例7-3。 |
|||
|---|---|---|---|---|
| + | - | |||
| 免疫荧光法 | + | 11(a) | 12(b) | |
| - | 2 (c) | 33(d) | ||
# 首先把数据放到一个矩阵中 dat <- matrix(c(11, 2, 12, 33), ncol=2, dimnames=list(免疫荧光法=c("+", "-"), 乳胶凝集法=c("+", "-"))) dat ## 乳胶凝集法 ## 免疫荧光法 + - ## + 11 12 ## - 2 33 # 这里期望的频数为(12+2)/2 = 7 > 5, 使用McNemar检验时有效的。 # McNemar检验默认correct = TRUE，这里与《医学统计学》教材上保持一致一致，因b+c<40而保持默认设置 mcnemar.test(dat, correct = TRUE) ## ## McNemar's Chi-squared test with continuity correction ## ## data: dat ## McNemar's chi-squared = 5.7857, df = 1, p-value = 0.01616
行×列表资料的卡方检验
这种情况也是用chisq.test函数。
| 例：测得某地5801人的ABO血型和MN血型结果下表,问两种血型系统之间是否有关联? | |||||
| MN血型 | 例题引自颜艳, 王彤, 刘洪波, 等. 医学统计学[G]. 5. 人民卫生出版社, 2020.第115页例7-8。 |
||||
|---|---|---|---|---|---|
| M | N | MN | |||
| ABO血型 | O | 431 | 490 | 902 | |
| A | 388 | 410 | 800 | ||
| B | 495 | 587 | 950 | ||
| AB | 137 | 179 | 32 | ||
dat <- " M N MN O 431 490 902 A 388 410 800 B 495 587 950 AB 137 179 32" dat <- read.table(textConnection(dat), row.names=1, header=TRUE) dat <- as.matrix(dat) chisq.test(dat) ## ## Pearson's Chi-squared test ## ## data: dat ## X-squared = 213.16, df = 6, p-value < 0.00000000000000022
显然，这两种血型间存在某种关联。
还有其它更多地使用场景，这里就不一一列举了，知道了函数怎么用，遇到具体问题地时候再根据自己地实际情况选择合适地函数使用适当地参数进行检验，有必要的时候现学也可以的，不要期望把所有的东西都学会了再开始干活，很多时候是再使用的过程中去学习。
费希尔精确检验 (Fisher's exact test)
前边再皮尔逊卡方检验时提到，在一些场合，使用费希尔精确检验更为合适。R里用fisher.test函数进行费希尔精确检验。一般说fisher.test适用于样本量比较小的情况，因为样本量比较大的时候，fisher.test计算起来会比较困难。但现在使用软件计算比较方便，对于四格表来说，非常大的样本用软件计算起来也很快，fisher.test也能用，fisher.test是一种精确检验，对于四格表资料来说，fisher.test优于chisq.test 9)。 但当表格行列超过2×2时，fisher.test的消耗计算资源就会变得不易接受，这时候chisq.test就是更好的选择。
需要注意，fisher.test并不适用于配对四格表资料，不能代替mcnemar.test。
fisher.test(x, y = NULL, workspace = 200000, hybrid = FALSE, hybridPars = c(expect = 5, percent = 80, Emin = 1), control = list(), or = 1, alternative = "two.sided", conf.int = TRUE, conf.level = 0.95, simulate.p.value = FALSE, B = 2000)
一般我们常用于四格表，只要构建四格表矩阵作为参数x即可。与chisq.test不同的时，fisher.test可以进行单尾检验，而chisq.test只能进行双尾检验。fisher.test默认进行双尾检验，即alternative = “two.sided”。
| 例（女士品茶实验）：据说在一次几位同事的座谈会上，Ronald Aylmer Fisher给藻类学家Muriel Bristol女士倒了杯茶，Bristol说自己更喜欢先加奶后加茶，大家并不相信她能够品出先奶后茶和先茶后奶的区别。一个同事William Roach就建议做个实验，于是Fisher快速地准备了8杯茶，4杯先奶后茶，4杯先茶后奶，仅加入奶和茶的顺序不一样，其余的都一样。然后把这八杯随机排序，让Bristol品尝。问Bristol真的能尝出先奶后茶还是先茶后奶么? | ||||
| 判断结果 | 此故事流传于网络，不同处说法又略有不同，这里综合了一下各方信息，不保证完全符合事实，但这并不影响它作为一个例子。 | |||
|---|---|---|---|---|
| 先奶后茶 | 先茶后奶 | |||
| 实际情况 | 先奶后茶 | 3 | 1 | |
| 先茶后奶 | 1 | 3 | ||
这里我们假设Bristol对各杯的品尝的结果之间不会相互影响。这里样本量比较小，构建四格列联表, 然后进行Fisher精确检验：
dat01 <- " 判断先奶后茶 判断先茶后奶 实际先奶后茶 3 1 实际先茶后奶 1 3" dat01 <- read.table(textConnection(dat01), row.names=1, header=TRUE) dat01 <- as.matrix(dat01) fisher.test(dat01, alternative="greater") ## ## Fisher's Exact Test for Count Data ## ## data: dat01 ## p-value = 0.2429 ## alternative hypothesis: true odds ratio is greater than 1 ## 95 percent confidence interval: ## 0.3135693 Inf ## sample estimates: ## odds ratio ## 6.408309
p-value = 0.2429 > 0.05，没有足够的证据表明Bristol女士能够尝出先茶后奶和先奶后茶的区别。这里由于Bristol女士有预先的宣称，采用alternative=“greater”, 我们只对她是否表现得比完全靠猜更好，也就是H0假设是女士并不能区分（完全靠猜或者表现比猜更差）。 如果要看这个结果是不是完全随机产生的，也就是Bristol女士完全可能是靠猜得到的，那么fisher.test(dat, alternative=“two.sided”)得到p-value = 0.4857，也就是无法拒接她可以完全靠猜得到这样的结果的假设。而如果想要看看Bristol女士是否可能把两种茶完全弄反而得到这样的结果，那么fisher.test(dat, alternative=“less”)得到p-value = 0.9857，也就这个结果也无证据表明Bristol女士会倾向于把两种茶弄混。
例：还有一种说法是当时Bristol完全答对了，这样检验如下：
dat02 <- " 判断先奶后茶 判断先茶后奶 实际先奶后茶 4 0 实际先茶后奶 0 4" dat02 <- read.table(textConnection(dat02), row.names=1, header=TRUE) dat02 <- as.matrix(dat02) fisher.test(dat02, alternative="greater") ## ## Fisher's Exact Test for Count Data ## ## data: dat02 ## p-value = 0.01429 ## alternative hypothesis: true odds ratio is greater than 1 ## 95 percent confidence interval: ## 2.003768 Inf ## sample estimates: ## odds ratio ## Inf
如果是这种情况，p-value = 0.01429，也就是说有足够的证据表明了Bristol女士能够尝出哪种茶是哪种茶。在这样的一个小样本量的情况下，只有完全判断正确才能有足够的证据支持Bristol能够尝出先茶后奶还是先奶后茶。
假如增加样本量10倍，而判断正确的比例依然还是各3/4，那么结果会怎样呢？
dat03 <- " 判断先奶后茶 判断先茶后奶 实际先奶后茶 30 10 实际先茶后奶 10 30" dat03 <- read.table(textConnection(dat03), row.names=1, header=TRUE) dat03 <- as.matrix(dat03) fisher.test(dat03, alternative="greater") ## ## Fisher's Exact Test for Count Data ## ## data: dat03 ## p-value = 0.000007437 ## alternative hypothesis: true odds ratio is greater than 1 ## 95 percent confidence interval: ## 3.423196 Inf ## sample estimates: ## odds ratio ## 8.703721
可以看出，样本量增加了之后，同样的正确率却有足够的证据表明她有区分的能力。
如果她刚好把所有的都搞反了呢？
dat04 <- " 判断先奶后茶 判断先茶后奶 实际先奶后茶 0 4 实际先茶后奶 4 0" dat04 <- read.table(textConnection(dat04), row.names=1, header=TRUE) dat04 <- as.matrix(dat04) fisher.test(dat04, alternative="greater")$p.value ## [1] 1 fisher.test(dat04, alternative="less")$p.value ## [1] 0.01428571 fisher.test(dat04, alternative="two.sided")$p.value ## [1] 0.02857143
比例比较（Comparation of Proportions）
R里用于比例比较的函数有prop.test、binom.test、chisq.test、……
prop.test用于检验一个比例是否等于特定值，或者两组或多组比例之间是否相等，其参数及默认值情况如下：
prop.test(x, n, p = NULL, alternative = c("two.sided", "less", "greater"), conf.level = 0.95, correct = TRUE)
x是向量时，必须指定n参数，n为测试的总次数，x为成功的次数，p为成功的概率；当x为矩阵或者表格时候，忽略n，x第一列为成功次数，第二列为失败次数。
binom.test(x, n, p = 0.5, alternative = c("two.sided", "less", "greater"), conf.level = 0.95)
x为向量，当x长度为1时，x为成功的次数，n为测试的总次数，p为预计的成功概率；当x长度为2时，忽略n，x第一个元素为成功次数，第二个元素为失败次数。
在进行单比例检验的时候，当样本量比较小时推荐使用binom.test，而当样本量比较大（n>30）时，可以用prop.test。prop.test进行单比例Z检验的时候要求成功频数及失败频数都要≥5。
二项检验（binom.test）
binom.test用于进行二项检验（Binomial Test），Binomial Test是一种精确检验，又叫Exact Binomial Test。
当一个实验有两种可能的结果（即，成功或失败），且预期有一个成功的概率，则可用二项检验来判断结果是否与预期有异。
例：商家宣传某种药物对某个病症的有效率超过60%，那你在人群中随机调查了40个案例，发现其中20个案例中这种药有效，你能接受商家的说法么？
binom.test(20, 40, p = 0.6, alternative = "greater") ## ## Exact binomial test ## ## data: 20 and 40 ## number of successes = 20, number of trials = 40, p-value = 0.9256 ## alternative hypothesis: true probability of success is greater than 0.6 ## 95 percent confidence interval: ## 0.3610917 1.0000000 ## sample estimates: ## probability of success ## 0.5
也就是说，这个调查结果并不支持商家宣称的超过60%的有效率。
那你能说这种药的有效率低于60%么？
binom.test(20, 40, p = 0.6, alternative = "less") ## ## Exact binomial test ## ## data: 20 and 40 ## number of successes = 20, number of trials = 40, p-value = 0.1298 ## alternative hypothesis: true probability of success is less than 0.6 ## 95 percent confidence interval: ## 0.0000000 0.6389083 ## sample estimates: ## probability of success ## 0.5
使用alternative=“less”，p-value = 0.1298 > 0.05, 也没有足够的证据表明说这种药的有效率低于60%。
binom.test(20, 40, p = 0.6, alternative = "two.sided") ## ## Exact binomial test ## ## data: 20 and 40 ## number of successes = 20, number of trials = 40, p-value = 0.2007 ## alternative hypothesis: true probability of success is not equal to 0.6 ## 95 percent confidence interval: ## 0.3380178 0.6619822 ## sample estimates: ## probability of success ## 0.5
这样的调查结果无法说明有效率与60%存在差异。
当然，如果你得到这个调查结果之后，发现了这种低于60%的可能性，于是你继续进行调查，然后发现80个人中只有40个是有效的，那么这时候binom.test(40, 80, p = 0.6, alternative = “less”)得到p-value = 0.0445 < 0.05，那么可以说有足够的证据表明这种药的有效率是低于60%的。
单比例Z检验（prop.test）
例：王老师R语言课程2024年秋选课学生中男生为45人，女生为33人，问选课学生中男生人数与女生人数是否有显著差异？ 也就是问男生所占的比例是否显著偏离0.5。这里用prop.test检验一个比例是否等于特定值。
prop.test(45, 45+33, 0.5) ## ## 1-sample proportions test with continuity correction ## ## data: 45 out of 45 + 33, null probability 0.5 ## X-squared = 1.5513, df = 1, p-value = 0.2129 ## alternative hypothesis: true p is not equal to 0.5 ## 95 percent confidence interval: ## 0.459909 0.686304 ## sample estimates: ## p ## 0.5769231
很显然，p > 0.05，不能推翻假设，也就是说并没有足够的证据表明选这门课的男生就多于女生，也就是说，下一年可能就是女生人数比男生多，这个是随机可能出现的。于是王老师查看了一下2023年秋的情况，男生37人，女生44人，prop.test(37, 37+44, 0.5)得到p-value = 0.505，也并没有证据表明男生比女生要少。也就是说，不存在某个因素（或者多个因素总体的效果）让王老师这门R语言课程男生比女生多或者少，看到的人数差异都在随机可能出现的范围。由于并不知道整体的研究生性别比例，所以就无法知道是不是男生或者女生更倾向于选这门课。
binom.test(45, 45+33, 0.5)得到p-value = 0.2127，与prop.test(45, 45+33, 0.5)得到的结果是一致的。
双比例Z检验（prop.test）
| **例：有研究对山西太原市396例肺癌新发病例及465名健康对照者进行了调查，其中部分结果如下表。 ** | ||||||
| 健康者 | 肺癌患者 | |||||
|---|---|---|---|---|---|---|
| 腺癌 | 鳞癌 | 小细胞癌 | 其它(大细胞癌、未分化和未定型肺癌) | 合计 | ||
| 吸烟 | 181 | 47 | 69 | 31 | 72 | 219 |
| 不吸烟 | 284 | 59 | 23 | 26 | 69 | 177 |
| 数据引自：穆丽娜, 苏佳, 俞顺章, 等. 控烟是预防肺癌的主要措施——记太原市肺癌病例对照流行病学调查[J]. 中国癌症杂志, 2021,31(04):335-343. | ||||||
看看吸烟者中患肺癌者的比例是不是比不吸烟者中患肺癌的比例要高，这里用prop.test检验两组比例是否相等，这里设置alternative=“greater”，因为这里有明确的方向性假设，也就是吸烟增加患癌风险：
prop.test(x = c(219, 177), n = c(219+181, 177+284), alternative = "greater") ## ## 2-sample test for equality of proportions with continuity correction ## ## data: c(219, 177) out of c(219 + 181, 177 + 284) ## X-squared = 22.41, df = 1, p-value = 0.000001101 ## alternative hypothesis: greater ## 95 percent confidence interval: ## 0.1058652 1.0000000 ## sample estimates: ## prop 1 prop 2 ## 0.5475000 0.3839479
结果显示差异是显著的，得到的p值和fisher.test的很接近。
换一种角度来问，患肺癌者中吸烟者的比例是不是比健康者中吸烟的比例要高，得到的p值是一样的：0.000001101。
prop.test(x = c(219, 181), n = c(219+177, 181+284), alternative = "greater") ## ## 2-sample test for equality of proportions with continuity correction ## ## data: c(219, 181) out of c(219 + 177, 181 + 284) ## X-squared = 22.41, df = 1, p-value = 0.000001101 ## alternative hypothesis: greater ## 95 percent confidence interval: ## 0.106019 1.000000 ## sample estimates: ## prop 1 prop 2 ## 0.5530303 0.3892473
这里得到的结果与卡方检验是一致的：
M <- " 患肺癌 没患肺癌 吸烟 219 181 不吸烟 177 284" M <- read.table(textConnection(M), row.names=1, header=TRUE) chisq.test(M) ## ## Pearson's Chi-squared test with Yates' continuity correction ## data: M ## X-squared = 22.41, df = 1, p-value = 0.000002202
p是单尾prop.test的两倍，相当于prop.test采用双尾检验。chiseq.test没有alternative 参数。
例：现有研究在广西16111例不孕不育患者中检出α地贫2322例10)，而先前有研究在广西13589例一般(未明确是否不孕不育)样本中检出α地贫1757例11)。现在问，不孕不育患者中是否有更高比例的患者具有α地贫？
这里实际上就是要检测这两个比例是否有显著性差异，这时候就可以使用比例差异检验。这里重在这个检验方法的介绍，忽略掉数据间的可比较性的问题，假设这几年期间发病率无变化，假设对α地贫的判断标准是一致的（实际上研究之间的方法的问题可能导致这里对两组数据的比较仅仅具有一定参考意义而不能作为结论性的，这个需要明确）。
prop.test(c(2322, 1757), c(16111, 13589)) ## ## 2-sample test for equality of proportions with continuity correction ## ## data: c(2322, 1757) out of c(16111, 13589) ## X-squared = 13.558, df = 1, p-value = 0.0002314 ## alternative hypothesis: two.sided ## 95 percent confidence interval: ## 0.006936168 0.022722588 ## sample estimates: ## prop 1 prop 2 ## 0.1441251 0.1292958
可以看到p < 0.0005，这两个比例是存在显著差异的，也就是说，α地贫对不孕不育的影响是显著的，α地贫患者比一般人更有可能不孕不育。
prop.test默认进行了连续矫正（continuity correction），如果不矫正，得到的P值与Z检验是一致的。使用pnorm函数实现z test, p值计算如下：
R1 <- 2322; R2 <- 1757; N1 <- 16111; N2 <- 13589 P1 <- R1/N1 ; P2 <- R2/N2 P <- (R1+R2)/(N1+N2) Z <- (P1-P2)/sqrt(P*(1-P)*(1/N1 + 1/N2)) p <- 2*(1-pnorm(Z))
得到的p值为0.0002164729，与prop.test(c(2322, 1757), c(16111, 13589), correct=FALSE)得到的结果是一致的。
多组比例比较（prop.test）
prop.test还能用于多组比例的比较，也就是各组的比例是否相等。
M <- " 吸烟 不吸烟 健康 181 284 腺癌 47 59 鳞癌 69 23 小细胞癌 31 26 其它 72 69 " M <- read.table(textConnection(Input), row.names=1, header=TRUE) M <- as.matrix(M)
prop.test(M) ## ## 5-sample test for equality of proportions without continuity correction ## ## data: M ## X-squared = 43.573, df = 4, p-value = 0.000000007869 ## alternative hypothesis: two.sided ## sample estimates: ## prop 1 prop 2 prop 3 prop 4 prop 5 ## 0.3892473 0.4433962 0.7500000 0.5438596 0.5106383
这里得到p < 0.05，也就是说存在一组或者某一组与其它的组比例存在差异，这时候可以用两两比较来找出与别的组存在差异的组，使用pairwise.prop.test函数。
两两比例比较（pairwise.prop.test）
pairwise.prop.test(x, n, p.adjust.method = p.adjust.methods, ...)
其中x和n的用法与
pairwise.prop.test(M) ## ## Pairwise comparisons using Pairwise comparison of proportions ## ## data: M ## ## 健康 腺癌 鳞癌 小细胞癌 ## 腺癌 1.00000 - - - ## 鳞癌 0.0000000043 0.00022 - - ## 小细胞癌 0.17844 1.00000 0.09646 - ## 其它 0.09646 1.00000 0.00350 1.00000 ## ## P value adjustment method: holm
可以看到，吸烟者的比例只有腺癌和健康者之间存在显著的差异。这或许意味着背后有一些机制。
拟合优度卡方检验（Chi-square Goodness of Fit Test）
拟合优度卡方检验用于比较多个观察比例与期望概率，这时候使用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是一个向量，而y没有给定时，进行goodness-of-fit test，检验比例是否与p指定的比例相等，如果p没有给定，则检验比例是否全部都相同。
例：孟德尔豌豆实验中，15株植物的556颗种子的形状和颜色分布为：圆且黄者315粒，皱且黄者101粒，圆且绿者108粒，皱且绿者32粒。按照预期，黄相对绿为显性，预计得到比例为3:1；圆相对皱为显性，预计得到比例也为3:1。两个性状自由组合，那么这四种种子的比例应该为9:3:3:1。问：这个数据符合这个比例么？也就是四种种子的比例是否依次为9/16, 3/16, 3/16, 1/16
chisq.test(x=c(315, 101, 108, 32), p = c(9/16, 3/16, 3/16, 1/16)) ## ## Chi-squared test for given probabilities ## ## data: c(315, 101, 108, 32) ## X-squared = 0.47002, df = 3, p-value = 0.9254
显然，没有足够的证据拒绝这个数据符合这个比例的假设。
这里同样要注意卡方检验的基本前提。
相关性分析
? cor
回归分析
逻辑回归分析
生存分析
install.packages("survival") vignette("survival")
其它
本网页已被访问15508次。