sars：拟合SAR模型的最新工具

Listenlii-生物信息知识分享

发布于 2020-05-29 11:34:38

1.1K0

发布于 2020-05-29 11:34:38

文章被收录于专栏：Listenlii的生物信息笔记

之前介绍过拟合种面积关系(species–arearelationship, SAR)工具：

R——mmSAR对种面积关系进行拟合

今年3月份又出现了一个更强大的工具：sars

近期研究表明只使用单一的模型不能很好地拟合所有SAR数据，多个模型叠加可能更有实际意义。因此作者开发了sars：

可以使用线性或非线性的回归拟合20个不同的模型(mmSAR只有8个模型)；

还可以计算多个模型的平均曲线；

可用bootstrapping的方法得到置信区间。

SAR研究中使用最广泛的是幂律模型(power model)。但是一些研究已经发现大尺度上的SAR符合的是S型曲线(反曲型)。针对SAR模型不统一的情况，目前有两种策略，一是多个模型进行拟合，根据一定的标准选出效果最优(如AIC最小)的模型；二是多个模型拟合，取平均曲线。但是目前没有R包能实现。之前的两个包：

BAT可拟合三种SAR模型：线性、幂律和对数模型。

mmSAR可拟合8种模型，但是相比于目前超过20种的模型也不够用。

Sars相比于mmSAR的优势在于：

绘图更友好。可绘制加权的多模型曲线；

确定模型的形状(如线性、上凸、下凸、S型)；

是否有渐近线；

利用物种-地点丰度矩阵对Coleman’s (1981)的随机布局模型进行拟合；

建立了岛屿生物地理学的一般动态模型(general dynamic model, GDM)

20个模型

>install.packages("sars")
>library(sars)

模型的拟合

非约束的Nelder–Mead最优化算法(R中的optim)最小化残差平方和来估计模型的参数。

##单个模型
>fit <- sar_loga(data = galap) 
>summary(fit) 
Model:
Logarithmic

Call: 
S == c + z * log(A)

Did the model converge: TRUE

Residuals:
      0%      25%      50%      75%     100% 
-164.000  -25.675   11.350   38.725  115.700 

Parameters:
     Estimate  Std. Error     t value    Pr(>|t|)        2.5%  97.5%
c   0.2915341  34.5985767   0.0084262   0.9933958 -68.9056193 69.489
z  30.2802630   8.3203931   3.6392827   0.0026812  13.6394768 46.921

R-squared: 0.49, Adjusted R-squared: 0.41
AIC: 141.78, AICc: 143.78, BIC: 144.1
Observed shape: convex up, Asymptote: FALSE


Warning: The fitted values of the model contain negative values (i.e. negative species richness values)

>plot(fit)
#多个模型
>fitC <- sar_multi(data = galap, obj = c("power", "loga", "monod"))
Now attempting to fit the 3 SAR models: 

--  multi_sars -------------------------------------------------- multi-model SAR --
> power : √
> loga  : √
> monod : √


>plot(fitC)
模型的一般性和同质性检验

>data(galap) 
>fit <- sar_linear(data = galap, normaTest ="lillie", homoTest = "cor.fitted") 
>summary(fit)
Model:
Linear model

Call: 
S == c + m*A

Did the model converge: TRUE

Residuals:
     0%     25%     50%     75%    100% 
-107.90  -52.10  -24.15   31.75  218.00 

Parameters:
   Estimate Std. Error   t value  Pr(>|t|)     2.5 %   97.5 %
c 70.314003  25.743130  2.731370  0.016228 15.100481 125.5275
m  0.185382   0.072884  2.543504  0.023410  0.029060   0.3417

R-squared: 0.32, Adjusted R-squared: 0.21
AIC: 146.36, AICc: 148.36, BIC: 148.68
Observed shape: linear, Asymptote: FALSE


Warning: The normality test selected indicated the model residuals are not normally distributed (i.e. P < 0.05)

>fit$normaTest
$test
[1] "lillie"

[[2]]

    Lilliefors (Kolmogorov-Smirnov) normality test

data:  res
D = 0.2146, p-value = 0.04725

多模型曲线

将每一个成功拟合模型的预测丰度值乘以模型的权重(AIC,AICC,BIC等)，然后对所有模型的结果值求和，单个模型的线性组合构建多模型平均曲线。

>data(niering) 
>fit <- sar_average(data= niering, obj =c("power","loga","koba","mmf","monod",
                                         "negexpo","chapman","weibull3","asymp"),
                   normaTest = "none", homoTest = "none", neg_check = FALSE, 
                   confInt = TRUE, ciN= 50) 

 Now attempting to fit the 9 SAR models: 

--  multi_sars -------------------------------------------------- multi-model SAR --
> power    : √
> loga     : √
> koba     : √
> mmf      : √
> monod    : √
> negexpo  : √
> chapman  : Warning: could not compute parameters statistics
> weibull3 : √
> asymp    : x

1 models could not be fitted and have been excluded  from the multi SAR


All models passed the model  validation checks

8 remaining models used to construct the multi  SAR:
 Power, Logarithmic, Kobayashi, MMF, Monod, Negative exponential, Chapman Richards, Cumulative Weibull 3 par. 
------------------------------------------------------------------------------------

Calculating sar_multi confidence intervals - this may take  some time:
  |==========================================================================| 100%

>par(mfrow = c(3,1))
#多个模型画图
>plot(fit, ModTitle = "a) Multimodel SAR")
#多模型平均曲线
>plot(fit, allCurves = FALSE, ModTitle =
       "c) Multimodel SAR with confidence intervals", confInt = TRUE)
#每个模型所占权重
>plot(fit, type = "bar", ModTitle = "b) Model weights", cex.lab = 1.3)