第五講多組平均水平的比較

1、Stata基本操作和數(shù)據(jù)分析入門第五講多組平均水平的比較一、復(fù)習(xí)和補(bǔ)充兩組比較的統(tǒng)計(jì)檢驗(yàn)1. 配對(duì)設(shè)計(jì)資料(又稱為 Dependent Samples)a)對(duì)于小樣本的情況下，如果配對(duì)的差值資料服從正態(tài)分布，用配對(duì)t 檢驗(yàn)(ttest差值變量=0)b)大樣本的情況下，可以用配對(duì) t 檢驗(yàn)c)小樣本的情況下，并且配對(duì)差值呈偏態(tài)分布，則用配對(duì)符號(hào)秩檢驗(yàn)(signrank 差值變量=0)2. 成組設(shè)計(jì)(Two Independent Samples)a)如果方差齊性并且大樣本情況下，可以用成組 t 檢驗(yàn)(ttest效應(yīng)指標(biāo)變量,by(分組變量)b)如果方差齊性并且兩組資料分別呈正態(tài)分布，可以用成組

2、t 檢驗(yàn)c)如果方差不齊，或者小樣本情況下偏態(tài)分布，則用秩和檢驗(yàn)(Ranksum test)51groupx07909309109209407709307409101015208307308801020900100081091083010608407808709501011101110011141107981110121194981110二、多組比較1. 完全隨機(jī)分組設(shè)計(jì)(要求各組資料之間相互)a)方差齊性并且以及每一組資料都服從正態(tài)分布(小樣本時(shí)要求)，則采用完全隨機(jī)設(shè)計(jì)的方差分析(即：?jiǎn)我蛩胤讲罘治?，One Way ANOVA)進(jìn)行分析。b)方差不齊或小樣本情況下資料偏態(tài)，則用 Krusk

3、al Wallis 檢驗(yàn)(H檢驗(yàn))例5.1 為研究胃癌與胃粘膜細(xì)胞中DNA含量（A.U）的，某醫(yī)師測(cè)得數(shù)據(jù)如下，試問四組人群的胃粘膜細(xì)胞中平均DNA含量是否相同？由于這四組對(duì)象的資料是相互的，因此屬于完全隨機(jī)分組類型四組 DNA 含量的平均水平相同嗎。如果每一組的。檢驗(yàn)是資料都正態(tài)分布并且方差齊性可以用 One way-ANOVA 進(jìn)行分析，反之用Kruskal Wallis 檢驗(yàn)。STATA 數(shù)據(jù)輸入格式53gx19.81112.73112.29112.53112.95組別groupDNA 含量（A.U）淺表型胃炎19.8112.7312.2912.5312.959.5312.68.912.

4、2714.2610.68腸化生214.6117.5415.11713.3915.3213.7418.2413.8112.6314.5316.17早期胃癌323.2620.820.623.517.8521.9122.1322.0419.5318.4121.4820.24晚期胃癌423.7319.4622.3919.5325.920.4320.7120.0523.4121.3421.3825.705419.53112.618.9112.27114.26110.68214.61217.54215.1217213.39215.32213.74218.24213.81212.63214.53216.17

5、323.26320.8320.6323.5317.85321.91322.13322.04319.53318.41321.48320.24423.73419.46422.39419.53425.9分組正態(tài)性檢驗(yàn)，a=0.05上述結(jié)果表明每一組資料都服從正態(tài)分布。55. sktest xif g=1Skewness/Kurtosis tests for Normalityjoint Variable |Pr(Skewness)Pr(Kurtosis) adj chi2(2)Prob>chi2+x |0.4910.4851.070.5861. sktest xif g=2Skewness/K

6、urtosis tests for Normalityjoint Variable |Pr(Skewness)Pr(Kurtosis) adj chi2(2)Prob>chi2+x |0.4820.5410.960.6201. sktest xif g=3Skewness/Kurtosis tests for Normalityjoint Variable |Pr(Skewness)Pr(Kurtosis) adj chi2(2)Prob>chi2+x |0.5270.7500.520.7704. sktest xif g=4Skewness/Kurtosis tests for

