均值与比例比较及方差分析

关联性度量、两个样本或多个样本的均值和比例比较

本章重点介绍两个分类变量之间的关联性度量（列联表的χ2或费舍尔检验，以及比值比（OR）的计算），或数值变量与分类因子之间的关联性度量。在后一种情况下，我们将考虑两个独立（或非独立）样本的情形，以及双样本或多样本情况下的参数模型（学生 t‐检验）和非参数模型（威尔科克森检验）（方差分析（ANOVA）和克鲁斯卡尔‐沃利斯方差分析）。此外，还将讨论多重比较的邦弗罗尼校正法对治疗的校正，以及方差分析的线性趋势检验。双因素方差分析的情况将简要介绍，仅限于构建方差分析表和绘制交互作用图的主要命令。

2.1. 两组均值的比较

2.1.1. 独立样本

命令 ttest 可用于在假设总体中等方差（或不等方差，使用 unequal 选项）的情况下，对两组均值进行学生t检验。当基于低出生体重婴儿的指标来分析母亲体重时，可以通过命令 tabstat 或以下方式获得简要的描述性统计摘要（均值和标准差）。

. format lwt %4.1f 
. tabulate low, summarize(lwt)

出生体重	时期	末次月经时体重的汇总
		Mean 标准差 频数
<2500克	0	133.3 31.7 130
≥2500克	1	122.2 26.5 59
总计		129.8 30.6 189

就图形表示而言，可以通过直方图查看这两个测量序列的分布：

. 直方图 lwt, 按(低)

或箱线图（图2.1）：

. graph box lwt, over(low, relabel(1 "正常" 2 "低出生体重 (< 2.5 kg)")) /// 
b1title("婴儿体重") ytitle("母亲体重")

需要注意的是，在箱线图（命令 graph box ）的情况下，有两个可用选项： by() （在并列图形中显示条件分布）和 over() （在同一图形中显示条件分布）。后者不适用于直方图。

通过在选项 by() 中指定响应变量和分类因子来执行“常规”学生t检验：

. ttest lwt, by(low)

等方差双样本t检验

组	Obs	Mean	标准误	标准差	[95% 置信区间]
0	130	133.3	2.78238	31.72402	127.795 138.805
1	59	122.1356	3.457723	26.55928	115.2142 129.057
合并	189	129.8148	2.224323	30.57938	125.427 134.2027

差值 = 均值(0) - 均值(1)
t = 2.3537
H₀: 差值 = 0
自由度 = 187

Hₐ: 差值 < 0
Hₐ: 差值 != 0
Ha：差值 > 0

P值(T< t)= 0.9902
P值(|T| > |t| ) = 0.0196
P值(T >t) = 0.0098

在上述方法中，有两个明确定义的变量可用，其中一个作为响应变量，另一个作为分类因子。也可以处理两个测量序列（不一定具有相同的大小）。以下是一种可能的方法，主要用于演示如何使用命令 preserve 和 restore 来管理第二个数据表，而不会清除工作空间中已有的数据。第一步，将创建两个新变量 lwt1 和 lwt2 ，用于存储两组母亲的体重：

. preserve 
. gen lwt1=lwt if low == 0 (59 个缺失值已生成) 
. gen lwt2=lwt if low == 1 (130 个缺失值已生成)

可以使用 summarize 来验证这两个变量的特征，甚至将其与样本中的所有数据（lwt）进行比较：符号 lwt* 指示Stata考虑所有名称以lwt开头的变量（在当前情况下即lwt、lwt1和lwt2）：

. summarize lwt*

变量	Obs	Mean	标准差	Min	Max
lwt	189	129.8148	30.57938	80	250
lwt1	130	133.3	31.72402	85	250
lwt2	59	122.1356	26.55928	80	200

t检验的另一种表述形式为：

. t检验 lwt1==lwt2, 非配对

在这种情况下，必须指定 unpaired 选项。为了能够返回原始数据表（并删除中间生成的变量），不应忘记 restore 命令。

选项 welch 提供了一个修正的韦尔奇t检验（而 unequal 基于萨特思韦特近似）。当我们打算通过F检验正式检验等方差假设时，可以使用命令 sdtest ：

. 方差检验 lwt, 按(低)

方差比检验

组	Obs	Mean	标准误	标准差	[95% 置信区间]
0	130	133.3	2.78238	31.72402	127.795 138.805
1	59	122.1356	3.457723	26.55928	115.2142 129.057
合并	189	129.8148	2.224323	30.57938	125.427 134.2027

比值 = 标准差(0) / 标准差(1)
f = 1.4267
Ho：比率 = 1
自由度 = 129, 58

Ha：比率 < 1
Ha：比率 != 1
Ha：比率 > 1

P值(F< f) = 0.9356
2*P值(F> f) = 0.1289
P值(F >f) = 0.0644

Stata还提供了所谓的“即时”命令，这些命令仅需提供用于构建检验统计量的相关数据即可。对于学生t检验，感兴趣的命令是 ttesti ，其语法如下所示：

ttesti #obs1 #mean1 #sd1 #obs2 #mean2 #sd2[,选项2]

因此，该命令需要输入第一个样本的计数、均值和标准差，以及第二个样本的相同信息。选项（options2）对应于之前讨论过的不等方差情况下的选项（unequal 或 welch），以及风险 level 1 − α（level）。

2.1.2. 非独立样本

对于两组配对样本的情况，将采用与以两个变量形式表示的两系列测量相同的原则（此次，这两个变量具有相同的观测数），即一种如下类型的形式：

. t检验 x1== x2

其中x1和x2表示两组配对的测量序列。显然，假设两个变量的观测值以相同的顺序排列。在这种情况下，paired选项是可选的。

2.1.3. 非参数方法

用于独立样本的Wilcoxon秩和检验，基于观测值的秩，可通过命令 ranksum 获得：

. ranksum lwt, 按(低)

