24小时热门版块排行榜    

查看: 1587  |  回复: 20
【奖励】 本帖被评价9次,作者学员762E0L增加金币 8
当前主题已经存档。

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

支持。祝新春快乐!
回复此楼
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
回帖支持 ( 显示支持度最高的前 50 名 )

zwdbordeaux

木虫 (正式写手)


首先祝大家新年快乐!

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

今天说的是如何在统计结果出来之后作表格的时候标注不同的显著水平。
为了明白起见,我从头开始。
通常情况下,如果我们测定了一组数据,如果他们不适合作图的话,我们就需要用表格来给出。给出数据的时候我们并不是直接给出原始的数据,而是对数据进行统计分析,计算实验处理的均值、标准误(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 ]
5楼2010-02-14 03:43:18
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
普通回帖

toltec

铜虫 (小有名气)


★★★★★ 五星级,优秀推荐

支持一下,R是很牛逼的!
2楼2010-02-13 19:19:55
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

estrid

金虫 (正式写手)


★★★★★ 五星级,优秀推荐

虽然不是学农业的,不过对数据处理感兴趣
支持lz
新年快乐
3楼2010-02-13 22:25:22
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zwdbordeaux

木虫 (正式写手)


引用回帖:
Originally posted by toltec at 2010-02-13 19:19:55:
支持一下,R是很牛逼的!

多谢,很高兴看到熟习R的人。

我本来还在合计是否介绍R会有些突兀呢!
4楼2010-02-14 03:28:17
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

★★★★★ 五星级,优秀推荐

好,我正想学习R分析。
6楼2010-02-14 08:06:46
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zwdbordeaux

木虫 (正式写手)


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

2)R包的下载安装(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) #察看这个数据。注意在R中,井号后面的东西均为自己的备忘和解释,并不会运行。
data(cabbages) #上载数据
cabbages #看数据里面的具体内容

#将数据写入到一个excel文件。在预定的文件夹中建立一个excel文件,cabbages.xls
#运用下面命令将数据写道粘贴板上
write.table(cabbages, "clipboard",sep="\t"
#在cab.xls中粘贴,同时把最上面一行向后移动一格,然后把第一列删掉,因为那时自动生成的行名
当然还有其他的写入方法,但是上面介绍的这种非常方便,当数据不是很大的时候可以应用。
请注意,上面的操作指示为了说明如何把数据从R输入到EXCEL中,通常情况下,我们在做实验的时候已经把数据保存在EXCEL中了,那么那个时候要做的事情便是从Excel中读入数据
# 将 excel.xls格式的文件另存为 .csv 格式(文件,另存为,下拉菜单中找逗号分隔符.csv)
#读入.csv格式的文件
cab<-read.csv(file.choose()) #这里给读入的数据命名为“cab”直接写你想用的名字就可以
#这个命令输入之后,会出现一个文件选择框,把你想选的文件点中就可。当然也可以直接指明路线,如果文件cabbages放在e盘data文件夹中,那么可以如下输入:
cab<-read.csv(“E\\data\\cabbages.csv”)
或者将双斜线\\换成单的反斜线/。请注意,R中区分大小写,还有上面命令中的引号,双斜线或单反斜线,以及文件的后缀.csv都不能有丝毫差错,不然就会出现错误。
#察看其中的变量和数据结构
str(cab)
# 可见卷心菜的品种,有两个 c39, 和 c52
#播种的时间,三个,d16,d20,d21
#headwt 卷心菜重量
#VC卷心菜维他命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 ]
7楼2010-02-17 04:32:07
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zwdbordeaux

木虫 (正式写手)


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()
8楼2010-02-17 04:37:53
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zwdbordeaux

木虫 (正式写手)


折线图的制作


前面对数据作了柱状图,柱状图适用于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)
9楼2010-02-17 04:46:25
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zwdbordeaux

木虫 (正式写手)


线性相关分析及作图
此处把两种相关性的分析和作图方法组合到一起,即运用所有观测值和运用均值。
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)
10楼2010-02-17 04:47:18
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zwdbordeaux

木虫 (正式写手)


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

请高手指出不足!
11楼2010-02-17 04:48:12
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zwdbordeaux

木虫 (正式写手)


boxplot 箱图

柱状图

折线图

简单的线性回归

[ Last edited by zwdbordeaux on 2010-2-17 at 05:06 ]
12楼2010-02-17 05:04:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

mythbird

木虫 (著名写手)


★★★★★ 五星级,优秀推荐

真是学术界的活雷锋啊,勇赞!
13楼2010-02-18 04:25:37
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

yangling_yzc

铁杆木虫 (小有名气)


