使用R来处理生态学数据越来越受到科研工作者的青睐，语义编程风格、漂亮的出图效果，能直接俘获众多用户。本文将生态学数据处理中经常会使用到的功能做个搜集整理。

本文假设读者有一些R的基础知识，对于R的编程规则有所了解。文章内容包括：

一、准备工作

二、数据分布形式

三、数据探索

四、均值比较

五、趋势检验

六、用回归分析不同因子的贡献大小

放一张R基础绘图函数绘制的图镇文

一、准备工作

# 清理内存中所有变量rm(list = ls())# 设置工作目录setwd("E:/BiostatisticsR/2016")# 读取工作目录中txt文件Data "temporary.txt", header=T) #列出表头(变量名)names(Data)#显示前几条记录head(Data)#散点图plot(Data$NMDS1,type="o",pch=1)#创建变量x 1:y 1:

最简单的绘图代码

#绘图符号及颜色plot(x=x,y=y,pch=c(1:25),col=c("red","blue","orange","green"),cex=1.5)points(x=(1:25),y=c(25:1),col="black")text(x=5,y=20,"look at me")#箱式图boxplot(Data$NMDS1~Data$Year)

二、 数据分布形式

#概率密度函数 dnorm(n, mean = 0, sd = 1)   #dbinom(n, mean = 0, sd = 1) #dnorm(1)#dnorm(1, mean=1, sd=1)#dnorm(c(1:10),mean=1,sd=1)#累积密度函数pnorm(n, mean = 0, sd = 1)#pnorm(c(1,2,3,4))#产生序列 seq#seq(min,max,length)#seq(1,10,length=1000)#产生随机数rnorm(n, mean = 0, sd = 1)rnorm(23, mean = 0, sd = 1)#取样函数 sample(x, size, replace = FALSE, prob = NULL)命令是从x中随机抽取size大小的样本for(i in c(1:10)){  print(sample(c(1:10),2))  }#二项式分布：dbinom(x, size, prob, log = FALSE)#已知某批鸡蛋的孵出率prob为0.9，抽取size为5个鸡蛋检查其孵化情况，发现最终x=3个鸡蛋孵化，求二项分布的概率。#R中的求解如下：#dbinom(3,5,0.9,log=F) for(i in c(1:9)){  if(i==1){    plot(c(1:1000),dbinom(c(1:1000),1000,0.25),type="l",col=1)    lines(c(1:1000),dbinom(c(1:1000),2000,0.25),type="l",col=2)  }  else  {lines(c(1:1000),dbinom(c(1:1000),1000,i/10),type="l",col=i)}}# 绘制正态分布图形plot(seq(-5,5,length=1000),     dnorm(seq(-5,5,length=1000),           mean = 0, sd = 0.5),      type="l",col=1)lines(seq(-5,5,length=1000),dnorm(seq(-5,5,length=1000), mean = 0, sd = 1), type="l",col=2)lines(seq(-5,5,length=1000),dnorm(seq(-5,5,length=1000), mean = 0, sd = 2), type="l",col=3)lines(seq(-5,5,length=1000),dnorm(seq(-5,5,length=1000), mean = 0, sd = 3), type="l",col=4)abline(v=-2,col="blue")lines(seq(-5,5,length=1000),dnorm(seq(-5,5,length=1000), mean = -2, sd = 3), type="l",col=4)

测试数据不同分布形式

m.pop 5                       sd.pop 2                      pop 3000,mean=m.pop, sd=sd.pop)  t "numeric",x2 "numeric",f "numeric",samples 4000  n 15          #随机取样for(i in c(1:samples))   {  y pop, n, replace = FALSE, prob = NULL)  mean.y y)  sd.y y)  t[i] # t=(样本均值-总体均值)/样本标准差    x2[i] 1)*var(}#画直方图及密度图par(mfrow=c(1,2))summary(t)hist(t,freq=F,breaks=seq(-1.5,1.5,by=0.005))lines(density(t),col=2) #density(x,bw)函数给出的是一个核密度估计(KDE) 百度百科：《核密度图详解》plot(density(t,bw=0.1))lines(density(t,bw=0.5),col=2)lines(density(t,bw=1),col=3)hist(x2,freq=F)lines(density(x2))#对t和x2做正态性监测shapiro.test(t)shapiro.test(x2)     #零假设：服从正态分布；备选假设：不服从正态分布