7、Normalityjoint Variable |Pr(Skewness)Pr(Kurtosis) adj chi2(2)Prob>chi2+x |0.2600.6161.750.4166420.43420.71420.05423.41421.34421.38425.7單因素方差分析的 STATA 命令：oneway效應(yīng)指標(biāo)變量 分組變量,tb其中t 表示計(jì)算每一組均數(shù)和標(biāo)準(zhǔn)差,b 表示采用Bonferroni 統(tǒng)計(jì)方法進(jìn)行兩兩比較。本例命令為 oneway x group,t b.onewayxg,tb方差齊性的檢驗(yàn)為：卡方=1.1354，度=3,P 值=0.769，因此可以認(rèn)為方差是

8、齊性的。H0：m1=m2=m3=m4 四組總體均數(shù)相同H1：m1，m2，m3，m4 不全相同a=0.05，相應(yīng)的統(tǒng)計(jì)量 F=77.87 以及相應(yīng)的度為 3 和 43，P 值<0.0001，因此 4 組均數(shù)的差別有統(tǒng)計(jì)學(xué)意義。56Comparison of x by g (Bonferroni)Row Mean-|Col Mean |123|Summary of xg |MeanStd. Dev.Freq.+1 |11.6863641.6884388112 |15.1733331.749173123 |20.9791671.7668279124 |22.00252.242908712+To

9、tal |17.5831914.608078947Analysis of VarianceSourceSSdfMSFProb > FBetween groups824.9425493274.9808577.870.0000Within groups151.839445433.53114987Total976.7819944621.2343912Bartlett's test for equal variances:chi2(3) =1.1354Prob>chi2 = 0.769上述輸出為兩兩比較的結(jié)果，在表格的每個(gè)單元中，第一行為兩組均數(shù)的差值，第二行為兩組均數(shù)比較檢驗(yàn)的

10、P 值。根據(jù)上述結(jié)果可以知道，第 2 組、第 3 組和第 4 組的 AU 均數(shù)均大于第 1 組的 AU 均數(shù)，并且差別有統(tǒng)計(jì)學(xué)意義。說明腸化生患者和胃癌患者的 DNA 的 AU 含量平均水平均高于正常人的 AU 平均水平，并且差別有統(tǒng)計(jì)學(xué)意義。第 3 組和第 4 組的 AU 均數(shù)也大于第 2 組的AU 平均水平，并且差別有統(tǒng)計(jì)學(xué)意義。說明胃癌患者的 DNA 的 AU 含量平均水平均高于腸化生患者的 AU 平均水平，并且差別有統(tǒng)計(jì)學(xué)意義。第 3 組和第 4 組兩組均數(shù)的差別沒有統(tǒng)計(jì)學(xué)意義，說明沒有足夠的證據(jù)可以 DNA 的AU 含量與的早期與晚期有。假如本例的資料不滿足方差分析的要求，則用 Kr

11、uskal Wallis 檢驗(yàn)，數(shù)據(jù)結(jié)構(gòu)同上。命令為：kwallis 效應(yīng)指標(biāo)變量, by(分組變量)本例令為 kwallisx,by(g)H0：4 組的AU 總體分布相同H1：4 組的AU 總體分布不全相同57+2 |3.48697(第 2 組樣本均數(shù)第 1 組樣本均數(shù))|0.000(H0:m1=m2 檢驗(yàn)的 P 值)|3 |9.29285.80583(第 3 組樣本均數(shù)第 2 組樣本均數(shù))|0.0000.000(H0:m3=m2 檢驗(yàn)的 P 值)|4 |10.31616.829171.02333(第 4 組樣本均數(shù)-第 3 組樣本均數(shù))|0.0000.0001.000(H0:m3=m4 檢

12、驗(yàn)的 P 值)a=0.05結(jié)果如下：說明：4 組 AU 的總體分布不全相同，然后秩和檢驗(yàn)，但a應(yīng)取小一些(多重比較時(shí)，會(huì)增大第一類錯(cuò)誤的概率)。根據(jù) Sidak 檢驗(yàn)的建議：a為多組比較總的檢驗(yàn)水平(一1a¢ = 1- (1-a ) ，其中 k 為要比較的次數(shù)，k般為 0.05)，a為兩兩比較時(shí)的檢驗(yàn)水平。1如本例：4 組兩兩比較共比C = 6 次，因此a¢ = 1- (0.95) = 0.0085 ，264對(duì)于比較第 1 組和第 2 組的 AU 分布差別的操作命令為：先計(jì)算中位數(shù)sort g組別變量排序by g:centilex,centile(50)計(jì)算各組中位數(shù)58-

