R 在植物生态生理学中的应用--从数据管理到图表制作

suncx

好久没了。。。。

很好学习中

qianshanyue

多谢分享　顶一个先

zwdbordeaux

不好意思，准备博后面试忙得不可开交，所以拖了。希望这周有时间更新一下。

piceayhb

期待中！

祝你顺利！

zwdbordeaux

呵呵，我都感觉自己有了弃坑的嫌疑，不过实际情况是没有。去丹麦面试之后得到了那个职位的offer,爽！

echo_july

恭喜楼主～！

用实际的例子来讲确实很好地说。。。不过，我找不到那个cabbages～～～文件

zwdbordeaux

[quote]引用第25楼echo_july于2008-12-06 18:48发表的“”:

恭喜楼主～！

用实际的例子来讲确实很好地说。。。不过，我找不到那个cabbages～～～文件[/quote]

你加载MASS这个包了吗？要首先加载这个包才可以。

library(MASS)

cab<-data.frame(cabbages)

cab

zwdbordeaux

察看了投稿截至日期为12月6日，周末不能上网，周五的时候本来要赶完的，结果有了突发状况，得到另外一个职位的offer，所以浪费了不少时间。不过现在已经完成了。

--------------------------------------------------------------------------------------------

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

-----------------------------------------------------------------------

表格格式有所变动，将就着看。

zwdbordeaux

5.2 柱状图的制作

##请注意，我在2009年4月9日发现了全新的作图方法，不要再浪费时间看着一部分了，请到http://cos.name/bbs/read.php?tid=14828看最简单的，最终极的方法！

R中的柱状图只能针对矩阵（matrix）对象才能进行，这点我很讨厌，很多时候都会忍不住把数据倒入excel，然后用sigmaplot去做。不过如果数据量非常大，也就值得花点时间在R中作一下。首先将前面的数据转换成矩阵格式，注意matrix(data)方程中的data需要是向量（vector）格式。这里把均值和误差值分别转换成矩阵。

matwt<-matrix(data=c(group$Meanwt),nrow=3,ncol=2,

dimnames=list(c("D16","D20","D22"),c("CV39","CV52")),byrow=FALSE)

matwtse<-matrix(data=c(group$SEwt),nrow=3,ncol=2,

dimnames=list(c("D16","D20","D22"),c("CV39","CV52")),byrow=FALSE)

matvc<-matrix(data=c(group$Meanvc),nrow=3,ncol=2,

dimnames=list(c("D16","D20","D22"),c("CV39","CV52")),byrow=FALSE)

matvcse<-matrix(data=c(group$SEvc),nrow=3,ncol=2,

dimnames=list(c("D16","D20","D22"),c("CV39","CV52")),byrow=FALSE)

数据准备好之后终于可以做柱状图了，这里运用gplots这个包里面的barplot2（）这个方程，因为这个方程可以直接把误差添加上，而R默认的barplot方程不具备这种功能。下面对多图组合也进行简要介绍。

#为了将多个图进行组合，往往需要在一张纸上画几个图，这些图很多时候可能共用一个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》一书第三章的解释。

library(gplots) #加载需要的包

par(mfcol=c(2,2),mai=c(0,0,0,0),omi=c(0.6,0.6,0.1,0.1))

#定义柱状图中的误差

ci.l<-matwt-matwtse

ci.u<-matwt+matwtse

#对重量做柱状图，详细的参数涵义参考： ?barplot2()，请特别注意xaxt的应用，它非常有用。

barplot2(matwt,beside=TRUE,col=c("black","gray","white"),legend=rownames(matwtse),plot.ci=TRUE,ci.l=ci.l,ci.u=ci.u,ylim=c(0,5.5),xaxt="n")

box()

#添加显著性差异的字母组合

text(x=c(1.5,2.5,3.5,5.5,6.5,7.5),y=c((group$Meanwt+group$SEwt)*1.1),group$Groupwt)

#添加小图的标志

legend("topleft","A",bty="n")

#对vc含量进行作图

ci.l<-matvc-matvcse

ci.u<-matvc+matvcse

barplot2(matvc,beside=TRUE,col=c("black","gray","white"),plot.ci=TRUE,ci.l=ci.l,ci.u=ci.u,ylim=c(0,79))

box()

text(x=c(1.5,2.5,3.5,5.5,6.5,7.5),y=c((group$Meanvc+group$SEvc)+3),group$Groupvc)

legend("topleft","C",bty="n")

#从图2可以看出，上面的这个柱状图的组合方法主要突出的是播种日期对重量和vc含量的影响。如果想强调品种之间的不同，可以通过对数据的重新组合画另外一种柱状图。