#画自由度分别为1:10，300的t分布曲线，和标准正态分布

for(i in c(1:10)){  if(i==1)  {plot(seq(-5,5,length=1000),dt(seq(-5,5,length=1000), df=i, ncp=0), type="l",col=i,ylim=c(0,0.45))}  else{    lines(seq(-5,5,length=1000),dt(seq(-5,5,length=1000), df=i, ncp=0), type="l",col=i)  }}lines(seq(-5,5,length=1000),dt(seq(-5,5,length=1000), df=300, ncp=0), type="l",col=i)lines(seq(-5,5,length=1000),dnorm(seq(-5,5,length=1000), mean = 0, sd = 1), type="l",col=1, lty=3,lwd=3)abline(v=0,col="lightgray")

画自由度分别为1:10，300的t分布曲线，和标准正态分布

#卡方分布rchisq(n,df,ncp=0 )for(i in c(1:10)){  if(i==1)  {plot(seq(0,40,length=1000),dchisq(seq(0,40,length=1000), df=i, ncp=0), type="l",col=i,ylim=c(0,0.5), xlim=c(0,40))}  else{    lines(seq(0,40,length=1000),dchisq(seq(0,40,length=1000), df=i, ncp=0), type="l",col=i)  }}lines(seq(0,40,length=1000),dchisq(seq(0,40,length=1000), df=20, ncp=0), type="l",col=1,lty=3,lwd=3)

#  F分布######随机抽样pop1 3000,mean=pop2 3000,mean=samples 3000   n1 15          n2 25          #随机取样for(i in c(1:samples))   {  y1   y2   var.y1   var.y2   f[i] = (var.y1/1.5)/(var.y2/2)}#画直方图及密度图hist(f,freq=F)lines(density(f))#对t和x2做正态性监测shapiro.test(f)plot(seq(0,4,length=1000),df(seq(0,4,length=1000), df1=4, df2=20, ncp=0), type="l",col=1,ylim=c(0,0.8), xlim=c(0,4))lines(seq(0,4,length=1000),df(seq(0,4,length=1000), df1=20, df2=4, ncp=0), type="l",col=2)n 100y 5, width 0.scal q(-2, 12, by = 0.01)ll length(scal)den_y y, n = ll, from = -plot(scal, den_y$y, type = "l")# use numerical integration# for example, you need pr(a < x < b)a 2b 4begin end prob # confirmtrue_prob 4, x=rnorm(100)plot(density(x))z=density(x)z$xz$y

三、 数据探索

# 读取工作目录中txt文件Data "temporary.txt", header=T) temp.data summary(temp.data)#均值，最小值，中位数，一分位数，三分位数，最大值mean(temp.data,na.rm=T,trim=0)#算术平均数range(temp.data)  #最大值最小值quantile(temp.data,probs=c(0,0.05,0.5,0.75,1))#分位数boxplot(temp.data) #箱式图var(temp.data)   #方差sd(temp.data)    #标准差   开根号plot(density(temp.data)) #画密度函数lines(seq(min(temp.data),max(temp.data),length=1000),#正态分布      dnorm(seq(min(temp.data),max(temp.data),length=1000),            mean = mean(temp.data), sd = sd(temp.data)),       type="l",col=2)rug(temp.data)    #添加样本#正态性检验shapiro.test(temp.data)#零假设：服从正态分布；备选假设：不服从正态分布qqnorm(temp.data)   #画QQ图x1,y1,abline(lm(y~x),col=2,lwd=2)temp.data   scale(Data$WT, center=T, scale=T)  )colnames(temp.data)"NMDS1")rownames(temp.data) 5:shapiro.test(log10(Data$Eutro))plot(density(Data$Eutro)) #画密度函数