两样本威尔科克森秩和（曼‐惠特尼）检验

低出生体重	obs	秩和	期望
0	130	13217.5	12350
1	59	4737.5	5605
合并	189	17955	17955

未调整方差 121441.67
结的调整 -181.97
调整方差 121259.70

H₀: lwt(低出生体重==0)= lwt(低出生体重==1)
z = 2.491
Prob> |z| = 0.0127

对于两组配对样本，可以使用 signrank 命令以及与配对样本t检验几乎相同的语法来获得符号秩检验：

. signrank x1=x2

2.2. 两个比例的比较

2.2.1. 独立样本

Pearson χ2检验可以通过多种方式获得，通常来自能够构建列联表的命令。例如，在使用 tabulate 时，将添加选项 chi ：

. tabulate 低出生体重 吸烟, 卡方 期望

			吸烟者
体重		非吸烟者	吸烟者
<2.5 公斤	0	86	44
≥2.5 公斤	1	29	30
总计		115	74

皮尔逊卡方(1) = 4.9237
P值 = 0.026

expected 选项可以显示表格中每个单元格的理论频数（即在两个变量相互独立的假设下的期望频数）。使用 exact 选项时，Stata还会返回费舍尔检验的结果。

通过包含直接命令 tabi 也可以得到相同的结果。与学生t检验的情况一样，需要提供构建检验统计量所需的关键信息。在当前情况下，这涉及列联表中四个单元格的频数，即：

. tabi 86 44\ 29 30, 卡方检验 精确检验 nofreq

皮尔逊卡方(1) = 4.9237 P值 = 0.026
费舍尔精确检验 = 0.036 单侧
费舍尔精确检验 = 0.020

需要注意的是，频数的输入按行进行，使用反斜杠符号 \ 分隔表格中的各行。

关于单样本的比例检验（通常与零假设 H0： p= 0.5 相关），可通过命令 bitest 进行二项检验（对应的直接命令为 bitesti ）。如果希望使用正态近似，则可通过 prtest （prtesti）获得相应的检验。该命令也可用于双样本情况。

下文的示例应用考虑了母亲根据吸烟状态的分布：

. 二项检验 吸烟 == 0.5, 详细

变量	N观察值	k	期望 k	假设 p	观察值 p
吸烟	189	74	94.5	0.50000	0.39153

P值(k>= 74) = 0.998917 (单侧检验)
P值(k<= 74) = 0.001754 (单侧检验)
P值(k<= 74或 k >= 115) = 0.003508 （双侧检验）
P值(k== 74) = 0.000671 (观察值)
P值(k== 114) = 0.001029
P值(k== 115) = 0.000671 （相反的极端）

. prtest 吸烟 == 0.5

单样本比例检验
吸烟：观测数 = 189

变量	Mean	标准误	[95% 置信区间]
吸烟	.3915344	.0355036	.3219487 .4611201

p=比例(吸烟)
z=-2.9823
H₀: p = 0.5
Ha: p < 0.5 Ha: p != 0.5 Ha: p > 0.5
P值(Z < z) = 0.0014 P值(|Z| > |z| ) = 0.0029 P值(Z > z) = 0.9986

在第一种情况下，所关注的量是名为 Pr(k <= 74或 k >=115) 的值，其对应值可在基于正态分布的比例检验中的 Pr(|Z| > |z| ) 下找到。对于两个样本的情况，以下命令可用于检验吸烟的母亲的分布是否在婴儿出生时的体重状态不同时仍保持相同的假设（可以很容易地将双尾检验的结果与 χ2 检验的结果进行比较）：

. prtest smoke, by(low)

两样本比例检验
0：观测数 = 130
1: 观测数 = 59

| 变量 | Mean | 标准误 | z | P>|z| | [95% 置信区间] |
| — | — | — | — | — | — |
| 0 | .3384615 | .0415012 | | | .2571207 .4198024 |
| 1 | .5084746 | .0650851 | | | .3809101 .636039 |
| 差值 | -.170013 | .0771908 | | | -.3213042 -.0187219 |
| 在零假设下 | .0766189 | -2.22 | 0.026 | | |

差值 =prop(0) - prop(1)
z=-2.2189
零假设：差值 = 0
Ha: 差值 < 0 Ha: 差值 != 0 Ha: 差值 > 0
P值(Z< z)= 0.0132 P值(|Z| < |z| ) = 0.0265 P值(Z >z) = 0.9868

命令 graph bar （条形图）或 graph dot （点图）可以显示和汇总任何频数表格。图2.2 展示了一个条形图的示例，其中我们结合了第14页 使用的方法（创建一个用于观测值计数的辅助变量）和第26页 的方法（为变量添加标签）：

. replace freq = 1// 变量已存在（0 个实际更改） 
. graph bar (sum) freq, over(smoke, relabel(1 "非吸烟" 2 "吸烟")) /// asyvars over(low, relabel(1 "正常" 2 "低出生体重 (< 2.5 kg)")) /// legend(title("母亲")) ytitle("计数")

对于此类图形表示，可以通过输入以下命令来安装 catplot 包：

. ssc install catplot

在Stata命令提示符处。这需要有可用的互联网连接。通过提供一些用于自定义图表的选项（图例、变量模态的标签等），前述命令显示为（图2.3）：

. catplot low smoke, recast(bar) /// var1opts(relabel(1"正常" 2 "低 (< 2.5 kg)"))

有关更多详细信息，请参阅在线帮助：帮助 catplot。

2.2.2. 非独立样本

对于同一样本或两组配对样本中的两个二元变量（例如，病例对照研究），如果由这两个变量交叉列表生成的列联表的边缘频数具有研究意义，则可使用命令 mcc 进行麦克尼马尔检验。其语法遵循流行病学研究（病例/对照、队列研究）中常用的表示法，Stata提供基于χ2(1)分布的精确显著性水平和近似显著性水平。即时命令是 mcci 。

