均数估计与t检验：SAS、R、Stata、SPSS实战指南

下午刚把第三章的作业看完十个里面至少有五个问我同一个问题均数估计和假设检验的公式明明会手算为什么换成SAS、R、Stata、SPSS就不知道点哪里、写什么命令。这个感觉我太懂了。学校课程往往只带一遍某一款软件真到了自己手里SAS输出一堆英文、R报错、Stata命令记混、SPSS菜单藏得深最后连结果对不对都不敢确认。这一章的核心其实不复杂先估计总体均数再检验假设真正让人头大的是把同样的例题在四款软件里都跑一遍。这篇文章我直接用一组常见例题把均数估计、单样本t检验、两独立样本t检验的手算结果算出来再分别给你SAS、R、Stata、SPSS的完整实现和输出解读。代码和菜单路径都已验证过你照着做就能出结果。1. 均数估计与假设检验先把逻辑想明白1.1 从点估计到置信区间均数估计要回答的问题是总体均值到底是多少样本均值是最自然的估计也就是点估计。问题是样本均值存在抽样误差换一拨样本数值就变光报一个点不够所以要有区间估计。经典的置信区间公式是[ \bar{x} \pm t_{1-\alpha/2,\ n-1} \times \frac{s}{\sqrt{n}} ]其中 (\bar{x}) 是样本均值(s) 是样本标准差(n) 是样本量(t) 是自由度为 (n-1) 的t分布临界值。为什么这里用t分布而不是标准正态分布的z值因为总体标准差通常未知只能拿样本标准差 (s) 去代替。样本量越小这个代替带来的不确定性越大t分布尾部更厚等于给结论多留一点余地。如果样本量到了三五十以上t分布的临界值已经和z值非常接近这时用哪种差别都不大。关于置信区间有一个常见误解95%置信区间不是“总体均值有95%的概率落在这个区间里”。正确理解是如果按同样的方式反复抽样100次大约有95个这样的区间能把总体均值罩住。软件给出的区间就是这一次抽样得到的那个区间。1.2 假设检验从原假设到P值假设检验解决的问题是有没有证据。最典型的两组设定原假设 (H_0)均数等于某个值或者两组均数相等备择假设 (H_1)均数不等于某个值或者两组均数不等。检验统计量就是信号和噪声的比值[ t \frac{\text{样本估计值} - \text{原假设中的值}}{\text{标准误}} ]这个比值是“相对差异”不是绝对差异。P值的意思是如果原假设成立出现当前这种样本结果甚至更极端结果的概率。P值越小说明样本数据越不像原假设所描述的世界于是拒绝原假设。P值不是“原假设成立的概率”这一点一定别搞混。至于显著性水平 (\alpha)一般取0.05就是给自己定一个犯错的容忍度。软件里最重要的几个数字就是t统计量、自由度、P值和置信区间。2. 一个例题两组数据手算结果先当标准答案2.1 例题设计与数据结构用一个医药领域最常见的场景比较两种降压药对收缩压下降值的影响。A药组和B药组各10名受试者每人记录治疗前后的收缩压下降值单位mmHg。数据结构采用长格式每一行是一个受试者。编号drugdecrease111221153184110511461971138111911610110112812210132614291527162517211182819272029这里drug1代表A药drug2代表B药。decrease是下降值数值越大说明药效越强。2.2 手算一组标准答案先算各组基本描述统计量组别n均值标准差标准误95%置信区间A药组1011.802.660.84(9.90, 13.70)B药组108.001.830.58(6.69, 9.31)A药组单样本t检验原假设为下降值均数等于0[ t \frac{11.80 - 0}{2.66 / \sqrt{10}} 14.04,\quad df9,\quad P0.001 ]两独立样本t检验比较A药和B药的下降值差异。两组样本方差分别是7.07和3.33合并方差[ s_p^2 \frac{9 \times 7.07 9 \times 3.33}{18} 5.20 ]合并标准误[ s_p \sqrt{5.20 \times \left(\frac{1}{10}\frac{1}{10}\right)} 1.02 ][ t \frac{11.80 - 8.00}{1.02} 3.73,\quad df18,\quad P0.01 ]两组均值差为3.80mmHg95%置信区间大概在(1.66, 5.94)区间不包含0说明A药对收缩压的下降效果在统计上比B药更明显。先把这个手算结果记下来后面四款软件跑出来的数字如果和它对不上说明操作有问题不是软件坏了。3. 用SAS实现均数置信区间与t检验3.1 数据步一次把数据集建好SAS里先要用数据步读入数据。因为drug是分组变量在INPUT语句里写清楚变量名即可。data hypertension; input drug decrease; datalines; 1 12 1 15 1 8 1 10 1 14 1 9 1 13 1 11 1 16 1 10 2 8 2 10 2 6 2 9 2 7 2 5 2 11 2 8 2 7 2 9 ; run;注意DATALINES后面直接放数据数据块结束的RUN;要放在最后。这个格式是SAS最基础也最实用的读入方式复制粘贴就能跑。3.2 用PROC MEANS输出均值和置信区间SAS中计算描述统计量用PROC MEANS加上CLM选项就能直接输出均值的置信区间。proc means datahypertension n mean std stderr clm alpha0.05; class drug; var decrease; run;输出中的N、Mean、Std Dev、Std Error、Lower 95% CL for Mean、Upper 95% CL for Mean分别对应手算表格里的样本量、均值、标准差、标准误和置信区间上下限。CLASS DRUG会自动按照分组输出不用提前排序比BY语句更方便。我习惯把CLM念成“Confidence Limits for Mean”它默认给的是95%区间。如果你想看90%或99%改ALPHA0.10或ALPHA0.01就行。3.3 用PROC TTEST完成单样本和两样本t检验单样本t检验检验A药组下降值是否不为0proc ttest datahypertension h00 plotsnone; var decrease; where drug1; run;WHERE DRUG1只保留A药组H00表示原假设的均数值是0。输出里要重点看Mean、t Value、Pr |t|和95% CL Mean。这里Pr |t|就是双侧P值。两独立样本t检验proc ttest datahypertension plotsnone; class drug; var decrease; run;CLASS DRUG指定分组变量。SAS的PROC TTEST会同时给出两组均值的置信区间以及两种方差假设下的t检验结果Pooled假设两组方差相等Satterthwaite不假设方差相等做了校正。旁边还有一个方差齐性的F检验。F检验P值大于0.05时看Pooled那行小于0.05时看Satterthwaite那行。这组例题F检验不显著所以主要看Pooled行t3.73P值约0.0016与手算一致。SAS默认还会画图和输出置信区间如果不想看图形PLOTSNONE能省掉一堆输出。4. 用R实现三分钟跑完分析4.1 数据准备与描述统计R最灵活的地方是向量和公式接口。先输入两组数据拼成长格式decrease_A - c(12, 15, 8, 10, 14, 9, 13, 11, 16, 10) decrease_B - c(8, 10, 6, 9, 7, 5, 11, 8, 7, 9) group - factor(rep(c(A, B), each 10)) decrease - c(decrease_A, decrease_B) df - data.frame(group group, decrease decrease)看分组均值、标准差和样本量aggregate(decrease ~ group, data df, FUN function(x) { c(n length(x), mean mean(x), sd sd(x)) })如果希望连置信区间一起算可以用t.test的单独输出也可以自己手写n - length(decrease_A) mean_A - mean(decrease_A) sd_A - sd(decrease_A) se_A - sd_A / sqrt(n) t_crit - qt(0.975, df n - 1) c(mean_A - t_crit * se_A, mean_A t_crit * se_A)这行代码就是置信区间公式的落地版建议每个人都亲手敲一遍。4.2 t.test函数单样本、独立样本、配对样本单样本t检验t.test(decrease_A, mu 0)MU0就是原假设的参考值。输出里有t14.037, df9, p-value很小95 percent confidence interval为9.90到13.70。两独立样本t检验t.test(decrease_A, decrease_B, var.equal TRUE)如果R没设VAR.EQUALTRUE默认用的是Welch检验也就是不假设方差相等。这组数据两组方差差别不大Welch结果也非常接近自由度会变成非整数t约等于3.70P值约0.002。手动做作业时为了和教材公式一致我会明确写var.equalTRUE实际研究我更推荐用Welch因为它少了一个容易翻车的假设。R也可以用公式写法t.test(decrease ~ group, data df, var.equal TRUE)如果处理的是同一批受试者治疗前后的数据用配对t检验before - c(...) after - c(...) t.test(after, before, paired TRUE)注意配对t检验本质上就是对差值做的单样本t检验。软件里如果给的是“差值”这一列直接跑单样本如果给的是“前测”“后测”两列就必须用PAIREDTRUE否则等于把两个变量当成独立样本会严重错。4.3 把R输出整理成可读表格R原样输出有点乱我习惯用broom包把t检验结果整理成数据框library(broom) tidy(t.test(decrease_A, mu 0))一行就得到estimate、statistic、p.value、conf.low、conf.high适合批量分析和写报告。这个方法在样本量多、检验做了几十个的时候特别好用就不用逐条复制控制台了。5. 用Stata实现命令比你想的更直白5.1 数据输入与均值置信区间Stata的命令风格比较口语化。先输入数据input drug decrease 1 12 1 15 1 8 1 10 1 14 1 9 1 13 1 11 1 16 1 10 2 8 2 10 2 6 2 9 2 7 2 5 2 11 2 8 2 7 2 9 end计算总体均值的置信区间ci mean decrease按分组计算均值和置信区间by drug, sort: ci mean decrease输出会按drug1和drug2各给一行包含Obs、Mean、Std. Err、Std. Dev和95% Conf. Interval。STATA里的Std. Err就是标准误Std. Dev是标准差别弄混。5.2 ttest命令与结果解读单样本t检验检验A药组下降值均数是否为0ttest decrease 0 if drug 1两独立样本t检验ttest decrease, by(drug)Stata输出有一个很有特色的地方除了双侧P值还会给单侧假设的P值。表格下方的Pr(T t) 是左侧单侧P值Pr(T t) 是右侧单侧P值Pr(|T| |t|) 是双侧P值。平时报告用的就是最后一个。Stata默认的BY形式t检验假设两组方差相等。如果不想做方差齐性假设加WELCH选项ttest decrease, by(drug) welch还可以用sdtest检查两组方差是否相等sdtest decrease, by(drug)这组数据SDTEST的P值不显著说明方差不齐问题不大标准t检验结果可信。Stata比较适合课堂演示因为命令短、输出规整几乎不需要额外调格式。缺点是某些默认假设不明显比如默认等方差容易让人忽略这个前提。6. 用SPSS实现点鼠标也能把逻辑讲清楚6.1 数据录入格式先把表格摆对SPSS是最容易上手的也是我见过出错最多的。常见问题是把两组独立样本的数据排成两列一行一个组别把重复测量数据也排成两列SPSS根本无法判断分析类型。本例子SPSS数据视图应该是三列编号、drug、decrease。每一行是其中一个受试者的数据drug用1和2编码。如果只有差值数据就一列decrease一行一个受试者。把drug定义成数值变量还是字符串变量无所谓但数值变量更省事。定义分组时直接输入1和2避免字符串引号问题。6.2 菜单路径均值置信区间与t检验如果想按分组看均值置信区间用Analyze Descriptive Statistics Explore把decrease选入Dependent Listdrug选入Factor List在Statistics里勾选DescriptivesPlots可以全部取消输出会给出每组的Mean、Std. Deviation、Std. Error和95% Confidence Interval for Mean。单样本t检验Analyze Compare Means One-Sample T Test把decrease选入Test VariableTest Value填0Options里确认Confidence Interval Percentage是95如果只检验A药组需要先用Data Select Cases选择drug1或者用Split File按drug分组后再跑。直接对全部20行做单样本t检验检验的是两组加在一起的总均值不是A药组这是新手最常踩的坑。两独立样本t检验Analyze Compare Means Independent-Samples T TestTest Variable选decreaseGrouping Variable选drug点击Define GroupsGroup 1填1Group 2填2输出里自动给Levene方差齐性检验和t检验结果SPSS输出比R更啰嗦但信息都在。关键看两行如果Levene检验P值大于0.05读“Equal variances assumed”那一行如果小于0.05读“Equal variances not assumed”那一行。本例Levene检验不显著读第一行t3.726df18Sig. (2-tailed)0.0016Mean Difference3.8095% Confidence Interval of the Difference约为1.66到5.94。6.3 SPSS Syntax如果不想一直点菜单点菜单当然直观但重复做多次时效率太低。SPSS也支持语法和菜单等价T-TEST /TESTVAL0 /VARIABLESdecrease /CRITERIACI(.95). T-TEST GROUPSdrug(1 2) /VARIABLESdecrease /CRITERIACI(.95).只要数据里变量名是decrease和drug这两行语法直接跑输出和菜单操作一模一样。建议把语法文件存起来下次换一批数据只改变量名就行。7. 实战中常见的坑从数据格式到结果解读7.1 独立样本和配对样本的数据格式千万别搞混这是四个软件里最容易翻车的问题。独立样本t检验要求长格式一个分组变量一个结果变量每行一个观测。配对样本t检验则要求两个结果变量分别放两列比如before和after。我还见过有同学把配对数据用SPSS的Independent-Samples T Test去跑直接把after和before当成两个组结果P值翻了一倍都不止。解决办法很简单先问自己“这两组数据是不是来自同一批人”。同一批人前后测用配对两拨人用独立。7.2 方差齐性到底该选哪种t检验传统教材会先做F检验或Levene检验再决定读哪一行。实际操作中很多统计课程老师会建议直接使用不假设方差相等的Welch检验因为它不需要经历“先检验方差再检验均值”的两步决策避免检方差这一步本身的不确定性。这组例题两种方式结论一致。但如果你自己跑数据时发现一方标准差明显大于另一方偏好几倍那就别犹豫了直接用Welch结果。在R里就是var.equalFALSE在SAS里看Satterthwaite行在Stata里加welch在SPSS里读“Equal variances not assumed”方法都是现成的。7.3 正态性假设t检验到底怕什么t检验要求样本均数的抽样分布近似正态。根据中心极限定理样本量较大时就算原始数据偏态也能忍。但当n很小比如每组小于15数据又明显偏态这时t检验的结论可能不稳。可以在正式分析前跑一下Shapiro-Wilk检验SAS用PROC UNIVARIATE NORMALR用shapiro.testStata用swilkSPSS用Analyze Descriptive Statistics Explore里的Normality Plots。如果正态性严重不满足考虑Wilcoxon秩和检验或符号秩检验。注意这里的“正态不强求”只针对均数比较不是哪都能放宽。7.4 P值大于0.05时别急着说“没有差异”假设检验的结论只有“拒绝原假设”和“不拒绝原假设”两种。P大于0.05说明现有数据不足以证明差异不代表两组真的完全一样。很多时候是样本量不够效应量又不小统计功效不足。我写结果时习惯把P值和置信区间一起报告比如“两组均值差为3.80mmHg95%CI为1.66到5.94P0.0016”。置信区间能同时看出效应的大小和方向比单看一个P值有信息量得多。如果区间跨过0无论P值看起来多接近0.05都说明差异证据不足。7.5 为什么四款软件输出有时候对不上不同软件版本、不同默认设置会让结果看起来略有出入。最常见的是三处自由度Welch检验自由度不是整数和常规t检验不同方差假设有的默认等方差有的默认不等方差单双侧所有软件默认都报告双侧P值但Stata额外列了单侧P值别读错。所以拿到软件输出先别急着贴到论文里用文章里的手算结果验一遍。手算、SAS、R、Stata、SPSS五套结果能对得上说明你真正理解了这一章。最后再分享一个我很推荐的验证方法把样本量、均值和标准差写进R或Stata试着手动构造t统计量再用软件自带命令跑一遍两边对比。这种“笨办法”能让你把公式和软件输出彻底打通。以后不管换什么软件、遇到多复杂的例子都知道自己要找哪个数字了。