支持一下，R是很牛逼的！

虽然不是学农业的，不过对数据处理感兴趣
支持lz
新年快乐

首先祝大家新年快乐！

呵呵，虽然设置了资源贴，但是没有上“硬菜”，呵呵。所以在这里先写点小的经验！

今天说的是如何在统计结果出来之后作表格的时候标注不同的显著水平。
为了明白起见，我从头开始。
通常情况下，如果我们测定了一组数据，如果他们不适合作图的话，我们就需要用表格来给出。给出数据的时候我们并不是直接给出原始的数据，而是对数据进行统计分析，计算实验处理的均值、标准误（SE），标准偏差（SD），或者变异系数（CV）等，之后对不同处理作方差分析，标注每个试验因子（factor）的显著性，对于有显著性的因子，可以做t检验或者多重比较，得出该因子下那些水平（level）之间会有所差异。

上面介绍了数据分析的大体流程，这里给出当作了方差分析，并且得出了结果之后，不同的因子可能在不同的显著水平（p value）上表现出显著差异。通常在sci发表的过程中，不同的水平会用不同的星号表示“
*** 表示 p<0.001 时显著差异，
** 表示p<0.01 时显著差异
* 表示 p<0.05时显著差异
NS 表示没有显著差异

很多计算程序给出的是p的值，这个时候如何快速的把方差分析得到的p值转变成 相应的符号呢？
我的方法是把p值放在excel的一列中，然后用if 函数搞定

假设p值放在 B 列，那么相应的if 函数是：

=IF(B1<0.001,"***",IF(B1<0.01,"**",IF(B1<0.05,"*","NS")))

之后拖拉一下，就可以把所有的结果弄出来

比如在B列中有数据（如下），之后在C列写入上面的公式，会得到如下的结果
p-value        sig
0.00079        ***
0.03        *
0.04        *
0.00012        ***
0.005        **
0.07        NS
0.4        NS
0.12        NS

[ Last edited by zwdbordeaux on 2010-2-14 at 03:44 ]

好，我正想学习R分析。

给对R感兴趣的同学推荐一个网站
cos:
http://cos.name/bbs/

我在上面曾经发过一个帖子，围绕R应用在生态生理学中的应用入门问题做了介绍。现在转帖到这个地方，后面我会运用同样的数据组，介绍如何在excel 和 sigmaplot中实现相关的功能。
请多批评指教，转帖请注明出处。
http://cos.name/bbs/read.php?tid=12277


开始用R是在2006年初，当时刚刚开始博士论文，因为课题涉及数学模型，导师推荐了R. 其实导师以前都是用S-Plus的，后来发现学生们学会了S-plus之后却买不起软件，所以便开始应用免费版的R了。我适逢其会，开始了R的学习和应用。现在博士即将毕业，对R的应用稍有心得，同时也深切感到这个软件的巨大潜力，实在是科研中的得力助手。但目前这个软件在植物生态生理学上的应用还有待相关学者的关注。写此小文以总结备忘自己的一些经验，同时希望能够吸引更多的研究者运用这个软件。

面对一个新的软件，初用者往往感到无从下手，读了很多说明，结果针对性太差，与真正应用到自己的数据管理，图表制作以及文章发表中距离较大。本文将以零起点用户为对象，运用一个实际例子，解释从数据管理到图表制作的过程，且保证这些得到的图表可以达到用来发表的标准。同时介绍R比其他统计软件SPSS和做图软件Sigmaplot的优缺点。下面将从几个方面逐步介绍：
1）R和Tinn-R的安装
2) 包的下载和安装
3）数据的管理（计算）及输入输出
4）统计分析（方差分析，多重比较，相关性分析）及结果输出
5）图的制作（折线图（点，线，误差），柱状图，回归线添加，方程添加，多图组合）
1）R 和Tinn-R 的安装
本文仅将R和Tinn-R 的应用和安装局限于windows操作系统下，其它操作系统可参考相关文献。