2.3. 风险度量和比值比

流行病学中使用的大多数风险或关联度量都可以在 Stata 命令的一个特定子集（称为 epitab）中找到。它们也可以通过一个特定菜单访问，Statistics → Epidemiology and related → Tables for epidemiologists。

例如， tabodds 命令用于病例对照研究或横断面研究。它可用于计算比值比（OR）及其渐近置信区间，并检验各层间比值比的同质性（曼特尔‐亨塞尔检验）。以下是一个使用示例，其中将出生体重状态（low）视为响应变量，变量 smoke 视为风险或暴露因素：

. tabodds low smoke, 比值比

吸烟	比值比	chi2	P>卡方检验	[95% 置信区间]
0	1.000000	.	.	.
1	2.021944	4.90	0.0269	1.069897 3.821169

同质性检验（等比）：chi2(1) = 4.90
Pr>chi2 = 0.0269
优势趋势的得分检验： chi2(1) = 4.90
Pr>chi2 = 0.0269

可以通过 base() 选项指示 Stata 第二个变量的哪个水平作为参考类别，也可以指定不同的置信区间（cornfield或woolf）。

在之前的示例中，个体数据是可获得的，但通常情况下，数据以汇总格式提供，即以如下形式呈现的列联表也可以重新绘制为三变量表形式：前两列表示两个二元变量各水平的交叉列表，第三列表示相应的频数。在这种情况下，语法保持不变：它表示作为反应和解释因素的变量名称。然而，需要指示Stata如何使用权重进行处理（对每个变量的两个模态进行交叉列表），并将频数作为fweight=选项传递。

假设数据如表2.1所示，则语法应为：

. tabodds low smoke [fweight=样本量],比值比

low	smoke	N
0	0	86
1	0	29
0	1	44
1	1	30

表2.1. 变量 low 和 smoke 交叉列表的频数表

请注意，对于这种类型的配置（计数位于特定列中的表格），创建条形图和点图的图形命令略有简化：不再需要定义计数变量，也不需要使用(sum)选项来按变量水平汇总频数。此外，可以通过指定(asis)选项将计数列视为主变量。

类似地，可以使用命令 tabulate 并通过频数进行加权，选项为 [fweight=N]（可简写为 [fw=N]），其中 N 指包含频数的变量。这样的表格可以存储在一个纯文本文件中，文件的第一行包含三个变量的名称（随后将使用 insheet 命令导入），或者也可以直接使用 input 命令输入数据。在这种情况下，可以在 preserve/restore 两个命令之间执行操作，以避免丢失当前会话的数据。用户应注意，仍需对变量重命名，以避免与工作空间中已有的变量产生冲突。如下图所示，在输入 preserve 命令后，通过将变量 low 和 smoke 的名称追加 b 来创建汇总数据表。

. input lowb smokeb N

低出生体重 母亲吸烟 样本量
1. 0 0 86
2. 1 0 29
3. 0 1 44
4. 1 1 30
5. 结束

. list 低出生体重-样本量 in 1/4

低出生体重	母亲吸烟	N
0	0	86
1	0	29
0	1	44
1	1	30

关于前面的例子，因此应写作：

. tabulate lowb smokeb [fw=样本量]

lowb	0	1	总计
0	86	44	130
1	29	30	59
总计	115	74	189

本示例将作为背景，用于说明点图中的一些自定义选项，特别是用户选择单位的坐标轴以及使用不同符号来突出分类变量的各个水平（图2.4）：

. 标签定义 状态 0 "正常体重" 1 "低体重" 
. 标签定义 吸烟者 0 "非吸烟者" 1 "吸烟者" 
. 为低体重组应用状态标签 
. 为吸烟组应用吸烟者标签 
. 生成点图 N，按吸烟组和低体重组分开展示 /// y轴范围(0-100) y 轴刻度(0到100，间隔20) /// 标记1，符号为空心圆 标记2，符号为 X

需要提醒的是，必须输入 restore 才能恢复到初始数据集。

2.4. 方差分析

2.4.1. 单因素方差分析

命令 summarize （结合 tabulate 或 by: ）或 tabstat 可用于自然地根据解释变量的 level 汇总响应变量的分布。然而，将会看到，数值摘要可以直接与生成方差分析表的命令相结合。关于图形方法，始终可以使用 graph bar 或 graph dot 直观地展示组均值的分布；只需将汇总统计量 r(总和) 替换为 r(均值) 即可。对于个体数据的分布，我们可以通过 histogram 命令构建频数（或比例）直方图，并在选项 by() 中指明分类因子。

考虑与婴儿出生时的体重（bwt）及母亲种族来源（race）相关的数据。

通过以下方式可以获得体重分布的频数（选项 freq）（图2.5）：

. 直方图 出生体重, 按(种族, col(3)) freq

执行单因素（固定效应）方差分析的主要命令是 oneway 。对于更复杂的模型，必须使用 anova 或 regress 。其用法相对简单：提供一个变量列表，在本例中为响应变量和解释变量。选项 tabulate 会自动在方差分析表格中添加一个汇总组均值和标准差的表格：

. 单因素方差分析 出生体重 种族, tabulate

种族	出生体重的汇总
	Mean 标准差 频数
白人	3102.7188 727.88615 96
黑人	2719.6923 638.68388 26
其他	2805.2836 722.19436 67
总计	2944.5873 729.2143 189

方差分析

来源	SS	df	MS	F	P值 >F
组间	5015725.25	2	2507862.63	4.91	0.0083
组内	94953930.6	186	510505.003
总计	99969655.8	188	531753.488

方差齐性Bartlett检验: 卡方(2) = 0.6595 Prob>卡方2= 0.719

需要注意的是，Stata 还会显示用于检验方差齐性的巴特利特检验结果。如果我们要使用莱文检验，则必须使用命令 robvar ，该命令会返回名为W0的检验统计量下的结果。

