醫(yī)用SAS統(tǒng)計分析（三）

醫(yī)用SAS統(tǒng)計分析第三講3/9/20231一、分類變量的統(tǒng)計推斷—卡方檢驗

在SAS系統(tǒng)中，對分類變量資料的基本統(tǒng)計分析方法主要通過FREQ過程實現(xiàn)的。FREQ過程的主要功能有：

1.產(chǎn)生一維或多維頻數(shù)表；2.計算各種表中格子的理論頻數(shù)、構(gòu)成比和各種率；3.對分類變量資料作相應(yīng)的假設(shè)檢驗。3/9/20232

FREQ過程的語句及說明

主要格式：

PROCFREQ[選擇項]；

TABLES表達(dá)式/[選擇項]；

WEIGHT<變量名>；說明：

PROC過程[選擇項]

data=數(shù)據(jù)集；規(guī)定PROCFREQ語句使用的數(shù)據(jù)集。

formchar(1,2,7)=‘|-+’;規(guī)定用來構(gòu)造列聯(lián)表單元的輪廓線和分隔線的字符（只有三個字符）。1為垂線，2為水平線，7為水平與垂直的交叉線。3/9/20233TABLES語句的表達(dá)式

表達(dá)式是要求FREQ過程分析處理的一維或多維表的清單。一維表有一個變量名表示，二維表由星號“*”聯(lián)接兩個變量名表示，如a*b表示變量a與b的二維表。三維表的形式為a*b*c。TABLES語句的[選擇項]1.普通選項out＝數(shù)據(jù)集

建立一個包含變量值和頻數(shù)計數(shù)的輸出數(shù)據(jù)集。如果TABLES語句中不止一個表達(dá)式，數(shù)據(jù)集的內(nèi)容相應(yīng)于TABLES語句中最后一個表達(dá)式的表格。3/9/202352.統(tǒng)計分析主要選項chisq

對每層作χ2檢驗，包括Pearsonχ2、似然比χ2

和Mantel-Haenszelχ2。此外還給出與χ2

檢驗有關(guān)的關(guān)聯(lián)指標(biāo)包括Phi系數(shù)、列聯(lián)系數(shù)和Cramer’sV。對于2×2表，給出Fisher精確概率。agree進(jìn)行配對χ2

檢驗(McNemar’s檢驗);一致性檢驗的Kappa值。exact

對大于2×2的列聯(lián)表計算Fisher精確概率。同時也給出CHISQ選項的全部統(tǒng)計量。cmh

給出Cochran-Mantel-Haenszel統(tǒng)計量，主要用于行×列表的統(tǒng)計分析。對于分層2×2表，cmh過程給出總體相對危險度估計及其可信區(qū)間，還給出各層關(guān)聯(lián)度指標(biāo)是否齊性的Breslow檢驗。3/9/20236measures對每層的二維表計算一系列關(guān)聯(lián)指標(biāo)及相應(yīng)的標(biāo)準(zhǔn)誤，包括Pearson和Spearman相關(guān)系數(shù)，以及Gamma和Kendall系數(shù)等。對于2×2表，還給出常用的危險度指標(biāo)及其可信區(qū)間。all給出chisq，measures，cmh所請求的全部統(tǒng)計量。Alpha=

給出α檢驗水準(zhǔn)。缺省為0.05.expected

給出期望頻數(shù)。3/9/20237四格表資料的卡方檢驗書中例10-1datachisq1;dor=1to2;doc=1to2;inputf@@;output;end;end;cards;63173168;procfreqFormchar(1,2,7)=‘|-+’;weightf;tablesr*c/chisq;run;3/9/20239卡方檢驗結(jié)果

TABLEOFRBYCRCFrequency|二維表每個格子的頻數(shù)Percent|每個格子的頻數(shù)在總頻數(shù)中的百分比。RowPct|行百分?jǐn)?shù)，每格子頻數(shù)占該行合計頻數(shù)的百分比。ColPct|*1|2|Total---------+--------+--------+1|63|17|80|35.20|9.50|44.69|78.75|21.25|*列百分?jǐn)?shù)，每格子頻數(shù)占|67.02|20.00|該列合計頻數(shù)的百分比。---------+--------+--------+2|31|68|99|17.32|37.99|55.31|31.31|68.69||32.98|80.00|---------+--------+--------+Total948517952.5147.49100.003/9/202310