13、> g = 1- Binom. Interp. - Variable |Obs PercentileCentile95% Conf. Interval+x |115012.299.72956412.7932-> g = 2- Binom. Interp. - Variable |Obs PercentileCentile95% Conf. Interval+x |125014.85513.7474516.91172-> g = 3- Binom. Interp. - Variable |Obs PercentileCentile95% Conf. Interval+Test:

14、 Equality of populations (Kruskal-Wallis test) g_Obs_RankSum11172.00212205.00312411.50412439.50chi-squared =37.814 with 3 bability =0.0001chi-squared with ties =37.816 with 3 d.f. probability =0.0001得到這 4 組中位數(shù)分別為：M1=12.29，M2=14.855，M3=21.14和M4=21.36ranksumxif g=1 | g=2,by(g)P 值<a，因此第 2 組 A

15、U 的平均水平要高于第1組的平均水平(M2>M1)，并且差別有統(tǒng)計(jì)學(xué)意義。第 1 組與第 3 組比較ranksum x if g=1| g=3,by(g)59Two-sample Wilcoxon rank-sum (Mann-Whitney) test g |obsrank sumexpected+1 |11661323 |12210144+combined |23276276unadjusted variance264.00adjustment for ties0.00Two-sample Wilcoxon rank-sum (Mann-Whitney) test g |obsran

16、k sumexpected+1 |11721322 |12204144+combined |23276276unadjusted variance264.00adjustment for ties0.00adjusted variance264.00Ho: x(g=1) = x(g=2)z = -3.693Prob > |z| =0.0002x |125021.1419.6055222.12043-> g = 4- Binom. Interp. - Variable |Obs PercentileCentile95% Conf. Interval+x |125021.3620.09

17、04223.69596P 值<a，因此第 3 組 AU 的平均水平要高于第 1 組的平均水平(M3>M1)，并且差別有統(tǒng)計(jì)學(xué)意義，其他比較類似進(jìn)行。要注意的：u在方差分析中，要求每一組資料服從正態(tài)分布(小樣本時(shí))，并不是要求各組資料服從一個(gè)正態(tài)分布(因?yàn)檫@就意味各組的總體均數(shù)相同，失去統(tǒng)計(jì)檢驗(yàn)的必要性)，所以不能把各組的資料合在一起作正態(tài)性檢驗(yàn)?？偟闹v，方差分析對(duì)正態(tài)性具有穩(wěn)健性，即：偏太大，故正態(tài)性檢驗(yàn)的a取 0.05態(tài)分布對(duì)方差分析的結(jié)果影響也就可以了。u樣本量較大時(shí)，方差分析對(duì)正態(tài)性要求大大降低(根據(jù)中心極限定理可知：樣本均數(shù)近似服從正態(tài)分布)。并且由于大多數(shù)情況下，樣本資料

18、只是近似服從正態(tài)分布而不是完全服從正態(tài)分布。由于在大樣本情況下，用正態(tài)性檢驗(yàn)就變?yōu)楹苊舾校瑢?duì)于不是完全服從正態(tài)分布的資料往往會(huì)拒絕正態(tài)性檢驗(yàn)的 H0：資料服從正態(tài)分布。因?yàn)檎龖B(tài)性檢驗(yàn)不能檢驗(yàn)資料是否近似服從正態(tài)分布，而是檢驗(yàn)是否服從正態(tài)分布。故在大樣本情況下，資料的近似正態(tài)性，應(yīng)用頻數(shù)圖進(jìn)行。u方差齊性對(duì)方差分析相對(duì)比較敏感，并且并不是隨著樣本量增大而方差齊性對(duì)方差分析減少影響的。但是當(dāng)各組樣本量接近相同或相同時(shí)，方差齊性對(duì)方差分析呈現(xiàn)某種穩(wěn)健性。即：只有60adjusted variance264.00Ho: x(g=1) = x(g=3)z = -4.062Prob > |z| =0

