第五讲

Size: px

Start display at page:

Download "第五讲"

缠镇祁
5 years ago
Views:

1 Stata 软件基本操作和数据分析入门第五讲多组平均水平的比较赵耐青一复习和补充两组比较的统计检验 1. 配对设计资料 ( 又称为 Dependent Samples) a) 对于小样本的情况下, 如果配对的差值资料服从正态分布, 用配对 t 检验 (ttest 差值变量 =0) b) 大样本的情况下, 可以用配对 t 检验 c) 小样本的情况下, 并且配对差值呈偏态分布, 则用配对符号秩检验 (signrank 差值变量 =0) 2. 成组设计 (Two Independent Samples) a) 如果方差齐性并且大样本情况下, 可以用成组 t 检验 (ttest 效应指标变量,by( 分组变量 )) b) 如果方差齐性并且两组资料分别呈正态分布, 可以用成组 t 检验 c) 如果方差不齐, 或者小样本情况下偏态分布, 则用秩和检验 (Ranksum test) group x

2 二多组比较 1. 完全随机分组设计 ( 要求各组资料之间相互独立 ) a) 方差齐性并且独立以及每一组资料都服从正态分布 ( 小样本时要求 ), 则采用完全随机设计的方差分析方法 ( 即 : 单因素方差分析,One Way ANOVA) 进行分析 b) 方差不齐或小样本情况下资料偏态, 则用 Kruskal Wallis 检验 (H 检验 ) 例 5.1 为研究胃癌与胃粘膜细胞中 DNA 含量 (A.U) 的关系, 某医师测得数据如下, 试问四组人群的胃粘膜细胞中平均 DNA 含量是否相同? 组别 group DNA 含量 (A.U) 浅表型胃炎肠化生早期胃癌

3 晚期胃癌由于这四组对象的资料是相互独立的, 因此属于完全随机分组类型的检验问题是考察四组 DNA 含量的平均水平相同吗如果每一组资料都正态分布并且方差齐性可以用 One way-anova 进行分析, 反之用 Kruskal Wallis 检验 STATA 数据输入格式 g x

4 分组正态性检验,α=0.05. sktest x if g== Skewness/Kurtosis tests for Normality joint Variable Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi x sktest x if g==2 Skewness/Kurtosis tests for Normality joint Variable Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi x sktest x if g==3 Skewness/Kurtosis tests for Normality joint Variable Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi x sktest x if g==4 Skewness/Kurtosis tests for Normality joint Variable Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi x 上述结果表明每一组资料都服从正态分布单因素方差分析的 STATA 命令 :oneway 效应指标变量分组变量,t b 其中 t 表示计算每一组均数和标准差,b 表示采用 Bonferroni 统计方法进行两

5 两比较本例命令为 oneway x group,t b. oneway x g,t b Summary of x g Mean Std. Dev. Freq Total Analysis of Variance Source SS df MS F Prob > F Between groups Within groups Total Bartlett's test for equal variances: chi2(3) = Prob>chi2 = 方差齐性的检验为 : 卡方 =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 ( 第 2 组样本均数 - 第 1 组样本均数 ) 0.000(H0:μ1=μ2 检验的 P 值 ) ( 第 3 组样本均数 - 第 2 组样本均数 ) (H0:μ3=μ2 检验的 P 值 ) ( 第 4 组样本均数 - 第 3 组样本均数 ) (H0:μ3=μ4 检验的 P 值 )\ 上述输出为两两比较的结果, 在表格的每个单元中, 第一行为两组均数的差值, 第二行为两组均数比较检验的 P 值根据上述结果可以知道, 第 2 组第 3 组和第 4 组的 AU 均数均大于第 1 组的 AU 均数, 并且差别有统计学意义说明肠化生患者和胃癌患者的 DNA 的 AU 含量平均水平均高于正常人的 AU 平均水平, 并且差别有统计学意义第 3 组和第 4 组的 AU 均数也大于第 2 组的 AU 平均水平, 并且差别有统计学意义说明胃