R的安装程序可以到 http://cran.r-project.org/bin/windows/base/ 下载，目前（2008年10月14日）最新的版本为 R-2.7.2-win32.exe, 下载后跟随安装向导，非常容易便可以安装。

Tinn-R是R的一个辅助软件，可以不用，但是建议应用。它可到 http://www.sciviews.org/Tinn-R/ 下载。之后的安装也很简单。

安装完R和 Tinn-R之后需要进行一点小的设置：
打开 R, 之后Edit (窗口上部工具条)--Gui Preference--Rgui Configuration Editor--Single or multiple windows(选中SDI)--Pager style (选中 single window)--Save。 这样R的设置就完成了。
打开Tinn-R, 在顶部窗口 R下拉菜单中找到 Hotkeys of R, 自定义一些热键.
首先在Assigned hotkes 中clear all
然后 Get hotkey method中定义你想要得键组合：比如Ctrl+Key:A
Tinn-R functional hotkeys中定义这个键组合要实现的功能：比如send: line
然后 add
注意底部option中选择 active

这个时候你可以尝试运用R了！

在Tinn-R 中file-New, 在新建的文件中输入： mean(c(4,5,6))
光标放在刚刚输入的那一行，然后Ctrol+A,这个简单的命令便自动输入到R中并且得到结果5。

２）Ｒ包的下载安装(install)，加载(load)及更新（update）

R的一个特点和优点便是它包括很多程序包（Package），这些包内包含各种数据和方程用来处理各种特殊的数据。在R安装的同时一些基本的包已经自己安装并且加载到工作环境中了，但是针对一些特殊的问题，我们难免需要寻求有针对性地程序包。为了获得这些程序包的内容，可以检索：http://pbil.univ-lyon1.fr/library/stats/html/00Index.html
比如要下载安装 plotrix 这个包，到 packages—install packages—select CRAN mirror (比如 Australia)—packages (选择 plotrix)。这样这个包便下载下来了。如果想用这个包，需要到 packages—load packages—plotrix 或者简单的输入 library (plotrix)。更新和加载方法一样。

3）数据的管理（计算）及输入输出
数据的管理及输入输出将主要结合R与Excel的互动，以一个虚拟研究为例，说明从数据录入，变量管理，从Excel到R的读入，以及从R到Excel的输出。同时介绍各种统计量的计算，如均值，标准差，变异系数等。

为简单起见，从MASS中找到了一个卷心菜的数据，这个例子JH Maindonald 的《Using R for data analysis and graphics 》中也用到过。这里将会更加详细的展开。
首先介绍一下这个“试验”：为了测定卷心菜品种和播种日期对其产量和品质（维C含量）的影响，试验采用了2个品种，在3个不同时间进行播种，之后测量了卷心菜的重量和VC的含量，每个组合10次重复。因为不同的播种时期涵盖了不同气候条件的影响，因此可以勉强算是一个植物生态生理学的例子吧。当这个实验出来之后，想要回答的几个问题：1）测定两个品种的产量和品质是否有所差别，若有哪个更好；2）播种日期对产量和品质是否有影响；3）播种日期若有影响，是否具有品种依赖性，也就说是否存在与品种之间的交互作用。

我深信100句解释不如一个例子的道理，所以才选择了这个非常容易便可以得到的数据，可以直接如下获得：
help(cabbages)　＃察看这个数据。注意在Ｒ中，井号后面的东西均为自己的备忘和解释，并不会运行。
data(cabbages)　＃上载数据
cabbages　＃看数据里面的具体内容