19、.0000當(dāng)各組樣本量相同時(shí)，方差齊性對(duì)方差分析結(jié)果的影響大大降低。這時(shí)隨著樣本量增大，影響會(huì)進(jìn)一步降低。相反，如果各組樣本量相差太大時(shí)，方差齊性對(duì)方差分析結(jié)果的影響很大。這時(shí)隨著樣本量增大，影響會(huì)進(jìn)一步加大。2. 隨機(jī)設(shè)計(jì)(處理組之間可能不)a)殘差(定義為： eij =- X j ，也就是隨機(jī)方差分析中i的誤差項(xiàng))的方差齊性且小樣本時(shí)正態(tài)分布，則用隨機(jī)的方差分析(無(wú)重復(fù)的兩因素方差分析,Two-way ANOVA)。b)不滿足方差齊性或小樣本時(shí)資料偏態(tài)，則對(duì)用秩變換后再用隨機(jī)的方差分析也可以直接用非參數(shù)隨機(jī)的秩和檢驗(yàn)Fredman test)。例2下表是某湖水中8個(gè)觀察地點(diǎn)不同季節(jié)取樣的氯

20、化物含量測(cè)定值，請(qǐng)問在不同季節(jié)該湖水中氯化物的含量有無(wú)差別？表2某湖水中不同季節(jié)的氯化物含量測(cè)定值（mg/L）locationno春夏秋冬1234567821.2822.7820.9019.9021.4922.3821.6722.0618.3319.8118.9321.2319.0917.9219.3919.6517.2716.5516.3617.8615.1116.5717.1916.5814.9114.8516.3015.7317.0514.3416.3114.33顯然同一地點(diǎn)不同季節(jié)的氯化物含量有一定的相關(guān)性，故不能采用完全隨機(jī)設(shè)計(jì)的方差分析對(duì)4個(gè)季節(jié)的氯化物含量進(jìn)行統(tǒng)計(jì)分61析。可以把

21、同一地點(diǎn)的4個(gè)季節(jié)氯化物含量視為一個(gè)，因此可以用隨機(jī)的方差分析進(jìn)行統(tǒng)計(jì)分析。設(shè)第8個(gè)地點(diǎn)在冬季的氯化物總體均數(shù)為m0，同樣在冬季，第i個(gè)地點(diǎn)的氯化物總體均數(shù)與第8個(gè)地點(diǎn)在冬季的氯化物總體均數(shù)相差bi，i=1，2，3，4，5，6，7。因此在冬季的這8個(gè)地點(diǎn)在冬季的氯化物總體均數(shù)可以表示為地點(diǎn)編號(hào)12345678m0+b1m0+b2m0+b3m0+b4m0+b5m0+b6m0+b7m0冬季氯化物均數(shù)假定在同一地區(qū)，春季的氯化物總體均數(shù)與冬季的氯化物總體均數(shù)相差a1，因此春節(jié)和冬季的氯化物總體均數(shù)可以表示為地點(diǎn)編號(hào)12345678m0+b1m0+b2m0+b3m0+b4m0+b5m0+b6m0+b7

22、m0冬季氯化物均數(shù)m0+a1+b1m0+a1+b2m0+a1+b3m0+a1+b4m0+a1+b5m0+a1+b6m0+a1+b7m0春季氯化物均數(shù)如果a10說明在同一地點(diǎn)，冬季和春季的氯化物總體均數(shù)相同；a1>0說明春季的氯化物含量平均高于冬季氯化物含量，反之a(chǎn)<0，說明春季氯化物含量均數(shù)低于冬季氯化物含量。同理假定在同一地區(qū)秋季的氯化物總體均數(shù)與冬季的氯化物總體均數(shù)分別相差a2和a3，則四個(gè)季節(jié)的氯化物總體均數(shù)可以表示為地點(diǎn)編號(hào)12345678m0+b1m0+b2m0+b3m0+b4m0+b5m0+b6m0+b7m0冬季氯化物均數(shù)春季氯化物均數(shù) m0+a1+b1m0+a1+b2