STATISTICSFORTABLEOFRBYCStatisticDFValueProb-------------------------------------------------------Chi-Square①139.9270.001LikelihoodRatioChi-Square②141.8600.001ContinuityAdj.Chi-Square③138.0470.001Mantel-HaenszelChi-Square④139.7040.001Fisher'sExactTest(Left)⑤1.000(Right)1.44E-10(2-Tail)1.64E-10PhiCoefficient⑥0.472ContingencyCoefficient⑦0.427Cramer'sV⑧0.472SampleSize=179①χ2值；②似然比χ2值；③連續(xù)性校正χ2值；④M-Hχ2值；⑤Fisher‘s

精確概率值；⑥Phi系數(shù)；⑦列聯(lián)系數(shù)；⑧可萊姆的V值3/9/202311程序datachi;dor=1to2;doc=1to2;inputf@@;output;end;end;cards;33263;procfreqformchar(1,2,7)='|-+';weightf;tablesr*c/chisqexpectednocolnorownopct;run;3/9/202313RCFrequency|Expected|1|2|Total---------+--------+--------+1|33|2|35|31.023|3.9773|---------+--------+--------+2|6|3|9|7.9773|1.0227|---------+--------+--------+Total39544

3/9/202314STATISTICSFORTABLEOFRBYCStatisticDFValueProb-----------------------------------------------------Chi-Square15.4220.020LikelihoodRatioChi-Square14.3670.037ContinuityAdj.Chi-Square13.0260.082Mantel-HaenszelChi-Square15.2990.021Fisher'sExactTest(Left)0.996(Right)0.050(2-Tail)0.050PhiCoefficient0.351ContingencyCoefficient0.331Cramer'sV0.351SampleSize=44WARNING:50%ofthecellshaveexpectedcountslessthan5.Chi-Squaremaynotbeavalidtest.3/9/202315關(guān)聯(lián)性檢驗(獨(dú)立性檢驗）

應(yīng)用一般的χ2檢驗H0：甲、乙兩試劑的檢驗結(jié)果無關(guān)系。

tables語句選擇\chisq;

如χ2檢驗有統(tǒng)計學(xué)意義，通過列聯(lián)系數(shù)（ContingencyCoefficient)說明兩者的關(guān)聯(lián)性。差別性檢驗（一致性檢驗）應(yīng)用McNemar'sTest

H0：兩總體的B=C，即兩試劑陽性率無差別。3/9/202317當(dāng)b+c<40時tables語句選擇\agree;分析結(jié)果中給出χ2檢驗統(tǒng)計量，Kappa值（一致性檢驗系數(shù)），ASE值（Kappa值的標(biāo)準(zhǔn)誤）及Kappa值的95%可信限。

Kappa值作為評價判斷或診斷一致性程度的指標(biāo)，

Kappa=(Pa-Pe)/(1-Pe)，Pa為實際觀察到的一致率；Pe為期望觀察一致率，即兩次化驗結(jié)果由于偶然機(jī)會所造成的一致率。本例Pa=(80+11)/132=0.6894,Pe=(T11+T22)/132=（75.68+6.68)/132=0.6239。Kappa=(0.6894-0.6239)/(1-0.6239)=0.1745。3/9/202318Kappa值應(yīng)在0到1之間，若Kappa值=1，說明兩次檢驗結(jié)果完全一致，若Kappa值=0，說明兩次檢驗結(jié)果完全沒有一致性。若Kappa值≥0.75，說明一致程度相當(dāng)滿意，若小于0.4，說明一致程度不夠理想。如95%可信限不包含負(fù)數(shù)，說明Kappa值有統(tǒng)計學(xué)意義，兩試劑結(jié)果（或兩次檢驗結(jié)果）有一致性。反之，無一致性。3/9/202319STATISTICSFORTABLEOFRBYCStatisticDFValueProb---------------------------------------------------------Chi-Square14.8670.027LikelihoodRatioChi-Square14.5810.032ContinuityAdj.Chi-Square13.8050.051Mantel-HaenszelChi-Square14.8300.028Fisher'sExactTest(Left)0.992(Right)0.028(2-Tail)0.040PhiCoefficient0.192ContingencyCoefficient0.189Cramer'sV0.192McNemar'sTest--------------Statistic=10.756DF=1Prob=0.001SimpleKappaCoefficient------------------------95%ConfidenceBoundsKappa=0.174ASE=0.0860.0060.342SampleSize=1323/9/202321