＃将数据写入到一个ｅｘｃｅｌ文件。在预定的文件夹中建立一个ｅｘｃｅｌ文件，cabbages．xls
#运用下面命令将数据写道粘贴板上
write.table(cabbages, "clipboard",sep="\t")
＃在ｃａｂ．ｘｌｓ中粘贴，同时把最上面一行向后移动一格，然后把第一列删掉，因为那时自动生成的行名
当然还有其他的写入方法，但是上面介绍的这种非常方便，当数据不是很大的时候可以应用。
请注意，上面的操作指示为了说明如何把数据从Ｒ输入到ＥＸＣＥＬ中，通常情况下，我们在做实验的时候已经把数据保存在ＥＸＣＥＬ中了，那么那个时候要做的事情便是从Excel中读入数据
# 将 excel.xls格式的文件另存为 .csv 格式（文件，另存为，下拉菜单中找逗号分隔符．ｃｓｖ）
#读入.csv格式的文件
cab<-read.csv(file.choose())　＃这里给读入的数据命名为“ｃａｂ”直接写你想用的名字就可以
＃这个命令输入之后，会出现一个文件选择框，把你想选的文件点中就可。当然也可以直接指明路线，如果文件ｃａｂｂａｇｅｓ放在ｅ盘ｄａｔａ文件夹中，那么可以如下输入：
cab<-read.csv(“E\\data\\cabbages.csv”)
或者将双斜线＼＼换成单的反斜线／。请注意，Ｒ中区分大小写，还有上面命令中的引号，双斜线或单反斜线，以及文件的后缀．ｃｓｖ都不能有丝毫差错，不然就会出现错误。
#察看其中的变量和数据结构
str(cab)
# 可见卷心菜的品种，有两个 c39, 和 c52
#播种的时间，三个，d16,d20,d21
#ｈｅａｄｗｔ 卷心菜重量
#VＣ卷心菜维他命C的含量
数据读入之后，便可以开始针对前面提出的科学问题开始选择必要的分析方法了，这里很明显需要进行方差分析。

4)方差分析：

写文章的时候最好把方差齐性以及数据正态性的检验写清楚。

R中 variance homogeneity的检验可以通过至少两种方法获得

1）Levene's test, 这个方法不需要数据正态分布作为前提，是一种robust的方法， R中的函数为levene.test(), 需要加载car包
2）Bartlett's test,这个方法是一个nonrobust的方法，需要数据正态分布， R中函数为bartlett.test(),不需要再额外加载包。

R中normality 的确定可以通过 Shapiro-Wilk test获得,R中函数为 shapiro.test()

下面以我以前运用的数据为例，进行示例:

## Levene's test, which is robust in the sense that it does not rely on the normality assumption.
##levene.test(car)
library(car)
library(MASS)
cab<-data.frame(cabbages)
cab$comtreat<-as.factor(paste(cab$Cult,cab$Date,sep=""))

with(cab,levene.test(HeadWt,comtreat))
# for HeadWt, p<-0.001, indicates that variance inhomogeneity!
## bartlett.test(), which is a nonrobust test with the normality assumption
##the null that the variances in each of the groups (samples) are the same
with(cab,bartlett.test(HeadWt,comtreat))
with(cab,bartlett.test(VitC,comtreat))
## resulted in the same conclusion as derived from the levene's test.
#### normality analysis
## Shapiro-Wilk test
#####shapiro.test(cab[cab$Cult=="c39","HeadWt"]) # normal distribution
with(cab,tapply(HeadWt,comtreat,shapiro.test))  
# at least one dataset is not a normal distribution

当不符合正态性或者方差不齐的时候需要采取特殊的分析方法如非参数估计等，但是这里依然假设数据符合上述假设，进行常规的方差分析和多重比较。

方差分析有两个方程 aov() 和 manova(), 前者适用于单变量数据，后者适用于多变量数据。在本例中观测的指标有两个，即重量和vc含量。因此可以应用manova(), 但此处也给出aov的应用说明。

datam<-cbind(HW=cab$HeadWt,VC=cab$VitC)
Cul<-as.factor(cab$Clut)
Date<-as.factor(cab$Date)
Datam<-cbind(cab$Headwt,cab$)
aov.1<-aov(cab$HeadWt~Cul*Date)
summary(aov.1)

