1. diversity
1.0 Read data of diversity and soil
################################################
setwd('F:\\wangjingSCI\\wangjing20161104\\data')
# dir()
#
require(ade4)
require(vegan)
require(PCNM)
#
dive0 <- read.csv("Alpha20161104.csv")
soil0 <- read.csv("Soil20161104.csv")
dive0 <- droplevels(dive0[dive0$样带!='灌乔带',])
soil0 <- droplevels(soil0[soil0$样带!='灌乔带',])
#
o1 <- order(dive0$编号)
o2 <- order(soil0$编号)
dive.scale <- as.data.frame(scale(dive0[o1,9:13]))
soil.scale <- as.data.frame(scale(soil0[o2,10:25]))
################################################
r.1 <- dive0$replication
dive.1 <- (dive.scale[r.1==1,] + dive.scale[r.1==2,])/2
soil.1 <- (soil.scale[r.1==1,] + soil.scale[r.1==2,])/2
xy0 <- dive0[, c("Dist_upper","Dist_water")]
names(xy0) <- c("x","y")
xy1 <- (xy0[r.1==1,] + xy0[r.1==2,])/2
add.y <- rep(c(-0.105,0,0.105),length=dim(xy1)[1])
xy1$y <- xy1$y+add.y
################################################
#
index <- ! names(soil.1) %in% c("C","C.N","铵态氮","硝态氮","粉粒","黏粒")
#
mite.env <- soil.1[, index]
mite.h <- dive.1
mite.xy <- xy1
#
index_1 <- mite.xy$x==1 & mite.xy$y<1.5
index_2 <- mite.xy$x==6 & mite.xy$y>2.5 & mite.xy$y<3.5
index_3 <- mite.xy$x==9 & mite.xy$y>2.5 & mite.xy$y<3.5
mite.h_add <- mite.h[index_1 | index_2 | index_3,]
mite.env_add <- mite.env[index_1 | index_2 | index_3,]
#
mite.xy_1 <- mite.xy[index_1,]
mite.xy_1$x <- mite.xy_1$x+1
mite.xy_2 <- mite.xy[index_2,]
mite.xy_2$x <- mite.xy_2$x-1
mite.xy_3 <- mite.xy[index_3,]
mite.xy_3$x <- mite.xy_3$x-1
#
mite.xy_add <- rbind(mite.xy_1, mite.xy_2, mite.xy_3)
#
mite.h <- rbind(mite.h, mite.h_add)
mite.env <- rbind(mite.env, mite.env_add)
mite.xy <- rbind(mite.xy, mite.xy_add)
###############################################
#
mite.h <- apply(mite.h, 2, as.numeric)
#
##################################################
#
# # head(mite.xy)
mite.PCNM.vegan <- pcnm(dist(mite.xy))
mite.PCNM <- as.data.frame(mite.PCNM.vegan$vectors)
nb.ev <- which(mite.PCNM.vegan$values > 0.0000001) #
mite.PCNM.pos <- mite.PCNM[, nb.ev]
#
##################################################
# #
# s.value(mite.xy, mite.xy[,"x"])
# s.value(mite.xy, mite.xy[,"y"])
# #
mite_x <- mite.xy
mite_x$y <- 0
PCNM_x.vegan <- pcnm(dist(mite_x))
PCNM_x <- as.data.frame(PCNM_x.vegan$vectors)
nb.ev <- which(PCNM_x.vegan$values > 0.0000001)
#
PCNM_x.pos <- PCNM_x[, nb.ev]
PCNM_x.pos <- apply(PCNM_x.pos, 2, as.numeric)
PCNM_x.pos <- as.matrix(PCNM_x.pos)
# #
# s.value(mite.xy, PCNM_x.pos[,2])
# s.value(mite.xy, PCNM_x[,6])
# #
###################################################
#
mite_y <- mite.xy
mite_y$x <- 0
#
PCNM_y.vegan <- pcnm(dist(mite_y))
PCNM_y <- as.data.frame(PCNM_y.vegan$vectors)
nb.ev <- which(PCNM_y.vegan$values > 0.0000001)
#
PCNM_y.pos <- PCNM_y[, nb.ev]
PCNM_y.pos <- apply(PCNM_y.pos, 2, as.numeric)
PCNM_y.pos <- as.matrix(PCNM_y.pos)
#
# #
# s.value(mite.xy, PCNM_y.pos[,2])
# s.value(mite.xy, PCNM_y.pos[,4])
##################################################
1.1 rda(mite.h, mite_Env.Pcnm)
###################################################
#
mite.h.1 <- as.data.frame(mite.h)
#
# #####################################################
# #
#
# mite.env.rda <- rda(mite.h.1, mite.env)
# #
# (mite.R2a <- RsquareAdj(mite.env.rda)$adj.r.squared)
# #
# (mite.Env.fwd <- forward.sel(mite.h.1, mite.env,
# nperm = 9999,
# adjR2thresh=mite.R2a,
# alpha=0.1))
# #
# (nb.sig.Env <- nrow(mite.Env.fwd)) # Number of
# (Env.sign <- sort(mite.Env.fwd[,2]))
# env.red <- mite.env[,c(Env.sign)]
# #
# # head(env.red)
# #
# ###################################################
#
mite.PCNM.rda <- rda(mite.h.1, mite.PCNM.pos)
#
(mite.R2a <- RsquareAdj(mite.PCNM.rda)$adj.r.squared)
## [1] 0.5820786
#
(mite.PCNM.fwd <- forward.sel(mite.h.1, mite.PCNM.pos,
nperm = 9999,
adjR2thresh=mite.R2a,
alpha=0.1))
## Testing variable 1
## Testing variable 2
## Testing variable 3
## Testing variable 4
## Testing variable 5
## Testing variable 6
## Testing variable 7
## Testing variable 8
## Testing variable 9
## Testing variable 10
## Testing variable 11
## Testing variable 12
## Testing variable 13
## Testing variable 14
## Testing variable 15
## Testing variable 16
## Testing variable 17
## Testing variable 18
## Testing variable 19
## Procedure stopped (adjR2thresh criteria) adjR2cum = 0.589052 with 19 variables (superior to 0.582079)
## variables order R2 R2Cum AdjR2Cum F pval
## 1 PCNM4 4 0.11128122 0.1112812 0.1037497 14.775410 0.0002
## 2 PCNM5 5 0.09635571 0.2076369 0.1940923 14.227844 0.0001
## 3 PCNM13 13 0.09108533 0.2987223 0.2805858 15.066639 0.0001
## 4 PCNM22 22 0.05495417 0.3536764 0.3311956 9.777966 0.0002
## 5 PCNM30 30 0.03429984 0.3879763 0.3611331 6.388938 0.0033
## 6 PCNM18 18 0.02955319 0.4175295 0.3866018 5.733356 0.0056
## 7 PCNM60 60 0.02863123 0.4461607 0.4115457 5.789943 0.0050
## 8 PCNM3 3 0.02827712 0.4744378 0.4365595 5.972196 0.0037
## 9 PCNM23 23 0.02418143 0.4986193 0.4575972 5.305264 0.0071
## 10 PCNM11 11 0.01990916 0.5185284 0.4743567 4.507219 0.0140
## 11 PCNM1 1 0.01920191 0.5377303 0.4906473 4.486139 0.0141
## 12 PCNM12 12 0.01901036 0.5567407 0.5070293 4.588981 0.0117
## 13 PCNM15 15 0.01850129 0.5752420 0.5231490 4.617067 0.0126
## 14 PCNM2 2 0.01810933 0.5933513 0.5391315 4.675975 0.0158
## 15 PCNM48 48 0.01320393 0.6065552 0.5498084 3.490219 0.0368
## 16 PCNM14 14 0.01259457 0.6191498 0.5599886 3.406170 0.0439
## 17 PCNM34 34 0.01256807 0.6317179 0.5703375 3.480873 0.0395
## 18 PCNM10 10 0.01188755 0.6436054 0.5800895 3.368857 0.0402
## 19 PCNM8 8 0.01106035 0.6546658 0.5890522 3.202796 0.0490
#
(nb.sig.PCNM <- nrow(mite.PCNM.fwd)) # Number of
## [1] 19
(PCNM.sign <- sort(mite.PCNM.fwd[,2]))
## [1] 1 2 3 4 5 8 10 11 12 13 14 15 18 22 23 30 34 48 60
PCNM.red <- mite.PCNM.pos[,c(PCNM.sign)]
# head(PCNM.red)
#########################
#
env.red <- mite.env
#
PCNM.env_red <- cbind(PCNM.red, env.red, mite.xy)
PCNM.env_red <- apply(PCNM.env_red, 2, as.numeric)
PCNM.env_red <- as.data.frame(PCNM.env_red )
#
# # head(PCNM.env_red)
#
mite.PCNM.rda2 <- vegan::rda(mite.h.1 ~., data=PCNM.env_red )
axes.test <- anova.cca(mite.PCNM.rda2, by="axis")
(nb.ax <- length(which(axes.test[,4] <= 0.05)))
## [1] 3
#
#########################
#
mod <- mite.PCNM.rda2
#
plot(mod, type="n")
# points(mod, pch=21, col="red", bg="yellow", cex=0.8)
text(mod, "species", col="blue", cex=0.8)
n1 = dim(PCNM.red)[2]
n2 = dim(PCNM.env_red)[2]-n1
text(mod, dis="cn",
col=c(rep(3,n1),rep(2,n2)),
cex=c(rep(0.5,n1),rep(1,n2)),
lty=c(rep(3,n1),rep(1,n2))
)
#
##########################
#
####################################################
1.1_xy
#
mite.PCNM.axes <- scores(mite.PCNM.rda2,
choices=c(1,2),
display="lc",
scaling=1)
#
PCNM.diversity <- cbind(mite.xy,mite.PCNM.axes[,1:2])
# write.csv(PCNM.diversity, "PCNM_axies_diversity.csv")
# ##
s.value(mite.xy, mite.PCNM.axes[,1]) # ade4 function:
#
shapiro.test(resid(lm(mite.PCNM.axes[,1] ~ ., data=mite.env))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 1] ~ ., data = mite.env))
## W = 0.98758, p-value = 0.3449
#
mite.PCNM.axis1.env <- lm(mite.PCNM.axes[,1]~., data=mite.env)
mite.PCNM.axis1.env <- stepAIC(mite.PCNM.axis1.env)
## Start: AIC=-329.83
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全磷 1 0.01391 6.4090 -331.57
## - 有效磷 1 0.02495 6.4201 -331.37
## - 全氮 1 0.02776 6.4229 -331.31
## - 有机质 1 0.03218 6.4273 -331.23
## <none> 6.3951 -329.83
## - 有效氮 1 0.18688 6.5820 -328.38
## - 全钾 1 0.33346 6.7286 -325.74
## - 有效钾 1 0.39094 6.7861 -324.71
## - 含水量 1 0.88867 7.2838 -316.22
## - pH 1 1.18837 7.5835 -311.38
## - 沙粒 1 1.38090 7.7760 -308.37
##
## Step: AIC=-331.57
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全钾 + 有效氮 +
## 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有效磷 1 0.01504 6.4241 -333.29
## - 全氮 1 0.02175 6.4308 -333.17
## - 有机质 1 0.04733 6.4564 -332.69
## <none> 6.4090 -331.57
## - 有效氮 1 0.18312 6.5922 -330.19
## - 全钾 1 0.33529 6.7443 -327.45
## - 有效钾 1 0.40937 6.8184 -326.14
## - 含水量 1 0.93428 7.3433 -317.24
## - pH 1 1.20997 7.6190 -312.82
## - 沙粒 1 1.37676 7.7858 -310.22
##
## Step: AIC=-333.29
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全钾 + 有效氮 +
## 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全氮 1 0.02666 6.4507 -334.80
## - 有机质 1 0.05670 6.4808 -334.24
## <none> 6.4241 -333.29
## - 有效氮 1 0.19686 6.6209 -331.67
## - 全钾 1 0.39191 6.8160 -328.19
## - 有效钾 1 0.44879 6.8729 -327.19
## - 含水量 1 0.94696 7.3710 -318.79
## - pH 1 1.32967 7.7537 -312.72
## - 沙粒 1 1.37069 7.7948 -312.08
##
## Step: AIC=-334.8
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全钾 + 有效氮 +
## 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有机质 1 0.03011 6.4808 -336.24
## <none> 6.4507 -334.80
## - 有效氮 1 0.18722 6.6380 -333.36
## - 全钾 1 0.37001 6.8208 -330.10
## - 有效钾 1 0.44920 6.8999 -328.72
## - 含水量 1 1.04273 7.4935 -318.82
## - 沙粒 1 1.35563 7.8064 -313.91
## - pH 1 1.42913 7.8799 -312.78
##
## Step: AIC=-336.24
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 全钾 + 有效氮 + 有效钾 +
## 沙粒
##
## Df Sum of Sq RSS AIC
## <none> 6.4808 -336.24
## - 有效氮 1 0.24035 6.7212 -333.87
## - 有效钾 1 0.42021 6.9011 -330.70
## - 全钾 1 0.58958 7.0704 -327.79
## - 含水量 1 1.07239 7.5532 -319.86
## - 沙粒 1 1.41995 7.9008 -314.46
## - pH 1 1.66194 8.1428 -310.84
(summ1_1 <- summary(mite.PCNM.axis1.env))
##
## Call:
## lm(formula = mite.PCNM.axes[, 1] ~ 含水量 + pH + 全钾 + 有效氮 +
## 有效钾 + 沙粒, data = mite.env)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50491 -0.17679 -0.01286 0.16602 0.60481
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.002237 0.021888 -0.102 0.91877
## 含水量 -0.113496 0.026247 -4.324 3.32e-05 ***
## pH -0.135918 0.025249 -5.383 4.02e-07 ***
## 全钾 -0.074893 0.023359 -3.206 0.00175 **
## 有效氮 0.062363 0.030464 2.047 0.04296 *
## 有效钾 -0.079051 0.029205 -2.707 0.00785 **
## 沙粒 -0.135862 0.027305 -4.976 2.34e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2395 on 113 degrees of freedom
## Multiple R-squared: 0.5305, Adjusted R-squared: 0.5056
## F-statistic: 21.28 on 6 and 113 DF, p-value: < 2.2e-16
# write.csv(summ1_1$coefficients, "axis.Dive1_env.csv")
# ##
#
shapiro.test(resid(lm(mite.PCNM.axes[,1] ~ ., data=mite.h.1))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 1] ~ ., data = mite.h.1))
## W = 0.99498, p-value = 0.9484
#
# # head(mite.h.1)
mite.PCNM.axis1.div <- lm(mite.PCNM.axes[,1] ~.,data= mite.h.1)
mite.PCNM.axis1.div <- stepAIC(mite.PCNM.axis1.div)
## Start: AIC=-426.29
## mite.PCNM.axes[, 1] ~ OTUs + chao1 + PD + shannon + invsimp
##
## Df Sum of Sq RSS AIC
## - shannon 1 0.00324 3.1147 -428.16
## - OTUs 1 0.01809 3.1296 -427.59
## - chao1 1 0.01811 3.1296 -427.59
## <none> 3.1115 -426.29
## - PD 1 0.21321 3.3247 -420.33
## - invsimp 1 0.37601 3.4875 -414.60
##
## Step: AIC=-428.16
## mite.PCNM.axes[, 1] ~ OTUs + chao1 + PD + invsimp
##
## Df Sum of Sq RSS AIC
## - chao1 1 0.01561 3.1303 -429.56
## - OTUs 1 0.02043 3.1351 -429.38
## <none> 3.1147 -428.16
## - PD 1 0.23425 3.3490 -421.46
## - invsimp 1 0.37383 3.4885 -416.56
##
## Step: AIC=-429.56
## mite.PCNM.axes[, 1] ~ OTUs + PD + invsimp
##
## Df Sum of Sq RSS AIC
## - OTUs 1 0.00959 3.1399 -431.20
## <none> 3.1303 -429.56
## - PD 1 0.25565 3.3860 -422.14
## - invsimp 1 0.54242 3.6727 -412.39
##
## Step: AIC=-431.2
## mite.PCNM.axes[, 1] ~ PD + invsimp
##
## Df Sum of Sq RSS AIC
## <none> 3.1399 -431.20
## - invsimp 1 0.76907 3.9090 -406.91
## - PD 1 1.49864 4.6385 -386.37
(summ1_2 <- summary(mite.PCNM.axis1.div))
##
## Call:
## lm(formula = mite.PCNM.axes[, 1] ~ PD + invsimp, data = mite.h.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.44335 -0.10209 -0.00713 0.10800 0.45789
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.005371 0.014961 0.359 0.72
## PD -0.188308 0.025199 -7.473 1.55e-11 ***
## invsimp -0.134774 0.025176 -5.353 4.37e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1638 on 117 degrees of freedom
## Multiple R-squared: 0.7725, Adjusted R-squared: 0.7686
## F-statistic: 198.7 on 2 and 117 DF, p-value: < 2.2e-16
# write.csv(summ1_2$coefficients, "axis.Dive1_div.csv")
# ##
# ##### ##############
s.value(mite.xy, mite.PCNM.axes[,2]) # ade4 function:
#
shapiro.test(resid(lm(mite.PCNM.axes[,2] ~ ., data=mite.env))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 2] ~ ., data = mite.env))
## W = 0.98632, p-value = 0.2693
#
mite.PCNM.axis2.env <- lm(mite.PCNM.axes[,2] ~ ., data=mite.env)
mite.PCNM.axis2.env <- step(mite.PCNM.axis2.env)
## Start: AIC=-495.62
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有效氮 1 0.00775 1.6142 -497.04
## - pH 1 0.01634 1.6228 -496.40
## - 有机质 1 0.02050 1.6269 -496.09
## - 全钾 1 0.02318 1.6296 -495.90
## <none> 1.6064 -495.62
## - 有效磷 1 0.02881 1.6353 -495.48
## - 全磷 1 0.08391 1.6904 -491.51
## - 沙粒 1 0.09142 1.6979 -490.97
## - 有效钾 1 0.15192 1.7584 -486.77
## - 全氮 1 0.18313 1.7896 -484.66
## - 含水量 1 0.33567 1.9421 -474.85
##
## Step: AIC=-497.04
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - pH 1 0.01837 1.6326 -497.68
## - 有效磷 1 0.02119 1.6354 -497.47
## - 全钾 1 0.02436 1.6385 -497.24
## - 有机质 1 0.02686 1.6410 -497.06
## <none> 1.6142 -497.04
## - 沙粒 1 0.08493 1.6991 -492.89
## - 全磷 1 0.08609 1.7003 -492.80
## - 有效钾 1 0.16798 1.7822 -487.16
## - 全氮 1 0.18274 1.7969 -486.17
## - 含水量 1 0.32801 1.9422 -476.84
##
## Step: AIC=-497.68
## mite.PCNM.axes[, 2] ~ 含水量 + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## <none> 1.6326 -497.68
## - 有效磷 1 0.03198 1.6645 -497.35
## - 全钾 1 0.03991 1.6725 -496.78
## - 有机质 1 0.04253 1.6751 -496.60
## - 沙粒 1 0.06784 1.7004 -494.80
## - 全磷 1 0.08159 1.7142 -493.83
## - 有效钾 1 0.16392 1.7965 -488.20
## - 全氮 1 0.17009 1.8027 -487.79
## - 含水量 1 0.33525 1.9678 -477.27
(summ2_1 <- summary(mite.PCNM.axis2.env))
##
## Call:
## lm(formula = mite.PCNM.axes[, 2] ~ 含水量 + 有机质 + 全氮 + 全磷 +
## 全钾 + 有效磷 + 有效钾 + 沙粒, data = mite.env)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30542 -0.07347 0.01390 0.08480 0.24245
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0007939 0.0110891 -0.072 0.943056
## 含水量 -0.0643870 0.0134862 -4.774 5.54e-06 ***
## 有机质 0.0330319 0.0194254 1.700 0.091847 .
## 全氮 0.0687633 0.0202205 3.401 0.000935 ***
## 全磷 -0.0370144 0.0157153 -2.355 0.020264 *
## 全钾 0.0253402 0.0153835 1.647 0.102339
## 有效磷 -0.0236239 0.0160197 -1.475 0.143130
## 有效钾 0.0507345 0.0151969 3.338 0.001147 **
## 沙粒 0.0280755 0.0130726 2.148 0.033915 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1213 on 111 degrees of freedom
## Multiple R-squared: 0.4063, Adjusted R-squared: 0.3636
## F-statistic: 9.497 on 8 and 111 DF, p-value: 6.208e-10
# ##
# write.csv(summ2_1$coefficients, "axis.Dive2_env.csv")
# # #####
#
shapiro.test(resid(lm(mite.PCNM.axes[,2] ~ ., data=mite.h.1))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 2] ~ ., data = mite.h.1))
## W = 0.99203, p-value = 0.7239
#
mite.PCNM.axis2.div <- lm(mite.PCNM.axes[,2]~., data=mite.h.1)
mite.PCNM.axis2.div <- step(mite.PCNM.axis2.div)
## Start: AIC=-582.55
## mite.PCNM.axes[, 2] ~ OTUs + chao1 + PD + shannon + invsimp
##
## Df Sum of Sq RSS AIC
## - shannon 1 0.00001 0.84613 -584.55
## - OTUs 1 0.00178 0.84790 -584.30
## - invsimp 1 0.00688 0.85300 -583.58
## - PD 1 0.01153 0.85766 -582.93
## <none> 0.84612 -582.55
## - chao1 1 0.93057 1.77669 -495.53
##
## Step: AIC=-584.55
## mite.PCNM.axes[, 2] ~ OTUs + chao1 + PD + invsimp
##
## Df Sum of Sq RSS AIC
## - OTUs 1 0.00392 0.85006 -585.99
## - invsimp 1 0.00695 0.85308 -585.57
## - PD 1 0.01220 0.85834 -584.83
## <none> 0.84613 -584.55
## - chao1 1 1.41276 2.25889 -468.71
##
## Step: AIC=-585.99
## mite.PCNM.axes[, 2] ~ chao1 + PD + invsimp
##
## Df Sum of Sq RSS AIC
## <none> 0.85006 -585.99
## - invsimp 1 0.02032 0.87038 -585.16
## - PD 1 0.08951 0.93956 -575.98
## - chao1 1 1.67243 2.52248 -457.47
(summ2_2 <- summary(mite.PCNM.axis2.div))
##
## Call:
## lm(formula = mite.PCNM.axes[, 2] ~ chao1 + PD + invsimp, data = mite.h.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.200429 -0.057235 -0.003718 0.053067 0.226659
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.004472 0.007822 -0.572 0.568579
## chao1 0.141460 0.009364 15.107 < 2e-16 ***
## PD -0.052387 0.014990 -3.495 0.000673 ***
## invsimp -0.022609 0.013575 -1.665 0.098532 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.0856 on 116 degrees of freedom
## Multiple R-squared: 0.6909, Adjusted R-squared: 0.6829
## F-statistic: 86.42 on 3 and 116 DF, p-value: < 2.2e-16
# ##
# write.csv(summ2_2$coefficients, "axis.Dive2_div.csv")
# #####
#
# windows(title="10 PCNM variables - mites")
# par(mfrow=c(2,5))
# for(i in mite.PCNM.fwd[,2][1:10] ){
# s.value(mite.xy, mite.PCNM.pos[,i], sub=i, csub=2)
# }
#############
1.2 rda(mite.h, mite_Env.Pcnm, mite_x.Pcnm)
###################################################
#
#
mite.h.1 <- residuals(lm(mite.h~PCNM_x.pos))
mite.h.1 <- as.data.frame(mite.h.1)
#
#####################################################
# #
# mite.env.rda <- rda(mite.h.1, mite.env)
# #
# (mite.R2a <- RsquareAdj(mite.env.rda)$adj.r.squared)
# #
# (mite.Env.fwd <- forward.sel(mite.h.1, mite.env,
# nperm = 9999,
# adjR2thresh=mite.R2a,
# alpha=0.1))
# #
# (nb.sig.Env <- nrow(mite.Env.fwd)) # Number of
# (Env.sign <- sort(mite.Env.fwd[,2]))
# env.red <- mite.env[,c(Env.sign)]
# #
# # head(env.red)
#
###################################################
#
mite.PCNM.rda <- rda(mite.h.1, mite.PCNM.pos)
#
(mite.R2a <- RsquareAdj(mite.PCNM.rda)$adj.r.squared)
## [1] 0.5001893
#
(mite.PCNM.fwd <- forward.sel(mite.h.1, mite.PCNM.pos,
nperm = 9999,
adjR2thresh=mite.R2a,
alpha=0.1))
## Testing variable 1
## Testing variable 2
## Testing variable 3
## Testing variable 4
## Testing variable 5
## Testing variable 6
## Testing variable 7
## Testing variable 8
## Testing variable 9
## Testing variable 10
## Testing variable 11
## Testing variable 12
## Testing variable 13
## Testing variable 14
## Procedure stopped (adjR2thresh criteria) adjR2cum = 0.504356 with 14 variables (superior to 0.500189)
## variables order R2 R2Cum AdjR2Cum F pval
## 1 PCNM5 5 0.10798513 0.1079851 0.1004257 14.284790 0.0001
## 2 PCNM13 13 0.10143048 0.2094156 0.1959013 15.010878 0.0002
## 3 PCNM22 22 0.06778928 0.2772049 0.2585119 10.879371 0.0003
## 4 PCNM30 30 0.04224932 0.3194542 0.2957830 7.139375 0.0030
## 5 PCNM18 18 0.04217015 0.3616244 0.3336254 7.530671 0.0019
## 6 PCNM60 60 0.03527422 0.3968986 0.3648755 6.609149 0.0032
## 7 PCNM23 23 0.02960883 0.4265074 0.3906641 5.782444 0.0057
## 8 PCNM12 12 0.02500497 0.4515124 0.4119817 5.060371 0.0110
## 9 PCNM11 11 0.02452846 0.4760408 0.4331714 5.149505 0.0093
## 10 PCNM15 15 0.02279394 0.4988348 0.4528563 4.957525 0.0104
## 11 PCNM24 24 0.01743199 0.5162668 0.4669976 3.891927 0.0267
## 12 PCNM48 48 0.01626749 0.5325342 0.4801082 3.723526 0.0286
## 13 PCNM34 34 0.01537783 0.5479121 0.4924673 3.605604 0.0353
## 14 PCNM10 10 0.01475502 0.5626671 0.5043560 3.542557 0.0350
#
(nb.sig.PCNM <- nrow(mite.PCNM.fwd)) # Number of
## [1] 14
(PCNM.sign <- sort(mite.PCNM.fwd[,2]))
## [1] 5 10 11 12 13 15 18 22 23 24 30 34 48 60
PCNM.red <- mite.PCNM.pos[,c(PCNM.sign)]
#
# # head(PCNM.red)
#########################
#
env.red <- mite.env
#
PCNM.env_red <- cbind(PCNM.red, env.red, mite.xy)
PCNM.env_red <- apply(PCNM.env_red , 2, as.numeric)
PCNM.env_red <- as.data.frame(PCNM.env_red )
#
# # head(PCNM.env_red)
#
mite.PCNM.rda2 <- vegan::rda(mite.h.1 ~., data=PCNM.env_red )
axes.test <- anova.cca(mite.PCNM.rda2, by="axis")
(nb.ax <- length(which(axes.test[,4] <= 0.05)))
## [1] 3
#
#########################
#
mod <- mite.PCNM.rda2
#
plot(mod, type="n")
# points(mod, pch=21, col="red", bg="yellow", cex=0.8)
text(mod, "species", col="blue", cex=0.8)
n1 = dim(PCNM.red)[2]
n2 = dim(PCNM.env_red)[2]-n1
text(mod, dis="cn",
col=c(rep(3,n1),rep(2,n2)),
cex=c(rep(0.5,n1),rep(1,n2)),
lty=c(rep(3,n1),rep(1,n2))
)
#
##########################
#
####################################################
1.2_xy
#################################################
#
mite.PCNM.axes <- scores(mite.PCNM.rda2, choices=c(1,2),
display="lc",
scaling=1)
PCNM.diversity <- cbind(mite.xy, mite.PCNM.axes[,1:2])
# write.csv(PCNM.diversity, "PCNM_axies_diversity.csv")
# ##
s.value(mite.xy, mite.PCNM.axes[,1]) # ade4 function:
###
#
shapiro.test(resid(lm(mite.PCNM.axes[,1] ~ ., data=mite.env)))
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 1] ~ ., data = mite.env))
## W = 0.95325, p-value = 0.0003773
#
mite.PCNM.axis1.env <- lm(mite.PCNM.axes[,1]~., data=mite.env)
mite.PCNM.axis1.env <- stepAIC(mite.PCNM.axis1.env)
## Start: AIC=-356.44
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全磷 1 0.01817 5.1416 -358.01
## - 有效磷 1 0.02166 5.1451 -357.93
## - 有机质 1 0.05580 5.1793 -357.14
## <none> 5.1235 -356.44
## - 有效氮 1 0.11204 5.2355 -355.84
## - 全氮 1 0.11711 5.2406 -355.73
## - 有效钾 1 0.11938 5.2428 -355.68
## - 全钾 1 0.59657 5.7200 -345.22
## - 含水量 1 0.85998 5.9834 -339.82
## - pH 1 1.07409 6.1975 -335.60
## - 沙粒 1 1.30609 6.4295 -331.19
##
## Step: AIC=-358.01
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全钾 + 有效氮 +
## 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有效磷 1 0.01085 5.1525 -359.76
## - 有机质 1 0.07893 5.2206 -358.19
## <none> 5.1416 -358.01
## - 全氮 1 0.10430 5.2459 -357.60
## - 有效氮 1 0.10864 5.2503 -357.51
## - 有效钾 1 0.12988 5.2715 -357.02
## - 全钾 1 0.61562 5.7572 -346.44
## - 含水量 1 0.90885 6.0505 -340.48
## - pH 1 1.09667 6.2383 -336.81
## - 沙粒 1 1.29365 6.4353 -333.08
##
## Step: AIC=-359.76
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全钾 + 有效氮 +
## 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## <none> 5.1525 -359.76
## - 有机质 1 0.08969 5.2422 -359.69
## - 有效氮 1 0.11184 5.2643 -359.18
## - 全氮 1 0.11412 5.2666 -359.13
## - 有效钾 1 0.14693 5.2994 -358.39
## - 全钾 1 0.69410 5.8466 -346.60
## - 含水量 1 0.91973 6.0722 -342.05
## - pH 1 1.19909 6.3516 -336.65
## - 沙粒 1 1.28875 6.4412 -334.97
(summ1_1 <- summary(mite.PCNM.axis1.env))
##
## Call:
## lm(formula = mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 +
## 全钾 + 有效氮 + 有效钾 + 沙粒, data = mite.env)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37543 -0.13561 -0.02085 0.11340 0.76231
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.001731 0.019696 -0.088 0.930135
## 含水量 -0.107220 0.024087 -4.451 2.04e-05 ***
## pH -0.121446 0.023895 -5.083 1.52e-06 ***
## 有机质 0.048410 0.034827 1.390 0.167307
## 全氮 -0.055535 0.035418 -1.568 0.119735
## 全钾 -0.097616 0.025244 -3.867 0.000186 ***
## 有效氮 0.044005 0.028350 1.552 0.123466
## 有效钾 -0.047813 0.026874 -1.779 0.077953 .
## 沙粒 -0.130418 0.024751 -5.269 6.80e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2155 on 111 degrees of freedom
## Multiple R-squared: 0.5706, Adjusted R-squared: 0.5396
## F-statistic: 18.44 on 8 and 111 DF, p-value: < 2.2e-16
# write.csv(summ1_1$coefficients, "axis.Dive1_env.csv")
# #####
#
shapiro.test(resid(lm(mite.PCNM.axes[,1] ~ ., data=mite.h.1))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 1] ~ ., data = mite.h.1))
## W = 0.99412, p-value = 0.899
# # head(mite.h.1)
mite.PCNM.axis1.div <- lm(mite.PCNM.axes[,1]~., data=mite.h.1)
mite.PCNM.axis1.div <- stepAIC(mite.PCNM.axis1.div)
## Start: AIC=-429.91
## mite.PCNM.axes[, 1] ~ OTUs + chao1 + PD + shannon + invsimp
##
## Df Sum of Sq RSS AIC
## - shannon 1 0.001476 3.0204 -431.85
## - OTUs 1 0.002154 3.0211 -431.82
## <none> 3.0190 -429.91
## - chao1 1 0.099763 3.1187 -428.01
## - invsimp 1 0.219968 3.2389 -423.47
## - PD 1 0.314572 3.3335 -420.01
##
## Step: AIC=-431.85
## mite.PCNM.axes[, 1] ~ OTUs + chao1 + PD + invsimp
##
## Df Sum of Sq RSS AIC
## - OTUs 1 0.01305 3.0335 -433.33
## <none> 3.0204 -431.85
## - chao1 1 0.14161 3.1621 -428.35
## - invsimp 1 0.23781 3.2583 -424.76
## - PD 1 0.31863 3.3391 -421.82
##
## Step: AIC=-433.33
## mite.PCNM.axes[, 1] ~ chao1 + PD + invsimp
##
## Df Sum of Sq RSS AIC
## <none> 3.0335 -433.33
## - chao1 1 0.13127 3.1648 -430.25
## - invsimp 1 0.47175 3.5052 -417.99
## - PD 1 1.51337 4.5469 -386.77
(summ1_2 <- summary(mite.PCNM.axis1.div))
##
## Call:
## lm(formula = mite.PCNM.axes[, 1] ~ chao1 + PD + invsimp, data = mite.h.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.47025 -0.09903 -0.01313 0.11006 0.45612
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.884e-17 1.476e-02 0.000 1.000
## chao1 4.878e-02 2.177e-02 2.240 0.027 *
## PD -2.218e-01 2.915e-02 -7.607 8.05e-12 ***
## invsimp -1.128e-01 2.655e-02 -4.247 4.39e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1617 on 116 degrees of freedom
## Multiple R-squared: 0.7472, Adjusted R-squared: 0.7406
## F-statistic: 114.3 on 3 and 116 DF, p-value: < 2.2e-16
# write.csv(summ1_2$coefficients, "axis.Dive1_div.csv")
# #####
# #####
s.value(mite.xy, mite.PCNM.axes[,2]) # ade4 function:
#
shapiro.test(resid(lm(mite.PCNM.axes[,2] ~ ., data=mite.env))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 2] ~ ., data = mite.env))
## W = 0.97971, p-value = 0.06694
#
mite.PCNM.axis2.env <- lm(mite.PCNM.axes[,2] ~ ., data=mite.env)
mite.PCNM.axis2.env <- step(mite.PCNM.axis2.env)
## Start: AIC=-596.29
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有机质 1 0.000619 0.69486 -598.18
## <none> 0.69424 -596.29
## - 有效氮 1 0.021990 0.71623 -594.55
## - pH 1 0.027890 0.72213 -593.57
## - 全磷 1 0.050438 0.74468 -589.88
## - 有效磷 1 0.057140 0.75138 -588.80
## - 全氮 1 0.077143 0.77138 -585.65
## - 沙粒 1 0.115923 0.81016 -579.76
## - 有效钾 1 0.174587 0.86883 -571.37
## - 全钾 1 0.208788 0.90303 -566.74
## - 含水量 1 0.268884 0.96312 -559.01
##
## Step: AIC=-598.18
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮 + 全磷 + 全钾 + 有效氮 +
## 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## <none> 0.69486 -598.18
## - 有效氮 1 0.021392 0.71625 -596.55
## - pH 1 0.027503 0.72236 -595.53
## - 全磷 1 0.057259 0.75212 -590.68
## - 有效磷 1 0.057462 0.75232 -590.65
## - 全氮 1 0.105085 0.79994 -583.28
## - 沙粒 1 0.118863 0.81372 -581.24
## - 有效钾 1 0.177347 0.87221 -572.91
## - 全钾 1 0.215991 0.91085 -567.70
## - 含水量 1 0.269815 0.96467 -560.81
(summ2_1 <- summary(mite.PCNM.axis2.env))
##
## Call:
## lm(formula = mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮 + 全磷 +
## 全钾 + 有效氮 + 有效磷 + 有效钾 + 沙粒, data = mite.env)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.145877 -0.063583 -0.000654 0.051326 0.172787
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0007526 0.0072670 -0.104 0.91771
## 含水量 -0.0581936 0.0089042 -6.536 2.05e-09 ***
## pH -0.0183292 0.0087842 -2.087 0.03923 *
## 全氮 0.0443373 0.0108705 4.079 8.60e-05 ***
## 全磷 -0.0300214 0.0099715 -3.011 0.00323 **
## 全钾 0.0606706 0.0103756 5.847 5.20e-08 ***
## 有效氮 0.0236627 0.0128585 1.840 0.06843 .
## 有效磷 -0.0385160 0.0127704 -3.016 0.00318 **
## 有效钾 0.0525032 0.0099089 5.299 6.06e-07 ***
## 沙粒 0.0400927 0.0092426 4.338 3.20e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07948 on 110 degrees of freedom
## Multiple R-squared: 0.5392, Adjusted R-squared: 0.5015
## F-statistic: 14.3 on 9 and 110 DF, p-value: 4.609e-15
# ##
# write.csv(summ2_1$coefficients, "axis.Dive2_env.csv")
# #####
shapiro.test(resid(lm(mite.PCNM.axes[,2] ~ ., data=mite.h.1))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 2] ~ ., data = mite.h.1))
## W = 0.98885, p-value = 0.4373
mite.PCNM.axis2.div <- lm(mite.PCNM.axes[,2]~., data=mite.h.1)
mite.PCNM.axis2.div <- step(mite.PCNM.axis2.div)
## Start: AIC=-591.91
## mite.PCNM.axes[, 2] ~ OTUs + chao1 + PD + shannon + invsimp
##
## Df Sum of Sq RSS AIC
## - PD 1 0.00023 0.78288 -593.87
## - shannon 1 0.00481 0.78746 -593.17
## - OTUs 1 0.00664 0.78929 -592.89
## <none> 0.78265 -591.91
## - invsimp 1 0.02013 0.80278 -590.86
## - chao1 1 0.41268 1.19533 -543.09
##
## Step: AIC=-593.87
## mite.PCNM.axes[, 2] ~ OTUs + chao1 + shannon + invsimp
##
## Df Sum of Sq RSS AIC
## - shannon 1 0.00458 0.78746 -595.17
## - OTUs 1 0.00973 0.79261 -594.39
## <none> 0.78288 -593.87
## - invsimp 1 0.02088 0.80376 -592.71
## - chao1 1 0.42535 1.20823 -543.80
##
## Step: AIC=-595.17
## mite.PCNM.axes[, 2] ~ OTUs + chao1 + invsimp
##
## Df Sum of Sq RSS AIC
## - OTUs 1 0.00652 0.79398 -596.18
## <none> 0.78746 -595.17
## - invsimp 1 0.01749 0.80495 -594.54
## - chao1 1 0.48072 1.26818 -539.99
##
## Step: AIC=-596.18
## mite.PCNM.axes[, 2] ~ chao1 + invsimp
##
## Df Sum of Sq RSS AIC
## <none> 0.79398 -596.18
## - invsimp 1 0.15161 0.94558 -577.21
## - chao1 1 0.66318 1.45716 -525.32
(summ2_2 <- summary(mite.PCNM.axis2.div))
##
## Call:
## lm(formula = mite.PCNM.axes[, 2] ~ chao1 + invsimp, data = mite.h.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.175489 -0.064017 -0.006106 0.060773 0.229363
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.419e-18 7.520e-03 0.000 1
## chao1 9.525e-02 9.635e-03 9.886 < 2e-16 ***
## invsimp -4.038e-02 8.543e-03 -4.727 6.42e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.08238 on 117 degrees of freedom
## Multiple R-squared: 0.4735, Adjusted R-squared: 0.4645
## F-statistic: 52.61 on 2 and 117 DF, p-value: < 2.2e-16
# ##
# write.csv(summ2_2$coefficients, "axis.Dive2_div.csv")
# #####
#
# windows(title="10 PCNM variables - mites")
# par(mfrow=c(2,5))
# for(i in mite.PCNM.fwd[,2][1:10] ){
# s.value(mite.xy, mite.PCNM.pos[,i], sub=i, csub=2)
# }
#############
1.3 rda(mite.h, mite_Env.Pcnm, mite_y.Pcnm)
###################################################
#
#
mite.h.1 <- residuals(lm(mite.h~PCNM_y.pos))
mite.h.1 <- as.data.frame(mite.h.1)
#
#####################################################
# #
# mite.env.rda <- rda(mite.h.1, mite.env)
# #
# (mite.R2a <- RsquareAdj(mite.env.rda)$adj.r.squared)
# #
# (mite.Env.fwd <- forward.sel(mite.h.1, mite.env,
# nperm = 9999,
# adjR2thresh=mite.R2a,
# alpha=0.1))
# #
# (nb.sig.Env <- nrow(mite.Env.fwd)) # Number of
# (Env.sign <- sort(mite.Env.fwd[,2]))
# env.red <- mite.env[,c(Env.sign)]
# #
# # head(env.red)
#
###################################################
#
mite.PCNM.rda <- rda(mite.h.1, mite.PCNM.pos)
#
(mite.R2a <- RsquareAdj(mite.PCNM.rda)$adj.r.squared)
## [1] 0.4121979
#
(mite.PCNM.fwd <- forward.sel(mite.h.1, mite.PCNM.pos,
nperm = 9999,
adjR2thresh=mite.R2a,
alpha=0.1))
## Testing variable 1
## Testing variable 2
## Testing variable 3
## Testing variable 4
## Testing variable 5
## Testing variable 6
## Testing variable 7
## Testing variable 8
## Testing variable 9
## Testing variable 10
## Testing variable 11
## Procedure stopped (adjR2thresh criteria) adjR2cum = 0.422462 with 11 variables (superior to 0.412198)
## variables order R2 R2Cum AdjR2Cum F pval
## 1 PCNM4 4 0.17869369 0.1786937 0.1717335 25.673558 0.0001
## 2 PCNM60 60 0.04099610 0.2196898 0.2063512 6.146970 0.0034
## 3 PCNM3 3 0.04048906 0.2601788 0.2410455 6.348468 0.0030
## 4 PCNM23 23 0.03988040 0.3000592 0.2757135 6.552334 0.0029
## 5 PCNM11 11 0.02850725 0.3285665 0.2991177 4.840132 0.0116
## 6 PCNM2 2 0.02796815 0.3565346 0.3223683 4.911532 0.0094
## 7 PCNM1 1 0.02749457 0.3840292 0.3455310 4.999249 0.0077
## 8 PCNM15 15 0.02649137 0.4105206 0.3680356 4.988372 0.0091
## 9 PCNM14 14 0.02565972 0.4361803 0.3900496 5.006155 0.0109
## 10 PCNM13 13 0.02076143 0.4569417 0.4071199 4.167133 0.0181
## 11 PCNM48 48 0.01890626 0.4758480 0.4224621 3.895580 0.0282
#
(nb.sig.PCNM <- nrow(mite.PCNM.fwd)) # Number of
## [1] 11
(PCNM.sign <- sort(mite.PCNM.fwd[,2]))
## [1] 1 2 3 4 11 13 14 15 23 48 60
PCNM.red <- mite.PCNM.pos[,c(PCNM.sign)]
#
# # head(PCNM.red)
#########################
#
env.red <- mite.env
#
PCNM.env_red <- cbind(PCNM.red, env.red, mite.xy)
PCNM.env_red <- apply(PCNM.env_red , 2, as.numeric)
PCNM.env_red <- as.data.frame(PCNM.env_red )
#
# # head(PCNM.env_red)
#
mite.PCNM.rda2 <- vegan::rda(mite.h.1 ~., data=PCNM.env_red )
axes.test <- anova.cca(mite.PCNM.rda2, by="axis")
(nb.ax <- length(which(axes.test[,4] <= 0.05)))
## [1] 2
#
#########################
#
mod <- mite.PCNM.rda2
#
plot(mod, type="n", xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5))
# points(mod, pch=21, col="red", bg="yellow", cex=0.8)
text(mod, "species", col="blue", cex=0.8)
n1 = dim(PCNM.red)[2]
n2 = dim(PCNM.env_red)[2]-n1
text(mod, dis="cn",
col=c(rep(3,n1),rep(2,n2)),
cex=c(rep(0.5,n1),rep(1,n2)),
lty=c(rep(3,n1),rep(1,n2))
)
#
##########################
#
####################################################
1.3_xy
#################################################
#
mite.PCNM.axes <- scores(mite.PCNM.rda2, choices=c(1,2),
display="lc",
scaling=1)
PCNM.diversity <- cbind(mite.xy, mite.PCNM.axes[,1:2])
# write.csv(PCNM.diversity, "PCNM_axies_diversity.csv")
# ##
s.value(mite.xy, mite.PCNM.axes[,1]) # ade4 function:
###
#
shapiro.test(resid(lm(mite.PCNM.axes[,1] ~ ., data=mite.env)))
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 1] ~ ., data = mite.env))
## W = 0.97117, p-value = 0.01113
#
mite.PCNM.axis1.env <- lm(mite.PCNM.axes[,1]~., data=mite.env)
mite.PCNM.axis1.env <- stepAIC(mite.PCNM.axis1.env)
## Start: AIC=-346.08
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全磷 1 0.00178 5.5871 -348.04
## - 有机质 1 0.00807 5.5934 -347.91
## - 沙粒 1 0.00899 5.5943 -347.89
## - 有效磷 1 0.01308 5.5984 -347.80
## - 有效氮 1 0.03057 5.6159 -347.43
## - 全钾 1 0.05327 5.6386 -346.94
## - 全氮 1 0.08267 5.6680 -346.32
## <none> 5.5854 -346.08
## - 有效钾 1 0.39059 5.9759 -339.97
## - pH 1 0.48803 6.0734 -338.03
## - 含水量 1 0.52261 6.1080 -337.35
##
## Step: AIC=-348.04
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全钾 + 有效氮 +
## 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有机质 1 0.00667 5.5938 -349.90
## - 沙粒 1 0.00785 5.5950 -349.87
## - 有效磷 1 0.01974 5.6069 -349.62
## - 有效氮 1 0.03003 5.6172 -349.40
## - 全钾 1 0.05458 5.6417 -348.88
## - 全氮 1 0.09022 5.6774 -348.12
## <none> 5.5871 -348.04
## - 有效钾 1 0.39956 5.9867 -341.75
## - pH 1 0.49387 6.0810 -339.88
## - 含水量 1 0.53999 6.1271 -338.97
##
## Step: AIC=-349.9
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 全氮 + 全钾 + 有效氮 + 有效磷 +
## 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 沙粒 1 0.00708 5.6009 -351.75
## - 有效磷 1 0.02376 5.6176 -351.39
## - 有效氮 1 0.02575 5.6196 -351.35
## - 全钾 1 0.05325 5.6471 -350.76
## <none> 5.5938 -349.90
## - 全氮 1 0.10857 5.7024 -349.59
## - 有效钾 1 0.43171 6.0255 -342.98
## - pH 1 0.49879 6.0926 -341.65
## - 含水量 1 0.53490 6.1287 -340.94
##
## Step: AIC=-351.75
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 全氮 + 全钾 + 有效氮 + 有效磷 +
## 有效钾
##
## Df Sum of Sq RSS AIC
## - 有效磷 1 0.02418 5.6251 -353.23
## - 有效氮 1 0.03281 5.6337 -353.05
## - 全钾 1 0.05666 5.6575 -352.54
## <none> 5.6009 -351.75
## - 全氮 1 0.11616 5.7170 -351.28
## - 有效钾 1 0.42469 6.0256 -344.98
## - 含水量 1 0.53019 6.1311 -342.89
## - pH 1 0.61759 6.2185 -341.20
##
## Step: AIC=-353.23
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 全氮 + 全钾 + 有效氮 + 有效钾
##
## Df Sum of Sq RSS AIC
## - 全钾 1 0.04252 5.6676 -354.33
## <none> 5.6251 -353.23
## - 全氮 1 0.11963 5.7447 -352.71
## - 有效氮 1 0.12042 5.7455 -352.69
## - 有效钾 1 0.40403 6.0291 -346.91
## - 含水量 1 0.51670 6.1418 -344.69
## - pH 1 0.59767 6.2227 -343.11
##
## Step: AIC=-354.33
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 全氮 + 有效氮 + 有效钾
##
## Df Sum of Sq RSS AIC
## <none> 5.6676 -354.33
## - 有效氮 1 0.09726 5.7648 -354.29
## - 全氮 1 0.31333 5.9809 -349.87
## - 有效钾 1 0.42892 6.0965 -347.57
## - 含水量 1 0.53973 6.2073 -345.41
## - pH 1 0.55516 6.2227 -345.11
(summ1_1 <- summary(mite.PCNM.axis1.env))
##
## Call:
## lm(formula = mite.PCNM.axes[, 1] ~ 含水量 + pH + 全氮 + 有效氮 +
## 有效钾, data = mite.env)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.39246 -0.16676 -0.02243 0.11753 0.58603
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.002369 0.020379 -0.116 0.90765
## 含水量 -0.080603 0.024463 -3.295 0.00131 **
## pH -0.070827 0.021195 -3.342 0.00113 **
## 全氮 0.056767 0.022612 2.510 0.01346 *
## 有效氮 0.038443 0.027485 1.399 0.16463
## 有效钾 -0.079824 0.027176 -2.937 0.00401 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.223 on 114 degrees of freedom
## Multiple R-squared: 0.2767, Adjusted R-squared: 0.245
## F-statistic: 8.723 on 5 and 114 DF, p-value: 4.973e-07
# write.csv(summ1_1$coefficients, "axis.Dive1_env.csv")
# #####
#
shapiro.test(resid(lm(mite.PCNM.axes[,1] ~ ., data=mite.h.1))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 1] ~ ., data = mite.h.1))
## W = 0.9919, p-value = 0.7112
# # head(mite.h.1)
mite.PCNM.axis1.div <- lm(mite.PCNM.axes[,1]~., data=mite.h.1)
mite.PCNM.axis1.div <- stepAIC(mite.PCNM.axis1.div)
## Start: AIC=-429.54
## mite.PCNM.axes[, 1] ~ OTUs + chao1 + PD + shannon + invsimp
##
## Df Sum of Sq RSS AIC
## - shannon 1 0.01856 3.0469 -430.81
## - OTUs 1 0.02344 3.0517 -430.61
## <none> 3.0283 -429.54
## - chao1 1 0.06470 3.0930 -429.00
## - PD 1 0.07909 3.1074 -428.44
## - invsimp 1 0.41431 3.4426 -416.15
##
## Step: AIC=-430.81
## mite.PCNM.axes[, 1] ~ OTUs + chao1 + PD + invsimp
##
## Df Sum of Sq RSS AIC
## - OTUs 1 0.00572 3.0526 -432.58
## - chao1 1 0.04634 3.0932 -430.99
## <none> 3.0469 -430.81
## - PD 1 0.10781 3.1547 -428.63
## - invsimp 1 0.39842 3.4453 -418.06
##
## Step: AIC=-432.58
## mite.PCNM.axes[, 1] ~ chao1 + PD + invsimp
##
## Df Sum of Sq RSS AIC
## - chao1 1 0.04088 3.0935 -432.98
## <none> 3.0526 -432.58
## - PD 1 0.56069 3.6133 -414.35
## - invsimp 1 0.61113 3.6637 -412.68
##
## Step: AIC=-432.98
## mite.PCNM.axes[, 1] ~ PD + invsimp
##
## Df Sum of Sq RSS AIC
## <none> 3.0935 -432.98
## - PD 1 0.55428 3.6477 -415.21
## - invsimp 1 0.74452 3.8380 -409.11
(summ1_2 <- summary(mite.PCNM.axis1.div))
##
## Call:
## lm(formula = mite.PCNM.axes[, 1] ~ PD + invsimp, data = mite.h.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42773 -0.09899 0.00971 0.11713 0.35608
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.271e-18 1.484e-02 0.000 1
## PD -1.321e-01 2.885e-02 -4.579 1.18e-05 ***
## invsimp -1.388e-01 2.615e-02 -5.307 5.37e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1626 on 117 degrees of freedom
## Multiple R-squared: 0.6052, Adjusted R-squared: 0.5985
## F-statistic: 89.69 on 2 and 117 DF, p-value: < 2.2e-16
# write.csv(summ1_2$coefficients, "axis.Dive1_div.csv")
# #####
# #####
s.value(mite.xy, mite.PCNM.axes[,2]) # ade4 function:
#
shapiro.test(resid(lm(mite.PCNM.axes[,2] ~ ., data=mite.env))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 2] ~ ., data = mite.env))
## W = 0.98499, p-value = 0.2051
#
mite.PCNM.axis2.env <- lm(mite.PCNM.axes[,2] ~ ., data=mite.env)
mite.PCNM.axis2.env <- step(mite.PCNM.axis2.env)
## Start: AIC=-502.17
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全钾 1 0.002745 1.5238 -503.95
## - 有效氮 1 0.002779 1.5238 -503.95
## - pH 1 0.012559 1.5336 -503.18
## - 有机质 1 0.013063 1.5341 -503.14
## - 有效磷 1 0.018099 1.5392 -502.75
## <none> 1.5211 -502.17
## - 全磷 1 0.029820 1.5509 -501.84
## - 沙粒 1 0.043621 1.5647 -500.78
## - 有效钾 1 0.132103 1.6532 -494.18
## - 全氮 1 0.189103 1.7102 -490.11
## - 含水量 1 0.264147 1.7852 -484.95
##
## Step: AIC=-503.95
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 有效氮 +
## 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有效氮 1 0.003022 1.5268 -505.72
## - pH 1 0.010375 1.5342 -505.14
## - 有机质 1 0.011691 1.5355 -505.04
## - 有效磷 1 0.017446 1.5413 -504.59
## <none> 1.5238 -503.95
## - 全磷 1 0.027429 1.5512 -503.81
## - 沙粒 1 0.043493 1.5673 -502.58
## - 有效钾 1 0.136750 1.6606 -495.64
## - 全氮 1 0.206174 1.7300 -490.73
## - 含水量 1 0.261419 1.7852 -486.95
##
## Step: AIC=-505.72
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 有效磷 +
## 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - pH 1 0.009407 1.5362 -506.98
## - 有机质 1 0.014672 1.5415 -506.57
## - 有效磷 1 0.015096 1.5419 -506.53
## <none> 1.5268 -505.72
## - 全磷 1 0.027905 1.5547 -505.54
## - 沙粒 1 0.040814 1.5676 -504.55
## - 有效钾 1 0.147964 1.6748 -496.62
## - 全氮 1 0.205058 1.7319 -492.59
## - 含水量 1 0.259522 1.7864 -488.88
##
## Step: AIC=-506.98
## mite.PCNM.axes[, 2] ~ 含水量 + 有机质 + 全氮 + 全磷 + 有效磷 +
## 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有机质 1 0.010230 1.5465 -508.18
## - 有效磷 1 0.011177 1.5474 -508.11
## <none> 1.5362 -506.98
## - 全磷 1 0.035782 1.5720 -506.22
## - 沙粒 1 0.065268 1.6015 -503.99
## - 有效钾 1 0.148945 1.6852 -497.87
## - 全氮 1 0.252382 1.7886 -490.73
## - 含水量 1 0.258043 1.7943 -490.35
##
## Step: AIC=-508.18
## mite.PCNM.axes[, 2] ~ 含水量 + 全氮 + 全磷 + 有效磷 + 有效钾 +
## 沙粒
##
## Df Sum of Sq RSS AIC
## - 有效磷 1 0.01702 1.5635 -508.87
## <none> 1.5465 -508.18
## - 全磷 1 0.02801 1.5745 -508.03
## - 沙粒 1 0.05650 1.6030 -505.88
## - 有效钾 1 0.17678 1.7232 -497.19
## - 含水量 1 0.27366 1.8201 -490.63
## - 全氮 1 0.64557 2.1920 -468.32
##
## Step: AIC=-508.87
## mite.PCNM.axes[, 2] ~ 含水量 + 全氮 + 全磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## <none> 1.5635 -508.87
## - 全磷 1 0.06579 1.6293 -505.92
## - 沙粒 1 0.06649 1.6300 -505.87
## - 有效钾 1 0.15998 1.7235 -499.18
## - 含水量 1 0.34137 1.9049 -487.17
## - 全氮 1 0.67898 2.2425 -467.59
(summ2_1 <- summary(mite.PCNM.axis2.env))
##
## Call:
## lm(formula = mite.PCNM.axes[, 2] ~ 含水量 + 全氮 + 全磷 + 有效钾 +
## 沙粒, data = mite.env)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.34953 -0.05244 -0.00449 0.07953 0.26744
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.001255 0.010701 -0.117 0.906829
## 含水量 -0.061600 0.012347 -4.989 2.19e-06 ***
## 全氮 0.085641 0.012172 7.036 1.57e-10 ***
## 全磷 -0.025457 0.011623 -2.190 0.030548 *
## 有效钾 0.045607 0.013353 3.415 0.000883 ***
## 沙粒 0.026465 0.012020 2.202 0.029694 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1171 on 114 degrees of freedom
## Multiple R-squared: 0.3987, Adjusted R-squared: 0.3723
## F-statistic: 15.12 on 5 and 114 DF, p-value: 2.237e-11
# ##
# write.csv(summ2_1$coefficients, "axis.Dive2_env.csv")
# #####
shapiro.test(resid(lm(mite.PCNM.axes[,2] ~ ., data=mite.h.1))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 2] ~ ., data = mite.h.1))
## W = 0.99291, p-value = 0.8038
mite.PCNM.axis2.div <- lm(mite.PCNM.axes[,2]~., data=mite.h.1)
mite.PCNM.axis2.div <- step(mite.PCNM.axis2.div)
## Start: AIC=-561.84
## mite.PCNM.axes[, 2] ~ OTUs + chao1 + PD + shannon + invsimp
##
## Df Sum of Sq RSS AIC
## - invsimp 1 0.00000 1.0055 -563.84
## - OTUs 1 0.00041 1.0059 -563.79
## - shannon 1 0.00643 1.0119 -563.08
## <none> 1.0055 -561.84
## - PD 1 0.02017 1.0257 -561.46
## - chao1 1 0.74626 1.7517 -497.23
##
## Step: AIC=-563.84
## mite.PCNM.axes[, 2] ~ OTUs + chao1 + PD + shannon
##
## Df Sum of Sq RSS AIC
## - OTUs 1 0.00042 1.0059 -565.79
## - shannon 1 0.00656 1.0121 -565.06
## <none> 1.0055 -563.84
## - PD 1 0.02025 1.0257 -563.45
## - chao1 1 0.80020 1.8057 -495.59
##
## Step: AIC=-565.79
## mite.PCNM.axes[, 2] ~ chao1 + PD + shannon
##
## Df Sum of Sq RSS AIC
## - shannon 1 0.01297 1.0189 -566.25
## <none> 1.0059 -565.79
## - PD 1 0.03488 1.0408 -563.70
## - chao1 1 1.40835 2.4143 -462.73
##
## Step: AIC=-566.25
## mite.PCNM.axes[, 2] ~ chao1 + PD
##
## Df Sum of Sq RSS AIC
## <none> 1.0189 -566.25
## - PD 1 0.32391 1.3428 -535.13
## - chao1 1 1.58134 2.6002 -455.83
(summ2_2 <- summary(mite.PCNM.axis2.div))
##
## Call:
## lm(formula = mite.PCNM.axes[, 2] ~ chao1 + PD, data = mite.h.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.272458 -0.061039 0.002137 0.064012 0.191726
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -7.313e-18 8.519e-03 0.000 1
## chao1 1.413e-01 1.048e-02 13.475 < 2e-16 ***
## PD -7.784e-02 1.276e-02 -6.099 1.42e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.09332 on 117 degrees of freedom
## Multiple R-squared: 0.6082, Adjusted R-squared: 0.6015
## F-statistic: 90.79 on 2 and 117 DF, p-value: < 2.2e-16
# ##
# write.csv(summ2_2$coefficients, "axis.Dive2_div.csv")
# #####
#
# windows(title="10 PCNM variables - mites")
# par(mfrow=c(2,5))
# for(i in mite.PCNM.fwd[,2][1:10] ){
# s.value(mite.xy, mite.PCNM.pos[,i], sub=i, csub=2)
# }
#############
2. L2
2.0 read data of L2.scale0
setwd('F:\\wangjingSCI\\wangjing20161104\\data')
dir()
## [1] "ade4_s.value" "ade4_s.value.r"
## [3] "Alpha20161104.csv" "axis.Dive1_div.csv"
## [5] "axis.Dive1_env.csv" "axis.Dive2_div.csv"
## [7] "axis.Dive2_env.csv" "axis.Phyla1_div.csv"
## [9] "axis.Phyla1_env.csv" "axis.Phyla2_div.csv"
## [11] "axis.Phyla2_env.csv" "chap7_web_PC20160905.R"
## [13] "L2.scale_11_dom.csv" "otu_table_sorted_L2.csv"
## [15] "PCNM_axies_diversity.csv" "PCNM_axies_phyla.csv"
## [17] "Soil20161104.csv"
require(vegan)
L2.scale0 <- read.csv("otu_table_sorted_L2.csv")
L2.scale0 <- droplevels(L2.scale0[L2.scale0$样带!='灌乔带',])
#
L2.scale <- as.data.frame(scale(L2.scale0[,9:19]))
soil.scale0 <- read.csv('Soil20161104.csv')
soil.scale0 <- droplevels(soil.scale0[soil.scale0$样带!='灌乔带',])
#
################################################
#
r.1 <- soil.scale0$replication
L2.1 <- (L2.scale[r.1==1,] + L2.scale[r.1==2,])/2
#
xy0 <- soil.scale0[, c("Dist_upper","Dist_water")]
names(xy0) <- c("x","y")
xy1 <- (xy0[r.1==1,] + xy0[r.1==2,])/2
add.y <- rep(c(-0.105, 0, 0.105),length=dim(xy1)[1])
xy1$y <- xy1$y+add.y
# plot(x=xy1$x, y=xy1$y)
###
#
mite.h <- L2.1
#
index_1 <- mite.xy$x==1 & mite.xy$y<1.5
index_2 <- mite.xy$x==4 & mite.xy$y>2.5 & mite.xy$y<3.5
index_3 <- mite.xy$x==7 & mite.xy$y>2.5 & mite.xy$y<3.5
#
mite.h_add <- mite.h[index_1 | index_2 | index_3,]
#
mite.h <- rbind(mite.h, mite.h_add)
#
mite.h <- apply(mite.h, 2, as.numeric)
#
############################################
2.1 rda(mite.h, mite_Env.Pcnm)
##################################################
#
mite.h.1 <- as.data.frame(mite.h)
#
dim(mite.h)
## [1] 120 11
#####################################################
# #
# # head(mite.env)
# dim(mite.env)
# mite.env.rda <- rda(mite.h.1, mite.env)
# #
# (mite.R2a <- RsquareAdj(mite.env.rda)$adj.r.squared)
# #
# (mite.Env.fwd <- forward.sel(mite.h.1, mite.env,
# nperm = 9999,
# adjR2thresh=mite.R2a,
# alpha = 0.1))
# #
# (nb.sig.Env <- nrow(mite.Env.fwd)) # Number of
# (Env.sign <- sort(mite.Env.fwd[,2]))
# env.red <- mite.env[,c(Env.sign)]
# #
# head(env.red)
#
###################################################
#
mite.PCNM.rda <- rda(mite.h.1, mite.PCNM.pos)
#
(mite.R2a <- RsquareAdj(mite.PCNM.rda)$adj.r.squared)
## [1] 0.3928379
#
(mite.PCNM.fwd <- forward.sel(mite.h.1, mite.PCNM.pos,
nperm = 9999,
adjR2thresh=mite.R2a,
alpha = 0.1))
## Testing variable 1
## Testing variable 2
## Testing variable 3
## Testing variable 4
## Testing variable 5
## Testing variable 6
## Testing variable 7
## Testing variable 8
## Testing variable 9
## Testing variable 10
## Testing variable 11
## Testing variable 12
## Testing variable 13
## Testing variable 14
## Testing variable 15
## Testing variable 16
## Testing variable 17
## Testing variable 18
## Testing variable 19
## Testing variable 20
## Testing variable 21
## Procedure stopped (adjR2thresh criteria) adjR2cum = 0.394281 with 21 variables (superior to 0.392838)
## variables order R2 R2Cum AdjR2Cum F pval
## 1 PCNM5 5 0.05378747 0.05378747 0.04576872 6.707713 0.0004
## 2 PCNM1 1 0.05224499 0.10603246 0.09075096 6.837679 0.0001
## 3 PCNM13 13 0.05141190 0.15744436 0.13565412 7.078203 0.0001
## 4 PCNM4 4 0.04209364 0.19953799 0.17169584 6.047468 0.0001
## 5 PCNM8 8 0.04065746 0.24019545 0.20687069 6.100188 0.0001
## 6 PCNM11 11 0.02966720 0.26986266 0.23109430 4.591457 0.0001
## 7 PCNM6 6 0.02284452 0.29270718 0.24850138 3.617436 0.0018
## 8 PCNM18 18 0.02072078 0.31342796 0.26394529 3.349986 0.0023
## 9 PCNM3 3 0.02069915 0.33412711 0.27964660 3.419432 0.0018
## 10 PCNM7 7 0.01965425 0.35378136 0.29449525 3.315153 0.0014
## 11 PCNM10 10 0.01776391 0.37154527 0.30753599 3.052729 0.0037
## 12 PCNM2 2 0.01538583 0.38693110 0.31817571 2.685317 0.0088
## 13 PCNM22 22 0.01503658 0.40196768 0.32862409 2.665202 0.0171
## 14 PCNM48 48 0.01474697 0.41671464 0.33894326 2.654672 0.0106
## 15 PCNM15 15 0.01343705 0.43015169 0.34796203 2.452325 0.0170
## 16 PCNM23 23 0.01266160 0.44281329 0.35626001 2.340588 0.0204
## 17 PCNM14 14 0.01236437 0.45517766 0.36437394 2.314820 0.0219
## 18 PCNM60 60 0.01229673 0.46747439 0.37256883 2.332225 0.0249
## 19 PCNM9 9 0.01167947 0.47915385 0.38019308 2.242402 0.0283
## 20 PCNM30 30 0.01110057 0.49025442 0.38727552 2.155893 0.0468
## 21 PCNM44 44 0.01091800 0.50117242 0.39428080 2.144958 0.0446
#
(nb.sig.PCNM <- nrow(mite.PCNM.fwd)) # Number of
## [1] 21
(PCNM.sign <- sort(mite.PCNM.fwd[,2]))
## [1] 1 2 3 4 5 6 7 8 9 10 11 13 14 15 18 22 23 30 44 48 60
PCNM.red <- mite.PCNM.pos[,c(PCNM.sign)]
# head(PCNM.red)
#########################
#
env.red <- mite.env
#
PCNM.env_red <- cbind(PCNM.red, env.red, mite.xy)
PCNM.env_red <- apply(PCNM.env_red , 2, as.numeric)
PCNM.env_red <- as.data.frame(PCNM.env_red )
#
# head(PCNM.env_red)
#
mite.PCNM.rda2 <- vegan::rda(mite.h.1 ~., data=PCNM.env_red )
axes.test <- anova.cca(mite.PCNM.rda2, by="axis")
(nb.ax <- length(which(axes.test[,4] <= 0.05)))
## [1] 8
#
#########################
#
#########################
#
mod <- mite.PCNM.rda2
#
plot(mod, scaling = "symmetric",type="n", xlim=c(-3,3), ylim=c(-3,3))
# points(mod, pch=21, col="red", bg="yellow", cex=0.8)
text(mod, "species", col="blue", cex=0.8)
n1 = dim(PCNM.red)[2]
n2 = dim(PCNM.env_red)[2]-n1
text(mod, dis="cn",
col=c(rep(3,n1),rep(2,n2)),
cex=c(rep(0.8,n1),rep(1,n2)),
lty=c(rep(3,n1),rep(1,n2))
)
#
##########################
#
####################################################
####################################################
#
2.1_xy
#################################################
#
mite.PCNM.axes <- scores(mite.PCNM.rda2, choices=c(1,2,3,4,5), display="lc",
scaling=1)
#
#####
PCNM_phyla <- cbind(mite.xy, mite.PCNM.axes[,1:2])
# write.csv(PCNM_phyla, "PCNM_axies_phyla.csv")
#####
#
s.value(mite.xy, mite.PCNM.axes[,1]) # ade4 function: s.value
#
#
shapiro.test(resid(lm(mite.PCNM.axes[,1] ~ ., data=mite.env))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 1] ~ ., data = mite.env))
## W = 0.9877, p-value = 0.3531
#
mite.PCNM.axis1.env <- lm(mite.PCNM.axes[,1]~., data=mite.env)
mite.PCNM.axis1.env <- step(mite.PCNM.axis1.env)
## Start: AIC=-369.32
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全钾 1 0.00017 4.6021 -371.32
## - 有机质 1 0.00108 4.6030 -371.29
## - 全磷 1 0.00417 4.6061 -371.21
## - 全氮 1 0.04032 4.6423 -370.27
## - 有效磷 1 0.07084 4.6728 -369.49
## <none> 4.6019 -369.32
## - 有效钾 1 0.24737 4.8493 -365.04
## - 有效氮 1 0.26844 4.8704 -364.52
## - 沙粒 1 0.43276 5.0347 -360.54
## - 含水量 1 0.64091 5.2429 -355.67
## - pH 1 0.93710 5.5390 -349.08
##
## Step: AIC=-371.32
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 有效氮 +
## 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有机质 1 0.00098 4.6031 -373.29
## - 全磷 1 0.00582 4.6079 -373.17
## - 全氮 1 0.05138 4.6535 -371.98
## - 有效磷 1 0.07136 4.6735 -371.47
## <none> 4.6021 -371.32
## - 有效钾 1 0.25065 4.8528 -366.95
## - 有效氮 1 0.26834 4.8705 -366.52
## - 沙粒 1 0.43268 5.0348 -362.53
## - 含水量 1 0.64863 5.2508 -357.49
## - pH 1 0.99996 5.6021 -349.72
##
## Step: AIC=-373.29
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 全氮 + 全磷 + 有效氮 + 有效磷 +
## 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全磷 1 0.00733 4.6104 -375.10
## - 有效磷 1 0.07149 4.6746 -373.44
## <none> 4.6031 -373.29
## - 全氮 1 0.08447 4.6876 -373.11
## - 有效钾 1 0.26570 4.8688 -368.56
## - 有效氮 1 0.27238 4.8755 -368.39
## - 沙粒 1 0.43296 5.0361 -364.50
## - 含水量 1 0.65418 5.2573 -359.35
## - pH 1 1.03274 5.6358 -351.00
##
## Step: AIC=-375.1
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 全氮 + 有效氮 + 有效磷 +
## 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全氮 1 0.07723 4.6877 -375.11
## <none> 4.6104 -375.10
## - 有效磷 1 0.10841 4.7188 -374.31
## - 有效钾 1 0.26361 4.8740 -370.43
## - 有效氮 1 0.27072 4.8811 -370.25
## - 沙粒 1 0.46869 5.0791 -365.48
## - 含水量 1 0.64757 5.2580 -361.33
## - pH 1 1.05683 5.6673 -352.33
##
## Step: AIC=-375.11
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有效氮 + 有效磷 + 有效钾 +
## 沙粒
##
## Df Sum of Sq RSS AIC
## <none> 4.6877 -375.11
## - 有效磷 1 0.13499 4.8226 -373.70
## - 有效钾 1 0.23201 4.9197 -371.31
## - 有效氮 1 0.33420 5.0219 -368.84
## - 沙粒 1 0.55853 5.2462 -363.60
## - 含水量 1 0.61685 5.3045 -362.27
## - pH 1 0.99793 5.6856 -353.95
(summ3_1 <- summary(mite.PCNM.axis1.env))
##
## Call:
## lm(formula = mite.PCNM.axes[, 1] ~ 含水量 + pH + 有效氮 + 有效磷 +
## 有效钾 + 沙粒, data = mite.env)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.46188 -0.15505 0.02212 0.15489 0.58708
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0009306 0.0186186 -0.050 0.960225
## 含水量 -0.0863161 0.0223842 -3.856 0.000192 ***
## pH -0.1026705 0.0209332 -4.905 3.16e-06 ***
## 有效氮 0.0920451 0.0324294 2.838 0.005378 **
## 有效磷 -0.0534611 0.0296369 -1.804 0.073916 .
## 有效钾 -0.0594017 0.0251179 -2.365 0.019739 *
## 沙粒 -0.0840697 0.0229115 -3.669 0.000373 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2037 on 113 degrees of freedom
## Multiple R-squared: 0.4748, Adjusted R-squared: 0.4469
## F-statistic: 17.03 on 6 and 113 DF, p-value: 6.223e-14
# ##
# write.csv(summ3_1$coefficients, "axis.Phyla1_env.csv")
# ##
##########
#
#
shapiro.test(resid(lm(mite.PCNM.axes[,1] ~ ., data=mite.h.1))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 1] ~ ., data = mite.h.1))
## W = 0.98885, p-value = 0.437
#
mite.PCNM.axis1.div <- lm(mite.PCNM.axes[,1]~., data=mite.h.1)
mite.PCNM.axis1.div <- step(mite.PCNM.axis1.div)
## Start: AIC=-503.05
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes + p__Cyanobacteria
##
## Df Sum of Sq RSS AIC
## - p__Cyanobacteria 1 0.002322 1.4873 -504.86
## - p__Chloroflexi 1 0.009430 1.4944 -504.29
## - p__Bacteroidetes 1 0.019285 1.5043 -503.50
## <none> 1.4850 -503.05
## - p__Nitrospirae 1 0.027705 1.5127 -502.83
## - p__Planctomycetes 1 0.063690 1.5487 -500.01
## - p__Actinobacteria 1 0.074884 1.5599 -499.15
## - p__Proteobacteria 1 0.086470 1.5714 -498.26
## - p__Gemmatimonadetes 1 0.099362 1.5843 -497.28
## - p__Verrucomicrobia 1 0.108595 1.5936 -496.58
## - p__Crenarchaeota 1 0.133218 1.6182 -494.74
## - p__Acidobacteria 1 0.152460 1.6374 -493.32
##
## Step: AIC=-504.86
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Chloroflexi 1 0.007220 1.4945 -506.28
## - p__Bacteroidetes 1 0.017510 1.5048 -505.46
## <none> 1.4873 -504.86
## - p__Nitrospirae 1 0.031025 1.5183 -504.39
## - p__Planctomycetes 1 0.072648 1.5599 -501.14
## - p__Actinobacteria 1 0.084766 1.5721 -500.21
## - p__Gemmatimonadetes 1 0.098699 1.5860 -499.15
## - p__Proteobacteria 1 0.107450 1.5947 -498.49
## - p__Verrucomicrobia 1 0.156834 1.6441 -494.83
## - p__Crenarchaeota 1 0.176977 1.6643 -493.37
## - p__Acidobacteria 1 0.185888 1.6732 -492.73
##
## Step: AIC=-506.28
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Bacteroidetes 1 0.01161 1.5061 -507.35
## <none> 1.4945 -506.28
## - p__Nitrospirae 1 0.03013 1.5247 -505.89
## - p__Planctomycetes 1 0.07101 1.5655 -502.71
## - p__Actinobacteria 1 0.08131 1.5758 -501.93
## - p__Gemmatimonadetes 1 0.09218 1.5867 -501.10
## - p__Verrucomicrobia 1 0.27054 1.7651 -488.32
## - p__Proteobacteria 1 0.27419 1.7687 -488.07
## - p__Acidobacteria 1 0.41790 1.9124 -478.70
## - p__Crenarchaeota 1 0.54405 2.0386 -471.03
##
## Step: AIC=-507.35
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Nitrospirae 1 0.01936 1.5255 -507.82
## <none> 1.5061 -507.35
## - p__Planctomycetes 1 0.06220 1.5683 -504.50
## - p__Actinobacteria 1 0.07101 1.5771 -503.83
## - p__Gemmatimonadetes 1 0.08261 1.5887 -502.95
## - p__Verrucomicrobia 1 0.25929 1.7654 -490.29
## - p__Proteobacteria 1 0.50925 2.0154 -474.40
## - p__Acidobacteria 1 0.66957 2.1757 -465.22
## - p__Crenarchaeota 1 1.41943 2.9255 -429.68
##
## Step: AIC=-507.82
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Verrucomicrobia + p__Planctomycetes +
## p__Actinobacteria + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## <none> 1.5255 -507.82
## - p__Actinobacteria 1 0.05219 1.5777 -505.79
## - p__Planctomycetes 1 0.07884 1.6043 -503.78
## - p__Gemmatimonadetes 1 0.08918 1.6147 -503.00
## - p__Verrucomicrobia 1 0.27537 1.8008 -489.91
## - p__Proteobacteria 1 0.74602 2.2715 -462.05
## - p__Acidobacteria 1 0.94061 2.4661 -452.18
## - p__Crenarchaeota 1 2.92327 4.4487 -381.38
(summ3_2 <- summary(mite.PCNM.axis1.div))
##
## Call:
## lm(formula = mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Verrucomicrobia + p__Planctomycetes +
## p__Actinobacteria + p__Gemmatimonadetes, data = mite.h.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.38634 -0.07942 0.00709 0.07780 0.24393
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.008832 0.010674 0.827 0.4097
## p__Proteobacteria 0.242652 0.032787 7.401 2.68e-11 ***
## p__Crenarchaeota 0.344404 0.023509 14.650 < 2e-16 ***
## p__Acidobacteria 0.165220 0.019882 8.310 2.50e-13 ***
## p__Verrucomicrobia 0.055512 0.012346 4.496 1.70e-05 ***
## p__Planctomycetes 0.042469 0.017652 2.406 0.0178 *
## p__Actinobacteria 0.026314 0.013442 1.958 0.0528 .
## p__Gemmatimonadetes 0.036372 0.014214 2.559 0.0118 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1167 on 112 degrees of freedom
## Multiple R-squared: 0.8291, Adjusted R-squared: 0.8184
## F-statistic: 77.61 on 7 and 112 DF, p-value: < 2.2e-16
# ##
# write.csv(summ3_2$coefficients, "axis.Phyla1_div.csv")
# ##
#
###################################
#
s.value(mite.xy, mite.PCNM.axes[,2]) # ade4 function: s.value
# ##
#
shapiro.test(resid(lm(mite.PCNM.axes[,2] ~ ., data=mite.env))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 2] ~ ., data = mite.env))
## W = 0.98775, p-value = 0.356
#
mite.PCNM.axis2.env <- lm(mite.PCNM.axes[,2] ~ ., data=mite.env)
mite.PCNM.axis2.env <- step(mite.PCNM.axis2.env)
## Start: AIC=-476.06
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有机质 1 0.002157 1.8929 -477.93
## - 有效氮 1 0.002474 1.8932 -477.91
## - 有效磷 1 0.004364 1.8951 -477.79
## - 全磷 1 0.005248 1.8960 -477.73
## - 全钾 1 0.005427 1.8962 -477.72
## - 全氮 1 0.023931 1.9146 -476.55
## <none> 1.8907 -476.06
## - pH 1 0.069249 1.9600 -473.75
## - 有效钾 1 0.114543 2.0053 -471.01
## - 含水量 1 0.185284 2.0760 -466.85
## - 沙粒 1 0.220592 2.1113 -464.82
##
## Step: AIC=-477.93
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮 + 全磷 + 全钾 + 有效氮 +
## 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有效氮 1 0.001710 1.8946 -479.82
## - 全钾 1 0.004615 1.8975 -479.64
## - 有效磷 1 0.006208 1.8991 -479.53
## - 全磷 1 0.007661 1.9005 -479.44
## - 全氮 1 0.024691 1.9176 -478.37
## <none> 1.8929 -477.93
## - pH 1 0.067367 1.9602 -475.73
## - 有效钾 1 0.112798 2.0057 -472.98
## - 含水量 1 0.193187 2.0861 -468.27
## - 沙粒 1 0.218437 2.1113 -466.82
##
## Step: AIC=-479.82
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮 + 全磷 + 全钾 + 有效磷 +
## 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全钾 1 0.004530 1.8991 -481.53
## - 全磷 1 0.007573 1.9022 -481.34
## - 有效磷 1 0.014899 1.9095 -480.88
## - 全氮 1 0.027050 1.9216 -480.12
## <none> 1.8946 -479.82
## - pH 1 0.071641 1.9662 -477.36
## - 有效钾 1 0.124328 2.0189 -474.19
## - 含水量 1 0.202939 2.0975 -469.61
## - 沙粒 1 0.235774 2.1304 -467.74
##
## Step: AIC=-481.53
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮 + 全磷 + 有效磷 + 有效钾 +
## 沙粒
##
## Df Sum of Sq RSS AIC
## - 有效磷 1 0.013289 1.9124 -482.70
## - 全磷 1 0.014955 1.9141 -482.59
## <none> 1.8991 -481.53
## - pH 1 0.067115 1.9662 -479.36
## - 全氮 1 0.068981 1.9681 -479.25
## - 有效钾 1 0.121585 2.0207 -476.09
## - 含水量 1 0.198632 2.0978 -471.60
## - 沙粒 1 0.236116 2.1352 -469.47
##
## Step: AIC=-482.7
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮 + 全磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全磷 1 0.005628 1.9180 -484.34
## <none> 1.9124 -482.70
## - pH 1 0.056669 1.9691 -481.19
## - 全氮 1 0.061444 1.9739 -480.90
## - 有效钾 1 0.170938 2.0833 -474.42
## - 含水量 1 0.243853 2.1563 -470.29
## - 沙粒 1 0.275220 2.1876 -468.56
##
## Step: AIC=-484.34
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## <none> 1.9180 -484.34
## - pH 1 0.052396 1.9704 -483.11
## - 全氮 1 0.056199 1.9742 -482.88
## - 有效钾 1 0.165310 2.0833 -476.42
## - 含水量 1 0.244408 2.1624 -471.95
## - 沙粒 1 0.283669 2.2017 -469.79
(summ4_1 <- summary(mite.PCNM.axis2.env))
##
## Call:
## lm(formula = mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮 + 有效钾 +
## 沙粒, data = mite.env)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.301121 -0.076500 0.004307 0.068307 0.294978
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.00437 0.01185 -0.369 0.713062
## 含水量 -0.05247 0.01377 -3.811 0.000225 ***
## pH 0.02325 0.01317 1.765 0.080291 .
## 全氮 -0.02426 0.01328 -1.828 0.070219 .
## 有效钾 -0.04564 0.01456 -3.135 0.002189 **
## 沙粒 0.05896 0.01436 4.106 7.6e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1297 on 114 degrees of freedom
## Multiple R-squared: 0.5349, Adjusted R-squared: 0.5145
## F-statistic: 26.22 on 5 and 114 DF, p-value: < 2.2e-16
# ##
# write.csv(summ4_1$coefficients, "axis.Phyla2_env.csv")
# ##
#
#
shapiro.test(resid(lm(mite.PCNM.axes[,2] ~ ., data=mite.h.1))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 2] ~ ., data = mite.h.1))
## W = 0.94641, p-value = 0.0001189
#
mite.PCNM.axis2.div <- lm(mite.PCNM.axes[,2]~., data=mite.h.1)
mite.PCNM.axis2.div <- step(mite.PCNM.axis2.div)
## Start: AIC=-504
## mite.PCNM.axes[, 2] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes + p__Cyanobacteria
##
## Df Sum of Sq RSS AIC
## - p__Bacteroidetes 1 0.000002 1.4732 -506.00
## - p__Proteobacteria 1 0.000007 1.4732 -506.00
## - p__Crenarchaeota 1 0.000511 1.4737 -505.96
## - p__Planctomycetes 1 0.002981 1.4762 -505.76
## - p__Chloroflexi 1 0.005158 1.4784 -505.59
## - p__Cyanobacteria 1 0.007991 1.4812 -505.36
## - p__Acidobacteria 1 0.014038 1.4873 -504.87
## - p__Verrucomicrobia 1 0.021664 1.4949 -504.25
## <none> 1.4732 -504.00
## - p__Nitrospirae 1 0.026015 1.4992 -503.90
## - p__Actinobacteria 1 0.051379 1.5246 -501.89
## - p__Gemmatimonadetes 1 0.072951 1.5462 -500.20
##
## Step: AIC=-506
## mite.PCNM.axes[, 2] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes + p__Cyanobacteria
##
## Df Sum of Sq RSS AIC
## - p__Proteobacteria 1 0.000147 1.4734 -507.99
## - p__Crenarchaeota 1 0.006214 1.4794 -507.50
## - p__Planctomycetes 1 0.007814 1.4810 -507.37
## - p__Cyanobacteria 1 0.010410 1.4836 -507.16
## - p__Chloroflexi 1 0.016812 1.4900 -506.64
## <none> 1.4732 -506.00
## - p__Verrucomicrobia 1 0.058851 1.5321 -503.30
## - p__Gemmatimonadetes 1 0.083863 1.5571 -501.36
## - p__Actinobacteria 1 0.086551 1.5598 -501.15
## - p__Nitrospirae 1 0.092675 1.5659 -500.68
## - p__Acidobacteria 1 0.101317 1.5745 -500.02
##
## Step: AIC=-507.99
## mite.PCNM.axes[, 2] ~ p__Crenarchaeota + p__Acidobacteria + p__Nitrospirae +
## p__Verrucomicrobia + p__Planctomycetes + p__Actinobacteria +
## p__Chloroflexi + p__Gemmatimonadetes + p__Cyanobacteria
##
## Df Sum of Sq RSS AIC
## - p__Planctomycetes 1 0.008476 1.4818 -509.30
## - p__Cyanobacteria 1 0.013887 1.4873 -508.87
## <none> 1.4734 -507.99
## - p__Chloroflexi 1 0.025006 1.4984 -507.97
## - p__Crenarchaeota 1 0.058041 1.5314 -505.36
## - p__Gemmatimonadetes 1 0.085047 1.5584 -503.26
## - p__Verrucomicrobia 1 0.112121 1.5855 -501.19
## - p__Actinobacteria 1 0.121407 1.5948 -500.49
## - p__Nitrospirae 1 0.262244 1.7356 -490.34
## - p__Acidobacteria 1 0.288249 1.7616 -488.55
##
## Step: AIC=-509.3
## mite.PCNM.axes[, 2] ~ p__Crenarchaeota + p__Acidobacteria + p__Nitrospirae +
## p__Verrucomicrobia + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes + p__Cyanobacteria
##
## Df Sum of Sq RSS AIC
## - p__Cyanobacteria 1 0.00976 1.4916 -510.52
## <none> 1.4818 -509.30
## - p__Chloroflexi 1 0.03322 1.5151 -508.64
## - p__Crenarchaeota 1 0.05209 1.5339 -507.16
## - p__Gemmatimonadetes 1 0.07990 1.5617 -505.00
## - p__Verrucomicrobia 1 0.10397 1.5858 -503.17
## - p__Actinobacteria 1 0.11793 1.5998 -502.12
## - p__Nitrospirae 1 0.31880 1.8007 -487.92
## - p__Acidobacteria 1 0.59557 2.0774 -470.76
##
## Step: AIC=-510.52
## mite.PCNM.axes[, 2] ~ p__Crenarchaeota + p__Acidobacteria + p__Nitrospirae +
## p__Verrucomicrobia + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## <none> 1.4916 -510.52
## - p__Chloroflexi 1 0.04813 1.5397 -508.70
## - p__Gemmatimonadetes 1 0.07651 1.5681 -506.51
## - p__Crenarchaeota 1 0.09312 1.5847 -505.25
## - p__Actinobacteria 1 0.11196 1.6036 -503.83
## - p__Verrucomicrobia 1 0.16573 1.6573 -499.87
## - p__Nitrospirae 1 0.32124 1.8129 -489.11
## - p__Acidobacteria 1 0.60174 2.0934 -471.85
(summ4_2 <- summary(mite.PCNM.axis2.div))
##
## Call:
## lm(formula = mite.PCNM.axes[, 2] ~ p__Crenarchaeota + p__Acidobacteria +
## p__Nitrospirae + p__Verrucomicrobia + p__Actinobacteria +
## p__Chloroflexi + p__Gemmatimonadetes, data = mite.h.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.42993 -0.06694 0.01804 0.07062 0.27938
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.004146 0.010572 -0.392 0.695662
## p__Crenarchaeota -0.035413 0.013392 -2.644 0.009362 **
## p__Acidobacteria 0.093301 0.013880 6.722 7.86e-10 ***
## p__Nitrospirae 0.054532 0.011103 4.911 3.11e-06 ***
## p__Verrucomicrobia -0.041215 0.011684 -3.528 0.000609 ***
## p__Actinobacteria 0.038916 0.013422 2.899 0.004499 **
## p__Chloroflexi 0.024785 0.013037 1.901 0.059857 .
## p__Gemmatimonadetes 0.033643 0.014036 2.397 0.018190 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1154 on 112 degrees of freedom
## Multiple R-squared: 0.6383, Adjusted R-squared: 0.6157
## F-statistic: 28.23 on 7 and 112 DF, p-value: < 2.2e-16
# ##
# write.csv(summ4_2$coefficients, "axis.Phyla2_div.csv")
# ##
#
# windows(title="10 PCNM variables - mites")
# par(mfrow=c(2,5))
# for(i in mite.PCNM.fwd[,2][1:10] ){
# s.value(mite.xy, mite.PCNM.pos[,i], sub=i, csub=2)
# }
2.2 rda(mite.h, mite_Env.Pcnm, mite_x.Pcnm)
###################################################
#
#
mite.h.1 <- residuals(lm(mite.h~PCNM_x.pos))
mite.h.1 <- as.data.frame(mite.h.1)
#
#####################################################
# #
# mite.env.rda <- rda(mite.h.1, mite.env)
# #
# (mite.R2a <- RsquareAdj(mite.env.rda)$adj.r.squared)
# #
# (mite.Env.fwd <- forward.sel(mite.h.1, mite.env,
# nperm = 9999,
# adjR2thresh=mite.R2a,
# alpha=0.1))
# #
# (nb.sig.Env <- nrow(mite.Env.fwd)) # Number of
# (Env.sign <- sort(mite.Env.fwd[,2]))
# env.red <- mite.env[,c(Env.sign)]
# #
# # head(env.red)
#
###################################################
#
mite.PCNM.rda <- rda(mite.h.1, mite.PCNM.pos)
#
(mite.R2a <- RsquareAdj(mite.PCNM.rda)$adj.r.squared)
## [1] 0.2893349
#
(mite.PCNM.fwd <- forward.sel(mite.h.1, mite.PCNM.pos,
nperm = 9999,
adjR2thresh=mite.R2a,
alpha=0.1))
## Testing variable 1
## Testing variable 2
## Testing variable 3
## Testing variable 4
## Testing variable 5
## Testing variable 6
## Testing variable 7
## Testing variable 8
## Testing variable 9
## Testing variable 10
## Testing variable 11
## Testing variable 12
## Testing variable 13
## Procedure stopped (adjR2thresh criteria) adjR2cum = 0.297046 with 13 variables (superior to 0.289335)
## variables order R2 R2Cum AdjR2Cum F pval
## 1 PCNM5 5 0.05985040 0.0598504 0.05188303 7.511940 0.0001
## 2 PCNM13 13 0.05693914 0.1167895 0.10169193 7.542799 0.0001
## 3 PCNM8 8 0.04948883 0.1662784 0.14471661 6.885637 0.0001
## 4 PCNM11 11 0.03508547 0.2013638 0.17358520 5.052150 0.0002
## 5 PCNM18 18 0.02582046 0.2271843 0.19328888 3.808841 0.0013
## 6 PCNM6 6 0.02496983 0.2521541 0.21244551 3.772958 0.0004
## 7 PCNM10 10 0.02025014 0.2724043 0.22692955 3.117137 0.0036
## 8 PCNM24 24 0.02024535 0.2926496 0.24166942 3.176974 0.0037
## 9 PCNM22 22 0.01778020 0.3104298 0.25401045 2.836291 0.0131
## 10 PCNM48 48 0.01744028 0.3278701 0.26620681 2.828308 0.0094
## 11 PCNM15 15 0.01589113 0.3437612 0.27692210 2.615270 0.0141
## 12 PCNM17 17 0.01517288 0.3589341 0.28703888 2.532499 0.0177
## 13 PCNM23 23 0.01490502 0.3738391 0.29704582 2.523205 0.0168
#
(nb.sig.PCNM <- nrow(mite.PCNM.fwd)) # Number of
## [1] 13
(PCNM.sign <- sort(mite.PCNM.fwd[,2]))
## [1] 5 6 8 10 11 13 15 17 18 22 23 24 48
PCNM.red <- mite.PCNM.pos[,c(PCNM.sign)]
#
# # head(PCNM.red)
#########################
#
env.red <- mite.env
#
PCNM.env_red <- cbind(PCNM.red, env.red, mite.xy)
PCNM.env_red <- apply(PCNM.env_red , 2, as.numeric)
PCNM.env_red <- as.data.frame(PCNM.env_red )
#
# # head(PCNM.env_red)
#
mite.PCNM.rda2 <- vegan::rda(mite.h.1 ~., data=PCNM.env_red )
axes.test <- anova.cca(mite.PCNM.rda2, by="axis")
(nb.ax <- length(which(axes.test[,4] <= 0.05)))
## [1] 6
#
#########################
#
mod <- mite.PCNM.rda2
#
plot(mod, type="n", xlim=c(-2.5, 2.5), ylim=c(-2.5, 2.5))
# points(mod, pch=21, col="red", bg="yellow", cex=0.8)
text(mod, "species", col="blue", cex=0.8)
n1 = dim(PCNM.red)[2]
n2 = dim(PCNM.env_red)[2]-n1
text(mod, dis="cn",
col=c(rep(3,n1),rep(2,n2)),
cex=c(rep(0.5,n1),rep(1,n2)),
lty=c(rep(3,n1),rep(1,n2))
)
#
##########################
#
####################################################
2.2_xy
##################################################
#
mite.PCNM.axes <- scores(mite.PCNM.rda2, choices=c(1,2,3,4,5), display="lc",
scaling=1)
#
#####
PCNM_phyla <- cbind(mite.xy, mite.PCNM.axes[,1:2])
# write.csv(PCNM_phyla, "PCNM_axies_phyla.csv")
#####
#
s.value(mite.xy, mite.PCNM.axes[,1]) # ade4 function: s.value
#
#
shapiro.test(resid(lm(mite.PCNM.axes[,1] ~ ., data=mite.env))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 1] ~ ., data = mite.env))
## W = 0.99015, p-value = 0.5476
#
mite.PCNM.axis1.env <- lm(mite.PCNM.axes[,1]~., data=mite.env)
mite.PCNM.axis1.env <- step(mite.PCNM.axis1.env)
## Start: AIC=-447.11
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全氮 1 0.00190 2.4085 -449.02
## - 全磷 1 0.00323 2.4099 -448.95
## - 全钾 1 0.01903 2.4257 -448.17
## - 有机质 1 0.02333 2.4300 -447.95
## - 有效钾 1 0.03766 2.4443 -447.25
## <none> 2.4066 -447.11
## - 有效磷 1 0.06887 2.4755 -445.73
## - 有效氮 1 0.18162 2.5882 -440.38
## - 含水量 1 0.53029 2.9369 -425.22
## - 沙粒 1 0.68856 3.0952 -418.92
## - pH 1 1.27164 3.6783 -398.21
##
## Step: AIC=-449.02
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全磷 + 全钾 + 有效氮 +
## 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全磷 1 0.00434 2.4129 -450.80
## - 有机质 1 0.02458 2.4331 -449.80
## - 全钾 1 0.02942 2.4379 -449.56
## - 有效钾 1 0.03766 2.4462 -449.16
## <none> 2.4085 -449.02
## - 有效磷 1 0.06807 2.4766 -447.67
## - 有效氮 1 0.18143 2.5900 -442.30
## - 含水量 1 0.53895 2.9475 -426.79
## - 沙粒 1 0.69819 3.1067 -420.47
## - pH 1 1.28335 3.6919 -399.76
##
## Step: AIC=-450.8
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全钾 + 有效氮 +
## 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有机质 1 0.02029 2.4332 -451.80
## - 全钾 1 0.02509 2.4380 -451.56
## <none> 2.4129 -450.80
## - 有效钾 1 0.04064 2.4535 -450.80
## - 有效磷 1 0.06524 2.4781 -449.60
## - 有效氮 1 0.17942 2.5923 -444.19
## - 含水量 1 0.55393 2.9668 -428.00
## - 沙粒 1 0.70015 3.1130 -422.23
## - pH 1 1.29127 3.7041 -401.36
##
## Step: AIC=-451.8
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 全钾 + 有效氮 + 有效磷 +
## 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全钾 1 0.01184 2.4450 -453.21
## <none> 2.4332 -451.80
## - 有效钾 1 0.05717 2.4903 -451.01
## - 有效磷 1 0.05928 2.4924 -450.91
## - 有效氮 1 0.16057 2.5937 -446.13
## - 含水量 1 0.54291 2.9761 -429.63
## - 沙粒 1 0.68200 3.1152 -424.15
## - pH 1 1.28846 3.7216 -402.80
##
## Step: AIC=-453.21
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有效氮 + 有效磷 + 有效钾 +
## 沙粒
##
## Df Sum of Sq RSS AIC
## <none> 2.4450 -453.21
## - 有效钾 1 0.05584 2.5008 -452.50
## - 有效磷 1 0.08064 2.5256 -451.32
## - 有效氮 1 0.16843 2.6134 -447.22
## - 含水量 1 0.54318 2.9882 -431.14
## - 沙粒 1 0.73148 3.1765 -423.81
## - pH 1 1.36396 3.8090 -402.02
(summ3_1 <- summary(mite.PCNM.axis1.env))
##
## Call:
## lm(formula = mite.PCNM.axes[, 1] ~ 含水量 + pH + 有效氮 + 有效磷 +
## 有效钾 + 沙粒, data = mite.env)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.39673 -0.07713 0.00013 0.08962 0.35103
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.403e-05 1.345e-02 0.007 0.99443
## 含水量 8.100e-02 1.617e-02 5.010 2.02e-06 ***
## pH 1.200e-01 1.512e-02 7.940 1.63e-12 ***
## 有效氮 -6.534e-02 2.342e-02 -2.790 0.00619 **
## 有效磷 4.132e-02 2.140e-02 1.931 0.05604 .
## 有效钾 2.914e-02 1.814e-02 1.607 0.11095
## 沙粒 9.621e-02 1.655e-02 5.814 5.74e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1471 on 113 degrees of freedom
## Multiple R-squared: 0.6553, Adjusted R-squared: 0.637
## F-statistic: 35.8 on 6 and 113 DF, p-value: < 2.2e-16
# ##
# write.csv(summ3_1$coefficients, "axis.Phyla1_env.csv")
# ##
##########
#
#
shapiro.test(resid(lm(mite.PCNM.axes[,1] ~ ., data=mite.h.1))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 1] ~ ., data = mite.h.1))
## W = 0.99155, p-value = 0.6786
#
mite.PCNM.axis1.div <- lm(mite.PCNM.axes[,1]~., data=mite.h.1)
mite.PCNM.axis1.div <- step(mite.PCNM.axis1.div)
## Start: AIC=-517.24
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes + p__Cyanobacteria
##
## Df Sum of Sq RSS AIC
## - p__Cyanobacteria 1 0.000971 1.3204 -519.15
## - p__Bacteroidetes 1 0.001285 1.3207 -519.12
## - p__Actinobacteria 1 0.017184 1.3366 -517.68
## - p__Nitrospirae 1 0.019700 1.3391 -517.46
## - p__Planctomycetes 1 0.021636 1.3411 -517.28
## <none> 1.3194 -517.24
## - p__Chloroflexi 1 0.024236 1.3437 -517.05
## - p__Proteobacteria 1 0.039909 1.3593 -515.66
## - p__Crenarchaeota 1 0.060762 1.3802 -513.83
## - p__Acidobacteria 1 0.090342 1.4098 -511.29
## - p__Verrucomicrobia 1 0.097354 1.4168 -510.69
## - p__Gemmatimonadetes 1 0.152605 1.4720 -506.10
##
## Step: AIC=-519.15
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Bacteroidetes 1 0.000530 1.3209 -521.10
## - p__Actinobacteria 1 0.018183 1.3386 -519.51
## <none> 1.3204 -519.15
## - p__Planctomycetes 1 0.024330 1.3447 -518.96
## - p__Nitrospirae 1 0.024631 1.3450 -518.93
## - p__Chloroflexi 1 0.033504 1.3539 -518.14
## - p__Proteobacteria 1 0.050945 1.3713 -516.60
## - p__Crenarchaeota 1 0.082705 1.4031 -513.86
## - p__Acidobacteria 1 0.114741 1.4351 -511.15
## - p__Verrucomicrobia 1 0.151999 1.4724 -508.07
## - p__Gemmatimonadetes 1 0.158459 1.4788 -507.55
##
## Step: AIC=-521.1
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Actinobacteria 1 0.02059 1.3415 -521.24
## <none> 1.3209 -521.10
## - p__Planctomycetes 1 0.03880 1.3597 -519.63
## - p__Nitrospirae 1 0.05172 1.3726 -518.49
## - p__Chloroflexi 1 0.10048 1.4214 -514.30
## - p__Gemmatimonadetes 1 0.16622 1.4871 -508.88
## - p__Verrucomicrobia 1 0.28386 1.6048 -499.74
## - p__Proteobacteria 1 0.28846 1.6094 -499.40
## - p__Acidobacteria 1 0.51537 1.8363 -483.57
## - p__Crenarchaeota 1 0.58735 1.9083 -478.96
##
## Step: AIC=-521.24
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Chloroflexi + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## <none> 1.3415 -521.24
## - p__Planctomycetes 1 0.03272 1.3742 -520.35
## - p__Nitrospirae 1 0.03333 1.3748 -520.30
## - p__Chloroflexi 1 0.11849 1.4600 -513.09
## - p__Gemmatimonadetes 1 0.25504 1.5965 -502.36
## - p__Verrucomicrobia 1 0.26332 1.6048 -501.74
## - p__Proteobacteria 1 0.27710 1.6186 -500.71
## - p__Acidobacteria 1 0.50210 1.8436 -485.09
## - p__Crenarchaeota 1 0.61843 1.9599 -477.75
(summ3_2 <- summary(mite.PCNM.axis1.div))
##
## Call:
## lm(formula = mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Chloroflexi + p__Gemmatimonadetes,
## data = mite.h.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.30828 -0.07067 -0.00210 0.05748 0.32472
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.325e-17 1.004e-02 0.000 1.00000
## p__Proteobacteria -2.453e-01 5.122e-02 -4.788 5.23e-06 ***
## p__Crenarchaeota -3.146e-01 4.398e-02 -7.153 9.61e-11 ***
## p__Acidobacteria -1.782e-01 2.765e-02 -6.446 3.07e-09 ***
## p__Nitrospirae -2.650e-02 1.596e-02 -1.661 0.09962 .
## p__Verrucomicrobia -7.507e-02 1.608e-02 -4.668 8.57e-06 ***
## p__Planctomycetes -3.015e-02 1.833e-02 -1.645 0.10271
## p__Chloroflexi 5.403e-02 1.726e-02 3.131 0.00223 **
## p__Gemmatimonadetes -6.013e-02 1.309e-02 -4.594 1.16e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1099 on 111 degrees of freedom
## Multiple R-squared: 0.8109, Adjusted R-squared: 0.7972
## F-statistic: 59.48 on 8 and 111 DF, p-value: < 2.2e-16
# ##
# write.csv(summ3_2$coefficients, "axis.Phyla1_div.csv")
# ##
#
###################################
#
s.value(mite.xy, mite.PCNM.axes[,2]) # ade4 function: s.value
# ##
#
shapiro.test(resid(lm(mite.PCNM.axes[,2] ~ ., data=mite.env))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 2] ~ ., data = mite.env))
## W = 0.98235, p-value = 0.1176
#
mite.PCNM.axis2.env <- lm(mite.PCNM.axes[,2] ~ ., data=mite.env)
mite.PCNM.axis2.env <- step(mite.PCNM.axis2.env)
## Start: AIC=-501.28
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有效磷 1 0.000064 1.5325 -503.27
## - 有机质 1 0.002885 1.5353 -503.05
## - 全氮 1 0.003161 1.5356 -503.03
## - 全钾 1 0.006406 1.5388 -502.78
## - 有效氮 1 0.022314 1.5547 -501.54
## - 全磷 1 0.022909 1.5553 -501.50
## <none> 1.5324 -501.28
## - pH 1 0.032019 1.5644 -500.80
## - 有效钾 1 0.083550 1.6160 -496.91
## - 沙粒 1 0.187151 1.7196 -489.45
## - 含水量 1 0.285786 1.8182 -482.76
##
## Step: AIC=-503.27
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全氮 1 0.003202 1.5357 -505.02
## - 有机质 1 0.003256 1.5357 -505.02
## - 全钾 1 0.006358 1.5388 -504.78
## <none> 1.5325 -503.27
## - 全磷 1 0.027743 1.5602 -503.12
## - 有效氮 1 0.031222 1.5637 -502.85
## - pH 1 0.033018 1.5655 -502.71
## - 有效钾 1 0.085811 1.6183 -498.73
## - 沙粒 1 0.188062 1.7206 -491.38
## - 含水量 1 0.287381 1.8199 -484.65
##
## Step: AIC=-505.02
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 有机质 + 全磷 + 全钾 + 有效氮 +
## 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有机质 1 0.000858 1.5366 -506.95
## - 全钾 1 0.012841 1.5485 -506.02
## <none> 1.5357 -505.02
## - pH 1 0.030655 1.5663 -504.65
## - 有效氮 1 0.030807 1.5665 -504.64
## - 全磷 1 0.033471 1.5692 -504.43
## - 有效钾 1 0.085555 1.6212 -500.52
## - 沙粒 1 0.185332 1.7210 -493.35
## - 含水量 1 0.287013 1.8227 -486.46
##
## Step: AIC=-506.95
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 全磷 + 全钾 + 有效氮 + 有效钾 +
## 沙粒
##
## Df Sum of Sq RSS AIC
## - 全钾 1 0.013431 1.5500 -507.91
## <none> 1.5366 -506.95
## - pH 1 0.030128 1.5667 -506.62
## - 有效氮 1 0.032519 1.5691 -506.44
## - 全磷 1 0.035024 1.5716 -506.25
## - 有效钾 1 0.092904 1.6295 -501.91
## - 沙粒 1 0.194585 1.7311 -494.65
## - 含水量 1 0.286248 1.8228 -488.45
##
## Step: AIC=-507.91
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 全磷 + 有效氮 + 有效钾 +
## 沙粒
##
## Df Sum of Sq RSS AIC
## - pH 1 0.021958 1.5719 -508.22
## <none> 1.5500 -507.91
## - 全磷 1 0.026650 1.5766 -507.86
## - 有效氮 1 0.033562 1.5835 -507.34
## - 有效钾 1 0.092875 1.6428 -502.93
## - 沙粒 1 0.184327 1.7343 -496.43
## - 含水量 1 0.292288 1.8423 -489.18
##
## Step: AIC=-508.22
## mite.PCNM.axes[, 2] ~ 含水量 + 全磷 + 有效氮 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## <none> 1.5719 -508.22
## - 全磷 1 0.03419 1.6061 -507.64
## - 有效氮 1 0.03492 1.6068 -507.59
## - 有效钾 1 0.09786 1.6698 -502.97
## - 沙粒 1 0.16238 1.7343 -498.43
## - 含水量 1 0.31644 1.8884 -488.21
(summ4_1 <- summary(mite.PCNM.axis2.env))
##
## Call:
## lm(formula = mite.PCNM.axes[, 2] ~ 含水量 + 全磷 + 有效氮 + 有效钾 +
## 沙粒, data = mite.env)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.27665 -0.06653 0.01080 0.07280 0.29903
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.00375 0.01073 -0.349 0.727426
## 含水量 -0.06138 0.01281 -4.791 5.05e-06 ***
## 全磷 0.01907 0.01211 1.575 0.118093
## 有效氮 -0.02507 0.01575 -1.591 0.114295
## 有效钾 -0.03817 0.01433 -2.664 0.008840 **
## 沙粒 0.04278 0.01247 3.432 0.000837 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1174 on 114 degrees of freedom
## Multiple R-squared: 0.5397, Adjusted R-squared: 0.5195
## F-statistic: 26.73 on 5 and 114 DF, p-value: < 2.2e-16
# ##
# write.csv(summ4_1$coefficients, "axis.Phyla2_env.csv")
# ##
#
#
shapiro.test(resid(lm(mite.PCNM.axes[,2] ~ ., data=mite.h.1))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 2] ~ ., data = mite.h.1))
## W = 0.96742, p-value = 0.00523
#
mite.PCNM.axis2.div <- lm(mite.PCNM.axes[,2]~., data=mite.h.1)
mite.PCNM.axis2.div <- step(mite.PCNM.axis2.div)
## Start: AIC=-531.71
## mite.PCNM.axes[, 2] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes + p__Cyanobacteria
##
## Df Sum of Sq RSS AIC
## - p__Cyanobacteria 1 0.000000 1.1695 -533.71
## - p__Planctomycetes 1 0.000115 1.1696 -533.70
## - p__Actinobacteria 1 0.000298 1.1698 -533.68
## - p__Proteobacteria 1 0.000579 1.1701 -533.65
## - p__Chloroflexi 1 0.001075 1.1705 -533.60
## - p__Bacteroidetes 1 0.001435 1.1709 -533.57
## - p__Crenarchaeota 1 0.002640 1.1721 -533.44
## - p__Acidobacteria 1 0.012279 1.1818 -532.46
## <none> 1.1695 -531.71
## - p__Verrucomicrobia 1 0.034719 1.2042 -530.20
## - p__Nitrospirae 1 0.040799 1.2103 -529.60
## - p__Gemmatimonadetes 1 0.080957 1.2504 -525.68
##
## Step: AIC=-533.71
## mite.PCNM.axes[, 2] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Planctomycetes 1 0.000176 1.1696 -535.70
## - p__Actinobacteria 1 0.000426 1.1699 -535.67
## - p__Proteobacteria 1 0.000926 1.1704 -535.62
## - p__Chloroflexi 1 0.001295 1.1708 -535.58
## - p__Bacteroidetes 1 0.001996 1.1715 -535.51
## - p__Crenarchaeota 1 0.004282 1.1738 -535.27
## - p__Acidobacteria 1 0.017423 1.1869 -533.94
## <none> 1.1695 -533.71
## - p__Verrucomicrobia 1 0.062460 1.2319 -529.47
## - p__Nitrospirae 1 0.069113 1.2386 -528.82
## - p__Gemmatimonadetes 1 0.087881 1.2573 -527.02
##
## Step: AIC=-535.7
## mite.PCNM.axes[, 2] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Actinobacteria + p__Chloroflexi + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Actinobacteria 1 0.000273 1.1699 -537.67
## - p__Proteobacteria 1 0.000846 1.1705 -537.61
## - p__Chloroflexi 1 0.001174 1.1708 -537.57
## - p__Bacteroidetes 1 0.002479 1.1721 -537.44
## - p__Crenarchaeota 1 0.005699 1.1753 -537.11
## <none> 1.1696 -535.70
## - p__Acidobacteria 1 0.027126 1.1968 -534.94
## - p__Verrucomicrobia 1 0.069242 1.2389 -530.79
## - p__Nitrospirae 1 0.092604 1.2623 -528.55
## - p__Gemmatimonadetes 1 0.094696 1.2643 -528.35
##
## Step: AIC=-537.67
## mite.PCNM.axes[, 2] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Chloroflexi + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Proteobacteria 1 0.000581 1.1705 -539.61
## - p__Chloroflexi 1 0.000932 1.1708 -539.57
## - p__Bacteroidetes 1 0.002207 1.1721 -539.44
## - p__Crenarchaeota 1 0.005964 1.1759 -539.06
## <none> 1.1699 -537.67
## - p__Acidobacteria 1 0.037475 1.2074 -535.88
## - p__Verrucomicrobia 1 0.080583 1.2505 -531.67
## - p__Gemmatimonadetes 1 0.102018 1.2719 -529.63
## - p__Nitrospirae 1 0.140885 1.3108 -526.02
##
## Step: AIC=-539.61
## mite.PCNM.axes[, 2] ~ p__Crenarchaeota + p__Acidobacteria + p__Bacteroidetes +
## p__Nitrospirae + p__Verrucomicrobia + p__Chloroflexi + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Chloroflexi 1 0.00035 1.1708 -541.57
## - p__Bacteroidetes 1 0.00250 1.1730 -541.35
## <none> 1.1705 -539.61
## - p__Crenarchaeota 1 0.08773 1.2582 -532.93
## - p__Gemmatimonadetes 1 0.12155 1.2920 -529.75
## - p__Verrucomicrobia 1 0.20747 1.3780 -522.03
## - p__Acidobacteria 1 0.55487 1.7254 -495.05
## - p__Nitrospirae 1 0.57031 1.7408 -493.98
##
## Step: AIC=-541.57
## mite.PCNM.axes[, 2] ~ p__Crenarchaeota + p__Acidobacteria + p__Bacteroidetes +
## p__Nitrospirae + p__Verrucomicrobia + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Bacteroidetes 1 0.00284 1.1737 -543.28
## <none> 1.1708 -541.57
## - p__Crenarchaeota 1 0.09167 1.2625 -534.53
## - p__Gemmatimonadetes 1 0.12611 1.2970 -531.30
## - p__Verrucomicrobia 1 0.20801 1.3789 -523.95
## - p__Acidobacteria 1 0.57171 1.7426 -495.86
## - p__Nitrospirae 1 0.59191 1.7628 -494.47
##
## Step: AIC=-543.28
## mite.PCNM.axes[, 2] ~ p__Crenarchaeota + p__Acidobacteria + p__Nitrospirae +
## p__Verrucomicrobia + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## <none> 1.1737 -543.28
## - p__Crenarchaeota 1 0.11112 1.2848 -534.43
## - p__Gemmatimonadetes 1 0.12578 1.2995 -533.06
## - p__Verrucomicrobia 1 0.21803 1.3917 -524.83
## - p__Nitrospirae 1 0.58912 1.7628 -496.47
## - p__Acidobacteria 1 0.63416 1.8078 -493.44
(summ4_2 <- summary(mite.PCNM.axis2.div))
##
## Call:
## lm(formula = mite.PCNM.axes[, 2] ~ p__Crenarchaeota + p__Acidobacteria +
## p__Nitrospirae + p__Verrucomicrobia + p__Gemmatimonadetes,
## data = mite.h.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.37396 -0.05211 0.01153 0.06765 0.23718
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.280e-18 9.263e-03 0.000 1.000000
## p__Crenarchaeota -4.001e-02 1.218e-02 -3.285 0.001354 **
## p__Acidobacteria 9.857e-02 1.256e-02 7.848 2.51e-12 ***
## p__Nitrospirae 7.690e-02 1.017e-02 7.564 1.08e-11 ***
## p__Verrucomicrobia -4.843e-02 1.052e-02 -4.602 1.09e-05 ***
## p__Gemmatimonadetes 4.011e-02 1.147e-02 3.495 0.000676 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1015 on 114 degrees of freedom
## Multiple R-squared: 0.6563, Adjusted R-squared: 0.6412
## F-statistic: 43.54 on 5 and 114 DF, p-value: < 2.2e-16
# ##
# write.csv(summ4_2$coefficients, "axis.Phyla2_div.csv")
# ##
#
# windows(title="10 PCNM variables - mites")
# par(mfrow=c(2,5))
# for(i in mite.PCNM.fwd[,2][1:10] ){
# s.value(mite.xy, mite.PCNM.pos[,i], sub=i, csub=2)
# }
2.3 rda(mite.h, mite_Env.Pcnm, mite_y.Pcnm)
###################################################
#
#
mite.h.1 <- residuals(lm(mite.h~PCNM_y.pos))
mite.h.1 <- as.data.frame(mite.h.1)
#
#####################################################
# #
# mite.env.rda <- rda(mite.h.1, mite.env)
# #
# (mite.R2a <- RsquareAdj(mite.env.rda)$adj.r.squared)
# #
# (mite.Env.fwd <- forward.sel(mite.h.1, mite.env,
# nperm = 9999,
# adjR2thresh=mite.R2a,
# alpha=0.1))
# #
# (nb.sig.Env <- nrow(mite.Env.fwd)) # Number of
# (Env.sign <- sort(mite.Env.fwd[,2]))
# env.red <- mite.env[,c(Env.sign)]
# #
# # head(env.red)
#
###################################################
#
mite.PCNM.rda <- rda(mite.h.1, mite.PCNM.pos)
#
(mite.R2a <- RsquareAdj(mite.PCNM.rda)$adj.r.squared)
## [1] 0.2801768
#
(mite.PCNM.fwd <- forward.sel(mite.h.1, mite.PCNM.pos,
nperm = 9999,
adjR2thresh=mite.R2a,
alpha=0.1))
## Testing variable 1
## Testing variable 2
## Testing variable 3
## Testing variable 4
## Testing variable 5
## Testing variable 6
## Testing variable 7
## Testing variable 8
## Testing variable 9
## Testing variable 10
## Testing variable 11
## Testing variable 12
## Testing variable 13
## Procedure stopped (adjR2thresh criteria) adjR2cum = 0.283494 with 13 variables (superior to 0.280177)
## variables order R2 R2Cum AdjR2Cum F pval
## 1 PCNM1 1 0.06298303 0.06298303 0.05504221 7.931551 0.0001
## 2 PCNM4 4 0.05702563 0.12000866 0.10496607 7.581891 0.0001
## 3 PCNM11 11 0.03576478 0.15577344 0.13394000 4.914219 0.0001
## 4 PCNM6 6 0.02753981 0.18331325 0.15490676 3.877960 0.0010
## 5 PCNM3 3 0.02495350 0.20826676 0.17354162 3.593002 0.0016
## 6 PCNM7 7 0.02369384 0.23196060 0.19117975 3.486024 0.0026
## 7 PCNM10 10 0.02141497 0.25337557 0.20671154 3.212427 0.0023
## 8 PCNM2 2 0.01955874 0.27293431 0.22053318 2.986003 0.0058
## 9 PCNM8 8 0.01929248 0.29222679 0.23431807 2.998379 0.0049
## 10 PCNM48 48 0.01777795 0.31000473 0.24670242 2.808420 0.0095
## 11 PCNM13 13 0.01777198 0.32777671 0.25930953 2.855263 0.0098
## 12 PCNM14 14 0.01705488 0.34483160 0.27135477 2.785349 0.0110
## 13 PCNM23 23 0.01693589 0.36176748 0.28349368 2.812774 0.0077
#
(nb.sig.PCNM <- nrow(mite.PCNM.fwd)) # Number of
## [1] 13
(PCNM.sign <- sort(mite.PCNM.fwd[,2]))
## [1] 1 2 3 4 6 7 8 10 11 13 14 23 48
PCNM.red <- mite.PCNM.pos[,c(PCNM.sign)]
#
# # head(PCNM.red)
#########################
#
env.red <- mite.env
#
PCNM.env_red <- cbind(PCNM.red, env.red, mite.xy)
PCNM.env_red <- apply(PCNM.env_red , 2, as.numeric)
PCNM.env_red <- as.data.frame(PCNM.env_red )
#
# # head(PCNM.env_red)
#
mite.PCNM.rda2 <- vegan::rda(mite.h.1 ~., data=PCNM.env_red )
axes.test <- anova.cca(mite.PCNM.rda2, by="axis")
(nb.ax <- length(which(axes.test[,4] <= 0.05)))
## [1] 7
#
#########################
#
mod <- mite.PCNM.rda2
#
plot(mod, type="n", xlim=c(-3,3), ylim=c(-3,3))
# points(mod, pch=21, col="red", bg="yellow", cex=0.8)
text(mod, "species", col="blue", cex=0.8)
n1 = dim(PCNM.red)[2]
n2 = dim(PCNM.env_red)[2]-n1
text(mod, dis="cn",
col=c(rep(3,n1),rep(2,n2)),
cex=c(rep(0.5,n1),rep(1,n2)),
lty=c(rep(3,n1),rep(1,n2))
)
#
##########################
#
####################################################
2.3_xy
##################################################
#
mite.PCNM.axes <- scores(mite.PCNM.rda2, choices=c(1,2,3,4,5), display="lc",
scaling=1)
#
#####
PCNM_phyla <- cbind(mite.xy, mite.PCNM.axes[,1:2])
# write.csv(PCNM_phyla, "PCNM_axies_phyla.csv")
#####
#
s.value(mite.xy, mite.PCNM.axes[,1]) # ade4 function: s.value
#
#
shapiro.test(resid(lm(mite.PCNM.axes[,1] ~ ., data=mite.env))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 1] ~ ., data = mite.env))
## W = 0.98952, p-value = 0.4927
#
mite.PCNM.axis1.env <- lm(mite.PCNM.axes[,1]~., data=mite.env)
mite.PCNM.axis1.env <- step(mite.PCNM.axis1.env)
## Start: AIC=-372.09
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 有效磷 1 0.00066 4.4977 -374.07
## - 全磷 1 0.00180 4.4988 -374.04
## - 沙粒 1 0.00777 4.5048 -373.88
## - 全钾 1 0.01643 4.5135 -373.65
## - 有机质 1 0.04816 4.5452 -372.81
## - 有效氮 1 0.06196 4.5590 -372.45
## <none> 4.4970 -372.09
## - 全氮 1 0.20219 4.6992 -368.81
## - pH 1 0.27604 4.7731 -366.94
## - 含水量 1 0.35716 4.8542 -364.92
## - 有效钾 1 0.42293 4.9200 -363.30
##
## Step: AIC=-374.07
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 全磷 1 0.00316 4.5009 -375.99
## - 沙粒 1 0.00825 4.5060 -375.85
## - 全钾 1 0.01615 4.5139 -375.64
## - 有机质 1 0.04813 4.5458 -374.79
## <none> 4.4977 -374.07
## - 有效氮 1 0.08167 4.5794 -373.91
## - 全氮 1 0.20159 4.6993 -370.81
## - pH 1 0.29679 4.7945 -368.40
## - 含水量 1 0.36296 4.8607 -366.76
## - 有效钾 1 0.44572 4.9434 -364.73
##
## Step: AIC=-375.99
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全钾 + 有效氮 +
## 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 沙粒 1 0.00686 4.5077 -377.80
## - 全钾 1 0.01299 4.5139 -377.64
## - 有机质 1 0.05505 4.5559 -376.53
## <none> 4.5009 -375.99
## - 有效氮 1 0.07862 4.5795 -375.91
## - 全氮 1 0.20012 4.7010 -372.77
## - pH 1 0.29823 4.7991 -370.29
## - 含水量 1 0.35986 4.8607 -368.76
## - 有效钾 1 0.44395 4.9448 -366.70
##
## Step: AIC=-377.8
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 全钾 + 有效氮 +
## 有效钾
##
## Df Sum of Sq RSS AIC
## - 全钾 1 0.01446 4.5222 -379.42
## - 有机质 1 0.05732 4.5650 -378.29
## - 有效氮 1 0.07192 4.5796 -377.90
## <none> 4.5077 -377.80
## - 全氮 1 0.19705 4.7048 -374.67
## - pH 1 0.30149 4.8092 -372.04
## - 含水量 1 0.36692 4.8746 -370.41
## - 有效钾 1 0.46178 4.9695 -368.10
##
## Step: AIC=-379.42
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 有机质 + 全氮 + 有效氮 +
## 有效钾
##
## Df Sum of Sq RSS AIC
## - 有机质 1 0.06232 4.5845 -379.78
## <none> 4.5222 -379.42
## - 有效氮 1 0.09014 4.6123 -379.05
## - 全氮 1 0.18476 4.7069 -376.61
## - 含水量 1 0.36016 4.8823 -372.22
## - pH 1 0.36579 4.8880 -372.09
## - 有效钾 1 0.45051 4.9727 -370.02
##
## Step: AIC=-379.78
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 全氮 + 有效氮 + 有效钾
##
## Df Sum of Sq RSS AIC
## - 有效氮 1 0.07175 4.6563 -379.91
## <none> 4.5845 -379.78
## - 全氮 1 0.13406 4.7186 -378.32
## - pH 1 0.30597 4.8905 -374.02
## - 含水量 1 0.31929 4.9038 -373.70
## - 有效钾 1 0.51778 5.1023 -368.94
##
## Step: AIC=-379.91
## mite.PCNM.axes[, 1] ~ 含水量 + pH + 全氮 + 有效钾
##
## Df Sum of Sq RSS AIC
## <none> 4.6563 -379.91
## - 全氮 1 0.17737 4.8336 -377.43
## - 含水量 1 0.26091 4.9172 -375.37
## - pH 1 0.35543 5.0117 -373.09
## - 有效钾 1 0.44992 5.1062 -370.84
(summ3_1 <- summary(mite.PCNM.axis1.env))
##
## Call:
## lm(formula = mite.PCNM.axes[, 1] ~ 含水量 + pH + 全氮 + 有效钾,
## data = mite.env)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.45532 -0.14062 0.00415 0.16177 0.48402
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.002076 0.018389 0.113 0.91032
## 含水量 0.053841 0.021210 2.538 0.01247 *
## pH 0.056121 0.018942 2.963 0.00371 **
## 全氮 -0.042015 0.020074 -2.093 0.03855 *
## 有效钾 0.073023 0.021906 3.333 0.00115 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2012 on 115 degrees of freedom
## Multiple R-squared: 0.2542, Adjusted R-squared: 0.2283
## F-statistic: 9.799 on 4 and 115 DF, p-value: 7.4e-07
# ##
# write.csv(summ3_1$coefficients, "axis.Phyla1_env.csv")
# ##
##########
#
#
shapiro.test(resid(lm(mite.PCNM.axes[,1] ~ ., data=mite.h.1))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 1] ~ ., data = mite.h.1))
## W = 0.99055, p-value = 0.584
#
mite.PCNM.axis1.div <- lm(mite.PCNM.axes[,1]~., data=mite.h.1)
mite.PCNM.axis1.div <- step(mite.PCNM.axis1.div)
## Start: AIC=-482.46
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes + p__Cyanobacteria
##
## Df Sum of Sq RSS AIC
## - p__Cyanobacteria 1 0.000031 1.7630 -484.46
## - p__Chloroflexi 1 0.017523 1.7805 -483.27
## - p__Bacteroidetes 1 0.022412 1.7854 -482.94
## - p__Nitrospirae 1 0.023525 1.7865 -482.87
## <none> 1.7629 -482.46
## - p__Verrucomicrobia 1 0.045433 1.8084 -481.41
## - p__Gemmatimonadetes 1 0.053594 1.8165 -480.87
## - p__Proteobacteria 1 0.077843 1.8408 -479.28
## - p__Actinobacteria 1 0.085619 1.8486 -478.77
## - p__Planctomycetes 1 0.090528 1.8535 -478.45
## - p__Crenarchaeota 1 0.111771 1.8747 -477.08
## - p__Acidobacteria 1 0.145815 1.9088 -474.92
##
## Step: AIC=-484.46
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Chloroflexi 1 0.019270 1.7822 -485.15
## - p__Bacteroidetes 1 0.027480 1.7905 -484.60
## <none> 1.7630 -484.46
## - p__Nitrospirae 1 0.036768 1.7997 -483.98
## - p__Gemmatimonadetes 1 0.057521 1.8205 -482.61
## - p__Verrucomicrobia 1 0.072754 1.8357 -481.61
## - p__Actinobacteria 1 0.107748 1.8707 -479.34
## - p__Proteobacteria 1 0.109943 1.8729 -479.20
## - p__Planctomycetes 1 0.122180 1.8852 -478.42
## - p__Crenarchaeota 1 0.164052 1.9270 -475.78
## - p__Acidobacteria 1 0.193942 1.9569 -473.93
##
## Step: AIC=-485.15
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Bacteroidetes 1 0.00827 1.7905 -486.60
## - p__Nitrospirae 1 0.01757 1.7998 -485.98
## <none> 1.7822 -485.15
## - p__Gemmatimonadetes 1 0.04855 1.8308 -483.93
## - p__Verrucomicrobia 1 0.05789 1.8401 -483.32
## - p__Actinobacteria 1 0.08880 1.8710 -481.32
## - p__Planctomycetes 1 0.10727 1.8895 -480.14
## - p__Proteobacteria 1 0.17219 1.9544 -476.09
## - p__Acidobacteria 1 0.31697 2.0992 -467.51
## - p__Crenarchaeota 1 0.31856 2.1008 -467.42
##
## Step: AIC=-486.6
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Nitrospirae 1 0.01033 1.8009 -487.91
## <none> 1.7905 -486.60
## - p__Gemmatimonadetes 1 0.04236 1.8329 -485.79
## - p__Verrucomicrobia 1 0.04977 1.8403 -485.31
## - p__Actinobacteria 1 0.08061 1.8711 -483.31
## - p__Planctomycetes 1 0.11411 1.9046 -481.18
## - p__Proteobacteria 1 0.32722 2.1177 -468.46
## - p__Acidobacteria 1 0.53826 2.3288 -457.06
## - p__Crenarchaeota 1 0.80587 2.5964 -444.00
##
## Step: AIC=-487.91
## mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Verrucomicrobia + p__Planctomycetes +
## p__Actinobacteria + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## <none> 1.8009 -487.91
## - p__Verrucomicrobia 1 0.03970 1.8405 -487.29
## - p__Gemmatimonadetes 1 0.04457 1.8454 -486.97
## - p__Actinobacteria 1 0.07032 1.8712 -485.31
## - p__Planctomycetes 1 0.13411 1.9350 -481.29
## - p__Proteobacteria 1 0.44896 2.2498 -463.20
## - p__Acidobacteria 1 0.80670 2.6076 -445.49
## - p__Crenarchaeota 1 1.44981 3.2507 -419.04
(summ3_2 <- summary(mite.PCNM.axis1.div))
##
## Call:
## lm(formula = mite.PCNM.axes[, 1] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Verrucomicrobia + p__Planctomycetes +
## p__Actinobacteria + p__Gemmatimonadetes, data = mite.h.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.25067 -0.08661 -0.00565 0.08366 0.38639
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.216e-17 1.158e-02 0.000 1.00000
## p__Proteobacteria -2.172e-01 4.111e-02 -5.284 6.29e-07 ***
## p__Crenarchaeota -2.968e-01 3.126e-02 -9.496 4.82e-16 ***
## p__Acidobacteria -1.608e-01 2.270e-02 -7.083 1.32e-10 ***
## p__Verrucomicrobia -2.582e-02 1.643e-02 -1.571 0.11893
## p__Planctomycetes -5.805e-02 2.010e-02 -2.888 0.00465 **
## p__Actinobacteria -3.392e-02 1.622e-02 -2.091 0.03876 *
## p__Gemmatimonadetes -2.673e-02 1.605e-02 -1.665 0.09874 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1268 on 112 degrees of freedom
## Multiple R-squared: 0.7116, Adjusted R-squared: 0.6935
## F-statistic: 39.47 on 7 and 112 DF, p-value: < 2.2e-16
# ##
# write.csv(summ3_2$coefficients, "axis.Phyla1_div.csv")
# ##
#
###################################
#
s.value(mite.xy, mite.PCNM.axes[,2]) # ade4 function: s.value
# ##
#
shapiro.test(resid(lm(mite.PCNM.axes[,2] ~ ., data=mite.env))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 2] ~ ., data = mite.env))
## W = 0.98681, p-value = 0.2964
#
mite.PCNM.axis2.env <- lm(mite.PCNM.axes[,2] ~ ., data=mite.env)
mite.PCNM.axis2.env <- step(mite.PCNM.axis2.env)
## Start: AIC=-442.49
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效磷 + 有效钾 + 沙粒
##
## Df Sum of Sq RSS AIC
## - 沙粒 1 0.000021 2.5011 -444.49
## - 全磷 1 0.000143 2.5012 -444.48
## - 有机质 1 0.011802 2.5129 -443.93
## - 有效氮 1 0.027003 2.5281 -443.20
## - 有效磷 1 0.027828 2.5289 -443.16
## - 有效钾 1 0.037896 2.5390 -442.69
## <none> 2.5011 -442.49
## - 全钾 1 0.053119 2.5542 -441.97
## - 全氮 1 0.072874 2.5740 -441.05
## - 含水量 1 0.118091 2.6192 -438.96
## - pH 1 0.241448 2.7425 -433.43
##
## Step: AIC=-444.49
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 有机质 + 全氮 + 全磷 + 全钾 +
## 有效氮 + 有效磷 + 有效钾
##
## Df Sum of Sq RSS AIC
## - 全磷 1 0.000128 2.5012 -446.48
## - 有机质 1 0.011822 2.5129 -445.92
## - 有效氮 1 0.027589 2.5287 -445.17
## - 有效磷 1 0.028198 2.5293 -445.15
## - 有效钾 1 0.038320 2.5394 -444.67
## <none> 2.5011 -444.49
## - 全钾 1 0.053108 2.5542 -443.97
## - 全氮 1 0.073495 2.5746 -443.02
## - 含水量 1 0.118309 2.6194 -440.94
## - pH 1 0.272420 2.7735 -434.08
##
## Step: AIC=-446.48
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 有机质 + 全氮 + 全钾 + 有效氮 +
## 有效磷 + 有效钾
##
## Df Sum of Sq RSS AIC
## - 有机质 1 0.011951 2.5132 -447.91
## - 有效氮 1 0.028025 2.5293 -447.15
## - 有效磷 1 0.034207 2.5354 -446.85
## - 有效钾 1 0.038322 2.5396 -446.66
## <none> 2.5012 -446.48
## - 全钾 1 0.061238 2.5625 -445.58
## - 全氮 1 0.076860 2.5781 -444.85
## - 含水量 1 0.119557 2.6208 -442.88
## - pH 1 0.272477 2.7737 -436.08
##
## Step: AIC=-447.91
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮 + 全钾 + 有效氮 + 有效磷 +
## 有效钾
##
## Df Sum of Sq RSS AIC
## - 有效钾 1 0.031912 2.5451 -448.40
## - 有效氮 1 0.037865 2.5511 -448.12
## - 有效磷 1 0.041229 2.5544 -447.96
## <none> 2.5132 -447.91
## - 全钾 1 0.059130 2.5723 -447.12
## - 全氮 1 0.075452 2.5886 -446.36
## - 含水量 1 0.135662 2.6489 -443.60
## - pH 1 0.262917 2.7761 -437.97
##
## Step: AIC=-448.4
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮 + 全钾 + 有效氮 + 有效磷
##
## Df Sum of Sq RSS AIC
## - 有效氮 1 0.024361 2.5695 -449.26
## <none> 2.5451 -448.40
## - 有效磷 1 0.053898 2.5990 -447.88
## - 全钾 1 0.055568 2.6007 -447.81
## - 全氮 1 0.093076 2.6382 -446.09
## - 含水量 1 0.176209 2.7213 -442.37
## - pH 1 0.268831 2.8139 -438.35
##
## Step: AIC=-449.26
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮 + 全钾 + 有效磷
##
## Df Sum of Sq RSS AIC
## - 有效磷 1 0.029537 2.5990 -449.88
## <none> 2.5695 -449.26
## - 全钾 1 0.056849 2.6263 -448.63
## - 全氮 1 0.076387 2.6458 -447.74
## - 含水量 1 0.155725 2.7252 -444.19
## - pH 1 0.245422 2.8149 -440.31
##
## Step: AIC=-449.88
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮 + 全钾
##
## Df Sum of Sq RSS AIC
## - 全钾 1 0.035874 2.6349 -450.24
## <none> 2.5990 -449.88
## - 全氮 1 0.116627 2.7156 -446.62
## - pH 1 0.224950 2.8239 -441.92
## - 含水量 1 0.240496 2.8395 -441.26
##
## Step: AIC=-450.24
## mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮
##
## Df Sum of Sq RSS AIC
## <none> 2.6349 -450.24
## - pH 1 0.19176 2.8266 -443.81
## - 含水量 1 0.21122 2.8461 -442.99
## - 全氮 1 0.26353 2.8984 -440.80
(summ4_1 <- summary(mite.PCNM.axis2.env))
##
## Call:
## lm(formula = mite.PCNM.axes[, 2] ~ 含水量 + pH + 全氮, data = mite.env)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35254 -0.09433 0.01132 0.10843 0.29286
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.00326 0.01377 -0.237 0.813295
## 含水量 -0.04364 0.01431 -3.049 0.002841 **
## pH 0.04108 0.01414 2.906 0.004392 **
## 全氮 -0.04965 0.01458 -3.406 0.000906 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1507 on 116 degrees of freedom
## Multiple R-squared: 0.2357, Adjusted R-squared: 0.216
## F-statistic: 11.93 on 3 and 116 DF, p-value: 7.284e-07
# ##
# write.csv(summ4_1$coefficients, "axis.Phyla2_env.csv")
# ##
#
#
shapiro.test(resid(lm(mite.PCNM.axes[,2] ~ ., data=mite.h.1))) # Normality test
##
## Shapiro-Wilk normality test
##
## data: resid(lm(mite.PCNM.axes[, 2] ~ ., data = mite.h.1))
## W = 0.9849, p-value = 0.2012
#
mite.PCNM.axis2.div <- lm(mite.PCNM.axes[,2]~., data=mite.h.1)
mite.PCNM.axis2.div <- step(mite.PCNM.axis2.div)
## Start: AIC=-509.9
## mite.PCNM.axes[, 2] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Bacteroidetes + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes + p__Cyanobacteria
##
## Df Sum of Sq RSS AIC
## - p__Bacteroidetes 1 0.000192 1.4028 -511.89
## - p__Proteobacteria 1 0.001554 1.4041 -511.77
## - p__Acidobacteria 1 0.002087 1.4046 -511.72
## - p__Actinobacteria 1 0.002254 1.4048 -511.71
## - p__Cyanobacteria 1 0.003665 1.4062 -511.59
## - p__Nitrospirae 1 0.007471 1.4100 -511.27
## - p__Crenarchaeota 1 0.008568 1.4111 -511.17
## - p__Chloroflexi 1 0.010055 1.4126 -511.05
## - p__Planctomycetes 1 0.010874 1.4134 -510.98
## - p__Verrucomicrobia 1 0.022514 1.4251 -509.99
## <none> 1.4026 -509.90
## - p__Gemmatimonadetes 1 0.122353 1.5249 -501.87
##
## Step: AIC=-511.89
## mite.PCNM.axes[, 2] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Actinobacteria + p__Chloroflexi +
## p__Gemmatimonadetes + p__Cyanobacteria
##
## Df Sum of Sq RSS AIC
## - p__Actinobacteria 1 0.005191 1.4080 -513.44
## - p__Cyanobacteria 1 0.005681 1.4084 -513.40
## - p__Proteobacteria 1 0.006702 1.4095 -513.31
## <none> 1.4028 -511.89
## - p__Acidobacteria 1 0.025302 1.4281 -511.74
## - p__Chloroflexi 1 0.025668 1.4284 -511.71
## - p__Nitrospirae 1 0.031982 1.4347 -511.18
## - p__Planctomycetes 1 0.036240 1.4390 -510.83
## - p__Verrucomicrobia 1 0.048201 1.4510 -509.83
## - p__Crenarchaeota 1 0.068352 1.4711 -508.18
## - p__Gemmatimonadetes 1 0.142804 1.5456 -502.25
##
## Step: AIC=-513.44
## mite.PCNM.axes[, 2] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Chloroflexi + p__Gemmatimonadetes +
## p__Cyanobacteria
##
## Df Sum of Sq RSS AIC
## - p__Cyanobacteria 1 0.003440 1.4114 -515.15
## - p__Proteobacteria 1 0.017572 1.4255 -513.95
## - p__Acidobacteria 1 0.020335 1.4283 -513.72
## <none> 1.4080 -513.44
## - p__Nitrospirae 1 0.026900 1.4348 -513.17
## - p__Chloroflexi 1 0.027388 1.4353 -513.13
## - p__Planctomycetes 1 0.034402 1.4424 -512.55
## - p__Verrucomicrobia 1 0.074762 1.4827 -509.23
## - p__Crenarchaeota 1 0.120294 1.5282 -505.60
## - p__Gemmatimonadetes 1 0.184042 1.5920 -500.70
##
## Step: AIC=-515.15
## mite.PCNM.axes[, 2] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Acidobacteria + p__Nitrospirae + p__Verrucomicrobia +
## p__Planctomycetes + p__Chloroflexi + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Acidobacteria 1 0.018112 1.4295 -515.62
## - p__Nitrospirae 1 0.023643 1.4350 -515.16
## <none> 1.4114 -515.15
## - p__Chloroflexi 1 0.025641 1.4370 -514.99
## - p__Proteobacteria 1 0.026425 1.4378 -514.92
## - p__Planctomycetes 1 0.031006 1.4424 -514.54
## - p__Verrucomicrobia 1 0.109359 1.5208 -508.19
## - p__Crenarchaeota 1 0.175464 1.5869 -503.09
## - p__Gemmatimonadetes 1 0.183846 1.5952 -502.46
##
## Step: AIC=-515.62
## mite.PCNM.axes[, 2] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Nitrospirae + p__Verrucomicrobia + p__Planctomycetes +
## p__Chloroflexi + p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## - p__Nitrospirae 1 0.00742 1.4369 -517.00
## <none> 1.4295 -515.62
## - p__Chloroflexi 1 0.04440 1.4739 -513.95
## - p__Planctomycetes 1 0.05465 1.4842 -513.12
## - p__Proteobacteria 1 0.16120 1.5907 -504.80
## - p__Gemmatimonadetes 1 0.18871 1.6182 -502.74
## - p__Verrucomicrobia 1 0.20427 1.6338 -501.59
## - p__Crenarchaeota 1 0.52326 1.9528 -480.19
##
## Step: AIC=-517
## mite.PCNM.axes[, 2] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Verrucomicrobia + p__Planctomycetes + p__Chloroflexi +
## p__Gemmatimonadetes
##
## Df Sum of Sq RSS AIC
## <none> 1.4369 -517.00
## - p__Chloroflexi 1 0.05056 1.4875 -514.85
## - p__Planctomycetes 1 0.05602 1.4929 -514.41
## - p__Gemmatimonadetes 1 0.18143 1.6184 -504.73
## - p__Proteobacteria 1 0.19649 1.6334 -503.62
## - p__Verrucomicrobia 1 0.24350 1.6804 -500.21
## - p__Crenarchaeota 1 0.72536 2.1623 -469.96
(summ4_2 <- summary(mite.PCNM.axis2.div))
##
## Call:
## lm(formula = mite.PCNM.axes[, 2] ~ p__Proteobacteria + p__Crenarchaeota +
## p__Verrucomicrobia + p__Planctomycetes + p__Chloroflexi +
## p__Gemmatimonadetes, data = mite.h.1)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.38949 -0.08337 0.00525 0.07469 0.26828
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.133e-17 1.029e-02 0.000 1.000000
## p__Proteobacteria -1.252e-01 3.185e-02 -3.931 0.000146 ***
## p__Crenarchaeota -2.020e-01 2.675e-02 -7.553 1.19e-11 ***
## p__Verrucomicrobia -6.215e-02 1.420e-02 -4.376 2.71e-05 ***
## p__Planctomycetes 3.632e-02 1.731e-02 2.099 0.038048 *
## p__Chloroflexi -3.227e-02 1.618e-02 -1.994 0.048566 *
## p__Gemmatimonadetes 4.882e-02 1.292e-02 3.777 0.000255 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1128 on 113 degrees of freedom
## Multiple R-squared: 0.5832, Adjusted R-squared: 0.5611
## F-statistic: 26.35 on 6 and 113 DF, p-value: < 2.2e-16
# ##
# write.csv(summ4_2$coefficients, "axis.Phyla2_div.csv")
# ##
#
# windows(title="10 PCNM variables - mites")
# par(mfrow=c(2,5))
# for(i in mite.PCNM.fwd[,2][1:10] ){
# s.value(mite.xy, mite.PCNM.pos[,i], sub=i, csub=2)
# }