#转置（transpose）前面的矩阵

tmatwt<-t(matwt)

tmatvc<-t(matvc)

#做新的柱状图

barplot2(tmatwt,beside=TRUE,ylim=c(0,5.5),col=c("black","gray"),legend=rownames(tmatwt),xaxt="n",yaxt="n")

box()

legend("topleft","B",bty="n")

barplot2(tmatvc,beside=TRUE,ylim=c(0,79),col=c("black","gray"),yaxt="n")

box()

legend("topleft","D",bty="n")

#添加x轴和y轴的标题。

mtext("品种",side=1,adj=0.5,outer=TRUE,line=2,at=0.25)

mtext("播种日期",side=1,adj=0.5,outer=TRUE,line=2,at=0.75)

mtext("重量(kg)",side=2,at=0.75,adj=0.5,outer=TRUE,line=2)

mtext("维他命含量(mg/kg)",side=2,at=0.25,adj=0.5,outer=TRUE,line=2)

zwdbordeaux

5.3 折线图的制作

##请注意，我在2009年4月9日发现了全新的作图方法，不要再浪费时间看着一部分了，请到http://cos.name/bbs/read.php?tid=14828看最简单的，最终极的方法！前面对数据作了柱状图，柱状图适用于x轴是分类数据的时候，比如本例的品种。当x轴是一系列数据的时候，比如说处理后的时间或者播种后的时间，这个时候更希望得到一种变化的趋势，所以折线图更加合适。为了显示如何做折线图，这个地方把播种日期作为数值来运用。

#在group数据中添加一列播种日期，注意另外一个播种日期列的属性(mode)为factor,也就说他们会被认为是字符。

group$date2<-rep(c(16,20,22),2)

#为了添加误差线，需要运用plotrix这个包中的plotCI这个方程，所以先加载这个包。

library(plotrix)

#计算重量和VC含量的变化范围，这样可以方便定义Y轴的刻度范围（ylim）

ran.wt<-range(cab$HeadWt)

ran.vc<-range(cab$VitC)

#同样的对图进行组合

par(mfrow=c(2,1),mai=c(0,0,0,0),omi=c(0.6,0.6,0.1,0.1))

#首先把图的框架画出来

plot(group$date2,group$Meanwt,type="n",ylim=ran.wt,xaxt="n")

#添加品种39的结果和它们的标准误差

plotCI(x=group[group$Cul=="39","date2"],y=group[group$Cul=="39","Meanwt"],liw=group[group$Cul=="39","SEwt"],uiw=group[group$Cul=="39","SEwt"],pch=2,add=TRUE)

#添加连接的折线

lines(x=group[group$Cul=="39","date2"],y=group[group$Cul=="39","Meanwt"],lty=1,lwd=2)

#添加品种52的数据点和标准误差

plotCI(x=group[group$Cul=="52","date2"],y=group[group$Cul=="52","Meanwt"],liw=group[group$Cul=="52","SEwt"],uiw=group[group$Cul=="52","SEwt"],pch=0,add=TRUE)

#添加折线

lines(x=group[group$Cul=="52","date2"],y=group[group$Cul=="52","Meanwt"],lty=2,lwd=2)

#添加小图标志

text(16,4,"A",font=2)

#添加图例

legend(18,4,legend=c("CV39","CV52"),lty=c(1,2),pch=c(1,0),bty="n")

#对vc含量进行类似操作。读者可以与之前比较看这些命令之间改动的地方，熟悉之后其实可以复制粘铁替换就轻松搞定了。

plot(group$date2,group$Meanvc,type="n",ylim=ran.vc)

plotCI(x=group[group$Cul=="39","date2"],y=group[group$Cul=="39","Meanvc"],liw=group[group$Cul=="39","SEvc"],uiw=group[group$Cul=="39","SEvc"],pch=2,add=TRUE)

lines(x=group[group$Cul=="39","date2"],y=group[group$Cul=="39","Meanvc"],lty=1,lwd=2)

plotCI(x=group[group$Cul=="52","date2"],y=group[group$Cul=="52","Meanvc"],liw=group[group$Cul=="52","SEvc"],uiw=group[group$Cul=="52","SEvc"],pch=0,add=TRUE)

lines(x=group[group$Cul=="52","date2"],y=group[group$Cul=="52","Meanvc"],lty=2,lwd=2)

text(16,80,"B",font=2)

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("维他命含量(mg/kg)",side=2,at=0.25,adj=0.5,outer=TRUE,line=2)

zwdbordeaux

5.3 线性相关分析及作图