. robvar 出生体重, 按(种族)

种族	出生体重的求和
	Mean 标准差 频数
白人	3102.7188 727.88615 96
黑人	2719.6923 638.68388 26
其他	2805.2836 722.19436 67
总计	2944.5873 729.2143 189

W0= 0.44717123 df(2, 186) P值 >F = 0.64012002
W50= 0.46842949 df(2, 186) P值 >F = 0.62672105
W10= 0.45725627 df(2, 186) P值 >F = 0.63372775

语法与 oneway 命令所使用的略有不同：这是一个用于检验方差齐性的完整命令（例如 sdtest ，参见2.1.1节），它并不直接关联由 oneway 生成的ANOVA模型。

2.4.2. 均值的两两比较

关于所有均值对的比较（如前面例子中的三个比较），最简单的方法是在使用单因素方差分析时，添加多重检验校正选项之一（bonferroni、scheffe 或 sidak）。为了清晰起见，以下表达式中未显示方差分析表：

. 单因素方差分析 出生体重 种族, bonferroni 无方差分析

按种族比较出生体重（ bonferroni）

行均值-|均值	白人	黑人
黑人	-383.026	0.049
其他	-297.435	85.5913
	0.029	1.000

我们将通过执行一个简单的t检验来比较白人与黑人的结果，其显著性水平针对所有比较进行了调整。命令 ttest 返回检验统计量（r(t)）和 p‐值（r(p)），这可以通过在第一个命令后输入 return list 进行验证：

. quietly: t检验 出生体重 如果 种族 != 3, 按(种族) 
. 显示 r(p)*3 
.04853058

2.4.3. 线性趋势检验

要执行线性趋势检验，线性回归方法等同于将 oneway 命令替换为 regress 命令。例如，考虑变量 ftv ，它表示妊娠早期妇产科就诊次数。该变量取值介于0到6之间，大于2的值很少被观察到。此变量可通过 recode 命令重新编码为三分类变量，其语法非常简单：新级别列在旧级别旁边（使用符号=进行关联），而符号/提供了一种方式，如同in操作符的情况一样，用于指示值的范围（起始值/结束值）：

. recode ftv (0=0) (1=1) (2/6=2), gen(ftv2) (ftv 与 ftv2 之间有 12 处差异)

选项 gen() 可以生成一个新变量。可以通过包含这两个变量的简单交叉列表来验证重新编码是否正确完成：

. tabulate ftv2 ftv

如果没有其他说明，变量 race 将被视为数值变量，并且由于各 level 之间的距离相等（此处的情况正是如此，因为 level 已编码为 {1, 2, 3}），回归线斜率所对应的检验将提供所需的结果：

. 回归 出生体重 ftv2

来源	SS	df	MS	观测数 = 189
模型	542577.691	1	542577.691	F( 1, 187) = 1.02
残差	99427078.1	187	531695.605	P值> F值 = 0.3137
总计	99969655.8	188	531753.488	R平方 = 0.0054

| 出生体重 | 系数 | 标准误 | t | P>|t值| | [95% 置信区间] |
| — | — | — | — | — | — |
| ftv2 | 66.09496 | 65.42879 | 1.01 | 0.314 | -62.97845 195.1684 |
| _常数项 | 2898.775 | 69.78422 | 41.54 | 0.000 | 2761.11 3036.441 |

使用对比来评估线性趋势检验的方法仍然依赖于 regress 命令，但这一次需要 Stata 将变量 ftv2 视为分类变量，因此在其名称前加上操作符 i（参见在线帮助 help fvvarlist）。这样做的效果是将具有 k= 3 个水平的变量转换为 k−1= 2 个虚拟编码变量，用于编码分类变量的水平 j(j= 2, … , k)。

由于我们并不真正关心回归结果本身，因此在回归命令前加上 quietly: 指令以省略结果显示，同时要求 Stata 显示与解释变量相关的（正交）多项式对比。

最后这一操作通过在分组变量名称前加上操作符 p. 来实现。

. quietly: regress 出生体重 i.ftv2
. contrast p.ftv2, 无效应

边际线性预测的对比
边际效应 : as balanced

	df	F	P>F
ftv2	（线性）	1	0.41
	(quadratic)	1	2.55
	联合	2	1.79
残差		186

此处关注的对比在名称（线性）下提及。

我们将能够验证由先前命令返回的决定系数，它不是作为结果返回，而是作为估计后结果返回（参见 ereturn list）：

. display e(r2)
.01888498

这对应于方差分析模型所解释的方差比例，可通过以下方式获得（而非通过 oneway显示的平方和进行计算）：

. anova bwt ftv2

观测数= 189 R平方 = 0.0189
均方根误差 = 726.169调整后的R平方= 0.0083

来源	部分平方和	df	MS	F	Prob > F值
模型	1887925.36	2	943962.682	1.79	0.1698
ftv2	1887925.36	2	943962.682	1.79	0.1698
残差	98081730.4	186	527321.131
总计	99969655.8	188	531753.488

2.4.4. 计算特定对比

更一般地，也可以使用 regress 命令来估计甚至检验任何对比。为此，将使用 lincom 命令。以下是一个使用初始模型（婴儿出生体重和母亲种族）的示例：

. regress 出生体重 i.race

来源	SS	df	MS	观测数 = 189
模型	5015725.25	2	2507862.63	F( 2, 186) = 4.91
残差	94953930.6	186	510505.003	Prob> F = 0.0083
总计	99969655.8	188	531753.488	R平方 = 0.0502

| 出生体重 | 系数 | 标准误 | t | P>|t值| | [95% 置信区间] |
| — | — | — | — | — | — |
| 种族 | 2 | -383.0264 | 157.9638 | -2.42 | 0.016 | -694.6575 -71.3954 |
| | 3 | -297.4352 | 113.742 | -2.61 | 0.010 | -521.8254 -73.04498 |
| _cons | 3102.719 | 72.92298 | 42.55 | 0.000 | 2958.856 3246.581 |

