Stata软件基本操作和数据分析入门第五讲 多组平均水平的比较赵耐青
复习和补充两组比较的统计检验配对设计资料(又称为Dependent Samples)
对于小样本的情况下，如果配对的差值资料服从正态分布，用配对t检验(ttest 差值变量=0)
大样本的情况下，可以用配对t检验小样本的情况下，并且配对差值呈偏态分布，则用配对符号秩检验(signrank 差值变量=0)
成组设计(Two Independent Samples)
如果方差齐性并且大样本情况下，可以用成组t检验(ttest 效应指标变量,by(分组变量))
如果方差齐性并且两组资料分别呈正态分布，可以用成组t检验如果方差不齐，或者小样本情况下偏态分布，则用秩和检验(Ranksum test)
group
x
0
79
0
93
0
91
0
92
0
94
0
77
0
93
0
74
0
91
0
101
0
83
0
73
0
88
0
102
0
90
0
100
0
81
0
91
0
83
0
106
0
84
0
78
0
87
0
95
0
101
1
101
1
100
1
114
1
86
1
106
1
107
1
107
1
94
1
89
1
104
1
98
1
110
1
89
1
103
1
89
1
121
1
94
1
95
1
92
1
109
1
98
1
98
1
120
1
104
1
110
多组比较完全随机分组设计(要求各组资料之间相互独立)
方差齐性并且独立以及每一组资料都服从正态分布(小样本时要求)，则采用完全随机设计的方差分析方法(即：单因素方差分析，One Way ANOVA)进行分析。
方差不齐或小样本情况下资料偏态，则用Kruskal Wallis 检验(H检验)
例5.1 为研究胃癌与胃粘膜细胞中DNA含量（A.U）的关系，某医师测得数据如下，试问四组人群的胃粘膜细胞中平均DNA含量是否相同？
组别
group
DNA含量（A.U）
浅表型胃炎
1
9.81
12.73
12.29
12.53
12.95
9.53
12.6
8.9
12.27
14.26
10.68
肠化生
2
14.61
17.54
15.1
17
13.39
15.32
13.74
18.24
13.81
12.63
14.53
16.17
早期胃癌
3
23.26
20.8
20.6
23.5
17.85
21.91
22.13
22.04
19.53
18.41
21.48
20.24
晚期胃癌
4
23.73
19.46
22.39
19.53
25.9
20.43
20.71
20.05
23.41
21.34
21.38
25.70
由于这四组对象的资料是相互独立的，因此属于完全随机分组类型的。检验问题是考察四组DNA含量的平均水平相同吗。如果每一组资料都正态分布并且方差齐性可以用One way-ANOVA进行分析，反之用Kruskal Wallis检验。
STATA数据输入格式
g
x
1
9.81
1
12.73
1
12.29
1
12.53
1
12.95
1
9.53
1
12.6
1
8.9
1
12.27
1
14.26
1
10.68
2
14.61
2
17.54
2
15.1
2
17
2
13.39
2
15.32
2
13.74
2
18.24
2
13.81
2
12.63
2
14.53
2
16.17
3
23.26
3
20.8
3
20.6
3
23.5
3
17.85
3
21.91
3
22.13
3
22.04
3
19.53
3
18.41
3
21.48
3
20.24
4
23.73
4
19.46
4
22.39
4
19.53
4
25.9
4
20.43
4
20.71
4
20.05
4
23.41
4
21.34
4
21.38
4
25.7
分组正态性检验，(=0.05
,sktest x if g==1
Skewness/Kurtosis tests for Normality
------- joint ------
Variable | Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi2
-------------+-------------------------------------------------------
x | 0.491 0.485 1.07 0.5861
,sktest x if g==2
Skewness/Kurtosis tests for Normality
------- joint ------
Variable | Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi2
-------------+-------------------------------------------------------
x | 0.482 0.541 0.96 0.6201
,sktest x if g==3
Skewness/Kurtosis tests for Normality
------- joint ------
Variable | Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi2
-------------+-------------------------------------------------------
x | 0.527 0.750 0.52 0.7704
,sktest x if g==4
Skewness/Kurtosis tests for Normality
------- joint ------
Variable | Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi2
-------------+-------------------------------------------------------
x | 0.260 0.616 1.75 0.4166
上述结果表明每一组资料都服从正态分布。
单因素方差分析的STATA命令：oneway 效应指标变量 分组变量,t b
其中t表示计算每一组均数和标准差,b表示采用Bonferroni统计方法进行两两比较。
本例命令为oneway x group,t b
,oneway x g,t b
| Summary of x
g | Mean Std,Dev,Freq.
------------+------------------------------------
1 | 11.686364 1.6884388 11
2 | 15.173333 1.749173 12
3 | 20.979167 1.7668279 12
4 | 22.0025 2.2429087 12
------------+------------------------------------
Total | 17.583191 4.6080789 47
Analysis of Variance
Source SS df MS F Prob > F
------------------------------------------------------------------------
Between groups 824.942549 3 274.98085 77.87 0.0000
Within groups 151.839445 43 3.53114987
------------------------------------------------------------------------
Total 976.781994 46 21.2343912
Bartlett's test for equal variances,chi2(3) = 1.1354 Prob>chi2 = 0.769
方差齐性的检验为：卡方=1.1354，自由度=3,P值=0.769，因此可以认为方差是齐性的。
H0：(1=(2=(3=(4 四组总体均数相同
H1：(1，(2，(3，(4不全相同
(=0.05，相应的统计量F=77.87以及相应的自由度为3和43，P值<0.0001，因此4组均数的差别有统计学意义。
Comparison of x by g
(Bonferroni)
Row Mean-|
Col Mean | 1 2 3
---------+---------------------------------
2 | 3.48697(第2组样本均数－第1组样本均数)
| 0.000(H0:(1=(2检验的P值)
|
3 | 9.2928 5.80583(第3组样本均数－第2组样本均数)
| 0.000 0.000(H0:(3=(2检验的P值)
|
4 | 10.3161 6.82917 1.02333(第4组样本均数-第3组样本均数)
| 0.000 0.000 1.000(H0:(3=(4检验的P值)\
上述输出为两两比较的结果，在表格的每个单元中，第一行为两组均数的差值，第二行为两组均数比较检验的P值。
根据上述结果可以知道，第2组、第3组和第4组的AU均数均大于第1组的AU均数，并且差别有统计学意义。说明肠化生患者和胃癌患者的DNA的AU含量平均水平均高于正常人的AU平均水平，并且差别有统计学意义。
第3组和第4组的AU均数也大于第2组的AU平均水平，并且差别有统计学意义。说明胃癌患者的DNA的AU含量平均水平均高于肠化生患者的AU平均水平，并且差别有统计学意义。
第3组和第4组两组均数的差别没有统计学意义，说明没有足够的证据可以DNA的AU含量与癌症的早期与晚期有关系。
假如本例的资料不满足方差分析的要求，则用Kruskal Wallis检验，数据结构同上。命令为：
kwallis 效应指标变量,by(分组变量)
本例的命令为 kwallis x,by(g)
H0：4组的AU总体分布相同
H1：4组的AU总体分布不全相同
(=0.05
结果如下：
Test,Equality of populations (Kruskal-Wallis test)
g _Obs _RankSum
1 11 72.00
2 12 205.00
3 12 411.50
4 12 439.50
chi-squared = 37.814 with 3 d.f.
probability = 0.0001
chi-squared with ties = 37.816 with 3 d.f.
probability = 0.0001
说明：4组AU的总体分布不全相同，然后秩和检验，但(应取小一些(多重比较时，会增大第一类错误的概率)。根据Sidak检验的建议：，其中k为要比较的次数，(为多组比较总的检验水平(一般为0.05)，(’为两两比较时的检验水平。
如本例：4组两两比较共比次，因此，
对于比较第1组和第2组的AU分布差别的操作命令为：
先计算中位数
sort g 组别变量排序
by g:centile x,centile(50) 计算各组中位数
-> g = 1
-- Binom,Interp,--
Variable | Obs Percentile Centile [95% Conf,Interval]
-------------+-------------------------------------------------------------
x | 11 50 12.29 9.729564 12.7932
-> g = 2
-- Binom,Interp,--
Variable | Obs Percentile Centile [95% Conf,Interval]
-------------+-------------------------------------------------------------
x | 12 50 14.855 13.74745 16.91172
-> g = 3
-- Binom,Interp,--
Variable | Obs Percentile Centile [95% Conf,Interval]
-------------+-------------------------------------------------------------
x | 12 50 21.14 19.60552 22.12043
-> g = 4
-- Binom,Interp,--
Variable | Obs Percentile Centile [95% Conf,Interval]
-------------+-------------------------------------------------------------
x | 12 50 21.36 20.09042 23.69596
得到这4组中位数分别为：M1=12.29，M2=14.855，M3=21.14和M4=21.36
ranksum x if g==1 | g==2,by(g)
Two-sample Wilcoxon rank-sum (Mann-Whitney) test
g | obs rank sum expected
-------------+---------------------------------
1 | 11 72 132
2 | 12 204 144
-------------+---------------------------------
combined | 23 276 276
unadjusted variance 264.00
adjustment for ties 0.00
----------
adjusted variance 264.00
Ho,x(g==1) = x(g==2)
z = -3.693
Prob > |z| = 0.0002
P值<(’，因此第2组AU的平均水平要高于第1组的平均水平(M2>M1)，并且差别有统计学意义。
第1组与第3组比较
ranksum x if g==1 | g==3,by(g)
Two-sample Wilcoxon rank-sum (Mann-Whitney) test
g | obs rank sum expected
-------------+---------------------------------
1 | 11 66 132
3 | 12 210 144
-------------+---------------------------------
combined | 23 276 276
unadjusted variance 264.00
adjustment for ties 0.00
----------
adjusted variance 264.00
Ho,x(g==1) = x(g==3)
z = -4.062
Prob > |z| = 0.0000
P值<(’，因此第3组AU的平均水平要高于第1组的平均水平(M3>M1)，并且差别有统计学意义，其他比较类似进行。
要注意的问题：
在方差分析中，要求每一组资料服从正态分布(小样本时)，并不是要求各组资料服从一个正态分布(因为这就意味各组的总体均数相同，失去统计检验的必要性)，所以不能把各组的资料合在一起作正态性检验。总的讲，方差分析对正态性具有稳健性，即：偏态分布对方差分析的结果影响不会太大，故正态性检验的(取0.05也就可以了。
样本量较大时，方差分析对正态性要求大大降低(根据中心极限定理可知：样本均数近似服从正态分布)。并且由于大多数情况下，样本资料只是近似服从正态分布而不是完全服从正态分布。由于在大样本情况下，用正态性检验就变为很敏感，对于不是完全服从正态分布的资料往往会拒绝正态性检验的H0：资料服从正态分布。因为正态性检验不能检验资料是否近似服从正态分布，而是检验是否服从正态分布。故在大样本情况下，考察资料的近似正态性，应用频数图进行考察。
方差齐性问题对方差分析相对比较敏感，并且并不是随着样本量增大而方差齐性对方差分析减少影响的。但是当各组样本量接近相同或相同时，方差齐性对方差分析呈现某种稳健性。即：只有当各组样本量相同时，方差齐性对方差分析结果的影响大大降低。这时随着样本量增大，影响会进一步降低。相反，如果各组样本量相差太大时，方差齐性对方差分析结果的影响很大。这时随着样本量增大，影响会进一步加大。
随机区组设计(处理组之间可能不独立)
残差(定义为：，也就是随机区组方差分析中的误差项)的方差齐性且小样本时正态分布，则用随机区组的方差分析(无重复的两因素方差分析,Two-way ANOVA)。
不满足方差齐性或小样本时资料偏态，则对用秩变换后再用随机区组的方差分析也可以直接用非参数随机区组的秩和检验Fredman test)。
例2下表是某湖水中8个观察地点不同季节取样的氯化物含量测定值，请问在不同季节该湖水中氯化物的含量有无差别？
表2 某湖水中不同季节的氯化物含量测定值（mg/L）
location no
春
夏
秋
冬
1
21.28
18.33
17.27
14.91
2
22.78
19.81
16.55
14.85
3
20.90
18.93
16.36
16.30
4
19.90
21.23
17.86
15.73
5
21.49
19.09
15.11
17.05
6
22.38
17.92
16.57
14.34
7
21.67
19.39
17.19
16.31
8
22.06
19.65
16.58
14.33
显然同一地点不同季节的氯化物含量有一定的相关性，故不能采用完全随机设计的方差分析方法对4个季节的氯化物含量进行统计分析。可以把同一地点的4个季节氯化物含量视为一个区组，因此可以用随机区组的方差分析进行统计分析。
设第8个地点在冬季的氯化物总体均数为(0，同样在冬季，第i个地点的氯化物总体均数与第8个地点在冬季的氯化物总体均数相差(i，i=1，2，3，4，5，6，7。因此在冬季的这8个地点在冬季的氯化物总体均数可以表示为地点编号
1
2
3
4
5
6
7
8
冬季氯化物均数
(0+(1
(0+(2
(0+(3
(0+(4
(0+(5
(0+(6
(0+(7
(0
假定在同一地区，春季的氯化物总体均数与冬季的氯化物总体均数相差(1，因此春节和冬季的氯化物总体均数可以表示为地点编号
1
2
3
4
5
6
7
8
冬季氯化物均数
(0+(1
(0+(2
(0+(3
(0+(4
(0+(5
(0+(6
(0+(7
(0
春季氯化物均数
(0+(1+(1
(0+(1+(2
(0+(1+(3
(0+(1+(4
(0+(1+(5
(0+(1+(6
(0+(1+(7
(0
如果(1＝0说明在同一地点，冬季和春季的氯化物总体均数相同；(1>0说明春季的氯化物含量平均高于冬季氯化物含量，反之(<0，说明春季氯化物含量均数低于冬季氯化物含量。
同理假定在同一地区，夏季和秋季的氯化物总体均数与冬季的氯化物总体均数分别相差(2和(3，则四个季节的氯化物总体均数可以表示为
地点编号
1
2
3
4
5
6
7
8
冬季氯化物均数
(0+(1
(0+(2
(0+(3
(0+(4
(0+(5
(0+(6
(0+(7
(0
春季氯化物均数
(0+(1+(1
(0+(1+(2
(0+(1+(3
(0+(1+(4
(0+(1+(5
(0+(1+(6
(0+(1+(7
(0
夏季氯化物均数
(0+(2+(1
(0+(2+(2
(0+(2+(3
(0+(2+(4
(0+(2+(5
(0+(2+(6
(0+(2+(7
(0
春季氯化物均数
(0+(3+(1
(0+(3+(2
(0+(3+(3
(0+(3+(4
(0+(3+(5
(0+(3+(6
(0+(3+(7
(0
根据上述总体均数表示，可以知道：在四个季节中的氯化物总体均数(同一地点)无变化就是H0：(1=(2=(3=0(在随机区组方差分析中称为无处理效应，但不能称4组的总体均数相同，因为在同一季节中不同地点的总体均数可能不同)。
H1：(1，(2，(3不全为0
Stata 数据输入格式
t
id
x
1
1
21.27589
1
2
22.77649
1
3
20.89943
1
4
19.9043
1
5
21.4929
1
6
22.38085
1
7
21.67344
1
8
22.06133
2
1
18.33405
2
2
19.80538
2
3
18.92919
2
4
21.22814
2
5
19.09215
2
6
17.9237
2
7
19.38569
2
8
19.64971
3
1
17.27141
3
2
16.54567
3
3
16.36019
3
4
17.85548
3
5
15.11296
3
6
16.56507
3
7
17.18734
3
8
16.58279
4
1
14.90559
4
2
14.85127
4
3
16.29782
4
4
15.7286
4
5
17.05169
4
6
14.34088
4
7
16.31367
4
8
14.33015
其中id表示观察地点编号，t＝1，2，3，4对应表示春节、夏季、秋季和冬季。
Stata操作命令：
anova x t id
,anova x t id
Number of obs = 32 R-squared = 0.8923
Root MSE = 1.01769 Adj R-squared = 0.8410
Source | Partial SS df MS F Prob > F
-----------+----------------------------------------------------
Model | 180.214326 10 18.0214326 17.40 0.0000
|
t | 177.344737 3 59.1149122 57.08 0.0000
id | 2.86958916 7,409941308 0.40 0.8942
|
Residual | 21.749618 21 1.0356961
-----------+----------------------------------------------------
Total | 201.963944 31 6.51496593
处理效应H0：(1＝(2＝(3＝0的检验对应的统计量
相应的P值<0.0001(计算机输出值是0.0000)，所以拒绝无效假设，可以认为4个季节的氯化物总体均数不全相同。
不同季节中的两两比较用LSD方法检验如下：
在输入anova x t id命令后，再输入regress命令便得到下列结果
Source | SS df MS Number of obs = 32
-------------+------------------------------ F( 10,21) = 17.40
Model | 180.214326 10 18.0214326 Prob > F = 0.0000
Residual | 21.749618 21 1.0356961 R-squared = 0.8923
-------------+------------------------------ Adj R-squared = 0.8410
Total | 201.963944 31 6.51496593 Root MSE = 1.0177
------------------------------------------------------------------------------
x Coef,Std,Err,t P>|t| [95% Conf,Interval]
------------------------------------------------------------------------------
_cons ((0) 15.37992,5966746 25.78 0.000 14.13906 16.62077
t
(1= 1 6.080619,5088458 11.95 0.000 5.022417 7.138822
(2= 2 3.816041,5088458 7.50 0.000 2.757838 4.874244
(3= 3 1.20765,5088458 2.37 0.027,1494472 2.265853
4 (dropped)
id
(1= 1 -.2092595,7196166 -0.29 0.774 -1.705784 1.287265
(2= 2,3387067,7196166 0.47 0.643 -1.157818 1.835231
(3= 3 -.034339,7196166 -0.05 0.962 -1.530864 1.462186
(4= 4,5231357,7196166 0.73 0.475 -.973389 2.01966
(5= 5,0314307,7196166 0.04 0.966 -1.465094 1.527955
(6= 6 -.353369,7196166 -0.49 0.628 -1.849894 1.143156
(7= 7,4840407,7196166 0.67 0.509 -1.012484 1.980565
8 (dropped)
其中，对应的假设检验H0：(1=0的统计量t=11.95，P值<0.001，95%可信区间为(5.022，7.139)，因此可以认为春季的氯化物平均高于冬季，差别有统计学意义。
，对应的假设检验H0：(2=0的统计量t=7.50，P值<0.001，95%可信区间为 (2.758，4.874)，因此可以认为夏季的氯化物平均高于冬季，差别有统计学意义。
，对应的假设检验H0：(3=0的统计量t=2.37，P值=0.027，95%可信区间为 (0.1494，2.266)，因此可以认为秋季的氯化物平均高于冬季，差别有统计学意义。
对于春季氯化物平均数((0+(1+(i)与夏季的氯化物平均数((0+(2+(i)比较对应为(1>(2、(1=(2和(1<(2的问题。因此需要检验H0：(1＝(2 vs H1：(1((2,相应的STATA命令(anova x t id命令和regress命令后)为test b[t[1]]=_b[t[2]]，得到下列结果
( 1) t[2] - t[3] = 0.0
F( 1,21) = 26.28
Prob > F = 0.0000
相应的统计量F=26.28，P值<0.0001，差别有统计学意义。由于(1的估计值>(2的估计值，所以可以认为春季氯化物平均高于夏季的氯化物含量。
同理检验H0：(1＝(3 vs H1：(1((3，只需输入命令test b[t[1]]=_b[t[3]]
检验H0：(2＝(3 vs H1：(2((3，只需输入命令test b[t[2]]=_b[t[3]]
此处不在详细叙述了。
由于随机区组方差分析要求残差()服从正态分布，再输入regress以后，只要输入predict 残差变量名,residual，就可以得到残差计算值。
本例用e表示残差变量名，因此输入predict e,residual
就可以得到残差计算值e，然后对残差进行正态性检验(sktest 残差变量名)
本例输入命令为,sktest e
结果如下（H0：残差服从正态分布 vs H1：残差偏态分布,(＝0.05）
Skewness/Kurtosis tests for Normality
------- joint ------
Variable | Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi2
-------------+-------------------------------------------------------
e | 0.699 0.586 0.46 0.7948
P值=0.93349>>(，因此可以认为资料近似服从正态分布。(大样本时，可以不考虑正态性问题)
如果资料呈偏态分布，可以对资料进行秩变换(Rank Transform)后，然后把变换后的秩视为原始数据进行随机区组的方差分析。
秩变换的STATA命令为 egen 秩变量名＝rank(观察变量名),by(区组变量)
为了说明上述操作分析的过程，故借用本例资料进行秩变换操作说明如下(本例资料正态分布，无需用秩变换，只是说明操作而言).
设用r表示秩变量名，则本例操作为
egen r=rank(x),by(id) 产生秩r
anova 命令anova r t id 结果如下
Number of obs = 32 R-squared = 0.9125
Root MSE =,408248 Adj R-squared = 0.8708
Source | Partial SS df MS F Prob > F
-----------+----------------------------------------------------
Model | 36.50 10 3.65 21.90 0.0000
|
t | 36.50 3 12.1666667 73.00 0.0000
id | 0.00 7 0.00 0.00 1.0000
|
Residual | 3.50 21,166666667
-----------+----------------------------------------------------
Total | 40.00 31 1.29032258
命令regress 结果如下
regress
Source | SS df MS Number of obs = 32
-------------+------------------------------ F( 10,21) = 21.90
Model | 36.50 10 3.65 Prob > F = 0.0000
Residual | 3.50 21,166666667 R-squared = 0.9125
-------------+------------------------------ Adj R-squared = 0.8708
Total | 40.00 31 1.29032258 Root MSE =,40825
------------------------------------------------------------------------------
r Coef,Std,Err,t P>|t| [95% Conf,Interval]
------------------------------------------------------------------------------
_cons 1.125,2393568 4.70 0.000,6272303 1.62277
t
1 2.75,2041241 13.47 0.000 2.325501 3.174499
2 2,2041241 9.80 0.000 1.575501 2.424499
3,75,2041241 3.67 0.001,3255006 1.174499
4 (dropped)
id
1 0,2886751 0.00 1.000 -.6003328,6003328
2 0,2886751 0.00 1.000 -.6003328,6003328
3 0,2886751 0.00 1.000 -.6003328,6003328
4 0,2886751 0.00 1.000 -.6003328,6003328
5 0,2886751 0.00 1.000 -.6003328,6003328
6 0,2886751 0.00 1.000 -.6003328,6003328
7 0,2886751 0.00 1.000 -.6003328,6003328
8 (dropped)
------------------------------------------------------------------------------
进一步两两比较
test _b[t[1]]=_b[t[2]]
( 1) t[1] - t[2] = 0.0
F( 1,21) = 13.50
Prob > F = 0.0014
,test _b[t[1]]=_b[t[3]]
( 1) t[1] - t[3] = 0.0
F( 1,21) = 96.00
Prob > F = 0.0000
,test _b[t[2]]=_b[t[3]]
( 1) t[2] - t[3] = 0.0
F( 1,21) = 37.50
Prob > F = 0.0000
解释如同上述，不再重复。
复习和补充两组比较的统计检验配对设计资料(又称为Dependent Samples)
对于小样本的情况下，如果配对的差值资料服从正态分布，用配对t检验(ttest 差值变量=0)
大样本的情况下，可以用配对t检验小样本的情况下，并且配对差值呈偏态分布，则用配对符号秩检验(signrank 差值变量=0)
成组设计(Two Independent Samples)
如果方差齐性并且大样本情况下，可以用成组t检验(ttest 效应指标变量,by(分组变量))
如果方差齐性并且两组资料分别呈正态分布，可以用成组t检验如果方差不齐，或者小样本情况下偏态分布，则用秩和检验(Ranksum test)
group
x
0
79
0
93
0
91
0
92
0
94
0
77
0
93
0
74
0
91
0
101
0
83
0
73
0
88
0
102
0
90
0
100
0
81
0
91
0
83
0
106
0
84
0
78
0
87
0
95
0
101
1
101
1
100
1
114
1
86
1
106
1
107
1
107
1
94
1
89
1
104
1
98
1
110
1
89
1
103
1
89
1
121
1
94
1
95
1
92
1
109
1
98
1
98
1
120
1
104
1
110
多组比较完全随机分组设计(要求各组资料之间相互独立)
方差齐性并且独立以及每一组资料都服从正态分布(小样本时要求)，则采用完全随机设计的方差分析方法(即：单因素方差分析，One Way ANOVA)进行分析。
方差不齐或小样本情况下资料偏态，则用Kruskal Wallis 检验(H检验)
例5.1 为研究胃癌与胃粘膜细胞中DNA含量（A.U）的关系，某医师测得数据如下，试问四组人群的胃粘膜细胞中平均DNA含量是否相同？
组别
group
DNA含量（A.U）
浅表型胃炎
1
9.81
12.73
12.29
12.53
12.95
9.53
12.6
8.9
12.27
14.26
10.68
肠化生
2
14.61
17.54
15.1
17
13.39
15.32
13.74
18.24
13.81
12.63
14.53
16.17
早期胃癌
3
23.26
20.8
20.6
23.5
17.85
21.91
22.13
22.04
19.53
18.41
21.48
20.24
晚期胃癌
4
23.73
19.46
22.39
19.53
25.9
20.43
20.71
20.05
23.41
21.34
21.38
25.70
由于这四组对象的资料是相互独立的，因此属于完全随机分组类型的。检验问题是考察四组DNA含量的平均水平相同吗。如果每一组资料都正态分布并且方差齐性可以用One way-ANOVA进行分析，反之用Kruskal Wallis检验。
STATA数据输入格式
g
x
1
9.81
1
12.73
1
12.29
1
12.53
1
12.95
1
9.53
1
12.6
1
8.9
1
12.27
1
14.26
1
10.68
2
14.61
2
17.54
2
15.1
2
17
2
13.39
2
15.32
2
13.74
2
18.24
2
13.81
2
12.63
2
14.53
2
16.17
3
23.26
3
20.8
3
20.6
3
23.5
3
17.85
3
21.91
3
22.13
3
22.04
3
19.53
3
18.41
3
21.48
3
20.24
4
23.73
4
19.46
4
22.39
4
19.53
4
25.9
4
20.43
4
20.71
4
20.05
4
23.41
4
21.34
4
21.38
4
25.7
分组正态性检验，(=0.05
,sktest x if g==1
Skewness/Kurtosis tests for Normality
------- joint ------
Variable | Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi2
-------------+-------------------------------------------------------
x | 0.491 0.485 1.07 0.5861
,sktest x if g==2
Skewness/Kurtosis tests for Normality
------- joint ------
Variable | Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi2
-------------+-------------------------------------------------------
x | 0.482 0.541 0.96 0.6201
,sktest x if g==3
Skewness/Kurtosis tests for Normality
------- joint ------
Variable | Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi2
-------------+-------------------------------------------------------
x | 0.527 0.750 0.52 0.7704
,sktest x if g==4
Skewness/Kurtosis tests for Normality
------- joint ------
Variable | Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi2
-------------+-------------------------------------------------------
x | 0.260 0.616 1.75 0.4166
上述结果表明每一组资料都服从正态分布。
单因素方差分析的STATA命令：oneway 效应指标变量 分组变量,t b
其中t表示计算每一组均数和标准差,b表示采用Bonferroni统计方法进行两两比较。
本例命令为oneway x group,t b
,oneway x g,t b
| Summary of x
g | Mean Std,Dev,Freq.
------------+------------------------------------
1 | 11.686364 1.6884388 11
2 | 15.173333 1.749173 12
3 | 20.979167 1.7668279 12
4 | 22.0025 2.2429087 12
------------+------------------------------------
Total | 17.583191 4.6080789 47
Analysis of Variance
Source SS df MS F Prob > F
------------------------------------------------------------------------
Between groups 824.942549 3 274.98085 77.87 0.0000
Within groups 151.839445 43 3.53114987
------------------------------------------------------------------------
Total 976.781994 46 21.2343912
Bartlett's test for equal variances,chi2(3) = 1.1354 Prob>chi2 = 0.769
方差齐性的检验为：卡方=1.1354，自由度=3,P值=0.769，因此可以认为方差是齐性的。
H0：(1=(2=(3=(4 四组总体均数相同
H1：(1，(2，(3，(4不全相同
(=0.05，相应的统计量F=77.87以及相应的自由度为3和43，P值<0.0001，因此4组均数的差别有统计学意义。
Comparison of x by g
(Bonferroni)
Row Mean-|
Col Mean | 1 2 3
---------+---------------------------------
2 | 3.48697(第2组样本均数－第1组样本均数)
| 0.000(H0:(1=(2检验的P值)
|
3 | 9.2928 5.80583(第3组样本均数－第2组样本均数)
| 0.000 0.000(H0:(3=(2检验的P值)
|
4 | 10.3161 6.82917 1.02333(第4组样本均数-第3组样本均数)
| 0.000 0.000 1.000(H0:(3=(4检验的P值)\
上述输出为两两比较的结果，在表格的每个单元中，第一行为两组均数的差值，第二行为两组均数比较检验的P值。
根据上述结果可以知道，第2组、第3组和第4组的AU均数均大于第1组的AU均数，并且差别有统计学意义。说明肠化生患者和胃癌患者的DNA的AU含量平均水平均高于正常人的AU平均水平，并且差别有统计学意义。
第3组和第4组的AU均数也大于第2组的AU平均水平，并且差别有统计学意义。说明胃癌患者的DNA的AU含量平均水平均高于肠化生患者的AU平均水平，并且差别有统计学意义。
第3组和第4组两组均数的差别没有统计学意义，说明没有足够的证据可以DNA的AU含量与癌症的早期与晚期有关系。
假如本例的资料不满足方差分析的要求，则用Kruskal Wallis检验，数据结构同上。命令为：
kwallis 效应指标变量,by(分组变量)
本例的命令为 kwallis x,by(g)
H0：4组的AU总体分布相同
H1：4组的AU总体分布不全相同
(=0.05
结果如下：
Test,Equality of populations (Kruskal-Wallis test)
g _Obs _RankSum
1 11 72.00
2 12 205.00
3 12 411.50
4 12 439.50
chi-squared = 37.814 with 3 d.f.
probability = 0.0001
chi-squared with ties = 37.816 with 3 d.f.
probability = 0.0001
说明：4组AU的总体分布不全相同，然后秩和检验，但(应取小一些(多重比较时，会增大第一类错误的概率)。根据Sidak检验的建议：，其中k为要比较的次数，(为多组比较总的检验水平(一般为0.05)，(’为两两比较时的检验水平。
如本例：4组两两比较共比次，因此，
对于比较第1组和第2组的AU分布差别的操作命令为：
先计算中位数
sort g 组别变量排序
by g:centile x,centile(50) 计算各组中位数
-> g = 1
-- Binom,Interp,--
Variable | Obs Percentile Centile [95% Conf,Interval]
-------------+-------------------------------------------------------------
x | 11 50 12.29 9.729564 12.7932
-> g = 2
-- Binom,Interp,--
Variable | Obs Percentile Centile [95% Conf,Interval]
-------------+-------------------------------------------------------------
x | 12 50 14.855 13.74745 16.91172
-> g = 3
-- Binom,Interp,--
Variable | Obs Percentile Centile [95% Conf,Interval]
-------------+-------------------------------------------------------------
x | 12 50 21.14 19.60552 22.12043
-> g = 4
-- Binom,Interp,--
Variable | Obs Percentile Centile [95% Conf,Interval]
-------------+-------------------------------------------------------------
x | 12 50 21.36 20.09042 23.69596
得到这4组中位数分别为：M1=12.29，M2=14.855，M3=21.14和M4=21.36
ranksum x if g==1 | g==2,by(g)
Two-sample Wilcoxon rank-sum (Mann-Whitney) test
g | obs rank sum expected
-------------+---------------------------------
1 | 11 72 132
2 | 12 204 144
-------------+---------------------------------
combined | 23 276 276
unadjusted variance 264.00
adjustment for ties 0.00
----------
adjusted variance 264.00
Ho,x(g==1) = x(g==2)
z = -3.693
Prob > |z| = 0.0002
P值<(’，因此第2组AU的平均水平要高于第1组的平均水平(M2>M1)，并且差别有统计学意义。
第1组与第3组比较
ranksum x if g==1 | g==3,by(g)
Two-sample Wilcoxon rank-sum (Mann-Whitney) test
g | obs rank sum expected
-------------+---------------------------------
1 | 11 66 132
3 | 12 210 144
-------------+---------------------------------
combined | 23 276 276
unadjusted variance 264.00
adjustment for ties 0.00
----------
adjusted variance 264.00
Ho,x(g==1) = x(g==3)
z = -4.062
Prob > |z| = 0.0000
P值<(’，因此第3组AU的平均水平要高于第1组的平均水平(M3>M1)，并且差别有统计学意义，其他比较类似进行。
要注意的问题：
在方差分析中，要求每一组资料服从正态分布(小样本时)，并不是要求各组资料服从一个正态分布(因为这就意味各组的总体均数相同，失去统计检验的必要性)，所以不能把各组的资料合在一起作正态性检验。总的讲，方差分析对正态性具有稳健性，即：偏态分布对方差分析的结果影响不会太大，故正态性检验的(取0.05也就可以了。
样本量较大时，方差分析对正态性要求大大降低(根据中心极限定理可知：样本均数近似服从正态分布)。并且由于大多数情况下，样本资料只是近似服从正态分布而不是完全服从正态分布。由于在大样本情况下，用正态性检验就变为很敏感，对于不是完全服从正态分布的资料往往会拒绝正态性检验的H0：资料服从正态分布。因为正态性检验不能检验资料是否近似服从正态分布，而是检验是否服从正态分布。故在大样本情况下，考察资料的近似正态性，应用频数图进行考察。
方差齐性问题对方差分析相对比较敏感，并且并不是随着样本量增大而方差齐性对方差分析减少影响的。但是当各组样本量接近相同或相同时，方差齐性对方差分析呈现某种稳健性。即：只有当各组样本量相同时，方差齐性对方差分析结果的影响大大降低。这时随着样本量增大，影响会进一步降低。相反，如果各组样本量相差太大时，方差齐性对方差分析结果的影响很大。这时随着样本量增大，影响会进一步加大。
随机区组设计(处理组之间可能不独立)
残差(定义为：，也就是随机区组方差分析中的误差项)的方差齐性且小样本时正态分布，则用随机区组的方差分析(无重复的两因素方差分析,Two-way ANOVA)。
不满足方差齐性或小样本时资料偏态，则对用秩变换后再用随机区组的方差分析也可以直接用非参数随机区组的秩和检验Fredman test)。
例2下表是某湖水中8个观察地点不同季节取样的氯化物含量测定值，请问在不同季节该湖水中氯化物的含量有无差别？
表2 某湖水中不同季节的氯化物含量测定值（mg/L）
location no
春
夏
秋
冬
1
21.28
18.33
17.27
14.91
2
22.78
19.81
16.55
14.85
3
20.90
18.93
16.36
16.30
4
19.90
21.23
17.86
15.73
5
21.49
19.09
15.11
17.05
6
22.38
17.92
16.57
14.34
7
21.67
19.39
17.19
16.31
8
22.06
19.65
16.58
14.33
显然同一地点不同季节的氯化物含量有一定的相关性，故不能采用完全随机设计的方差分析方法对4个季节的氯化物含量进行统计分析。可以把同一地点的4个季节氯化物含量视为一个区组，因此可以用随机区组的方差分析进行统计分析。
设第8个地点在冬季的氯化物总体均数为(0，同样在冬季，第i个地点的氯化物总体均数与第8个地点在冬季的氯化物总体均数相差(i，i=1，2，3，4，5，6，7。因此在冬季的这8个地点在冬季的氯化物总体均数可以表示为地点编号
1
2
3
4
5
6
7
8
冬季氯化物均数
(0+(1
(0+(2
(0+(3
(0+(4
(0+(5
(0+(6
(0+(7
(0
假定在同一地区，春季的氯化物总体均数与冬季的氯化物总体均数相差(1，因此春节和冬季的氯化物总体均数可以表示为地点编号
1
2
3
4
5
6
7
8
冬季氯化物均数
(0+(1
(0+(2
(0+(3
(0+(4
(0+(5
(0+(6
(0+(7
(0
春季氯化物均数
(0+(1+(1
(0+(1+(2
(0+(1+(3
(0+(1+(4
(0+(1+(5
(0+(1+(6
(0+(1+(7
(0
如果(1＝0说明在同一地点，冬季和春季的氯化物总体均数相同；(1>0说明春季的氯化物含量平均高于冬季氯化物含量，反之(<0，说明春季氯化物含量均数低于冬季氯化物含量。
同理假定在同一地区，夏季和秋季的氯化物总体均数与冬季的氯化物总体均数分别相差(2和(3，则四个季节的氯化物总体均数可以表示为
地点编号
1
2
3
4
5
6
7
8
冬季氯化物均数
(0+(1
(0+(2
(0+(3
(0+(4
(0+(5
(0+(6
(0+(7
(0
春季氯化物均数
(0+(1+(1
(0+(1+(2
(0+(1+(3
(0+(1+(4
(0+(1+(5
(0+(1+(6
(0+(1+(7
(0
夏季氯化物均数
(0+(2+(1
(0+(2+(2
(0+(2+(3
(0+(2+(4
(0+(2+(5
(0+(2+(6
(0+(2+(7
(0
春季氯化物均数
(0+(3+(1
(0+(3+(2
(0+(3+(3
(0+(3+(4
(0+(3+(5
(0+(3+(6
(0+(3+(7
(0
根据上述总体均数表示，可以知道：在四个季节中的氯化物总体均数(同一地点)无变化就是H0：(1=(2=(3=0(在随机区组方差分析中称为无处理效应，但不能称4组的总体均数相同，因为在同一季节中不同地点的总体均数可能不同)。
H1：(1，(2，(3不全为0
Stata 数据输入格式
t
id
x
1
1
21.27589
1
2
22.77649
1
3
20.89943
1
4
19.9043
1
5
21.4929
1
6
22.38085
1
7
21.67344
1
8
22.06133
2
1
18.33405
2
2
19.80538
2
3
18.92919
2
4
21.22814
2
5
19.09215
2
6
17.9237
2
7
19.38569
2
8
19.64971
3
1
17.27141
3
2
16.54567
3
3
16.36019
3
4
17.85548
3
5
15.11296
3
6
16.56507
3
7
17.18734
3
8
16.58279
4
1
14.90559
4
2
14.85127
4
3
16.29782
4
4
15.7286
4
5
17.05169
4
6
14.34088
4
7
16.31367
4
8
14.33015
其中id表示观察地点编号，t＝1，2，3，4对应表示春节、夏季、秋季和冬季。
Stata操作命令：
anova x t id
,anova x t id
Number of obs = 32 R-squared = 0.8923
Root MSE = 1.01769 Adj R-squared = 0.8410
Source | Partial SS df MS F Prob > F
-----------+----------------------------------------------------
Model | 180.214326 10 18.0214326 17.40 0.0000
|
t | 177.344737 3 59.1149122 57.08 0.0000
id | 2.86958916 7,409941308 0.40 0.8942
|
Residual | 21.749618 21 1.0356961
-----------+----------------------------------------------------
Total | 201.963944 31 6.51496593
处理效应H0：(1＝(2＝(3＝0的检验对应的统计量
相应的P值<0.0001(计算机输出值是0.0000)，所以拒绝无效假设，可以认为4个季节的氯化物总体均数不全相同。
不同季节中的两两比较用LSD方法检验如下：
在输入anova x t id命令后，再输入regress命令便得到下列结果
Source | SS df MS Number of obs = 32
-------------+------------------------------ F( 10,21) = 17.40
Model | 180.214326 10 18.0214326 Prob > F = 0.0000
Residual | 21.749618 21 1.0356961 R-squared = 0.8923
-------------+------------------------------ Adj R-squared = 0.8410
Total | 201.963944 31 6.51496593 Root MSE = 1.0177
------------------------------------------------------------------------------
x Coef,Std,Err,t P>|t| [95% Conf,Interval]
------------------------------------------------------------------------------
_cons ((0) 15.37992,5966746 25.78 0.000 14.13906 16.62077
t
(1= 1 6.080619,5088458 11.95 0.000 5.022417 7.138822
(2= 2 3.816041,5088458 7.50 0.000 2.757838 4.874244
(3= 3 1.20765,5088458 2.37 0.027,1494472 2.265853
4 (dropped)
id
(1= 1 -.2092595,7196166 -0.29 0.774 -1.705784 1.287265
(2= 2,3387067,7196166 0.47 0.643 -1.157818 1.835231
(3= 3 -.034339,7196166 -0.05 0.962 -1.530864 1.462186
(4= 4,5231357,7196166 0.73 0.475 -.973389 2.01966
(5= 5,0314307,7196166 0.04 0.966 -1.465094 1.527955
(6= 6 -.353369,7196166 -0.49 0.628 -1.849894 1.143156
(7= 7,4840407,7196166 0.67 0.509 -1.012484 1.980565
8 (dropped)
其中，对应的假设检验H0：(1=0的统计量t=11.95，P值<0.001，95%可信区间为(5.022，7.139)，因此可以认为春季的氯化物平均高于冬季，差别有统计学意义。
，对应的假设检验H0：(2=0的统计量t=7.50，P值<0.001，95%可信区间为 (2.758，4.874)，因此可以认为夏季的氯化物平均高于冬季，差别有统计学意义。
，对应的假设检验H0：(3=0的统计量t=2.37，P值=0.027，95%可信区间为 (0.1494，2.266)，因此可以认为秋季的氯化物平均高于冬季，差别有统计学意义。
对于春季氯化物平均数((0+(1+(i)与夏季的氯化物平均数((0+(2+(i)比较对应为(1>(2、(1=(2和(1<(2的问题。因此需要检验H0：(1＝(2 vs H1：(1((2,相应的STATA命令(anova x t id命令和regress命令后)为test b[t[1]]=_b[t[2]]，得到下列结果
( 1) t[2] - t[3] = 0.0
F( 1,21) = 26.28
Prob > F = 0.0000
相应的统计量F=26.28，P值<0.0001，差别有统计学意义。由于(1的估计值>(2的估计值，所以可以认为春季氯化物平均高于夏季的氯化物含量。
同理检验H0：(1＝(3 vs H1：(1((3，只需输入命令test b[t[1]]=_b[t[3]]
检验H0：(2＝(3 vs H1：(2((3，只需输入命令test b[t[2]]=_b[t[3]]
此处不在详细叙述了。
由于随机区组方差分析要求残差()服从正态分布，再输入regress以后，只要输入predict 残差变量名,residual，就可以得到残差计算值。
本例用e表示残差变量名，因此输入predict e,residual
就可以得到残差计算值e，然后对残差进行正态性检验(sktest 残差变量名)
本例输入命令为,sktest e
结果如下（H0：残差服从正态分布 vs H1：残差偏态分布,(＝0.05）
Skewness/Kurtosis tests for Normality
------- joint ------
Variable | Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi2
-------------+-------------------------------------------------------
e | 0.699 0.586 0.46 0.7948
P值=0.93349>>(，因此可以认为资料近似服从正态分布。(大样本时，可以不考虑正态性问题)
如果资料呈偏态分布，可以对资料进行秩变换(Rank Transform)后，然后把变换后的秩视为原始数据进行随机区组的方差分析。
秩变换的STATA命令为 egen 秩变量名＝rank(观察变量名),by(区组变量)
为了说明上述操作分析的过程，故借用本例资料进行秩变换操作说明如下(本例资料正态分布，无需用秩变换，只是说明操作而言).
设用r表示秩变量名，则本例操作为
egen r=rank(x),by(id) 产生秩r
anova 命令anova r t id 结果如下
Number of obs = 32 R-squared = 0.9125
Root MSE =,408248 Adj R-squared = 0.8708
Source | Partial SS df MS F Prob > F
-----------+----------------------------------------------------
Model | 36.50 10 3.65 21.90 0.0000
|
t | 36.50 3 12.1666667 73.00 0.0000
id | 0.00 7 0.00 0.00 1.0000
|
Residual | 3.50 21,166666667
-----------+----------------------------------------------------
Total | 40.00 31 1.29032258
命令regress 结果如下
regress
Source | SS df MS Number of obs = 32
-------------+------------------------------ F( 10,21) = 21.90
Model | 36.50 10 3.65 Prob > F = 0.0000
Residual | 3.50 21,166666667 R-squared = 0.9125
-------------+------------------------------ Adj R-squared = 0.8708
Total | 40.00 31 1.29032258 Root MSE =,40825
------------------------------------------------------------------------------
r Coef,Std,Err,t P>|t| [95% Conf,Interval]
------------------------------------------------------------------------------
_cons 1.125,2393568 4.70 0.000,6272303 1.62277
t
1 2.75,2041241 13.47 0.000 2.325501 3.174499
2 2,2041241 9.80 0.000 1.575501 2.424499
3,75,2041241 3.67 0.001,3255006 1.174499
4 (dropped)
id
1 0,2886751 0.00 1.000 -.6003328,6003328
2 0,2886751 0.00 1.000 -.6003328,6003328
3 0,2886751 0.00 1.000 -.6003328,6003328
4 0,2886751 0.00 1.000 -.6003328,6003328
5 0,2886751 0.00 1.000 -.6003328,6003328
6 0,2886751 0.00 1.000 -.6003328,6003328
7 0,2886751 0.00 1.000 -.6003328,6003328
8 (dropped)
------------------------------------------------------------------------------
进一步两两比较
test _b[t[1]]=_b[t[2]]
( 1) t[1] - t[2] = 0.0
F( 1,21) = 13.50
Prob > F = 0.0014
,test _b[t[1]]=_b[t[3]]
( 1) t[1] - t[3] = 0.0
F( 1,21) = 96.00
Prob > F = 0.0000
,test _b[t[2]]=_b[t[3]]
( 1) t[2] - t[3] = 0.0
F( 1,21) = 37.50
Prob > F = 0.0000
解释如同上述，不再重复。