23、m0+a1+b3m0+a1+b4m0+a1+b5m0+a1+b6m0+a1+b7m0m0+a2+b1m0+a2+b2m0+a2+b3m0+a2+b4m0+a2+b5m0+a2+b6m0+a2+b7m0夏季氯化物均數(shù)m0+a3+b1m0+a3+b2m0+a3+b3m0+a3+b4m0+a3+b5m0+a3+b6m0+a3+b7m0秋季氯化物均數(shù)根據(jù)上述總體均數(shù)表示，可以知道：在四個(gè)季節(jié)中的氯化物總體均數(shù)(同一地點(diǎn))無(wú)變化就是H0：a1=a2=a3=0(在隨機(jī)方差分析中稱為無(wú)處理效應(yīng)，但不能稱4組的總體均數(shù)相同，因?yàn)樵谕患竟?jié)中不同地62點(diǎn)的總體均數(shù)可能不同)。H1：a1，a2，a3不全為0Sta

24、ta數(shù)據(jù)輸入格式其中id 表示觀察地點(diǎn)編號(hào)，t1，2，3，4 對(duì)應(yīng)表示春節(jié)、夏季、秋63tidx1121.275891222.776491320.899431419.90431521.49291622.380851721.673441822.061332118.334052219.805382318.929192421.228142519.092152617.92372719.385692819.649713117.271413216.545673316.360193417.855483515.112963616.565073717.187343816.582794114.905594214.

25、851274316.297824415.72864517.051694614.340884716.313674814.33015季和冬季。Stata 操作命令：anova x t id處 理 效 應(yīng) H0 ： a1 a2 a3 0的 檢 驗(yàn) 對(duì) 應(yīng) 的 統(tǒng) 計(jì) 量F = MS處理= 18.021 = 57.08MS誤差1.036相應(yīng)的P 值<0.0001(計(jì)算機(jī)輸出值是 0.0000)，所以拒絕無(wú)效假設(shè)，可以認(rèn)為 4 個(gè)季節(jié)的氯化物總體均數(shù)不全相同。不同季節(jié)中的兩兩比較用 LSD檢驗(yàn)如下：在輸入 anovaxtid 命令后，再輸入 regress 命令便得到下列結(jié)果64Number of

26、 obs =32R-squared=0.8923 Root MSE= 1.01769Adj R-squared =0.8410Source |Partial SSdfMSFProb > F+M|180.2143261018.021432617.400.0000|t |177.344737359.114912257.080.0000id |2.869589167 .4099413080.400.8942|Residual |21.749618211.0356961+Total |201.96394431 6.5149659365Source |SSdfMSNumber of obs =32

27、+F( 10,21) =17.40M|180.21432610 18.0214326Prob > F=0.0000Residual |21.749618211.0356961R-squared=0.8923-+-Adj R-squared =0.8410Total |201.96394431 6.51496593Root MSE=1.0177其中a1 = 6.081，對(duì)應(yīng)的假設(shè)檢驗(yàn) H0：a1=0 的統(tǒng)計(jì)量 t=11.95，P 值<0.001，95%區(qū)間為(5.022，7.139)，因此可以認(rèn)為春季的氯化物平均高于冬季，差別有統(tǒng)計(jì)學(xué)意義。a2 = 3.816 ，對(duì)應(yīng)的假設(shè)檢驗(yàn) H0