自定义开N次方函数

sqrtn function(x, n){  x^(1/n) #开N次方函数}sqrtn(2,0.89)

NMDS.T Data, center=T, scale=T)) NMDS.T Data+NMDS.T Data+plot(NMDS.T[,2])cwt Datacwt 1,cwt 2]colnames(cwt)[7] "WT"cwt[,7] Datacwt.mean cwt.mean^(1/2)plot(NMDS.T,type="l")Mann_Kendall(Data$CWT)cor(Data$CWT,Data$WT)   #相关系数

四、均值比较

# t检验#   1、单样本t检验(one sample t-test)#      已知总体平均值μ，但总体标准差σ未知# 检查前提条件# 正态分布#已知夏天的温度是22度。#人家都说江苏没有春天，只有冬天和夏天#现有太湖1992-2012的春天温度，请验证太湖春天温度和夏天是否一致。。。。#解析：总体均值已知，但方差未知，使用单样本t检验比较均值#准备数据temp.data colnames(temp.data) "WT"test.data mean.1 ## 温度平均值sd.1   ## 温度的方差max.1 min.1 shapiro.test(test.data)     #零假设：服从正态分布；备选假设：不服从正态分布#检验正态分布的直观表现#思想：先画密度柱形图#      后画密度曲线#      添加一条相应的正态分布密度曲线hist(test.data,freq=F) lines(density(test.data),lty=3,lwd=2,col="red")          #画密度函数lines(  seq(min.1,max.1,length=10),                      #x轴  dnorm(                                             #求正态分布下的密度    seq(min.1,max.1,length=10),    mean = mean.1,    sd=sd.1  ),  type="l",  lwd=2,  col="blue"      )t.test(test.data, mu=22, var.equal=FALSE) #零假设：两个平均值相等；备选假设：两个均值不相等。t.test(test.data, mu=30, var.equal=FALSE) #零假设：两个平均值相等；备选假设：两个均值不相等。#mu表示平均值，var.equal = FALSE表示方差未知rm(test.data,max.1, min.1,mean.1,sd.1)#  2、两样本t检验(two-sample t-test)# 检查前提条件# 正态性检验# Bartlett Test方差齐性检验(参数)# 要将两个欲对比数据放在一列。。。。。。。。。。。# 方差齐性分析对离群点非常敏感。NMDS1 NMDS2 NMDS1 1]NMDS2 1]shapiro.test(NMDS1)     #零假设：服从正态分布；备选假设：不服从正态分布shapiro.test(NMDS2)     #零假设：服从正态分布；备选假设：不服从正态分布temp.data #第一个变量，#第二个变量                        group=factor(                                                             c(rep(1,length(NMDS1)),                                  rep(2,20)                                  )                              )                      )bartlett.test(value.test~group,data=temp.data)  #零假设：方差齐。#备选假设：方差不齐rm(temp.data)#独立样本t检验t.test(NMDS1,NMDS2,            #零假设：两个平均值相等；备选假设：两个均值不相等。       var.equal=F,                      #方差齐性检验         alternative="two.sided"           #检验方式：双尾、大于、小于      ) #  3、成对样本t检验#  条件如前t.test(NMDS1,NMDS2, paired=TRUE)    #零假设：两个平均值相等；备选假设：两个均值不相等。#其他案例#成对t检验有意思的地方serial1 1000)hist(serial1,freq=F) lines(density(serial1),lty=3,lwd=2,col="red")          #画密度函数lines(  seq(-3,4,length=1000),                             #x轴  dnorm(                                             #求正态分布下的密度    seq(-3,4,length=1000),    mean = 0,    sd=1  ),  type="l",  lwd=2,  col="blue"      )serial1 serial2 t.test(serial1,serial2,            #零假设：两个平均值相等；备选假设：两个均值不相等。       var.equal=T,                      #方差齐性检验         alternative="two.sided"           #检验方式：双尾、大于、小于    ) t.test(serial1, serial2, var.equal=T, paired=TRUE)t.test(serial1, serial2+2, var.equal=T, paired=TRUE)#  4、非参数t检验# 在无法满足t检验的前提下# 优点：不需要满足任何分布# 缺点：检验功效不强# Wilcoxon秩和检验(Mann-Whitney U检验)wilcox.test(Data$NMDS1,Data$NMDS2, paired=T)            #零假设：两个平均值相等；备选假设：两个均值不相等。