行×列表資料的卡方檢驗

行×列表資料可分為三種清況：①雙向無序行×列表，即行變量和列變量均為名義變量（如民族、職業(yè)、血型等），包括多個樣本率的比較。②單向有序行×列表，即行變量或列變量中有一個是名義變量，另一個為有序變量（如檢驗結(jié)果的-，+，++，+++；治療效果的痊愈，良好，有效，無效等）；③雙向有序行×列表，即行變量和列變量均為有序變量。在SAS中，三種情況的統(tǒng)計量是通過TABLES語句中選擇CMH（Cochren-Mantel-Haenzsel）項實現(xiàn)的。CMH檢驗主要出現(xiàn)三個統(tǒng)計量：

1.NonzeroCorrelation(非零相關(guān)）：針對雙向有序行×列表。

2.RowMeanScoresDiffer（行平均分不同）：針對單向有序行×列表。3/9/202322

3.GeneralAssociation

(一般關(guān)聯(lián)）：針對雙向無序行×列表。對雙向無序行×列表在TABLES語句中選擇EXACT，可計算行×列表的精確概率值。雙向無序行×列表

例10-3dataa;dor=1to2;doc=1to4;inputf@@;output;end;end;cards;584959184327338;procfreq;weightf;tablesr*c/exactcmhnoprint;run;3/9/202323列變量為有序變量的行×列表例10-4

datab;dor=1to2;doc=1to3;inputf@@;output;end;end;cards;738414206;procfreq;weightf;tablesr*c/cmhnoprint;run;3/9/202325SUMMARYSTATISTICSFORRBYCCochran-Mantel-HaenszelStatistics(BasedonTableScores)StatisticAlternativeHypothesisDFValueProb--------------------------------------------------------------1NonzeroCorrelation15.9010.0152RowMeanScoresDiffer15.9010.0153GeneralAssociation27.0240.030TotalSampleSize=1163/9/202326

雙向有序且屬性相同的列聯(lián)表的一致性檢驗例兩名放射科醫(yī)師對200名棉屑沉著病可疑患者的診斷結(jié)果見下表，問兩醫(yī)生診斷結(jié)果否一致。

表200例棉屑沉著病可疑患者的診斷結(jié)果第一人檢查第二人檢查合計正常Ⅰ型Ⅱ型正常Ⅰ型Ⅱ型合計786055610013328375428471452003/9/202329程序dataagree1;dor=1to3;doc=1to3;inputf@@;output;end;end;cards;78506561301032;procfreqformchar(1,2,7)='|-+';weightf;tablesr*c/agreenocolnorownopct;run;3/9/202330

TABLEOFRBYCRCFrequency|1|2|3|Total---------+--------+--------+--------+1|78|5|0|83---------+--------+--------+--------+2|6|56|13|75---------+--------+--------+--------+3|0|10|32|42---------+--------+--------+--------+Total847145200

STATISTICSFORTABLEOFRBYCTestofSymmetry----------------Statistic=0.482DF=3Prob=0.923KappaCoefficientsStatisticValueASE95%ConfidenceBounds------------------------------------------------------SimpleKappa0.7370.0410.6570.816WeightedKappa0.7940.0320.7310.858SampleSize=2003/9/202331多層行×列表資料的卡方檢驗

CMH統(tǒng)計量也適用于多層行×列表資料的的卡方檢驗，即按一個或多個因素分層后，研究行變量和列變量間的聯(lián)系。其可通過控制分層變量的影響后，檢驗行變量和列變量的關(guān)系。

多層的2×2表

表男女使用別嘌呤醇引發(fā)皮疹比較性別藥物發(fā)病未發(fā)病男使用未使用使用未使用53610583364519518女3/9/202332

程序datachisq1;dosex=1to2;/*1為男，2為女*/dodrug=1to2;/*1為使用，2為未使用*/docase=1to2;/*1為發(fā)病，2為未發(fā)病*/inputf@@;output;end;end;end;cards;53336645101958518;procfreqformchar(1,2,7)='|-+';weightf;tablessex*drug*case/cmhmeasuresnopctnorownocol;run;3/9/202333TABLE1OFDRUGBYCASECONTROLLINGFORSEX=1DRUGCASEFrequency|1|2|Total---------+--------+--------+1|5|33|38---------+--------+--------+2|36|645|681---------+--------+--------+Total41678719

EstimatesoftheRelativeRisk(Row1/Row2)95%TypeofStudyValueConfidenceBounds--------------------------------------------------Case-Control2.7151.0007.369Cohort(Col1Risk)2.4891.0365.980Cohort(Col2Risk)0.9170.8091.039SampleSize=719

3/9/202334TABLE2OFDRUGBYCASECONTROLLINGFORSEX=2DRUGCASEFrequency|1|2|Total---------+--------+--------+1|10|19|29---------+--------+--------+2|58|518|576---------+--------+--------+Total68537605

EstimatesoftheRelativeRisk(Row1/Row2)95%TypeofStudyValueConfidenceBounds-----------------------------------------------------Case-Control4.7012.08610.592Cohort(Col1Risk)3.4241.9605.983Cohort(Col2Risk)0.7290.5590.950SampleSize=6053/9/202335

（1）SUMMARYSTATISTICSFORCASEBYDRUGCONTROLLINGFORSEXCochran-Mantel-HaenszelStatistics(BasedonTableScores)StatisticAlternativeHypothesisDFValueProb--------------------------------------------------------------1NonzeroCorrelation119.5130.0012RowMeanScoresDiffer119.5130.0013GeneralAssociation119.5130.001（1）由于行變量和列變量都只有2個，可以看作雙向無序資料。三種結(jié)果是一致的，結(jié)果為：χ2mch=19.513，所對應(yīng)的p=0.001，說明使用別嘌呤醇與發(fā)生藥物性皮疹是有關(guān)系的。3/9/202336

（2）EstimatesoftheCommonRelativeRisk(Row1/Row2)95%TypeofStudyMethodValueConfidenceBounds--------------------------------------------------------------Case-ControlMantel-Haenszel3.7562.0886.757(OddsRatio)Logit3.7772.0117.094CohortMantel-Haenszel3.0441.8584.989(Col1Risk)Logit3.1241.9515.001CohortMantel-Haenszel0.8380.7740.906(Col2Risk)Logit0.8790.7850.985TheconfidenceboundsfortheM-Hestimatesaretest-based.

（3）Breslow-DayTestforHomogeneityoftheOddsRatiosChi-Square=0.703DF=1Prob=0.402TotalSampleSize=13243/9/202337（2）本部分有3個指標(biāo)：

第1個指標(biāo)是病例對照研究（Case-Control）的比數(shù)比(OR)值及可信區(qū)間。一個是用Mantel-Haenszel方法計算的值，另一個是用Logit方法計算的值。本例說明使用別嘌呤醇發(fā)生藥物性皮疹的危險性是不用該藥危險性的3.756倍或3.777倍。第2個和第3個指標(biāo)是隊列研究（Cohort）的相對危險度(RR)值及可信區(qū)間。第2個指標(biāo)是使用藥物發(fā)病率是不使用藥物發(fā)病率的倍數(shù)；第3個指標(biāo)是使用藥物未發(fā)病率是不使用藥物未發(fā)病率的倍數(shù)。（3）Breslow-Day檢驗

是檢驗各層間危險度的方向是否一致，P>0.05說明總體方向是一致的，P≤0.05說明總體方向是不一致的。如方向一致可以用總的OR值表示，否則不能用總的OR值表示。3/9/202338

多層的行×列表程序10-6datachisq2;dohospital=1to2;dotrt=1to2;doeffect=1to3;inputf@@;output;end;end;end;cards;237220135186213132;procfreqformchar(1,2,7)='|-+';weightf;tableshospital*trt*effect/cmhnopctnocol;run;3/9/202339二、二項分布與Poisson分布二項分布二項分布的統(tǒng)計說明二項分布的應(yīng)用條件

1.兩分類對立結(jié)果資料；2。試驗結(jié)果是相互獨(dú)立，互不影響的。

二項分布的概率函數(shù)X=0,1,2,…,n。

二項分布的分布函數(shù)至多有x例陽性的概率為，即下側(cè)累積概率：至少有x例陽性的概率為，即上側(cè)累積概率：3/9/202340SAS中二項分布函數(shù)SAS中二項分布函數(shù)表達(dá)式為：probbnml(π,n,x);式中π

為事件發(fā)生的概率，n為樣本含量，x為陽性事件個數(shù)。表達(dá)式計算的結(jié)果為發(fā)生陽性事件數(shù)≤x的概率。例12-1datajin;p=probbnml(0.5,4,2);q=1-p;procprint;run;3/9/202341結(jié)果OBSPQ10.68750.3125表明死亡≤2的概率（0+1+2的概率）為0.6875；死亡>2的概率（3+4的概率）為0.3124。用二項分布檢驗治療效果

例12-3根據(jù)以往的經(jīng)驗用一般療法治療某病，其病死率為40%，治愈率為60%。今用某種新藥治療該病人5名，這5名均治愈了。問該項新藥是否比一般療法為優(yōu)。

當(dāng)實際率大于理論率時，計算上則累積概率。本例應(yīng)計算：

p=probbnml(0.6,5,5)-probbnml(0.6,5,4)3/9/202342程序12-3databinom2;p=probbnml(0.6,5,5)-probbnml(0.6,5,4);procprint;run;OBSP10.07776假如10名病人，結(jié)果1名死亡，9名治愈，則應(yīng)計算P(x≥9）的概率。程序12-4databinom3;p=probbnml(0.6,10,10)-probbnml(0.6,10,8);Procprint;Run;OBSP10.0463573/9/202343求二項分布的概率分布

例12-5

設(shè)π=0.5，n=4，求x=0，1，2，3，4的概率databinom4;dor=0to4;p=probbnml(0.5,4,r);q=1-p;ifr=0thend=p;elsed=probbnml(0.5,4,r)-probbnml(0.5,4,r-1);output;end;procprint;run;OBSRPQD100.06250.93750.0625210.31250.68750.2500320.68750.31250.3750430.93750.06250.2500541.00000.00000.06253/9/202344

Poisson分布

Poisson分布的統(tǒng)計說明Poisson分布的應(yīng)用條件同二項分布，Poisson分布是二項分布的特例，當(dāng)π很小，n很大時，Poisson分布非常接近二項分布。

Poisson分布概率函數(shù)（X=0,1,2,…）λ=π×n

Poisson分布的分布函數(shù)

下側(cè)累積概率：上側(cè)累積概率：3/9/202345SAS中Poisson分布函數(shù)SAS中Poisson分布函數(shù)表達(dá)式為：poisson(λx);式中λ為均數(shù)（λ

=n×π），x為事件發(fā)生個數(shù)。表達(dá)式計算的結(jié)果為發(fā)生事件數(shù)≤x的概率。例12-6有人觀察紅細(xì)胞計數(shù)池中400個小格，數(shù)出每小格中紅細(xì)胞數(shù)，其均數(shù)為3.6175，試計算每格中恰有細(xì)胞數(shù)5個的概率及小格數(shù)。datapoisson1;p=poisson(3.6175,5)-poisson(3.6175,4)d=p*400;procprint;run;OBSPD10.1386155.44483/9/202346樣本計數(shù)與總體均數(shù)差別的統(tǒng)計意義檢驗

例12-8已知在一培養(yǎng)液中，有細(xì)菌數(shù)為每毫升3個，今采集放在5℃冰箱的1ml培養(yǎng)液的細(xì)菌數(shù)5個，能否說明培養(yǎng)液中細(xì)菌數(shù)有增加？本例λ=3，x=5>3，計算x≥5的上則累積概率。

可計算p=1-p(n≤4）的概率。Datapoisson3;P=1-poisson(3,4);Procprint;Run;OBSP10.18474因P>0.05，尚不能認(rèn)為培養(yǎng)液中細(xì)菌數(shù)有增長。3/9/202347

例

某地區(qū)以往胃癌發(fā)病率為萬分之一?，F(xiàn)調(diào)查10萬人，發(fā)現(xiàn)3例胃癌病人。試分析該地區(qū)現(xiàn)在的胃癌發(fā)病率是否低于以往的發(fā)病率。本例π=0.001，n=100000，x=3，λ=100000×0.0001=10，現(xiàn)X<λ，應(yīng)計算≤x的下則累積概率。dataa;p=poisson(10,3);procprint;run;OBSP10.010336因P<0.05，可以認(rèn)為該地區(qū)現(xiàn)在的胃癌發(fā)病率低于以往的發(fā)病率。3/9/202348三、非參數(shù)統(tǒng)計

在SAS中，非參數(shù)統(tǒng)計主要由univariate過程和npar1way過程來實現(xiàn)，前者在前面的第一講中已經(jīng)介紹，它可以進(jìn)行配對設(shè)計差值的符號秩和檢驗（Wilcoxon配對法）；后者是一個單因素的非參數(shù)方差分析過程，可進(jìn)行成組設(shè)計的兩樣本（Wilcoxon法）或多樣本比（Kruskal-Wallis法即H檢驗）的秩和檢驗。

注意：①npar1way過程不能處理按頻數(shù)輸入的資料。對單項有序資料的頻數(shù)表，可用前面講到的cmh過程的RowMeanScoresDiffer（行平均分不同）來檢驗。

②SAS不提供非參檢驗兩兩比較的方法，其余統(tǒng)計軟件里也只有PEMS提供這一功能。

3/9/202349③npar1way含意：“npar”是“非參”的英文縮寫，“1way”代表一個因素，合起來就是“單因素的非參數(shù)檢驗”。千萬不要寫成“NPAPLWAY”!NPAP1WAY過程的語法格式PROCNPAR1WAY[DATA=<數(shù)據(jù)集名>[選項]];

CLASS<分組變量名>;

必需，指定要分析的分組變量

EXACT<關(guān)鍵字>;

要求程序在必要時計算確切概率

OUTPUT<OUT=數(shù)據(jù)集名><選項>;

指定統(tǒng)計結(jié)果的輸出數(shù)據(jù)集

VAR<結(jié)果變量名>;

指定要分析的因變量

BY<變量名列>;

統(tǒng)計按指定的變量分組進(jìn)行，要求已排序。3/9/202350PROC過程的[選項]MISSING

將缺失值也用于統(tǒng)計分析

ANOVA

同時進(jìn)行方差分析

MEDIAN

要求進(jìn)行中位數(shù)檢驗

NOPRINT

禁止統(tǒng)計結(jié)果在OUTPUT視窗內(nèi)輸出

SAVAGE

要求對樣本進(jìn)行SAVAGE得分分析

WILCOXON

要求進(jìn)行WILCOXON秩和檢驗，通常必選。

在省略所有選項的情況下，SAS系統(tǒng)默認(rèn)輸出所有的統(tǒng)計結(jié)果，這恰恰說明了非參數(shù)檢驗方法的不完善。

3/9/202351

配對資料的符號秩和檢驗

和配對t檢驗過程一樣，用Procunivariate過程，結(jié)果觀察SgnRank(即T值）Pr>=|s|(即P值）。如例13-1T值=-8.5，P=0.3594。（SAS程序略）

兩樣本成組比較程序13-2datanpar1;doc=1to2;inputn;doI=1ton;inputx@@;output;end;end;cards;1213414610411912416110783113129971237701181018510713294;procnpar1waywilcoxon;classc;Varx;run;3/9/