★★★★★ 五星级,优秀推荐

高手高手,学习了
14楼2010-02-20 00:10:40
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

cato8163

银虫 (正式写手)


★★★★★ 五星级,优秀推荐

小可,学过SAS,是用SAS分析完后,再用Excel制图。
小可,想问楼主,与R软件的区别,各自的优缺点,或者R软件有什么特殊之处?
16楼2010-03-03 09:50:10
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zwdbordeaux

木虫 (正式写手)


引用回帖:
Originally posted by cato8163 at 2010-03-03 09:50:10:
小可,学过SAS,是用SAS分析完后,再用Excel制图。
小可,想问楼主,与R软件的区别,各自的优缺点,或者R软件有什么特殊之处?

R是免费的,功能非常强大,简单易学,有一个非常大的团队在不断的扩展.
任何用户都可把自己的东西以包的形式与同行共享.
所有的处理都可以在R中搞定,你不需要再换(从sas--excel)软件了。


当然还有别的区别,我只是通俗的说一下。
17楼2010-03-03 23:26:22
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zwdbordeaux

木虫 (正式写手)


不好意思,最近因为小木虫的故障加上我自己更换了城市,所以没有及时更新.

今天上来冒个泡, 等一切安顿下来再接续我的更新!
18楼2010-03-03 23:28:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

土匪乙

金虫 (小有名气)


★★★ 三星级,支持鼓励

晕乎了……
19楼2010-03-05 09:35:15
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

zwdbordeaux

木虫 (正式写手)


A Quick Guide to Teaching R Programming to Computational Biology Students

推荐一篇文章:

Eglen SJ (2009) A Quick Guide to Teaching R Programming to Computational Biology Students. PLoS Comput Biol 5(8): e1000482. doi:10.1371/journal.pcbi.1000482
免费下载的文章。

下面是我的读书笔记:
引言:为什么要在计算生物学中运用R?
计算生物学家运用R的一个主要原因是因为 Bioconductor 项目 (http://www.bioconductor.org).该项目包含一系列可用于分析基因组数据的R包。这些包已经被广泛用来阐述期刊文章中运用的分析算法。很多计算生物学家将R和Bioconductor作为自己研究过程中的基础工具。
R是一种现代的,功能强大的编程语言,可便于将新的思路与面向对象的特性相接合用于迅速开发软件。
R中丰富的内置函数让它成为大高通量数据分析和统计模拟的理想工具,它所采用的R包系统让他人提供的代码易于分享。
最后R可生成高质量的图表,确保整个研究的过程从模型模拟数据分析到论文发表均可在R中完成。

之后作者开始介绍:
如何教学生运用R
1)课程材料
他们提供了一种课程要点笔记,这个笔记除了让学生入门之外,可以在随后的运用过程中成为一种参考工具。
2)课堂试验
通过解决不同的问题,学习这个软件。
3)用R作图(没有具体的内容,只是给出了一些建议。)
4)可重复性研究
强调了LaTex和Sewave的运用,给出了他们结合使用的一个例子。
5)参考书籍(列出了一些参考的数据,并根据他们的内容作了简要评述)
6)有用网站 (给出了三个有用的网站)
6.1 http://www.rseek.org,这个大家应该都知道了,r搜索引擎
6.2 http://germain.its.maine.edu/~hiebeler/comp/matlabR.html
一个介绍Matlab与R之间函数关系的网站。如果有人已经用matlab了,可以在这个网站找到R中对应的函数。
6.3 http://addictedtor.free.fr/graphiques 图片集合
7)学习R中通常遇到的问题
7.1 语法错误和入门
7.2预先赋值的变量(pre-allocation of variable,不确定是否翻译正确!)
7.3矢量化(vectorization,介绍如何避免循环,有点意思!)
7.4矢量化例子(可以看看)
8)数据类型(作者提到在运用或者选择list, dataframe, matrix时可能的困惑,但是提到只有大量练习才能真正掌握)

最后就是小结了。等到有一定基础之后,推荐去cran和Bioconductor网站上去寻宝。

------------------------
总体感觉:
具体的和细节的东西不是很多,但是把各种材料和资源作了一个梳理,读了之后知道起码该如何下嘴了。值得一读。即使不是生物学的也可以看看。
20楼2010-03-14 20:23:57
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
简单回复
jinyixiao15楼
2010-02-20 22:45   回复  
  谢谢分享!!!!!!!
相关版块跳转 我要订阅楼主 学员762E0L 的主题更新
☆ 无星级 ★ 一星级 ★★★ 三星级 ★★★★★ 五星级
普通表情 高级回复 (可上传附件)
信息提示
请填处理意见