方差分析(ANOVA)

#1、单因素方差分析#用于比较多组均值之间的差异显著性检验#思想：将变异分解为组间和组内。#要求：1、各组观测值服从正态分布或者近似正态分布#      2、各组之间方差具有齐性。#      3、各组之间相互独立。#gl(n,k,length=n*k,labels=1:nordered=FALSE),n 为水平数，k为重复次数.#我们的数据是1992-2012年四个采样点的数据，请比较每五年温度变化是否显著。#解析：每5年一组就有4个处理，1992-96、1997-01、2002-07、2008-2012。t检验已经不适合Data$treat 4,#gl(n, k, length = n*k, labels = seq_len(n), ordered = FALSE)#n处理数#k重复数Data$sites 1,table(Data$treat)                  #计算不同处理的个数tapply(Data$WT,list(Data$treat),mean)tapply(Data$WT,list(Data$treat),sd)bartlett.test(WT~treat, data=Data) #零假设：方差齐qqPlot(lm(WT~treat,data=Data))     #正态性检验：要求所有点均落在95%置信区间内。需要car包支持#否则需要一个一个验证kruskal.test(WT~treat,data=Data)   #方差不齐时使用，秩和检验 null: 有区别  Alternative：没区别m1 data=summary(m1)outlierTest(m1)                    #离群点分析(需要car包支持)#多重比较TukeyHSD(m1)                       #Tukey多重比较是比较多的一种plot(TukeyHSD(m1))                 #假如置信区间包括0，则说明二者差异不显著。boxplot(WT~treat,data=Data)tuk "Tukey"))plot(cld(tuk, level=0.05),col="lightgrey",xlab="处理组", ylab="温度(OC)")#相同字母表示差异不显著t.test(Data$WT[Data$treat=="B"],Data$WT[Data$treat=="C"])#单因素协方差分析#ANCOVA#不同采样点时间不一样，是否会影响到结果呢？#采样时间不一样，属于非人为能控制，因此考虑作为协变量##实现方式1：定量m1data=summary(m1)#实现方式2：直观ancova(WT~sites+treat,data=Data)         #零假设：没有差异；备选假设：有差异。HH包#如果斜率相等，说明该协变量对总体没有影响library(effects)effect("treat",m1)#双因素方差分析#目的：不同点位间、不同年份之间的温度是否有差异？#解析：不同处理是一个因素，不同点位是一个因素#采用双因素方差分析tapply(Data$WT, list(Data$treat, Data$sites), FUN=mean) #求均值tapply(Data$WT, list(Data$treat, Data$sites), FUN=sd)   #求方差bartlett.test(WT~sites, data=Data)                     #零假设：方差齐bartlett.test(WT~treat, data=Data)                     #零假设：方差齐qqPlot(lm(WT~treat*sites,data=Data))                   #正态性检验：要求所有点均落在95%置信区间内m1 data=summary(m1)#差异可视化interaction.plot(Data$treat, Data$sites,                  Data$WT, type="b",                 col=c("red","blue","black","gray"),                 pch=c(16,18,12,5),                 main="比较"                 )interaction2wt(WT~treat*sites,data=Data)  #需要HH包支持#以上方差分析基本基于aov函数#线性回归做方差分析m1 data=anova(m1)

五、趋势检验