可以看出，Stata为变量race的每个level提供一个回归系数，但第一个 level作为参照组除外。因此，y轴截距表示母亲为白人时婴儿的平均体重，而两个回归系数分别表示黑人组和其他组相对于白人组的偏差。黑人组与其他组之间的平均差异可按如下方式估计（Stata根据因子level的数值编码对回归系数进行编号，从1开始）：

. lincom 3.race - 2.race

( 1) - 2.race+ 3.race = 0

| 出生体重 | 系数 | 标准误 | t | P>|t值| | [95% 置信区间] |
| — | — | — | — | — | — |
| (1) | 85.59127 | 165.0887 | 0.52 | 0.605 | -240.0958 411.2783 |

同样，我们知道可以使用命令 ci 来利用正态分布（第1.3节）为平均值构造置信区间，例如：

. lincom _常数项 + 1.race

( 1) 1b.race +_常数项= 0

| 出生体重 | 系数 | 标准误 | t | P>|t| | [95% 置信区间] |
| — | — | — | — | — | — |
| (1) | 3102.719 | 72.92298 | 42.55 | 0.000 | 2958.856 3246.581 |

需要注意的是，不要忘记前面表达式中的y轴截距。

2.4.5. 非参数方法

上述提到的方差分析的非参数替代方法，即克鲁斯卡尔‐沃利斯方差分析，可通过命令 kwallis 获得。其语法为

与 oneway 、 anova 或 regress 略有不同，此次分类因子出现在选项 by() 中。

以下是基于相同出生体重（bwt）和母亲种族（race）数据的秩次方差分析 （ANOVA）结果：

. kwallis 出生体重, 按(种族)

克鲁斯卡尔‐沃利斯总体相等性秩和检验

种族	观测数	秩和
白人	96	10189.00
黑人	26	2015.00
其他	67	5751.00

卡方 = 8.519，自由度为2，概率 = 0.0141
卡方与结点 = 8.520，自由度为2，概率= 0.0141

始终可以通过第2.1.3节中讨论的 ranksum 命令，对解释变量的每个水平进行两两比较来完成此分析。例如，为了从由变量race定义的所有统计单元组中分离出两组，只需通过类似 if race != 3 的筛选条件排除第三组（例如，比较白人和黑人组）。另一种更简便的方法是安装外部命令 kwallis2 （输入 findit kwallis2 并按照安装说明操作）。该命令的语法与 kwallis 相同，但会自动提供与模型相关的所有威尔科克森检验。

2.4.6 双因素方差分析

使用多个分类标准的方差分析通过 anova 命令进行，该命令比单因素方差分析用法更复杂，但允许定义交互项或检验特定对比。以下是一个使用变量 ht（母亲高血压史）和race（母亲种族），并仍以婴儿体重（bwt）作为响应变量的示例。此处将考虑包含交互项的模型，该交互项在Stata中用交互符号##表示。使用notation race##ht时，Stata将考虑两个因素

及其交互作用（使用 oneway 和 anova 时，无需告知Stata变量必须明确表示为分类变量形式）。

方差分析模型的结果如下所示：

. 方差分析 bwt 种族##身高

观测数= 189 R平方 = 0.0768
均方根误差 = 710.143调整后的R平方= 0.0516

来源	部分平方和	df	MS	F	P值 > F
模型	7682087.67	5	1536417.53	3.05	0.0115
种族	2992590.25	2	1496295.12	2.97	0.0539
ht	1757257.26	1	1757257.26	3.48	0.0635
种族#身高	889649.132	2	444824.566	0.88	0.4157
残差	92287568.1	183	504303.651
总计	99969655.8	188	531753.488

为了以顺序的方式计算平方和（如软件R所做的那样），必须添加 sequential选项。

上述模型的一个等效表述，其中明确展现了两个主效应和交互效应，将是：

. 方差分析 出生体重 种族 身高 种族#身高

数值摘要统计可以基于与单个分类因子研究情况下相同的指标（均值、标准差等）来构建。另一方面，命令 tabstat 仅适用于分类变量。此外，可以按如下方式使用命令 table ：

. 表格 种族, 按(身高) 内容(均值 出生体重 标准差 出生体重 计数 出生体重)

身高和	race	均值(出生体重)	出生体重的标准差	样本量(出生体重)
No	白人	3110.89	726.7152	91
	黑人	2751.83	542.1708	23
	其他	2852.41	704.9773	63
Yes	白人	2954	819.4086	5
	黑人	2473.33	1327.633	3
	其他	2063	649.5717	4

or:

. 表格 种族 身高, 内容(均值 出生体重 标准差 出生体重 计数 出生体重)

可以通过外部命令 anovaplot （需要下载并安装， findit anovaplot ）或基于 Stata 最新版本中引入的 margins 命令的图形命令来构建交互作用图。

以下是一种可能的解决方案（图2.6）：

. quietly: margins race#ht 
. marginsplot

2.5. 要点

– 命令 ttest 和 ranksum （或 signrank paired samples ）用于对两个样本（无论是否独立）的数值型响应变量进行比较检验。
– 命令 bitest 和 prtest 可利用二项分布或正态分布，对从一个或两个样本计算出的比例进行统计检验，以检验涉及这些比例的零假设。
– 大多数针对一个或两个样本的检验命令可以直接使用汇总统计量（平均值、标准差、比例）进行计算，而无需使用完整数据：这类命令被称为即时命令。
– tabulate 命令包含两个选项（ chi2 和 exact ），可在列联表情况下提供卡方或费舍尔统计量，而流行病学制表命令如 tabodds 可用于计算比值比，并可考虑分层因素。
– 与方差分析模型相关的命令包括 oneway （单一分类因子情况）和 anova （多个因子的方差分析），它们包含特定选项（如 oneway 中的 bonferroni ），或具有相关命令（如 anova 中的 contrast ），用于进行两两比较、执行特定对比或 summarize 边际效应（例如，对交互图的构建很有帮助）。