6 癌患者的 DNA 的 AU 含量平均水平均高于肠化生患者的 AU 平均水平, 并且差别有统计学意义第 3 组和第 4 组两组均数的差别没有统计学意义, 说明没有足够的证据可以 DNA 的 AU 含量与癌症的早期与晚期有关系假如本例的资料不满足方差分析的要求, 则用 Kruskal Wallis 检验, 数据结构同上命令为 : kwallis 效应指标变量, by( 分组变量 ) 本例的命令为 kwallis x,by(g) H 0 :4 组的 AU 总体分布相同 H 1 :4 组的 AU 总体分布不全相同 α=0.05 结果如下 : Test: Equality of populations (Kruskal-Wallis test) g _Obs _RankSum chi-squared = with 3 d.f. probability = chi-squared with ties = with 3 d.f. probability = 说明 :4 组 AU 的总体分布不全相同, 然后秩和检验, 但 α 应取小一些 ( 多重比较时, 会增大第一类错误的概率 ) 根据 Sidak 检验的建议 : α 1 =, 其中 k 为要比较的次数,α 1 (1 α) k 为多组比较总的检验水平 ( 一般为 0.05),α 为两两比较时的检验水平 2 如本例 :4 组两两比较共比 C = 次, 因此 a = 1 (0.95) 6 = , 对于比较第 1 组和第 2 组的 AU 分布差别的操作命令为 : 先计算中位数 sort g 组别变量排序 by g:centile x,centile(50) 计算各组中位数 -> g = 1 -- Binom. Interp. -- Variable Obs Percentile Centile [95% Conf. Interval] x > g = 2 -- Binom. Interp. -- Variable Obs Percentile Centile [95% Conf. Interval] x > g = 3 -- Binom. Interp. -- Variable Obs Percentile Centile [95% Conf. Interval]

7 x > g = 4 -- Binom. Interp. -- Variable Obs Percentile Centile [95% Conf. Interval] x 得到这 4 组中位数分别为 :M 1 =12.29,M 2 =14.855,M 3 =21.14 和 M 4 =21.36 ranksum x if g==1 g==2,by(g) Two-sample Wilcoxon rank-sum (Mann-Whitney) test g obs rank sum expected combined unadjusted variance adjustment for ties adjusted variance Ho: x(g==1) = x(g==2) z = Prob > z = P 值 <α, 因此第 2 组 AU 的平均水平要高于第 1 组的平均水平 (M 2 >M 1 ), 并且差别有统计学意义第 1 组与第 3 组比较 ranksum x if g==1 g==3,by(g) Two-sample Wilcoxon rank-sum (Mann-Whitney) test g obs rank sum expected combined unadjusted variance adjustment for ties adjusted variance Ho: x(g==1) = x(g==3)

8 z = Prob > z = P 值 <α, 因此第 3 组 AU 的平均水平要高于第 1 组的平均水平 (M 3 >M 1 ), 并且差别有统计学意义, 其他比较类似进行要注意的问题 : 在方差分析中, 要求每一组资料服从正态分布 ( 小样本时 ), 并不是要求各组资料服从一个正态分布 ( 因为这就意味各组的总体均数相同, 失去统计检验的必要性 ), 所以不能把各组的资料合在一起作正态性检验总的讲, 方差分析对正态性具有稳健性, 即 : 偏态分布对方差分析的结果影响不会太大, 故正态性检验的 α 取 0.05 也就可以了样本量较大时, 方差分析对正态性要求大大降低 ( 根据中心极限定理可知 : 样本均数近似服从正态分布 ) 并且由于大多数情况下, 样本资料只是近似服从正态分布而不是完全服从正态分布由于在大样本情况下, 用正态性检验就变为很敏感, 对于不是完全服从正态分布的资料往往会拒绝正态性检验的 H 0 : 资料服从正态分布因为正态性检验不能检验资料是否近似服从正态分布, 而是检验是否服从正态分布故在大样本情况下, 考察资料的近似正态性, 应用频数图进行考察方差齐性问题对方差分析相对比较敏感, 并且并不是随着样本量增大而方差齐性对方差分析减少影响的但是当各组样本量接近相同或相同时, 方差齐性对方差分析呈现某种稳健性即 : 只有当各组样本量相同时, 方差齐性对方差分析结果的影响大大降低这时随着样本量增大, 影响会进一步降低相反, 如果各组样本量相差太大时, 方差齐性对方差分析结果的影响很大这时随着样本量增大, 影响会进一步加大 2. 随机区组设计 ( 处理组之间可能不独立 ) a) 残差 ( 定义为 : eij = Xij + X Xi. X. j, 也就是随机区组方差分析中的误差项 ) 的方差齐性且小样本时正态分布, 则用随机区组的方差分析 ( 无重复的两因素方差分析,Two-way ANOVA) b) 不满足方差齐性或小样本时资料偏态, 则对用秩变换后再用随机区组的方差分析也可以直接用非参数随机区组的秩和检验 Fredman test) 例 2 下表是某湖水中 8 个观察地点不同季节取样的氯化物含量测定值, 请问在不同季节该湖水中氯化物的含量有无差别? 表 2 某湖水中不同季节的氯化物含量测定值 (mg/l) location no 春夏秋冬显然同一地点不同季节的氯化物含量有一定的相关性, 故不能采用完全随机设计的方差分析方法对 4 个季节的氯化物含量进行统计分析可以把同一地点的 4 个季节氯化物含量视为一个区组, 因此可以用随机区组的方差分析进行统计分析