# 读取工作目录中txt文件Data "temporary.txt", header=T) temp.data as.data.frame(tapply(Data$WT, colnames(temp.data) "WT")temp.data$Year as.numeric(rownames(temp.data))#M-K检验：Mann-Kendall趋势检验#需要样本遵循一定的分布、不受少数异常值的干扰#适用于水文、气象等非正态分布的数据。MannKendall(temp.data$WT)    #零假设：没有趋势；备择假设：有趋势plot(temp.data$WT,type = "1")#线性趋势#拟合直线plot(temp.data$Year,temp.data$WT,type="o")m1 summary(m1)#平滑函数#不能定量，只能定性plot(temp.data$Year, temp.data$WT, type="o")lines(lowess(x=temp.data$Year, temp.data$WT), lty=4, lwd=2,col="red")#分段回归library("segmented")x 1:y 1992:m1 plot(y,x)M.break list(y=c(plot(M.break,add=TRUE,link=FALSE,lwd=2, xaxt="n")summary(M.break)#降季节趋势Data "a1.txt", header=T) Cal 12)stl1 "per")plot(stl1)names(stl1)d names(d)d as.data.frame(stl1$time.series)$trendplot(d)stl1$time.series#M-K突变检验#根据魏凤英老师书编写代码#20151012根据网友提示修改代码#Qingyi1207@126.com#Mann_Kendall(序列数据,开始年份)Mann_Kendall function(timeserial,start){  Mann_Kendall_sub function(timeserial){    r     s     U         for(i in 2:length(timeserial))    {r[i] 0        for(j in 1:i)    {      if(timeserial[i]>timeserial[j]){r[i]=r[i]+1}    }            s[i] 0    for (ii in 2:i){       s[i]     }            U[i] 0    U[i] 1)/            }        r[1] 0        s[1] 0        U[1] 0            LST list(r = r, s = s, U = U)        return (LST)  }    timeserial_rev     y1   y2     y2$U     LST list(UF=y1,UB=y2)  plot(x=c(start:(start-1+length(timeserial))),y=y1$U,type="l",lwd=2, ylim=c(min(y1$U,y2$U,-1.96),max(y1$U,y2$U,1.96)),xlab="序列",ylab="UF=黑；UB=红")      points(x=c(start:(start-1+length(timeserial))),y=y2$U,type="l",lty=3,col="red",lwd=2)    abline(h=1.969,lty=4,lwd=0.5)  abline(h=-1.96,lty=4,lwd=0.5)      abline(h=0,col="gray",lwd=0.5)  return(LST)  }temperature d1900)

六、用回归分析不同因子的贡献大小