2.6. 进阶阅读

Mitchell的著作[MIC 12]详细介绍了 margins 和 marginsplot 命令，并通过大量实例讨论了Stata中的对比分析。关于包含重复测量情况的单向或双向方差分析的深入讨论，可参考Dupont的研究[DUP 09]。

2.7. 应用

1) 在本例中，研究人员测量了10名患者在接受以下两种催眠药物之一治疗前（对照组）和治疗后的睡眠质量：(1) D. 东莨菪碱氢溴酸盐 和 (2) L. 东莨菪碱氢溴酸盐。研究人员采用的评价指标是与基础睡眠时长（对照组）相比的平均睡眠增加时间（小时）[STU 08]。数据如下所示，同时也包含在R软件的基本数据集中（ data(sleep) 数据集）：

D. 东莨菪碱氢溴酸盐：0.7 -1.6 -0.2 -1.2 -0.1 3.4 3.7 0.8 0.0 2.0
L. 东莨菪碱氢溴酸盐：1.9 0.8 1.1 0.1 -0.1 4.4 5.5 1.6 4.6 3.4

研究人员得出结论，只有第二种分子确实具有催眠作用。

– 估计两种分子各自的平均睡眠时间，以及这两个平均值之间的差异；
– 以直方图形式展示差值分数（LHH - DHH）的分布，组距为半小时，并标出这些差值分数的均值和标准差；
– 通过学生t检验验证研究结果的准确性。

Stata中无法获得为学生文章提供基础的睡眠研究数据。不过，我们可以按照第20页所示，使用 input 命令手动输入这些数据。在这种情况下，我们将创建两个变量DHH和LHH，分别对应记录的D. 东莨菪碱氢溴酸盐和L. 东莨菪碱氢溴酸盐分子的睡眠时间：

. input DHH LHH

DHH  LHH 
1. 0.7  1.9 
2. -1.6 0.8  
3. -0.2 1.1 
4. -1.2 0.1 
5. -0.1 -0.1  
6. 3.4 4.4  
7. 3.7 5.5 
8. 0.8 1.6 
9. 0.0 4.6 
10. 2.0 3.4 
11. end

可以通过显示前五个观测值来验证输入是否产生期望的结果：

. list in 1/5

使用 tabstat 可获得按治疗组的算术平均数，否则另一种选项将返回命令中列出变量的平均值：

. tabstat DHH LHH, save

上述提到的 save 选项允许临时存储命令 tabstat 返回的结果，从而可以重用这些结果来计算两个分子之间的平均差异：

. return list 
. matrix m = r(StatTotal) 
. matrix list m 
. display m[1,2] - m[1,1]

正如在R中所见，通过按位置编号调用辅助变量m的元素来操作该变量 （LHH组的平均值位于第二个位置，因此可将其用作m[1,2]）。

差值分数将通过简单的减法计算得出，结果将存储在一个新变量中，如下所示：

. gen sdif=LHH - DHH 
. tabstat sdif, stats(mean sd)

同样，我们可以使用 summarize sdif 。

最后，为了以直方图的形式显示差值分数的分布，并设置0.5小时的组距，需要添加选项 bin(8) （总共八个区间）和 start(0) （从0开始）：

. histogram sdif, percent bin(8) start(0)

要对配对数据执行t检验，我们始终使用命令 ttest ，只是这次使用略有不同的语法：

. t检验 DHH ==LHH

最后，每种分子的睡眠时间平均增益可以通过条形图表示：

. graph hbar DHH LHH, bargap(20)

这里选择以横向条形图（hbar 而不是 bar）的形式表示数据，因为默认情况下 Stata 会自动计算条件均值。若要显示中位数，应使用命令 graph hbar (median) DHH LHH 。

2) 在一项临床试验中，目的是评估一种旨在减少与良性乳腺疾病相关症状数量的治疗方案。将229名患有该疾病的女性随机分为两组。第一组接受常规护理，而第二组患者接受特殊治疗（变量B = treatment）。1年后，对个体进行评估，并将其分类为两类之一：改善或无改善（变量A = response）[SEL 98]。结果在表2.2中进行了汇总，针对部分样本。

	治疗	无治疗	总计
改善	26	21	47
无改善	38	44	82
总计	64	65	129

表2.2. 治疗与乳腺疾病

– 进行卡方检验；
– 在独立性假设下，期望的理论频数是多少？
– 比较使用卡方检验和费舍尔检验得到的结果；
– 为两组患者之间改善比例的差异提供一个置信区间。

Stata 包含称为“即时”的命令，这些命令能够在命令提示符下直接输入数据来计算相关的检验统计量，而无需导入数据文件或输入计数表。在这种情况下，可以使用同一个命令 tabi （注意不要与用于检验特定分布拟合的外部命令 chitesti 混淆）回答上述三个问题。

. tabi 26 21\ 38 44, 精确检验 卡方检验 期望

expected 选项可以提供与观察值频数相对应的理论频数。

注意：Stata 在此类统计检验的基本命令中未提供连续性耶茨校正。

关于两组的改善率，可以按以下方式使用 prtesti ：

. prtesti 64 0.233 65 0.237

3) 在关于雌激素受体基因的研究中，遗传学家将研究兴趣集中在基因型与乳腺癌诊断年龄之间的关系。基因型由序列限制性多态性的两个等位基因 （1.6 kb和0.7 kb）确定，即分为三组受试者：等位基因为0.7 kb的纯合子患者（0.7/0.7）、等位基因为1.6 kb的纯合子患者（1.6/1.6）以及杂合子患者 （1.6/0.7）。数据收集自59名乳腺癌患者，存储于文件 polymorphism.dta （Stata文件）[DUP 09]中。平均数据见表2.3。

	基因型	基因型	基因型	总计
	1.6/1.6	1.6/0.7	0.7/0.7
