赖江山的博客分享 http://blog.sciencenet.cn/u/laijiangshan 生态、统计与R语言

博文

phylolm.hp包发布及案例代码

已有 366 次阅读 2024-12-7 17:13 |个人分类:rdacca.hp|系统分类:科研笔记

11月13日,我们发布新包phylolm.hp:phylolm.hp包在CRAN上发布,用于量化谱系、环境因子对性状变异的影响. 此包基于“平均共享方差理论”,研究者可以对phylolm包中phylolm()和phyloglm()输出的模型进行 R² 分解,不仅能获得环境因子、谱系对性状变异的单个R²(individual R²)(各个individual R²总和刚好等于全模型的总 R²);还能够清晰展现环境因子与谱系之间的共享R², 从而量化研究谱系和环境因子对性状的影响。当环境变量少于 3 个时,phylolm.hp 还可生成更直观的变差分解韦恩图。

包发布之后,受到广泛关注,公号内容已经累计浏览2300多次。有学者使用之后给我提了要求,是否能让预测变量分组(比如气候组、土壤组、性状组),并且评估预测变量组之间、谱系与预测变量之间的共同解释比例,以及独立解释比例和单个 R²。我重新升级一下包,以实现对预测变量进行分组再进行分析。目前 phylolm.hp第3版现已上传至CRAN平台,我们诚邀各位用户从CRAN更新包,并尝试新的分组分解功能。若有任何问题或建议,欢迎通过电子邮件与我联系:lai@njfu.edu.cn

目前phyloglm.hp()函数可以有3个输入参数:mod是phylolm中包中phyloglm()和phylolm()函数输出的模型对象、以及一个逻辑参数,用于控制是否输出commonality分析结果。另外加了一个iv参数,用一个列表来对模型中预测变量进行变量组合,变量名务必是在mod中预测变量的名称,建议自己取好组的名字。

以下是用新模拟数据形成的案例,期望大家用真实数据进行测试:

library(phylolm.hp)

set.seed(231)

tre <- rcoal(60)

taxa <- sort(tre$tip.label) 

b0 <- 0      

b1 <- 0.3    

b2 <- 0.5 

b3 <- 0.4

x <- rTrait(n=1, phy=tre, model="lambda", parameters=list(ancestral.state=0, sigma2=15, lambda=0.9))          

x2 <- rTrait(n=1, phy=tre, model="lambda",  

parameters=list(ancestral.state=0, sigma2=10, lambda=0.9))  

x3 <- rTrait(n=1, phy=tre, model="lambda",  

parameters=list(ancestral.state=0, sigma2=13, lambda=0.9)) 

y <- b0 + b1 * x + b2 * x2 + b3*x3+ rTrait(n=1, phy=tre, model="lambda",

parameters=list(ancestral.state=0, sigma2=5, lambda=0.9))            

dat <- data.frame(trait=y[taxa], pred=x[taxa], pred2=x2[taxa],pred3=x3[taxa])

fit <- phylolm(trait ~ pred + pred2 + pred3, data=dat, phy=tre, model="lambda")

phyloglm.hp(fit,commonality=TRUE)

plot(phyloglm.hp(fit,commonality=TRUE),commonality=TRUE)

iv=list(env1="pred",env2=c("pred2","pred3"))

phyloglm.hp(fit,iv)

plot(phyloglm.hp(fit,iv))

plot(phyloglm.hp(fit,iv,commonality=TRUE),commonality=TRUE)

set.seed(123456)

tre <- rtree(50)

x1 <- rTrait(n=1, phy=tre)  

x2 <- rTrait(n=1, phy=tre)

x3 <- rTrait(n=1, phy=tre)

X <- cbind(rep(1, 50), x1, x2, x3)

y <- rbinTrait(n=1, phy=tre, beta=c(-1, 0.9, 0.9, 0.5), alpha=1, X=X)

dat <- data.frame(trait01=y, predictor1=x1, predictor2=x2, predictor3=x3)

fit <- phyloglm(trait01 ~ predictor1 + predictor2 + predictor3, phy=tre, data=dat)

phyloglm.hp(fit)

phyloglm.hp(fit,commonality=TRUE)

plot(phyloglm.hp(fit,commonality=TRUE),commonality=TRUE)

iv=list(env1="predictor1",env2=c("predictor2","predictor3"))

phyloglm.hp(fit,iv)

plot(phyloglm.hp(fit,iv))

plot(phyloglm.hp(fit,iv,commonality=TRUE),commonality=TRUE)



https://blog.sciencenet.cn/blog-267448-1463334.html

上一篇:MEE的文章google学术引用超过600次,mark一下!
下一篇:单因素方差分析的多重比较如何获取每次比较的effect sizes和confidence intervals
收藏 IP: 218.2.103.*| 热度|

1 杨正瓴

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2024-12-22 16:50

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部