vlambda博客
学习文章列表

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

#绘图符号及颜色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 <- rnorm(3000,mean=m.pop, sd=sd.pop) #生成平均值为5，方差为2的总体

t <- vector(mode="numeric",length=0) #存放t值x2 <- vector(mode="numeric",length=0) #存放卡方值f <- vector(mode="numeric",length=0) #存放F值

samples <- 4000 #随机取样次数n <- 15 #每个样本含量

#随机取样for(i in c(1:samples)) { y <- sample(pop, n, replace = FALSE, prob = NULL) mean.y <- mean(y) sd.y <- sd(y) t[i] <- (mean.y-m.pop)/sd.y # t=(样本均值-总体均值)/样本标准差  x2[i] <- ((n-1)*var(y))/(sd.pop^2) # x2 = (n-1)样本方差/总体方差}

#画直方图及密度图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) #零假设：服从正态分布；备选假设：不服从正态分布

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")

#卡方分布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 <- rnorm(3000,mean=5, sd=1.5)pop2 <- rnorm(3000,mean=7, sd=2)samples <- 3000 #随机取样次数n1 <- 15 #每个样本含量n2 <- 25 #每个样本含量#随机取样for(i in c(1:samples)) { y1 <- sample(pop1, n1, replace = FALSE, prob = NULL)  y2 <- sample(pop2, n2, replace = FALSE, prob = NULL)   var.y1 <- var(y1)  var.y2 <- 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 <- rnorm(n, 5, 1)

width <- 0.01scal <- seq(-2, 12, by = 0.01)

ll <- length(scal)
den_y <- density(y, n = ll, from = -2, to = 12)  # here defines the range and width
plot(scal, den_y$y, type = "l")

# use numerical integration# for example, you need pr(a < x < b)a <- 2b <- 4
begin <- which.max(den_y$x[den_y$x < a])end <- which.max(den_y$x[den_y$x < b])
prob <- sum(den_y$y[begin:end]*width)
# confirmtrue_prob <- pnorm(4, 5, 1) - pnorm(2, 5, 1) # increase n to 100000, true_prob = prob. In small sample, they are different
x=rnorm(100)plot(density(x))z=density(x)z$xz$y

# 读取工作目录中txt文件Data <- read.table("temporary.txt", header=T) #header=T表示第一行是表头

temp.data <- Data$NMDS1

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图x<-c(1,2,3,4,5)y<-c(1,2,3,4,5)abline(lm(y~x),col=2,lwd=2)

temp.data <- as.data.frame( scale(Data$WT, center=T, scale=T) )colnames(temp.data)<-c("NMDS1")rownames(temp.data) <- c(5:84)

shapiro.test(log10(Data$Eutro))plot(density(Data$Eutro)) #画密度函数

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

sqrtn(2,0.89)

NMDS.T <- as.data.frame(scale(Data, center=T, scale=T)) #中心标准化NMDS.T <- log(Data+100,base=10) #对数转换NMDS.T <- sqrt(Data+100) #平方根plot(NMDS.T[,2])



cwt <- Data$CWTcwt <- Data[,c(1,2,5,6,7,10)]cwt <- cwt[,-2]colnames(cwt)[7] <- "WT"

cwt[,7] <- Data$WTcwt.mean <- mean(cwt)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 <- as.data.frame(tapply(Data$WT,list(Data$Year),mean))colnames(temp.data) <- "WT"
test.data <- temp.data$WTmean.1 <- mean(test.data) ## 温度平均值sd.1 <- sd(test.data) ## 温度的方差max.1 <- max(test.data)min.1 <- min(test.data)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 <- as.data.frame(tapply(Data$NMDS1,list(Data$Year),mean))NMDS2 <- as.data.frame(tapply(Data$NMDS2,list(Data$Year),mean))NMDS1 <- NMDS1[,1]NMDS2 <- NMDS2[,1]

shapiro.test(NMDS1) #零假设：服从正态分布；备选假设：不服从正态分布shapiro.test(NMDS2) #零假设：服从正态分布；备选假设：不服从正态分布