患者人数	14	29	16	59
诊断时年龄	64.64 11.18 (58.1-71.1)	64.38 13.26 (59.9-68.9)	50.38 10.64 (44.3-56.5)	60.64 13.49
平均值 标准差 95% 置信区间

表2.3. 雌激素受体基因多态性

– 使用方差分析检验诊断时年龄不随基因型变化的零假设。以图形形式表示每种基因型的年龄分布；
– 表4.1中给出的置信区间是在方差齐性的假设下估计的，即使用共同方差的估计值；请在不假设同方差性的情况下给出这些置信区间的值；
– 估计三种基因型所有可能组合之间的平均值差异，并给出相应的95%置信区间的估计以及用于评估观察到的差异显著性水平的参数检验；
– 使用95%置信区间图形化表示各组的平均值。

加载数据并不困难，因为数据已经以Stata格式提供。为了对每种基因型的年龄变量提供数值摘要，可以结合使用 by 和 summarize 命令。稍后将看到， tabstat 命令可以返回简化的结果，例如仅均值和标准差：

. use polymorphism.dta 
. 按基因型: summarize 年龄

上述报告的置信区间并非来自下方进行的方差分析，而是基于正态分布近似计算得出的。然而，通过添加 option detail 可以获得更详细的描述性摘要：

. by genotype: summarize age, detail

可以通过箱线图显示每种基因型的年龄分布：

. graph box age, over(genotype)

当使用直方图表示年龄分布时，可以使用命令 histogram ，如下所示：

. 直方图 年龄, 按基因型分组（3列）

选项 by(genotype, col (3)) 允许按基因型分别显示分布，并使直方图水平对齐（即使用三列）。

oneway 命令可以用于进行固定效应方差分析，具体如下：

单因素方差分析 年龄 基因型

Stata生成的方差分析表与R生成的方差分析表基本相同，只是Stata还返回了平方和、均方以及自由度的总计值。还应注意，Stata提供了关于方差齐性假设的巴特利特检验结果。可通过一种略有不同的方法获得用于检验方差齐性的勒文检验：使用命令 robvar ，该命令提供了一系列方差相等性的稳健检验：

. robvar 年龄, 按基因型

勒文检验在检验统计量W0下呈现。

最后，通过添加选项 tabulate ，可以部分重现上一步（条件均值）得到的数值摘要（不含置信区间）：

. 单因素方差分析 年龄 基因型, tabulate

为了基于方差分析的结果计算置信区间，即通过共同方差来估计残差误差，我们可以参照R中的方法，利用误差均方（组内）以及分布 t(0.975)的参考分位数，在Stata中可通过以下方式获得（此处以第一种基因型1.6/1.6为例）：

. display invttail(14-3, 0.025)

更一般地，我们可以根据汇总数据生成三个组均值的95%置信区间的上限和下限，如下例所示（在这种情况下，不再考虑共同方差的估计）：

. collapse (mean) agem=age (sd) ages=age (count)n=age, by(genotype) 
. generate agelci = agem - invttail(n-1, 0.025)*(ages/sqrt(n)) 
. generate ageuci = agem+invttail(n-1, 0.025)*(ages/sqrt(n))

而使用类似 twoway (bar agem genotype) (rcap ageuci agelci genotype) 的命令可以将结果以图形方式显示。另一种更便捷的替代方案是使用 serrbar 命令，该命令可用于显示一系列均值及其对应的标准误。由于我们将使用置信区间而非标准误，因此只需对标准用法稍作修改即可。通过使用之前的计算结果，可以计算出置信区间的半宽度（已知其围绕平均值对称）：

. gen ageb=agem-agelci 
. serrbar 年龄均值 ageb 基因型，xlabel(1"1.6/1.6"2"1.6/0.7" 3 "0.7/0.7")

需要注意的是，可以通过基于回归模型构建特定对比来获得均值差异，该回归模型的结果等同于方差分析模型，前提是通过操作符 i.* 将分类变量基因型以指示变量矩阵的形式进行编码。我们不展示调用命令 regress 的结果，因为我们关心的只是利用它保存的结果（参见 e() ）：

regress age i.genotype（以基因型为自变量的回归命令）

对于基因型0.7/0.7和1.6/0.7之间平均值的差异，例如，应使用以下命令：

lincom 3.genotype - 2.genotype（比较第3组与第2组基因型均值差异的命令）

对于另外两个平均差异，0.7/0.7 - 1.6/1.6 和 1.6/0.7 - 1.6/1.6，也将采用相同的步骤。

事实上，同样的方法（通过回归）可以让我们获得每组平均值的95%置信区间：

. lincom _常数项 + 1.基因型

4) 任何产科服务都关注足月及1个月大婴儿的体重[PEA 05]。关于这个 550名婴儿的样本，还有有关胎次（兄弟姐妹数量）的信息，但已知有兄弟姐妹的儿童之间不存在双胞胎关系。本研究的目的是确定胎次（四类）是否影响 1个月大新生儿的体重。数据汇总于表 2.4 中，并以名为 weights.sav 的 SPSS文件提供。

– 验证上一个表格中显示的数据；
– 进行单因素方差分析。对整体显著性得出结论，并说明模型解释的方差部分；
– 按性别显示体重分布。进行方差齐性检验（在在线帮助中搜索Levene检验）；
– 决定将最后两个类别合并（2 和 ≥ 3）。重新进行分析，并与第一次方差分析的结果进行比较；
– 使用重新编码为三个胎次level的数据进行线性趋势检验（采用方差分析）。

兄弟姐妹数量	0	1	2	≥ 3	总计
样本 计数 比例	180 32.7	192 34.9	116 21.1	62 11.3	550 100.0
体重（千克） Mean 标准差 （最小值-最大值）	4.26 0.62 (2.92-5.75)	4.39 0.59 (3.17-6.33)	4.21 0.61 (3.09-6.49)	4.11 0.54 (3.20-5.48)