在做方差分析的时候首先确定作用因子，将他们的属性定义为factor, 然后定义作用因子之间的相互关系，a+b表示两个因子单独效应，a*b包括了单独效应和互作效应。从上面的结果我们可以看到，品种和播种日期对重量具有显著影响，且二者有显著互作效应。运用下面的方法可以一次得出两个以上变量的方差分析结果

aov.2<-manova(datam~Cul*Date)
summary.aov(aov.2)

由结果可见，品种和播种日期对VC含量也有显著影响但是二者之间没有互作效应。
5）多重比较
当我们发现某个影响因子对研究指标具有显著影响之后，随之而来的问题便是这种影响是怎样的趋势，那个品种和播种日期组合下的重量最高，那个品种的VC含量更高? 为了回答这些问题，便需要多重比较了。在这里不求过多纠缠多重比较背后的统计学原理，而是仅给出他们常用的方法。

在做多重比较之前，首先看一下整体的结果，这里我给数据添加了一个分类组，就是把品种和播种日期组合在一起。然后对数据作boxplot, 看它们中值，和变化范围。

cab$comtreat<-as.factor(paste(cab$Cult,cab$Date,sep=""))
par(mfrow=c(2,1),mai=c(0,0,0,0),omi=c(0.6,0.6,0.1,0.1))
plot(cab$comtreat,cab$HeadWt,xaxt="n")
legend("topright","A",bty="n")
plot(cab$comtreat,cab$VitC)
legend("topright","B",bty="n")
mtext("品种和播种日期",side=1,adj=0.5,outer=TRUE,line=2)
mtext("重量(kg)",side=2,at=0.75,adj=0.5,outer=TRUE,line=2)
mtext("维他命含量(g/kg)",side=2,at=0.25,adj=0.5,outer=TRUE,line=2)


# 看一下重量和VC含量之间的相关性
plot(cab$HeadWt,cab$VitC)
lm1<-lm(VitC~HeadWt,data=cab)
summary(lm1)
abline(lm1)

Tukey多重比较
虽然R自己有一个Tukey多重比较的方程，但是我发现 agricolae这个包中的HSD.test方程非常好用。这个包是专门针对农业田间试验和植物育种试验而开发的，实用性比较高。
library(agricolae) #加载包
DFerror<-df.residual(aov.1)
MSerror<-deviance(aov.1)/DFerror
cmp.1<-HSD.test(cab$HeadWt,Cul,DFerror,MSerror,group=TRUE)

# 结果表明品种C39的重量显著高于另一品种
#为了得到所有组合之间的比较，我们进行另外一种方差分析：
aov.3<-aov(HeadWt~comtreat,data=cab)
summary(aov.3)
DFerror<-df.residual(aov.3)
MSerror<-deviance(aov.3)/DFerror
cmp.2<-HSD.test(cab$HeadWt,cab$comtreat,DFerror,MSerror,group=TRUE)

#从中提取处理，均值，误差和分组字母，构成一个矩阵一便后面柱状图的操作。
group1<-matrix(c(cmp.2$means,cmp.2$M,cmp.2$std.err),
nrow=6,ncol=3,byrow=TRUE,dimnames=list(cmp.2$trt,c("mean","significant","se")))

barplot(group1[1:6,1],beside=TRUE,angle=c(0,45,90,180,0,0),density=c(-1,20,20,20,-1,-1),ylim=c(0,4.5),col=c(1,2,3,4,5,6))
box()

5. 图表制作

