r语言生物统计.docx

r语言生物统计.docx

《r语言生物统计.docx》由会员分享，可在线阅读，更多相关《r语言生物统计.docx（33页珍藏版）》请在冰豆网上搜索。

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