temp.data <- data.frame(value.test=c(NMDS1,NMDS2), #第一个变量，#第二个变量 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 <- rnorm(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 <- serial1[order(serial1)]serial2 <- rev(serial1)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) #零假设：两个平均值相等；备选假设：两个均值不相等。

#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 <- gl(4,20,80,labels=c("A", "B", "C","D"))#gl(n, k, length = n*k, labels = seq_len(n), ordered = FALSE)#n处理数#k重复数

Data$sites <- rep(c(1,3,4,5),20)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 <- aov(WT~treat,data=Data) #零假设：没有差异；备选假设：有差异summary(m1)outlierTest(m1) #离群点分析（需要car包支持）#多重比较TukeyHSD(m1) #Tukey多重比较是比较多的一种plot(TukeyHSD(m1)) #假如置信区间包括0，则说明二者差异不显著。

boxplot(WT~treat,data=Data)tuk <- glht(m1, linfct=mcp(treat="Tukey"))#需要multcomp软件包支持 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：定量m1<-aov(WT~sites+treat,data=Data) #零假设：没有差异；备选假设：有差异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 <- aov(WT~treat*sites, data=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 <- lm(WT~treat*sites,data=Data)anova(m1)

# 读取工作目录中txt文件Data <- read.table("temporary.txt", header=T) #header=T表示第一行是表头



temp.data <- as.data.frame(tapply(Data$WT, list(Data$Year), FUN=mean)) #求均值colnames(temp.data) <- c("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 <- lm(WT~Year, data=temp.data)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 <- c(1:10,10:1,1:10)y <- c(1992:2021)m1 <- lm(x~y)plot(y,x)M.break <- segmented(m1, seg.Z=~y, psi=list(y=c(1993,2000))) plot(M.break,add=TRUE,link=FALSE,lwd=2, xaxt="n")summary(M.break)

#降季节趋势Data <- read.table("a1.txt", header=T) #header=T表示第一行是表头Cal <- ts(Data$SunShine, frequency=12)stl1 <- stl(Cal, "per")plot(stl1)names(stl1)d <- summary(stl1)names(d)d <- as.data.frame(stl1$time.series)$trendplot(d)stl1$time.series#M-K突变检验#根据魏凤英老师书编写代码#20151012根据网友提示修改代码#[email protected]#Mann_Kendall(序列数据,开始年份)Mann_Kendall <- function(timeserial,start){#参数传递为时间序列和开始年份 Mann_Kendall_sub <- function(timeserial){#计算秩 r <- c() s <- c() U <- c()  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] <- s[i]+r[ii] }   U[i] <- 0 U[i] <- (s[i]- (i*(i+1)/4))/sqrt(i*(i-1)*(2*i+5)/72)   }  r[1] <- 0  s[1] <- 0  U[1] <- 0   LST <- list(r = r, s = s, U = U)  return (LST) }  timeserial_rev <- rev(timeserial)  y1 <- Mann_Kendall_sub(timeserial) y2 <- Mann_Kendall_sub(timeserial_rev)  y2$U <- -(rev(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 <- temp.data$WTd<-Mann_Kendall(temperature,1900)

#问题的提出：请分析水温、电导率、营养水平对第一轴的贡献大小#考虑通过多元分析来解决#1、先分析模型，检查共线性#2、检验模型显著性及残差分布状态#3、对比不同变量的贡献
# 方差膨胀因子（VIF）计算函数# Written by Qingyi @20150605# 方差膨胀因子计算公式及定义见：http://en.wikipedia.org/wiki/Variance_inflation_factor# 需要传入拟比较的变量和所在的数据集VIF <- function(variable,data){  var.nums <- length(variable) #获取变量个数 vif.table <- matrix(,nrow=var.nums,ncol=2) #声明一个矩阵，用于返回计算结果   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 <- lm(",text,sep="") text=paste(text,")",sep="") #字符串拼接结束。  eval(parse(text=text)) #执行lm() temp.d <- summary(temp.lm) #取得r2所在数据集  vif.table[i,1] <- variable[i] #标明变量 vif.table[i,2] <- round(1/(1-temp.d$r.squared),2) #给出具体的数字，显示两位小数  }  colnames(vif.table) <- c("Var.", "VIF") #给表格Header vif.table <- as.data.frame(vif.table) #转化为数据框 return(vif.table) #返回值}names(Data)VIF(c("WT","Cond","Eutro"),"Data") #以10为阈值

l.model <- lm(NMDS1~WT+Cond+Eutro,data=Data)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 <- regsubsets(NMDS1~WT+Cond+Eutro,data=Data,nbest=3)plot(leaps.Raw,scale="adjr2")#根据图中的显示选择R最大的组合







# 相对重要性计算函数# 20151007摘自《R语言实战》# 可参考：Johnson(2000, Multivariate Behavioral Research, 35, 1-19)# @Qingyirelweights <- function(fit, ...){ R <- cor(fit$model) nvar <- ncol(R) rxx <- R[2:nvar, 2:nvar] rxy <- R[2:nvar, 1] svd <- eigen(rxx) evec <- svd$vectors ev <- svd$values delta <- diag(sqrt(ev)) lambda <- evec %*% delta %*% t(evec) lambdasq <- lambda^2 beta <- solve(lambda) %*% rxy rsquare <- colSums(beta^2) rawwgt <- lambdasq %*% beta ^2 import <- (rawwgt/rsquare)*100 lbls <- names(fit$model[2:nvar]) rownames(import) <- lbls 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 <- rstudent(fit) 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 <- Data[,c(2,12,15,16)]plot(s.cor)







#第一种方法：通过标准化的系数比较s.cor.T <- as.data.frame(scale(s.cor, center=T, scale=T)) l.model <- lm(NMDS1~WT+Cond+Eutro,data=s.cor.T)summary(l.model)

#第二种方法：通过比较偏相关系数#需要corpcor包支持library("corpcor")xcor <- cor(s.cor)xpcor <- cor2pcor(xcor)xpcor

#第三种方法：通过相对重要指数#相对权重（relative weight）衡量各个因子的重要性###############对所有可能子模型添加一个预测变量引起R平方平均增加量的近似值。#其他方法可以根据相关系数、但是适用于预测变量之间不相关的情况

#还可以通过标准化数据之后回归系数比较relweights(l.model,col="red")#数据标准化之后结果与前面一致



library("relaimpo")calc.relimp(l.model, type="lmg") #需要relaimpo包支持



#广义可加模型（GAM）#对数据分布不做特别要求#使用数据说话#缺点：#晦涩难懂library("mgcv")

gamc <- gam(NMDS1~s(WT)+s(Cond)+s(Eutro),data=Data)summary(gamc)gam.check(gamc)AIC(gamc) #用于选择最优模型plot(gamc,page=1)



d <- as.data.frame(predict(gamc,type="terms"))plot(d$`s(WT)`,type="l",ylim = c(-1,1))abline(h=0,col="gray")lines(d$`s(Eutro)`,type = "1",col="red")