现在多重比较的结果也出来了，剩下的便是怎么把结果呈现出来。对这个数据组，这里尽量运用多种方式，包括三线表、柱状图和折线图。
5.1 结果输出及表的制作
#从中提取重量的均值、误差和分组字母，然后把结果输出到Excel中进行表的制作。
#把重量多重比较的结果提取出来，包括的信息有处理、品均值、标准误差和分组字母。
groupwt<-data.frame(Treat=cmp.2$trt,Cul=substr(cmp.2$trt,2,3),Date=substr(cmp.2$trt,5,6), Meanwt=cmp.2$means,Groupwt=cmp.2$M,SEwt=cmp.2$std.err)
#同样对vc含量作类似操作
aov.4<-aov(VitC~comtreat,data=cab)
summary(aov.4)
DFerror<-df.residual(aov.4)
MSerror<-deviance(aov.4)/DFerror
cmp.3<-HSD.test(cab$VitC,cab$comtreat,DFerror,MSerror,group=TRUE)
groupvc<-data.frame(Treat=cmp.3$trt,Cul=substr(cmp.3$trt,2,3),Date=substr(cmp.3$trt,5,6),
Meanvc=cmp.3$means,Groupvc=cmp.3$M,SEvc=cmp.3$std.err)
#将两组数据合并到一个数据中，当然这个并非必需的，只是展示一下相应的数据操作。
group<-merge(groupwt,groupvc,by=c("Cul","Date"))
#将结果输出到Excel中
write.table(group,"clipboard",sep="\t")
之后经过很容易的操作便可以得到表1的结果。
表1. 品种和播种日期对重量和Vc含量的方差分析。不同字母表示在P<0.05水平上显著差异; ***表示P<0.001; **表示P<0.01; ns表示非显著差异。
------------------------------------------------------------------------------
品种和播种日期   重量   Vc含量
        kg     mg/kg
--------------------------------------------------------------------------------
c39d16   3.18   a   50.3   cd
c39d20   2.8   ab   49.4   d
c39d21   2.74   ab   54.8   bcd
c52d16   2.26   bc   62.5   b
c52d20   3.11   ab   58.9   bc
c52d21   1.47   c   71.8   a
品种   ***   ***
播种日期   ***   ***
品种x播种日期   **   ns
-----------------------------------------------------------------------
表格格式有所变动，将就着看。

[ Last edited by zwdbordeaux on 2010-2-17 at 04:44 ]，

R作图概述

R的作图功能非常强大，且在不断拓展更新之中。但是与sigmaplot或origin等作图软件相比R作图对初学者来说可能比较繁琐，这主要是因为R作图过程中对图中的每个构成要素均需通过特定的函数或者参数来设定，而对函数功能和参数设定的熟悉、理解及至融会贯通需要一定的时间。不过，一旦掌握了其中的窍门，将会大大提供作图的效率，因此值得拿出时间和精力学习一下。


如在其它作图软件中一样，R中的图也是由很多构成元素组成，比如图的外框（box）、图题（title）、图例（legend）、坐标轴（axes)、 坐标轴名称（label)、图内的点（points）、线（line）、文字（text）等等。对于每个元素，又有不同的属性，比如线的线型（实线、虚线、点划线等）。这些元素的刻画和各元素属性的调控可以通过不同的函数和函数中的参数来实现。R中提供了两类作图函数，一类是高水平作图函数，如plot()、matplot()，可以一次性画出前面提到的所有构成要素；另一类为低水平作图函数，如points()、lines()，专门画出或调控图中的某一特定构成元素。

比如我们经常会用的一种作图方法就将数据由符号（symbol）表示，然后用线串联起来。

plot(x=c(1:3),y=c(4,7,6),xlab="x轴名称"，ylab="y轴名称"，main="图题"，type="b")

上面这个命令中，我们运用一个高水平作图函数，通过相应参数对图的x坐标、y坐标以及坐标轴的名称和图的名称作了设定。

如果想在同一个图中画两条以上的曲线，那么可以在前面运行基础上增加：

points(x=c(2,2.5,3),y=c(4,5,6.5),pch=3,col=2)
lines(x=c(2,2.5,3),y=c(4,5,6.5),lty=2,col=4)

