r语言生物统计.docx
《r语言生物统计.docx》由会员分享,可在线阅读,更多相关《r语言生物统计.docx(33页珍藏版)》请在冰豆网上搜索。
![r语言生物统计.docx](https://file1.bdocx.com/fileroot1/2023-1/6/ad3b61da-9a0f-4202-ae43-1503c6aa11fc/ad3b61da-9a0f-4202-ae43-1503c6aa11fc1.gif)
r语言生物统计
R语言概述
R是属于GNU系统的一个自由、免费、源代码开放的软件,它是一个用于统计计算和统计制图的优秀工具。
R是S语言的一种实现。
S语言是由AT&T贝尔实验室开发的一种用来进行数据探索、统计分析、作图的解释型语言。
最初S语言的实现版本主要是S-PLUS。
S-PLUS是一个商业软件,它基于S语言,并由MathSoft公司的统计科学部进一步完善。
后来Auckland大学的RobertGentleman和RossIhaka及其他志愿人员开发了一个R系统。
R的使用与S-PLUS有很多类似之处,两个软件有一定的兼容性。
S-PLUS的使用手册,只要经过不多的修改就能成为R的使用手册。
所以有人说:
R,是S-PLUS的一个“克隆”。
但是请不要忘了:
Risfree。
R是一套完整的数据处理、计算和制图软件系统。
其功能包括:
数据存储和处理系统;数组运算工具(其向量、矩阵运算方面功能尤其强大);完整连贯的统计分析工具;优秀的统计制图功能;简便而强大的编程语言:
可操纵数据的输入和输入,可实现分支、循环,用户可自定义功能。
与其说R是一种统计软件,还不如说R是一种数学计算的环境,因为R并不是仅仅提供若干统计程序、使用者只需指定数据库和若干参数便可进行一个统计分析。
R的思想是:
它可以提供一些集成的统计工具,但更大量的是它提供各种数学计算、统计计算的函数,从而使使用者能灵活机动的进行数据分析,甚至创造出符合需要的新的统计计算方法。
R是一个免费的自由软件,它有UNIX、LINUX、MacOS和WINDOWS版本,都是可以免费下载和使用的。
R的主要是www.r-project.org。
在那儿可以下载到R的安装程序、各种外挂程序和文档。
在R的安装程序中只包含了8个基础模块,其他外在模块可以通过CRAN获得(cran.r-project.org)。
以下简述RFORWINDOWS的安装和使用:
在网址:
cran.r-project.org下可以找到R的各个版本的安装程序和源代码。
点击进入:
Windows(95andlater),再点击:
base,下载SetupR.exe,约18兆,此便是RFORWINDOWS的安装程序。
双击SetupR.exe,按照提示一步步安装即可。
安装完成后,程序会创建R程序组并在桌面上创建R主程序的快捷方式(也可以在安装过程中选择不要创建)。
通过快捷方式运行R,便可调出R的主窗口(如下图1-1)。
类似于许多以编程方式为主要工作方式的软件,R的界面简单而朴素,只有不多的几个菜单和快捷按钮。
快捷按钮下面的窗口便是命令输入窗口,它也是部分运算结果的输出窗口,有些运算结果则会输出在新建的窗口中。
主窗口上方的一些文字是刚运行R时出现的一些说明和指引。
文字下的:
>符号便是R的命令提示符,在其后可输出命令;>后的矩形是光标。
R一般是采用交互方式工作的,在命令提示符后输入命令,回车后便会输出结果。
第一节 数据的输入
注:
在R中,大小写是敏感的,即大小写是不同的。
一、变量数据直接输入
1、数值变量 如果有一个变量,变量名为x,其数据如下:
12,15,46,23,15。
命令语句如下:
x<-c(12,15,46,23,15) 或 c(12,15,46,23,15)->x
其中:
x为变量名;<-与->为赋值符;c( )为向量建立函数,表示把括号中的数据建立为一个向量。
以上命令语句建立了一个数值变量x。
2、字符变量 字符变量的建立与数值变量一致。
字符串使用引号(单、双均可),如:
y<-c(“er”,”sdf”,”eir”,”jk”,”dim”) 或 c(“er”,”sdf”,”eir”, ”jk”,”dim”)->y
将建立字符变量y。
3、逻辑变量 逻辑变量中的元素是:
TRUE(或简写为T)、FALSE(F)、NA(表示缺省)。
请注意必须都是大写。
逻辑变量可以直接输入,如:
z<-c(TRUE,FALSE,F,T,NA) 可建立逻辑变量z。
1可以表示T,0可以表示F。
所以可以直接写成:
z<-c(1,0,0,1,NA) 在计算时,T和F也是当作1和0来用的。
4、其他输入方式
1)利用“:
”
“:
”可以产生连续的整数序列,如 > x<-1:
10 ↙ > x ↙
[1] 1 2 3 4 5 6 7 8 9 10
2)利用seq函数
产生有规律的数列的更一般方法是采用seq函数,它的一般用法是:
seq(下界,上界,间距),如:
> x=seq(1,10,0.5) ↙
> x ↙
[1] 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 7.5 8.0 8.5 9.0 9.5 10.0
3)利用rep函数 rep函数重复第一个变量若干次
> y<-rep(1,5) ↙
> y ↙
[1] 1 1 1 1 1
二、建立数据集
在上面,我们建立起了三个变量:
x,y,z。
如果要建立起一个三变量的数据集,即含有x、y、z的二维表,则可使用data.frame()函数。
当然,在这种情况下,要求变量的长度是一样的。
dd<-data.frame(x,y,z)
便可建立数据集dd。
edit(dd)
可调用R中的数据编辑器显示、编辑数据集dd。
如下图。
如果eidt( )编辑的是一维的向量,如:
edit(x),那么它便会调用记事本来编辑,而不会调用数据编辑器。
如下图。
三、从其他文件中读取数据
在应用统计学中,数据量一般是比较大的,变量也很多。
如果用上述方法来建立数据集,好像辛苦了一些。
上述方法适用于少量数据、少量变量的分析。
对于大量数据和变量,一般应在其他软件中输好,再读R中处理。
1、读入输好的纯文本文件
若数据已经输入一个纯文本文件“c:
/test.txt”。
其格式如下:
Price Floor Area Rooms Age Cent.heat
01 52.00 111.0 830 5 6.2 no
02 54.75 128.0 710 5 7.5 no
03 57.50 101.0 1000 5 4.2 no
04 57.50 131.0 690 6 8.8 no
05 59.75 93.0 900 5 1.9 yes
其中第一行为变量名,第一列为记录序号。
则利用read.table( )函数可读入数据,如:
read.table(“c:
/test.txt”)->test
便把数据集读入,并命名为test。
注:
文件的后缀不必一定要.txt,关键文件要为纯文本,里面不能有特殊格式符。
如果数据文件中没有第一列记录序号,如:
PriceFloorAreaRoomsAgeCent.heat
52.00111.083056.2no
54.75128.071057.5no
57.50101.0100054.2no
57.50131.069068.8no
59.7593.090051.9yes
则命令语句为:
read.table(“c:
/test.txt”,header=TRUE)->test
系统根据每个变量第一个值的类型,自动识别变量类型,如以上数据集中,除最后一个“Cent.hear”是字符变量,其他均为数值变量。
2、读入其他格式的数据库要读入其他格式数据库,必须先安装“foreign”模块。
它不属于R的8个在模块,需要在使用前安装。
安装方法很简便,只需键入命令:
library(foreign)即可。
四、向量运算
1、加减乘除对向量进行加减乘除运算,其含义是对向量的每一个元素进行运算:
例:
>x<-c(1,4,7,8,10)↙
>y<-x+2↙>y↙
[1]3691012
向量之间的运算,则是对应元素之间的运算。
例:
>x<-c(1,4,7,8,10)↙>y<-x+2↙>z<-x+y↙>z↙[1]410161822
2、添加数据与替换数据
添加函数的句法是:
append(向量,添加部分,after=起点)。
例:
>x<-seq(1,15,2)↙
>x↙
[1]13579111315
>append(x,20:
30,after=4)↙
[1]135720212223242526272829309111315
如果after参数缺省,则默认在向量的最后插入。
>x<-seq(1,15,2)↙
>append(x,20:
30)↙
[1]135791113152021222324252627282930
替换数据的句法是:
replace(向量,替换位置,替换值)。
例:
>x<-seq(1,15,2)↙
>replace(x,c(1,4,7),-1)↙
[1]-135-1911-115
第二节统计描述
一、数值变量的统计描述
在R中没有一个直接的统计描述命令,可以把常见的指标都一起算出来。
如概述中所说,R往往是提供一些比较基础的统计命令,需要使用者自己编程组合。
1、和:
sum(x),积:
prod()
例:
>x<-c(12,14,16,18)↙(建立一个数据向量x)
>sum(x)↙(求x的和)
[1]60。
>prod(x)↙
[1]48384
2、算术均数:
mean(x)
例:
>x<-c(12,14,16,18)↙(建立一个数据向量x)
>mean(x)↙(求x的算术平均数)
将输出x的均数15。
但是,算术平均数有一个缺点就是对outlier非常敏感,如果少数几个值非常大,那么整个算术平均数也随之增加。
因此,为了避免这种情况,可以使用trimmedmean——也就是说,先对最大和最小的数剔除了(一般会按照百分比剔除,比如剔除最大和最小的a%的个案剔除出去),然后再计算算术平均值。
mean命令也可以计算trimmdmean,只要增加一个trim参数就可以了。
>c(1,2,3,4,5,6,7,8,9,12,1,31,4,15,16,17)->z1↙
>mean(z1)↙
[1]8.8125>mean(z1,trim=0.1)↙(剔除最大和最小10%的数据)
[1]7.785714
有时候,计算mean时的权重希望不是一样的,那么可以使用weighted.mean命令计算。
比如在评委评分的时候,可能资深评委的评委更重要,给予他们更大的权重就是这种情况。
比如有五位评委,他们的某一次评分分别是2,3,4,5,4,而他们的权重是0.1,0.2,0.2,0.3,0.6那么可以使用以下的命令求出加权平均数。
>weighted.mean(x=c(2,3,4,5,4),w=c(0.1,0.2,0.2,0.3,0.6))↙
[1]3.928571
3、中位数:
median(x)
例:
>y<-c(11,13,14,12,17,10)↙
>median(y)↙
将输出中位数12.5
4、极差:
max(x)-min(x)
例:
>y<-c(11,13,14,12,17,10)↙
>y=max(y)-min(y)↙
将输出结果7。
还有一个命令是range命令,他同时返回最大值和最小值。
>y<-c(11,13,14,12,17,10)↙
>range(y)↙
[1]1017
5、样本方差:
var(x)
例:
>x<-c(9,8,7,10,12,10,11,14,8,9)↙
>var(x)↙
将输出结果4.4。
6、样本标准差:
sd(x)
例:
>x<-c(9,8,7,10,12,10,11,14,8,9)↙
>sd(x)↙
将输出结果2.097618。
请注意都是“样本”。
7、变异系数:
(sd(x)/mean(x))*100
8、四分位差:
IQR(x)
>x<-c(9,8,7,10,12,10,11,14,8,9)↙
>IQR(x)↙
输出结果2.5。
注:
四分位差表示的,是75%分位值和25%分位值的差通过以上命令和程序,可以求出想要的描述指标。
二、正态性检验与方差齐次性检验
1.正态性检验可采用Shapiro-Wilk检验。
程序如下:
>x<-c(9,8,7,10,12,10,11,14,8,9)↙
>shapiro.test(x)↙
输出结果如下:
Shapiro-Wilknormalitytest
data:
x
W=0.9513,p-value=0.6843
2.方差齐次性分析
例:
>y<-c(10,12,14,15,17,29)↙
>x<-c(9,8,7,10,12,10,11,14,8,9)↙
>var.test(x,y)↙
结果如下:
Ftesttocomparetwovariances
data:
xandy
F=0.097,numdf=9,denomdf=5,p-value=0.003206
alternativehypothesis:
trueratioofvariancesisnotequalto1
95percentconfidenceinterval:
0.014516800.43493188
sampleestimates:
ratioofvariances
0.0969875
第三节t检验
一、样本均数与总体均数的假设检验
例:
随机抽测了10只兔的直肠温度,其数据为:
38.7、39.0、38.9、39.6、39.1、39.8、38.5、39.7、39.2、38.4(℃),已知该品种兔直肠温度的总体平均数0µ=39.5(℃),试检验该样本平均温度与0µ是否存在显著差异?
R程序如下:
>x<-c(38.7,39,38.9,39.6,39.1,39.8,38.5,39.7,39.2,38.4)↙(建立变量x)
>t.test(x,mu=39.5,alternative="two.sided",conf.level=0.95)↙
注:
以上程序中,x为数据变量;mu=39.5为总体均数;alternative=“two.sided”表示行双侧检验,如果alternative=“less”,为μ1<μ2的单侧检验,反之,alternative=“greater”,为μ1>μ2的单侧检验;conf.level=0.95指检验水准为0.05。
输出结果如下:
OneSamplet-test
data:
x
t=-2.641,df=9,p-value=0.02687
alternativehypothesis:
truemeanisnotequalto39.5
95percentconfidenceinterval:
38.7388239.44118
sampleestimates:
meanofx
39.09
输出假设检验结果、备择假设、总体均数95%可信区间和样本均数。
从结果中可以看出,t=-2.641,自由度为9,样本均数95%的可信区间为[38.73882,39.44118],假设检验的概率为0.02687,因此,0.01
二、配对资料T检验
例:
11只60日龄的雄鼠在x射线照射前后之体重数据见下表(单位:
g):
检
验雄鼠在照射x射线前后体重差异是否显著?
编号1234567891011
照射前25.724.421.125.226.423.821.522.923.125.129.5
照射后22.523.220.623.425.420.420.621.922.623.524.3
R程序如下:
>x<-c(25.7,24.4,21.1,25.2,26.4,23.8,21.5,22.9,23.1,25.1,29.5)↙
>y<-c(22.5,23.2,20.6,23.4,25.4,20.4,20.6,21.9,22.6,23.5,24.3)↙
>t.test(x,y,alternative="two.sided",paired=TRUE,conf.level=0.95)↙
结果如下:
Pairedt-testdata:
xandyt=4.1334,df=10,p-value=0.002033
alternativehypothesis:
truedifferenceinmeansisnotequalto0
95percentconfidenceinterval:
(注:
95%置信度样本差值的可信区间)
0.8506492.840260
sampleestimates:
meanofthedifferences
1.845455
结果大家可以自己分析。
例2:
某猪场从10窝大白猪的仔猪中,每窝抽出性别相同、体重接近的仔猪2头,将每窝两头仔猪随机地分配到两个饲料组,进行饲料对比试验,试验时间30天,增重结果见下表。
试检验两种饲料喂饲的仔猪平均增重差异是否显著?
窝号12345678910
饲料Ⅰ10.011.212.110.511.19.810.812.512.09.9
饲料Ⅱ9.510.511.89.512.08.89.711.211.09.0
R程序如下:
>x<-c(10,11.2,12.1,10.5,11.1,9.8,10.8,12.5,12.0,9.9)↙
>y<-c(9.5,10.5,11.8,9.5,12,8.8,9.7,11.2,11,9)↙
>t.test(x,y,alternative="two.sided",paired=TRUE,,conf.level=0.95)↙
结果如下:
Pairedt-test
data:
xandy
t=3.4553,df=9,p-value=0.007214
alternativehypothesis:
truedifferenceinmeansisnotequalto0
95percentconfidenceinterval:
0.23826031.1417397
sampleestimates:
meanofthedifferences
0.69
三、两个样本的T检验
例1:
分别测定了10只大耳白家兔、11只青紫蓝家兔在停食18小时后正常血
糖值如下,问该两个品种家兔的正常血糖值是否有显著差异?
(单位:
kg)
大耳白57120101137119117104735368
青紫蓝8936825039325782963188
R程序如下:
>x<-c(57,120,101,137,119,117,104,73,53,68)↙
>y<-c(89,36,82,50,39,32,57,82,96,31,88)↙
首先进行方差齐次性分析。
>var.test(x,y)↙
Ftesttocomparetwovariances
data:
xandy
F=1.3454,numdf=9,denomdf=10,p-value=0.6479
alternativehypothesis:
trueratioofvariancesisnotequalto1
95percentconfidenceinterval:
0.35602665.3330292
sampleestimates:
ratioofvariances
1.345411
从结果中可以看出,方差可以认为是同质的。
>t.test(x,y,alternative="two.sided",paired=FALSE,var.equal=TRUE,conf.level=0.95)↙TwoSamplet-test
data:
xandy
t=2.7179,df=19,p-value=0.01365
alternativehypothesis:
truedifferenceinmeansisnotequalto0
95percentconfidenceinterval:
7.56431458.235686
sampleestimates:
meanofxmeanofy
94.962.0
例2:
有人曾对公雏鸡作了性激素效应试验。
将22只公雏鸡完全随机地分为两组,每组11只。
一组接受性激素A(睾丸激素)处理;另一组接受激素C(雄甾烯醇酮)处理。
在第15天取它们的鸡冠个别称重,所得数据如下:
激素鸡冠重量(mg)
A57120101137119117104735368118
C8930825039225732963188
R程序如下:
>x<-c(57,120,101,137,119,117,104,73,53,68,118)↙
>y<-c(89,30,82,50,39,22,57,32,96,31,88)↙
>var.test(x,y)↙
Ftesttocomparetwovariances
data:
xandy
F=1.0934,numdf=10,denomdf=10,p-value=0.8904
alternativehypothesis:
trueratioofvariancesisnotequalto1
95percentconfidenceinterval:
0.29419014.0641018
sampleestimates:
ratioofvari