envfit结果如何创建？

Posted 2023-01-27

我对如何从envfit()程序包中的vegan函数重新创建结果有疑问。

这里是envfit()与排序和环境矢量一起使用的示例。

data(varespec)
data(varechem)
ord <- metaMDS(varespec)
chem.envfit <- envfit(ord, varechem, choices = c(1,2), permutations = 999)
chem.scores.envfit <- as.data.frame(scores(chem.envfit, display = "vectors"))
chem.scores.envfit

“表中显示的值是线性回归的标准化系数，该线性回归用于将向量投影到排序中。这些是单位长度箭头的方向。”-Plotted envfit vectors not matching NMDS scores

中的注释

也来自?envfit：

连续变量（向量）的打印输出给出了方向余弦，是单元头的坐标长度向量。在图中，它们通过它们的相关性（平方r2）的根，因此弱预测变量的箭头较短比强大的预测指标您可以使用以下方式查看缩放后的相对长度命令分数。

有人可以明确地告诉我正在运行什么线性模型，正在使用什么标准化系数，以及在哪里使用余弦来创建这些值？

答案

我可能不应该在该答案中说“标准化”。

对于varechem中的每一列（变量）和排序的前两个轴（choices = 1:2），线性模型为：

\hat(env_j) = \beta_1 * scr1 + \beta_2 * scr2

其中env_j是varechem中的$ j $ th变量，scr1和scr2是所考虑的第一和第二条轴上的轴分数（即，由choices = 1:2定义的平面，但是到更高的维度），而\beta是这对轴分数的回归系数。

[我们在模型中没有权重，因为我们（加权）将varechem和轴分数中的所有变量居中，权重实际上仅涉及CCA，capscale()和DCA方法，因为它们本身就是加权模型。

[由轴分数跨越的空间中箭头的头部是该模型的系数-我们实际上进行了归一化（在其他答复中我误称为“标准化”），以使箭头具有单位长度。这些值（NMDS1输出中的NMDS2和envfit列）在https://en.wikipedia.org/wiki/Direction_cosine的意义上是直接余弦。

这里是不涉及权重并且env中所有变量都是数字的情况下的简化操作，如您的示例。 （请注意，出于效率原因，我们实际上并不这样做：如果您确实需要详细信息，请参阅vectorfit()后面的代码以了解QR分解。）

## extract the axis scores for the axes we want, 1 and 2
scrs <- scores(ord, choices = c(1,2))

## centre the scores (note not standardising them)
scrs <- as.data.frame(scale(scrs, scale = FALSE, center = TRUE))

## centre the environmental variables - keep as matrix
env <- scale(varechem, scale = FALSE, center = TRUE)

## fit the linear models with no intercept
mod <- lm(env ~ NMDS1 + NMDS2 - 1, data = scrs)

## extract the coefficients from the models
betas <- coef(mod)

## normalize coefs to unit length
##   i.e. betas for a  particular env var have sum of squares = 1
t(sweep(betas, 2L, sqrt(colSums(betas^2)), "/"))
最后一行给出：

> t(sweep(betas, 2L, sqrt(colSums(betas^2)), "/"))
               NMDS1      NMDS2
N        -0.05731557 -0.9983561
P         0.61972792  0.7848167
K         0.76646744  0.6422832
Ca        0.68520442  0.7283508
Mg        0.63252973  0.7745361
S         0.19139498  0.9815131
Al       -0.87159427  0.4902279
Fe       -0.93600826  0.3519780
Mn        0.79870870 -0.6017179
Zn        0.61755690  0.7865262
Mo       -0.90308490  0.4294621
Baresoil  0.92487118 -0.3802806
Humdepth  0.93282052 -0.3603413
pH       -0.64797447  0.7616621
在这种情况下，它复制（除了显示更多的符号数字）envfit()返回的值：

> chem.envfit

***VECTORS

            NMDS1    NMDS2     r2 Pr(>r)    
N        -0.05732 -0.99836 0.2536  0.045 *  
P         0.61973  0.78482 0.1938  0.099 .  
K         0.76647  0.64228 0.1809  0.095 .  
Ca        0.68520  0.72835 0.4119  0.006 ** 
Mg        0.63253  0.77454 0.4270  0.003 ** 
S         0.19139  0.98151 0.1752  0.109    
Al       -0.87159  0.49023 0.5269  0.002 ** 
Fe       -0.93601  0.35198 0.4450  0.002 ** 
Mn        0.79871 -0.60172 0.5231  0.002 ** 
Zn        0.61756  0.78653 0.1879  0.100 .  
Mo       -0.90308  0.42946 0.0609  0.545    
Baresoil  0.92487 -0.38028 0.2508  0.061 .  
Humdepth  0.93282 -0.36034 0.5201  0.001 ***
pH       -0.64797  0.76166 0.2308  0.067 .  
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Permutation: free
Number of permutations: 999

另一答案

根据as.mlm.cca here的文档，vegan使用其自己的算法cca。