请注意，我们在运用各种高级或低级作图函数的时候，运用了一些参数如pch控制符号形状（symbol），lty控制线型、xlab控制x轴的名称等。通常情况下高级和低级作图函数中的参数是通用的，比如pch、lty、lwd等，不过也有特殊的情况。为了弄清各个函数中包含的参数以及各个参数的含义，需要学会善于运用R中的帮助文件。在函数名称前面加上问号，就可以轻松得到函数功能的详细介绍和相关参数含义及取值范围的解释。比如我们在R中输入?plot()，就可以看到关于plot()这个函数的所有信息。此外，在《R Graphics(Murrell,2006)一书中，作者系统解释了各种函数的功能和相应参数的含义及取值范围，所以我们强烈推荐感兴趣的读者详细阅读该书。

对R作图有了基本了解之后，我们正式转入论文发表过程中的运用。显然前面给出的例子距离论文发表过程中的运用还有较大差距。论文发表过程中，我们往往需要在数据点（symbol）上添加数值的观测误差（se）或者标准偏差
（sd），同时还需要将类似的图归纳整合在一起。这就需要我们运用一些特殊的作图函数，并对图中的各个元素进行精细调控。

#为了将多个图进行组合，往往需要在一张纸上画几个图，这些图很多时候可能共用一个x轴或者一个y轴（见图2）。对此R中有一个par()方程进行设定。这个方程中有很多参数，本文中将要重复应用的参数包括mfcol=c(a,b), 定义组合图中有a行b列共a*b个图; mai=c(a,b,c,d)和 omi=c(a,b,c,d)分别定义图的边（margin）和图的外边空间大小。这个空间所对应的图上的区域是以后图形组合的关键，强烈建议参考《R Graphics》一书第三章的解释。
当我们知道了控制图形布局和每个图中元素的调控方法之后，后面表示决定具体该如何进行设置。在此我们建议先把想要最后得到的构架以草图的用笔画出，做到心中有图。比如本文的图2,看到这个图的框架之后，我们便可以知道需要对各个部分作如下设定：

1）图板上需要安排2行2列，共4个图，可通过par（mfrow=c(2,2)）得到；
2）图之间没有间隔，也就说所有图都紧密联接，可通过par(mai=c(rep(0,4)))得到；
3）整个画板上只需要底部和左边留出空间以便显示坐标轴和坐标名称，而顶部和右边不需要留有空间，可通过par(omi=c(0.6,0.6,0,0))获得，注意参数单位为inch。
4）图的类型为柱状图，且带有误差线，可通过sciplot包中的bargraph.CI()得到；
5）A和B两个图不需要x轴，而B和D两个图不需要y轴，可通过高级作图函数中（此处为bargraph.CI()）参数xaxt="n"和/或yaxt="n"得到；
6）x轴和y轴的名称有两个，建议通过mtext(outer=TRUE)加上par(omi=c(0.6,0.6,0,0))获得；
7）每个次级图上均有一个字母（ABCD）且每个柱上有多重比较的显著性字母，可通过text（）得到；
8）只有A和B中具有图例，可以通过legend()得到；

当然还有一些其它参数的精细调控, 在此不再一一解释，不过此处涉及到的参数均反复需要且非常实用，请读者参考R的帮助文件，实践运行下面命令，体会其中的思路。


柱状图的制作

#加载数据
library(MASS)
cab<-data.frame(cabbages)
cab[1,]

##加载这个终极包，人家叫scientific plot，真是没有白叫！注意，没有下载的请先下载！
library(sciplot)
bargraph.CI(Cult, HeadWt, group =Date , data =cab,
           xlab = NA, ylab = NA, cex.lab = 1.5, x.leg = 1,
           col = "black", angle = 45, cex.names = 1.25,
           density = c(0,20,100), legend = TRUE,ylim=c(0,5))
box()

bargraph.CI(Date, HeadWt, group =Cult , data =cab,
           xlab = NA, ylab = NA, cex.lab = 1.5, x.leg = 1,
           col = "black", angle = 45, cex.names = 1.25,
           density = c(0,20), legend = TRUE,ylim=c(0,5))
box()