#问题的提出：请分析水温、电导率、营养水平对第一轴的贡献大小#考虑通过多元分析来解决#1、先分析模型，检查共线性#2、检验模型显著性及残差分布状态#3、对比不同变量的贡献# 方差膨胀因子(VIF)计算函数# Written by Qingyi @20150605# 方差膨胀因子计算公式及定义见：http://en.wikipedia.org/wiki/Variance_inflation_factor# 需要传入拟比较的变量和所在的数据集VIF function(variable,data){    var.nums #获取变量个数  vif.table var.nums,ncol=    for(i in 1:var.nums)                        #循环计算各个变量的VIF  {    #以下部分用于对字符串进行拼接，用于执行lm()函数。代码并不是最优。        text=""    for(ii in variable[-i])    {      text=paste(text,ii,sep="+")    }         text=paste(text, ",data=",sep="")    text=paste(text,data,sep="")    text=paste(variable[i],substr(text, 2, nchar(text)),sep="~")          text=paste("temp.lm ,text,sep="")    text=paste(text,")",sep="")    #字符串拼接结束。      eval(parse(text=text))                #执行lm()    temp.d #取得r2所在数据集        vif.table[i,1] #标明变量    vif.table[i,2] 1/(  }    colnames(vif.table) "Var.",   vif.table as.data.frame(vif.table)       return(vif.table)                         #返回值}names(Data)VIF(c("WT","Cond","Eutro"),"Data")  #以10为阈值l.model residplot(l.model)summary(l.model)par(mfrow=c(2,2))plot(l.model)names(l.model)shapiro.test(as.data.frame(l.model$residuals)[,1])     #零假设：服从正态分布；备选假设：不服从正态分布#W = 0.97704, p-value = 0.1593#符合正态分布#全子集回归分析#选择最佳模型#主要是通过调整后的R平方来衡量#需要leaps包支持#全子集回归可用leaps包中的regsubsets()函数实现。你能通过R平方、调整R平方或#Mallows Cp统计量等准则来选择“最佳”模型library(leaps)par(mfrow=c(1,1),oma=c(0.5,0.5,0.5,0.5), mar=c(0.5,0,0,0), cex=0.8, mgp=c(2.3,0.8,0),tck=-0.01, family = 'serif')leaps.Raw 3)plot(leaps.Raw,scale="adjr2")#根据图中的显示选择R最大的组合# 相对重要性计算函数# 20151007摘自《R语言实战》# 可参考：Johnson(2000, Multivariate Behavioral Research, 35,  1-19)# @Qingyirelweights function(fit, ...){  R   nvar   rxx 2:nvar,   rxy 2:nvar,   svd   evec   ev   delta   lambda   lambdasq 2  beta   rsquare 2)  rawwgt 2  import 100  lbls 2:nvar])  rownames(import)   colnames(import) "Weights"  barplot(t(import), names.arg=lbls, ylab="% of R-Square", xlab="Predictor Variables", main="Relative Importance of Predictor Variables", sub=paste("R-Square", round(rsquare, digits=5)),...)  return(import)}names(Data)VIF(c("CWT","WT","ss","TP","Cond","Eutro"),"Data")VIF(c("Cloud","Global","Diffuse","Global"),"Data")##绘制残差函数##摘录自《R语言实战》##Qingyi@20151007residplot function(fit, nbreaks=10){  z   hist(z, breaks=nbreaks, freq=F, xlab="Studentized Residual", main="Distribution of Errors")  rug(jitter(z), col="brown")  curve(dnorm(x,mean=mean(z),sd=sd(z)), add=T, col="blue", lwd=2)    lines(density(z)$x, density(z)$y, col="red", lwd=2, lty=2)  legend("topright", legend=c("Normal Curve", "Kernel Density Curve"), lty=1:2, col=c("blue", "red"), cex=0.7)  }s.cor   2,plot(s.cor)#第一种方法：通过标准化的系数比较s.cor.T as.data.frame(scale(s.cor, center=T, scale=T)) l.model summary(l.model)#第二种方法：通过比较偏相关系数#需要corpcor包支持library("corpcor")xcor xpcor xpcor#第三种方法：通过相对重要指数#相对权重(relative weight)衡量各个因子的重要性###############对所有可能子模型添加一个预测变量引起R平方平均增加量的近似值。#其他方法可以根据相关系数、但是适用于预测变量之间不相关的情况#还可以通过标准化数据之后回归系数比较relweights(l.model,col="red")#数据标准化之后结果与前面一致library("relaimpo")calc.relimp(l.model, type="lmg") #需要relaimpo包支持#广义可加模型(GAM)#对数据分布不做特别要求#使用数据说话#缺点：#晦涩难懂library("mgcv")gamc summary(gamc)gam.check(gamc)AIC(gamc)       #用于选择最优模型plot(gamc,page=1)d as.data.frame(predict(gamc,type=plot(d$`s(WT)`,type="l",ylim = c(-1,1))abline(h=0,col="gray")lines(d$`s(Eutro)`,type = "1",col="red")

#*****************************************************

#    中科院南京地理与湖泊研究所 生物统计与R运用          *

#   Writen By Deng Jianming  qingyi1207@126.com    *

#                  Wrote on 20160907                                   *

#*****************************************************