本文通过多个实例展示了使用R、SAS、STATA和SPSS进行数据读取、描述统计、正态性检验、置信区间计算以及t检验。涉及的内容包括夏皮罗正态性检验、样本均数的置信区间计算、单样本和配对样本t检验，以及两独立样本t检验。案例涵盖了医学研究、统计分析等多个场景。
例3-1
R代码及结果
#install.packages("readxl") # 第一运行该文件时加载readxl包用于读取Excel输入数据，将本行开始的#号删除即可
library( readxl ) # 调用外部数据程序包readxl
# 读取Data文件夹下的03_01例题的Excel数据文件（“../Data/E03_01.xls”代表上一级文件夹Data中的E03_01.xls）
E03_01 <- read_excel( "../Data/E03_01.xls" ) # 符号“../” 表示回溯到上一级文件夹
#################### 描述读取的例题数据 ###########################
ls() # 查看工作数据文件名称列表
str(E03_01) # 查看数据结构：no = 样本号，mean=样本均数，sd=标准差
head(E03_01) # 查看数据前几例
summary(E03_01) # 各列的统计量（未分组）
typeof(E03_01) # E03_01的数据类型
#################### 下面为例题的统计分析 ###########################
shapiro.test(E03_01$mean) # 夏皮罗正态性检验-样本均数
help( shapiro.test ) # 夏皮罗正态性检验函数shapiro.test的详细说明
结果表示如下：
library( readxl ) # 调用外部数据程序包readxl
> library(readxl)
> E03_01 <- read_excel("F:/Data/E03_01.xls")
> View(E03_01)
> ls() # 查看工作数据文件名称列表
[1] "E03_01"
> str(E03_01) # 查看数据结构：no = 样本号，mean=样本均数，sd=标准差
tibble[,3] [100 x 3] (S3: tbl_df/tbl/data.frame)
$ no : num [1:100] 1 2 3 4 5 6 7 8 9 10 ...
$ mean: num [1:100] 124 121 125 123 124 ...
$ sd : num [1:100] 5.94 6.28 3.93 5.17 5.14 3.94 6.02 5.17 3.42 5.12 ...
> head(E03_01) # 查看数据前几例
# A tibble: 6 x 3
no mean sd
<dbl> <dbl> <dbl>
1 1 124. 5.94
2 2 121. 6.28
3 3 125. 3.93
4 4 123. 5.17
5 5 124. 5.14
6 6 124. 3.94
> summary(E03_01) # 各列的统计量（未分组）
no mean sd
Min. : 1.00 Min. :119.4 Min. :2.000
1st Qu.: 25.75 1st Qu.:122.8 1st Qu.:4.117
Median : 50.50 Median :123.8 Median :4.865
Mean : 50.50 Mean :123.7 Mean :4.906
3rd Qu.: 75.25 3rd Qu.:124.7 3rd Qu.:5.758
Max. :100.00 Max. :127.1 Max. :8.320
> typeof(E03_01) # E03_01的数据类型
[1] "list"
> shapiro.test(E03_01$mean) # 夏皮罗正态性检验-样本均数
Shapiro-Wilk normality test
data: E03_01$mean
W = 0.99336, p-value = 0.9086 # p>0.5 尚不能认为这些样本的总体不服从正态分布
> help( shapiro.test ) # 夏皮罗正态性检验函数shapiro.test的详细说明
例3-2 在例3-1中抽得第26号样本的均数为122.96（cm)，标准差为4.77（cm），求其总体均数的95%置信区间。
用下面公式计算：
SAS代码及结果
data ex3_2;
n=10;
mean=122.96;
std=4.77;
t=tinv(0.975,n-1);*t分布的分位数*
pts=t*std/sqrt(n);
lclm=mean-pts;
uclm=mean+pts;
proc print;
var lclm uclm;
run;
STATA代码及结果
cii mean 10 122.96 4.77
例3-3 某地抽取正常成年人200名，测得其血清胆固醇的均数为3.64mmol/L,标准差为1.20mmol/L,试估计该地正常成年人血清胆固醇均数的95%置信区间。
本例n>60,可以用正态近似的方法计算置信区间。
SAS代码及结果
data ex3_3;
n=200;
mean=3.64;
std=1.20;
u=PROBIT(0.975); /* 正态分布的分位数*/
pts=u*std/sqrt(n);
lclm=mean-pts;
uclm=mean+pts;
proc print;
var lclm uclm;
run;
例3-5 某医生测量了36名从事铅作业男性工人的血红蛋白含量，算得其均数为130.83g/L，标准差为25.74g/L。问从事铅作业工人的血红蛋白是否不同于正常成年男性平均值140g/L？
SAS代码1及结果
data ex3_5;
n=36;
s_m=130.83;
std=25.74;
p_m=140;
df=n-1;
t=(s_m-p_m)/(std/sqrt(n));
p=(1-probt(abs(t),df))*2;
proc print;
var t p;
run;
SAS代码2及结果
data ex3_5;
input x@@;
cards;
112 137 129 126 88 90 105 178 130 128 126 103
172 116 125 90 96 162 157 151 135 113 175 129
165 171 128 128 160 110 140 163 100 129 116 127
;
run;
proc ttest h0=140;
var x;
run;
R代码及结果
#install.packages("readxl") # 第一运行该文件时加载readxl包用于读取Excel输入数据，将本行开始的#号删除即可
library( readxl ) # 调用外部数据程序包readxl
# 读取Data文件夹下的03_05例题的Excel数据文件（“../Data/E03_05.xls”代表上一级文件夹Data中的E03_05.xls）
E03_05 <- read_excel( "../Data/E03_05.xls" ) # 符号“../” 表示回溯到上一级文件夹
#################### 描述读取的例题数据 ###########################
ls() # 查看工作数据文件名称列表
str(E03_05) # 查看数据结构：hb=血红蛋白含量g/L
head(E03_05) # 查看数据前几例
summary(E03_05) # 各列的统计量（未分组）
typeof(E03_05) # E03_05的数据类型
#################### 下面为例题的统计分析 ###########################
#单样本t检验，原假设H0：μ=140 备择假设H1:μ！=140
t.test(E03_05$hb,alternative="two.sided",mu=140)
help( t.test ) # 单样本t检验函数t.test的详细说明
结果：
> E03_05 <- read_excel("F:/Data/E03_05.xls")
> View(E03_05)
> knitr::opts_chunk$set(echo = TRUE)
> ls() # 查看工作数据文件名称列表
[1]"E03_05"
> str(E03_05) # 查看数据结构：hb=血红蛋白含量g/L
tibble[,1] [36 x 1] (S3: tbl_df/tbl/data.frame)
$ hb: num [1:36] 112 137 129 126 88 90 105 178 130 128 ...
> head(E03_05) # 查看数据前几例
> summary(E03_05) # 各列的统计量（未分组）
hb
Min. : 88.0
1st Qu.:112.8
Median :128.0
Mean :130.8
3rd Qu.:152.5
Max. :178.0
> typeof(E03_05) # E03_05的数据类型
[1] "list"
> #单样本t检验，原假设H0：μ=140 备择假设H1:μ！=140
> t.test(E03_05$hb,alternative="two.sided",mu=140)
One Sample t-test
data: E03_05$hb
t = -2.1367, df = 35, p-value = 0.03969
alternative hypothesis: true mean is not equal to 140
95 percent confidence interval:
122.1238 139.5428
sample estimates:
mean of x
130.8333
STATA代码及结果
ttest x=140 /*x是变量，140是总体均数*/
ttesti 36 130.83 25.74 140 /*已知样本量、均数、标准差和总体均数的简洁命令*/
结果如上图
SPSS代码及结果
鼠标点击：分析—比较平均值—单样本T检验
T-TEST
/TESTVAL = 140
/MISSING = ANALYSIS
/VARIABLES = hb
/CRITERIA = CI(.95).
例3-6 为比较两种方法对乳酸饮料中脂肪含量测定结果是否不同，随机抽取了10份乳酸饮料制品，分别用脂肪酸水解法和哥特里－罗紫法测定其结果如表3-6第(1)~(3)栏。问两法测定结果是否不同？
笔算如下：
SAS代码及结果
MEANS或UNIVARIATE或TTEST，三种过程均对配对资料做两样本均数比较的t检验
/*用MEANS作配对资料两个样本均数比较的t检验*/
data ex3_6;
input x1 x2 @@;
d=x1-x2;
cards;
0.840 0.580
0.591 0.509
0.674 0.500
0.632 0.316
0.687 0.337
0.978 0.517
0.750 0.454
0.730 0.512
1.200 0.997
0.870 0.506
;
proc means t prt;
var d;
run;
/*用UNIVARIATE过程作配对资料两样本均数比较的t检验*/
proc univariate data=ex3_6;
var d;
run;
/*用TTEST过程作配对资料两样本均数比较的t检验*/
proc ttest data=ex3_6;
var d;
run;
R代码及结果
# install.packages("readxl") # 第一运行该文件时加载readxl包用于读取Excel输入数据，将本行开始的#号删除即可
library( readxl ) # 调用外部数据程序包readxl
# 读取Data文件夹下的03_06例题的Excel数据文件（“../Data/E03_06.xls”代表上一级文件夹Data中的E03_06.xls）
E03_06 <- read_excel( "../Data/E03_06.xls" ) # 符号“../” 表示回溯到上一级文件夹
#################### 描述读取的例题数据 ###########################
ls() # 查看工作数据文件名称列表
str(E03_06) # 查看数据结构：no=编号，X1=哥特里－罗紫法，X2=脂肪酸水解法
head(E03_06) # 查看数据前几例
summary(E03_06) # 各列的统计量（未分组）
typeof(E03_06) # E03_05的数据类型
#################### 下面为例题的统计分析 ###########################
#配对样本t检验-检验两种方法的检验结果是否不同
#原假设H0：μd=0
t.test(E03_06$x1-E03_06$x2,mu=0)
t.test(E03_06$x1, E03_06$x2, paired=TRUE)
help( t.test )
结果如下：
> library(readxl)
> E03_06 <- read_excel("F:/Data/E03_06.xls")
> View(E03_06)
> knitr::opts_chunk$set(echo = TRUE)
> ls() # 查看工作数据文件名称列表
[1] "E03_06"
> str(E03_06) # 查看数据结构：no=编号，X1=哥特里－罗紫法，X2=脂肪酸水解法
tibble[,3] [10 x 3] (S3: tbl_df/tbl/data.frame)
$ no: num [1:10] 1 2 3 4 5 6 7 8 9 10
$ x1: num [1:10] 0.84 0.591 0.674 0.632 0.687 0.978 0.75 0.73 1.2 0.87
$ x2: num [1:10] 0.58 0.509 0.5 0.316 0.337 0.517 0.454 0.512 0.997 0.506
> head(E03_06) # 查看数据前几例
> summary(E03_06) # 各列的统计量（未分组）
no x1 x2
Min. : 1.00 Min. :0.5910 Min. :0.3160
1st Qu.: 3.25 1st Qu.:0.6773 1st Qu.:0.4655
Median : 5.50 Median :0.7400 Median :0.5075
Mean : 5.50 Mean :0.7952 Mean :0.5228
3rd Qu.: 7.75 3rd Qu.:0.8625 3rd Qu.:0.5158
Max. :10.00 Max. :1.2000 Max. :0.9970
> typeof(E03_06) # E03_05的数据类型
[1] "list"
> t.test(E03_06$x1-E03_06$x2,mu=0)
One Sample t-test
data: E03_06$x1 - E03_06$x2
t = 7.926, df = 9, p-value = 2.384e-05
alternative hypothesis: true mean is not equal to 0
95 percent confidence interval:
0.1946542 0.3501458
sample estimates:
mean of x
0.2724
> t.test(E03_06$x1, E03_06$x2, paired=TRUE)
Paired t-test
data: E03_06$x1 and E03_06$x2
t = 7.926, df = 9, p-value = 2.384e-05
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.1946542 0.3501458
sample estimates:
mean of the differences
0.2724
STATA代码及结果
ttest x1=x2
SPSS代码及结果
分析—比较平均值----成对样本t检验
T-TEST
PAIRS = x1 WITH x2 (PAIRED)
/CRITERIA = CI(.95)
/MISSING = ANALYSIS.
例3-7 为探讨青藤碱抗兔动脉粥样硬化作用，采用单纯高脂饲料12周喂养方法建立动脉粥样硬化家兔模型。将12周造模成功家兔16只随机等分为2组，模型组给予高饲料100g/d；青藤碱组除给予高饲料100g/d，在饲料中添加青藤碱109mg/(kg/d)。连续喂养3周，实验结束后心脏采血测定心脏血液中高密度脂蛋白，能否认为模型组与青藤碱组中高密度脂蛋白含量的总体均数不同？
本题采用两样本T检验
SAS代码及结果
data ex3_7;
input x @@;
if _n_<9 then c=1;
else c=2;
cards;
0.66 0.76 0.79 0.88 0.78 0.66 0.75 0.88
0.58 0.69 0.59 0.70 0.69 0.68 0.58 0.60
;
proc ttest;
var x;
class c;
run;
R代码及结果
library( readxl ) # 调用外部数据程序包readxl
# 读取Data文件夹下的03_07例题的Excel数据文件（“../Data/E03_07.xls”代表上一级文件夹Data中的E03_07.xls）
E03_07 <- read_excel( "../Data/E03_07.xls" ) # 符号“../” 表示回溯到上一级文件夹
#################### 描述读取的例题数据 ###########################
ls() # 查看工作数据文件名称列表
str(E03_07) # 查看数据结构：x=HDL含量（mmol/L),Group=分组（1=青藤碱组；2=模型组）
head(E03_07) # 查看数据前几例
summary(E03_07) # 各列的统计量（未分组）
typeof(E03_07) # E03_07的数据类型
#################### 下面为例题的统计分析 ###########################
#两独立样本t检验——检验模型组和青藤碱组的HDL含量是否不同
t.test(x~Group,E03_07,var.equal=TRUE)
help( t.test ) # 单样本t检验函数t.test的详细说明
结果如下：
> library(readxl)
> E03_07 <- read_excel("F:/Data/E03_07.xls")
> View(E03_07)
> knitr::opts_chunk$set(echo = TRUE)
> ls() # 查看工作数据文件名称列表
[1] "E03_06" "E03_07"
> str(E03_07) # 查看数据结构：x=HDL含量（mmol/L),Group=分组（1=青藤碱组；2=模型组）
tibble[,2] [16 x 2] (S3: tbl_df/tbl/data.frame)
$ Group: num [1:16] 1 1 1 1 1 1 1 1 2 2 ...
$ x : num [1:16] 0.66 0.76 0.79 0.88 0.78 0.66 0.75 0.88 0.58 0.69 ...
> head(E03_07) # 查看数据前几例
> summary(E03_07) # 各列的统计量（未分组）
Group x
Min. :1.0 Min. :0.5800
1st Qu.:1.0 1st Qu.:0.6450
Median :1.5 Median :0.6900
Mean :1.5 Mean :0.7044
3rd Qu.:2.0 3rd Qu.:0.7650
Max. :2.0 Max. :0.8800
> typeof(E03_07) # E03_07的数据类型
[1] "list"
>
> #################### 下面为例题的统计分析 ###########################
> #两独立样本t检验——检验模型组和青藤碱组的HDL含量是否不同
>
> t.test(x~Group,E03_07,var.equal=TRUE)
Two Sample t-test
data: x by Group
t = 3.689, df = 14, p-value = 0.00243
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
0.05494115 0.20755885
sample estimates:
mean in group 1 mean in group 2
0.77000 0.63875
STATA代码及结果
ttest x,by(group)
代码说明：
ttest x1=x2,unpaird /*原始数据以两个变量x1,x2分别录入*/
ttest x,by(group) /*原始数据以一个变量录入，用变量group指示分组*/
ttesti 8 0.77 0.08 8 0.64 0.06 /*已知两组样本量、均数、标准差时候的简洁命令*/
sdtest x,by(group) /*方差齐性检验*/
ttest x1=x2,unpune[w] /*方差不齐时候的近似t检验*/
SPSS代码及结果
分析–比较平均值–独立样本t检验
T-TEST
GROUPS = grp(1 2)
/MISSING = ANALYSIS
/VARIABLES = x
/CRITERIA = CI(.95).
扫一扫
3290
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