28、：a2=0 的統(tǒng)計(jì)量 t=7.50，P 值<0.001，95%區(qū)間為 (2.758，4.874)，因此可以認(rèn)為夏季的氯化物平均高于冬季，差別有統(tǒng)計(jì)學(xué)意義。a3 = 1.208 ，對(duì)應(yīng)的假設(shè)檢驗(yàn) H0：a3=0 的統(tǒng)計(jì)量 t=2.37，P 值=0.027，95%區(qū)間為 (0.1494，2.266)，因此可以認(rèn)為秋季的氯化物平均高于冬季，差別有統(tǒng)計(jì)學(xué)意義。對(duì)于春季氯化物平均數(shù)(m0+a1+bi)與夏季的氯化物平均數(shù)(m0+a2+bi)比較對(duì)應(yīng)為a1>a2、a1=a2 和a1<a2 的。因此需要檢驗(yàn) H0：a1a2 vsH1：a1¹a2 ，相應(yīng)的 STATA 命令(ano

29、va x t id 命令和 regress 命令66xCoef.Std. Err.tP>|t|95% Conf. Interval_cons (m0)15.37992.596674625.780.00014.1390616.62077ta1=16.080619.508845811.950.0005.0224177.138822a2=23.816041.50884587.500.0002.7578384.874244a3=31.20765.50884582.370.027.14944722.2658534(dropped)idb1=1-.2092595.7196166-0.290.774-

30、1.7057841.287265b2=2.3387067.71961660.470.643-1.1578181.835231b3=3-.034339.7196166-0.050.962-1.5308641.462186b4=4.5231357.71961660.730.475-.9733892.01966b5=5.0314307.71961660.040.966-1.4650941.527955b6=6-.353369.7196166-0.490.628-1.8498941.143156b7=7.4840407.71961660.670.509-1.0124841.9805658(droppe

31、d)后)為 test _bt1=_bt2，得到下列結(jié)果相應(yīng)的統(tǒng)計(jì)量 F=26.28，P 值<0.0001，差別有統(tǒng)計(jì)學(xué)意義。由于a1的估計(jì)值>a2 的估計(jì)值，所以可以認(rèn)為春季氯化物平均高于夏季的氯化物含量。同 理 檢 驗(yàn) H0 ： a1 a3vs H1 ： a1¹a3 ， 只 需 輸 入 命 令 testbt1=_bt3檢驗(yàn)H0：a2a3 vs H1：a2¹a3，只需輸入命令 test bt2=_bt3此處不在詳細(xì)敘述了。ij + X - Xi - X j )服從正態(tài)分布，由于隨機(jī)方差分析要求殘差( e再輸入regress 以后，只要輸入 predict 殘差變

32、量名,residual，就可以得到殘差計(jì)算值。本例用 e 表示殘差變量名，因此輸入 predicte,residual就可以得到殘差計(jì)算值 e，然后對(duì)殘差進(jìn)行正態(tài)性檢驗(yàn)(sktest 殘差變量名)本例輸入命令為:sktest eH1：殘差偏態(tài)分布,a0.05）結(jié)果如下（H0：殘差服從正態(tài)分布vsP 值=0.93349>>a，因此可以認(rèn)為資料近似服從正態(tài)分布。(大樣本時(shí)，67Skewness/Kurtosis tests for Normalityjoint Variable |Pr(Skewness)Pr(Kurtosis) adj chi2(2)Prob>chi2+e |0

33、.6990.5860.460.7948( 1)t2 - t3 = 0.0F( 1,21) =26.28 Prob > F =0.0000)可以不考慮正態(tài)性如果資料呈偏態(tài)分布，可以對(duì)資料進(jìn)行秩變換(Rank Transform)后，然后把變換后的秩視為原始數(shù)據(jù)進(jìn)行隨機(jī)的方差分析。秩變換的STATA 命令為 egen秩變量名rank(觀察變量名),by(變量)為了說明上述操作分析的過程，故借用本例資料進(jìn)行秩變換操作說明如下(本例資料正態(tài)分布，無(wú)需用秩變換，只是說明操作而言).設(shè)用r 表示秩變量名，則本例操作為egen r=rank(x) ,by(id) 產(chǎn)生秩 ranova 命令 anovartid結(jié)果如下命令regress結(jié)果如下68regressSource