统计学中r检验的基本步骤：共线性检验和CAP分析
统计学中r检验的基本步骤：共线性检验和CAP分析此时我们能看到存在一个Imaginary的解释量，这意味着可能存在负引擎，为了避免这种情况我们可以把负值去除。代码如下：vare_cap <- capscale(varespec ~. Condition(Humdepth Baresoil ) env2 dist="bray") vare_cap结果表明了约束部分解释了39.21%的物种变化，控制部分解释了25.09%的物种变化，未被解释的部分有41.39%。这里经过检验我们发现铁和铝、钙和镁、钾和硫这6个变量两两间存在显著的共线性，因此我们只取共线性变量的其中之一用于接下来的CAP分析。这里我们选择铝、钙和钾这三个变量与其它变量组成环境因子矩阵。#select env %>% data.frame() %>% dplyr::select(-Fe -Mg -S) ->env2 env22.
前言：在进行物种多样性分析的过程中，排序分析有时候能够用来表征样方之间物种的β多样性（即群落的组成变化）。排序分析可以分为间接梯度排序(indirect gradient analysis)和直接梯度排序(direct gradient analysis)。间接梯度排序又叫非约束性排序，该分析是为了寻求潜在的或在间接的环境梯度下解释物种数据的变化包括PCA，PCoA，NDMS；直接排序又叫约束性排序。它是指在特定的梯度上（环境轴）上探讨物种的变化情况，包括 RDA CCA CAP或dbRDA。排序的目的是在一个可视化的低维空间或平面重新排列目标样本，使样本之间的距离最大程度地反映出平面散点图内样本之间的关系信息。
1. 共线性检验在进行排序分析之前尤其是在直接梯度排序分析前，一般需要检验使用的环境因子之间是否存在共线性问题。因为，环境因子（相当于自变量）之间如果存在严重的共线性会使物种矩阵被过度解读，从而使排序结果不准确。R语言中Hmisc包提供的varclus函数可以基于聚类分析的算法评估目标变量的共线性和冗余性。
代码如下：
setwd(choose.dir())
rm(list = ls())
library(vegan)
library(tidyverse)
library(Hmisc)
#import data
data(varespec)# species data
data(varechem)# env data
#Co-linearity analysis
env<-as.matrix(varechem)
co_linear <- varclus(env similarity="spear") # spearman is the default
co_linear
plot(co_linear)#
一般来说当spearman ρ2 > 0.7时则表明分支上的变量可能存在严重的共线性，此时应该根据研究的目的去除其中一个变量。
这里经过检验我们发现铁和铝、钙和镁、钾和硫这6个变量两两间存在显著的共线性，因此我们只取共线性变量的其中之一用于接下来的CAP分析。这里我们选择铝、钙和钾这三个变量与其它变量组成环境因子矩阵。
#select
env %>%
data.frame() %>%
dplyr::select(-Fe -Mg -S) ->env2
env2
CAP分析使用的是vegan包中的capscale函数。该函数最初是作为约束性近似分析的一个变体而开发的（Anderson & Willis 2003），但这些发展使它与dbRDA相似。然而，它抛弃了具有负特征值的虚构维度，排序和显著性检验区域只基于真实维度和正特征值。
与RDA和CCA分析不同的是CAP分析必须控制（约束）某些环境因子，评估未被约束的环境因子对物种矩阵的解释量。这里我们控制Baresoil和Humdepth两个变量。
vare_cap <- capscale(varespec ~. Condition(Humdepth Baresoil ) env2
dist="bray")
vare_cap
结果表明了约束部分解释了39.21%的物种变化，控制部分解释了25.09%的物种变化，未被解释的部分有41.39%。
此时我们能看到存在一个Imaginary的解释量，这意味着可能存在负引擎，为了避免这种情况我们可以把负值去除。代码如下：
#remove negative eigenvalues with additive constant
vare_cap2 <- capscale(varespec ~. Condition(Humdepth Baresoil ) env2
dist="bray" add =TRUE)
vare_cap2
此时的结果表明了约束部分解释了37.61%的物种变化，控制部分解释了20.49%的物种变化，未被解释的部分有41.91%。并且没有了Imaginary的解释量。
进一步查看每个轴的重要性，代码如下：
summary(vare_cap2)
结果表明CAP1解释了15.02%，CAP2解释了12.82%。
plot(vare_cap2)
1）模型整体的显著性检验
anova(vare_cap2)
结果表明我们用来做例子构建的模型并不显著
2）单一变量的显著性
anova(vare_cap2 by = "term")
结果表明这些变量中只有N和Al显著影响物种矩阵。
3）轴的显著性检验
anova(vare_cap2 by = "axis")
轴检验表明了所有的轴对物种矩阵的变化解释量都不显著。
4. 画图一般来说anova检验不显著就没必要继续画图了，应该考虑选择其它的模型，这里为了演示，我们对以上数据进行画图。
代码如下：
#-----------------plot---------------------
scrs<-scores(vare_cap2 display=c("sp" "wa" "lc" "bp" "cn"))
names(scrs)
#计算箭头变化倍数并计算箭头的长度
multiplier <- vegan:::ordiArrowMul(scrs$biplot)
multiplier
df_arrows<- scrs$biplot*multiplier#这一步是计算箭头的长度
df_arrows
colnames(df_arrows)
rownames(df_arrows)
#添加样方坐标
site<-scrs$sites %>%
as_tibble(rownames = "sample")
site
group <- env2 %>%
as_tibble(rownames = "sample") %>%
mutate(group = if_else(Humdepth>2 "deep" "superficial")) %>%
select(sample group)
site2 <- site %>%
left_join(group by = "sample")
df_arrows=as.data.frame(df_arrows)
ggplot(data=site2 aes(CAP1 CAP2))
geom_hline(aes(yintercept=0) colour="#d8d6d6" linetype=5)
geom_vline(aes(xintercept=0) colour="#d8d6d6" linetype=5)
geom_point(aes(color=group) size=3.5 shape=16)
geom_segment(data=df_arrows aes(x = 0 y = 0 xend = CAP1 yend = CAP2)
arrow = arrow(angle= 15 length = unit(0.25 "cm")) color="#000000" alpha=0.5) #画箭头
geom_text(data=as.data.frame(df_arrows*1.1) aes(CAP1 CAP2 label = rownames(df_arrows))
color="#000000" alpha=0.7) #箭头对应的信息
scale_color_manual(values = c("#2692d6" "#8cc1e3")) #点的颜色
theme_bw()
theme(panel.grid=element_blank())
labs(x = "CAP1 (15.02%)" y = "CAP2 (12.82%)")
