24小时热门版块排行榜

返回列表

【奖励】本帖被评价9次，作者学员762E0L增加金币 8 个

当前主题已经存档。

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

[资源] ★★★★★ 五星级,优秀推荐

支持。祝新春快乐！

回复此楼

» 猜你喜欢

1楼 2010-02-12 10:05:02

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

zwdbordeaux

木虫 (正式写手)

应助: 0 (幼儿园)
贵宾: 0.086
金币: 1718.4
帖子: 770
在线: 52.6小时
虫号: 493623

给对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 ]