9 设第 8 个地点在冬季的氯化物总体均数为 μ 0, 同样在冬季, 第 i 个地点的氯化物总体均数与第 8 个地点在冬季的氯化物总体均数相差 β i,i=1,2,3,4,5,6,7 因此在冬季的这 8 个地点在冬季的氯化物总体均数可以表示为地点编号冬季氯化物均数 μ 0 +β 1 μ 0 +β 2 μ 0 +β 3 μ 0 +β 4 μ 0 +β 5 μ 0 +β 6 μ 0 +β 7 μ 0 假定在同一地区, 春季的氯化物总体均数与冬季的氯化物总体均数相差 α 1, 因此春节和冬季的氯化物总体均数可以表示为地点编号冬季氯化物均数 μ 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, 则四个季节的氯化物总体均数可以表示为地点编号冬季氯化物均数 μ 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 根据上述总体均数表示, 可以知道 : 在四个季节中的氯化物总体均数 ( 同一地点 ) 无变化就是 H 0 :α 1 =α 2 =α 3 =0( 在随机区组方差分析中称为无处理效应, 但不能称 4 组的总体均数相同, 因为在同一季节中不同地点的总体均数可能不同 ) H 1 :α 1,α 2,α 3 不全为 0 Stata 数据输入格式 t id x

10 其中 id 表示观察地点编号,t=1,2,3,4 对应表示春节夏季秋季和冬季 Stata 操作命令 : anova x t id. anova x t id Number of obs = 32 R-squared = Root MSE = Adj R-squared = Source Partial SS df MS F Prob > F Model t id Residual Total MS处理处理效应 H 0 :α 1 =α 2 =α 3 =0 的检验对应的统计量 F = = = MS 相应的 P 值 <0.0001( 计算机输出值是 ), 所以拒绝无效假设, 可以认为 4 个季节的氯化物总体均数不全相同不同季节中的两两比较用 LSD 方法检验如下 : 在输入 anova x t id 命令后, 再输入 regress 命令便得到下列结果误差

11 Source SS df MS Number of obs = F( 10, 21) = Model Prob > F = Residual R-squared = Adj R-squared = Total Root MSE =

12 x Coef. Std. Err. t P> t [95% Conf. Interval] _cons (μ0) t α1= α2= α3= id 4 (dropped) β1= β2= β3= β4= β5= β6= β7= (dropped) 其中 α 1 = 6.081, 对应的假设检验 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 = (2.758,4.874), 因此可以认为夏季的氯化物平均高于冬季, 差别有统计学意义, 对应的假设检验 H0:α3=0 的统计量 t=2.37,p 值 =0.027,95% 可信区间为 α 3 = (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) = Prob > F = 相应的统计量 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]] 此处不在详细叙述了由于随机区组方差分析要求残差 ( eij = Xij + X Xi. X. j ) 服从正态分布, 再输入 regress 以后, 只要输入 predict 残差变量名,residual, 就可以得到残差计算值本例用 e 表示残差变量名, 因此输入 predict e,residual

13 就可以得到残差计算值 e, 然后对残差进行正态性检验 (sktest 残差变量名 ) 本例输入命令为 : sktest e 结果如下 (H 0 : 残差服从正态分布 vs H 1 : 残差偏态分布, α=0.05) Skewness/Kurtosis tests for Normality joint Variable Pr(Skewness) Pr(Kurtosis) adj chi2(2) Prob>chi e P 值 = >>α, 因此可以认为资料近似服从正态分布 ( 大样本时, 可以不考虑正态性问题 ) 如果资料呈偏态分布, 可以对资料进行秩变换 (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 = Root MSE = Adj R-squared = Source Partial SS df MS F Prob > F Model t id Residual Total

14 命令 regress regress 结果如下 Source SS df MS Number of obs = F( 10, 21) = Model Prob > F = Residual R-squared = Adj R-squared = Total Root MSE = r Coef. Std. Err. t P> t [95% Conf. Interval] _cons t (dropped) id (dropped) 进一步两两比较 test _b[t[1]]=_b[t[2]] ( 1) t[1] - t[2] = 0.0 F( 1, 21) = Prob > F = test _b[t[1]]=_b[t[3]] ( 1) t[1] - t[3] = 0.0 F( 1, 21) = Prob > F = test _b[t[2]]=_b[t[3]] ( 1) t[2] - t[3] = 0.0 F( 1, 21) = Prob > F = 解释如同上述, 不再重复

第一讲 Stata操作入门

第一讲 Stata操作入门 Stata 软件基本操作和数据分析入门赵耐青张文彤第一讲 Stata 操作入门第一节概况 Stata 最初由美国计算机资源中心 (Computer Resource Center) 研制, 现在为 Stata 公司的产品, 其最新版本为 7.0 版它操作灵活简单易学易用, 是一个非常有特色的统计分析软件, 现在已越来越受到人们的重视和欢迎, 并且和 SAS SPSS 一起, 被称为新的三大权威统计软件