此处把两种相关性的分析和作图方法组合到一起，即运用所有观测值和运用均值。

par(mfrow=c(2,1),mai=c(0,0,0,0),omi=c(0.6,0.6,0.1,0.1))

#运用所有观测值作图和分析

plot(cab$HeadWt,cab$VitC,type="n",xaxt="n")

points(cab[cab$Cul=="c39","HeadWt"],cab[cab$Cul=="c39","VitC"],pch=1)

points(cab[cab$Cul=="c52","HeadWt"],cab[cab$Cul=="c52","VitC"],pch=0)

#线性模型的分析

lm1<-lm(VitC~HeadWt,data=cab)

#结果的总结

summary(lm1)

#在图中添加拟合的直线

abline(lm1)

添加方程

text(1.7,47, "y=77.6-7.57x\n(P<0.0001,n=58)")

legend("topright",legend=c("CV39","CV52"),pch=c(1,0),title="A")

#运用均值作图

plot(cab$HeadWt,cab$VitC,type="n")

#添加x轴方向的误差线

plotCI(x=group[group$Cul=="39","Meanwt"],y=group[group$Cul=="39","Meanvc"],err="x",liw=group[group$Cul=="39","SEwt"],uiw=group[group$Cul=="39","SEwt"],pch=2,add=TRUE)

plotCI(x=group[group$Cul=="52","Meanwt"],y=group[group$Cul=="52","Meanvc"],err="x",liw=group[group$Cul=="52","SEwt"],uiw=group[group$Cul=="52","SEwt"],pch=0,add=TRUE)

#添加y轴方向的误差线

plotCI(x=group[group$Cul=="39","Meanwt"],y=group[group$Cul=="39","Meanvc"],err="y",liw=group[group$Cul=="39","SEvc"],uiw=group[group$Cul=="39","SEvc"],pch=2,add=TRUE)

plotCI(x=group[group$Cul=="52","Meanwt"],y=group[group$Cul=="52","Meanvc"],err="y",liw=group[group$Cul=="52","SEvc"],uiw=group[group$Cul=="52","SEvc"],pch=0,add=TRUE)

#对均值作线性拟合

lm2<-lm(Meanvc~Meanwt,data=group)

summary(lm2)

abline(lm2)

text(1.7,47, "y=87.5-11.4x\n(P<0.05,n=58)")

legend("topright",legend="B",bty="n")

mtext("重量(kg)",side=1,at=0.5,adj=0.5,outer=TRUE,line=2)

mtext("维他命含量(mg/kg)",side=2,at=0.5,adj=0.5,outer=TRUE,line=2)

zwdbordeaux

6. 小结

R的功能非常强大，这里仅仅给出了最经常运用到的一些操作。不过相信如果能把这里给出的代码全部弄懂会用，生态生理中的基本问题就都可以轻松解决了。当然，这些给出的代码可能不是最简单和最有效的，给出的另外一个目的也是在于让大家批评指正，共同提高。最后将所有代码整理放在一个.txt文件中，将它们复制到R中或者运用Tinn-R运行之后便可以得到给出的图。

-----------------------------------

需要完整文本和全部代码者请留下邮件。

请高手指出不足！

danielbaggio

楼主辛苦了。

p535fans

看来R还真是不错啊

unangu

小错误：

Cul<-as.factor(cab$Clut)

应为 Cul <- as.factor(cab$Cult)

看到了好文章，忙着试试，发现了点小错误，希望也能做点贡献。

unangu

仔细一看，16楼的东西，好像还有问题。

unangu

很有收获，很想再研究一下，谢谢了。

unangu@gmail.com

cloud_wei

[quote]引用第34楼unangu于2008-12-20 16:43发表的“”:

小错误：

Cul<-as.factor(cab$Clut)

应为 Cul <- as.factor(cab$Cult)

看到了好文章，忙着试试，发现了点小错误，希望也能做点贡献。[/quote]

这个我已经在正式pdf中已经改了，thanks

cloud_wei

[quote]引用第35楼unangu于2008-12-20 16:52发表的“”:

仔细一看，16楼的东西，好像还有问题。[/quote]

还请赐教，我是本文的pdf编辑

zwdbordeaux

[quote]引用第34楼unangu于2008-12-20 16:43发表的“”:

小错误：

Cul<-as.factor(cab$Clut)

应为 Cul <- as.factor(cab$Cult)

看到了好文章，忙着试试，发现了点小错误，希望也能做点贡献。[/quote]

多谢!确实如你指出的!

后来的运行命令文本中已经作了修改,这个地方望了修改了!