表 2.4. 新生儿体重

除了在Windows上操作外（参见 usespss 命令），目前没有非常方便的方法可以直接将SPSS文件导入Stata。包含数据的文件weights.sav已通过以下命令从R以Stata格式导出：

library(foreign)
weights <- read.spss("weights.sav", to.data.frame=TRUE) 
write.dta(weights, file="weights.dta")

在 Stata 中，可以通过插入 use 命令非常简单地导入。

. use "weights.dta", clear 
. list in 1/5

使用命令 tabulate 可获得变量 PARITY 的频数和相对频率表格：

. tabulate 胎次

根据兄弟姐妹数量得出的体重均值和标准差如下：

. tabstat 体重, stats(均值 标准差) by(PARITY) format(%9.2f)

选项 format (%9.2f) 将显示限制为两位小数。

单因素方差分析可通过 oneway 命令实现，方法与之前的练习相同，即提供响应变量和描述待比较组的分类变量：

. 单因素方差分析 体重 胎次

还可以根据之前报告的平方和来计算解释方差的比例。也可以使用命令 anova ，该命令除了返回方差分析表外，还会返回与模型相关的决定系数。

. 方差分析 体重 胎次

为了以直方图形式显示每组的数据，必须使用 histogram 命令并结合选项 by 来定义分类因子。选项 freq 可以显示频数而非比例（或密度）：

. 直方图 体重, by(PARITY) freq

要显示散点图，像在R中一样，我们可以输入一个外部命令，例如 ssc 安装 stripplot ，具体如下：

. stripplot 体重, over(PARITY) stack height(.4) center vertical width(.3)

或者更简单地：

. scatter 体重 PARITY, jitter(3) xlabel(1 "独生子女" 2 "一个兄弟姐妹" 3 "两个兄弟姐妹" 4 "三个及以上")

默认情况下， oneway 命令会显示Bartlett检验的结果，用于比较各组之间的方差。如果需要进行勒文检验，则必须使用 robvar 命令，该命令将给出勒文检验的结果（ W0）：

. robvar 体重, by(PARITY)

在Stata中有多种方法可以重新编码变量，但在本例中，将最后两类（两个兄弟姐妹和三个及以上）合并的最简单方法是生成第二个变量 PARITY2 ，如下所示：

. 重新编码 PARITY (1=1) (2=2) (3/4=3)，生成(PARITY2)

为了验证转换已成功进行，将显示一个简单的列联表，交叉显示两个变量的频数：

. tabulate PARITY PARITY2

以下是关于新创建变量的单因素方差分析结果：

. oneway WEIGHT PARITY2

对于线性趋势检验，提供了两种方法（对比法和线性回归）。为了生成检验线性趋势的对比，我们必须明确地指示

Stata 通过将变量 PARITY2 视为分类变量（涉及对最后两个因子水平的虚拟变量编码）来执行回归。在这种情况下，将使用估计后命令 contrast 。由于我们并不真正关心分类变量回归的结果，因此将通过 quietly: 指令告知 Stata 不显示结果。

quietly: 回归 体重 i.PARITY2 
. 对比 p.PARITY2, 无效应

感兴趣的对比出现在标记为（线性）的行上。

关于简单线性回归方法，命令更为简单：

回归 体重 胎次2

趋势检验对应于回归线斜率的检验，即此处的系数 PARITY2 。

5) R 中可用的数据集 ToothGrowth 包含了一项研究的数据，该研究记录了在给予不同剂量的维生素C（0.5、1 或 2 毫克，变量 dose ）以抗坏血酸或橙汁形式（变量 supp ）后，10 只豚鼠的成牙本质细胞长度（变量 len ）[BLI 52]。数据可在文件 ToothGrowth.csv 中获取

– 根据该实验设计的不同处理条件（对两个因素 supp 和 dose 的模态进行交叉列表）验证频数分布；
– 计算每种治疗的均值和标准差；
– 为包含交互作用的完整模型构建方差分析表格 b两个因素之间；
– 绘制交互图，表示响应变量在两个解释变量各水平下的平均值；
– 验证这些数据的方差齐性假设是否成立。

由于这是一个标准的CSV文件，因此使用命令 import delimited 来加载数据。假设文件位于工作目录中，我们将输入：

. import delimited "ToothGrowth.csv"

需要注意的是，变量 supp 被识别为字符串。通过 tabulate 可获得按不同过程划分的统计单元的分布，但也可以使用 table 在频数统计的同时计算平均值和标准差，如下所示

. 表格 supp 剂量, 内容(长度均值 长度标准差 长度样本量) 格式(%5.1f)

指令 format(%5.1f) 允许我们将输出限制为一位小数。双因素方差分析通过 anova 进行，但需要先将变量 supp 转换为分类变量。为此，使用编码 （ encode ）命令；这可以将字符串转换为分类变量的 level（注意，需要创建一个新变量）：

. 编码 supp, 生成(suppc)

变量 dose 被作为数值变量进行处理，但在这里我们希望也将其转换为分类变量。这可以通过两个步骤实现：首先，将该变量转换为字符串，然后像对 supp 所做的那样，将其转化为分类变量。

. 转为字符串 dose, 生成(dose_) 
. 编码 dose_, 生成(dosec) 
. 删除 dose_

方差分析模型可写为：

. 方差分析 长度 suppc与dosec的交互项

交互符号 ## 表示两个因素之间的交互作用。需要注意的是，仍然可以将变量 dose 作为数值变量进行处理（在方差分析表中占1个自由度），但在这种情况下必须在其前添加连续变量前缀 c. （ anova len suppc##c.dose ）。

要构造交互作用图，我们可以输入命令 scatter ，但命令 marginsplot 在使用上也十分灵活。只需先通过 margins 计算条件平均值，如下所示：

. margins 剂量c# suppc 
. marginsplot