之前的推文分享过一篇
跟着Nature Ecology&Evolution学数据分析：R语言做方差分解的一个简单小例子
内容是重复了一篇 Nature Ecology&Evolution期刊论文的方差分解过程，虽然对运行过程还是不太理解，但是能得到结果。今天的推文我们利用得到的结果复现一下论文中的Figure4a 的右半部分，左侧的堆积柱形图以及如何把两个图组合到一起争取再单独出一期推文介绍（其实是最近有点创作乏力，一篇推文的内容水成两篇 哈哈哈哈）
image.png
首先是运行之前推文的代码得到画图数据
代码语言：javascript复制datatotal<-read.table("datasetmultifunctionality.txt", header=T, sep="t")
colnames(datatotal)
#####logtransformation moments
datatotal[,c(12,13,16,17)]<-log(datatotal[,c(12,13,16,17)])
datatotal[,14]<-log(datatotal[,14]-min(datatotal[,14]) 1)
datatotal[,15]<-log(datatotal[,15]-min(datatotal[,15]) 1)
datatotal[,18]<-log(datatotal[,18]-min(datatotal[,18]) 1)
datatotal[,19]<-log(datatotal[,19]-min(datatotal[,19]) 1)
#####Zscorring environmental variables
datatotal$ELEVATION<-(datatotal$ELEVATION-mean(datatotal$ELEVATION))/sd(datatotal$ELEVATION)
datatotal$LAT<-(datatotal$LAT-mean(datatotal$LAT))/sd(datatotal$LAT)
datatotal$SINLONG<-(datatotal$SINLONG-mean(datatotal$SINLONG))/sd(datatotal$SINLONG)
datatotal$COSLONG<-(datatotal$COSLONG-mean(datatotal$COSLONG))/sd(datatotal$COSLONG)
datatotal$SLO<-(datatotal$SLO-mean(datatotal$SLO))/sd(datatotal$SLO)
datatotal$ARIDITY<-(datatotal$ARIDITY-mean(datatotal$ARIDITY))/sd(datatotal$ARIDITY)
datatotal$SAND<-(datatotal$SAND-mean(datatotal$SAND))/sd(datatotal$SAND)
datatotal$PH<-(datatotal$PH-mean(datatotal$PH))/sd(datatotal$PH)
datatotal$SR<-(datatotal$SR-mean(datatotal$SR))/sd(datatotal$SR)
#####Zscorring moments
datatotal$CWM_logH<-(datatotal$CWM_logH-mean(datatotal$CWM_logH))/sd(datatotal$CWM_logH)
datatotal$CWV_logH<-(datatotal$CWV_logH-mean(datatotal$CWV_logH))/sd(datatotal$CWV_logH)
datatotal$CWS_logH<-(datatotal$CWS_logH-mean(datatotal$CWS_logH))/sd(datatotal$CWS_logH)
datatotal$CWK_logH<-(datatotal$CWK_logH-mean(datatotal$CWK_logH))/sd(datatotal$CWK_logH)
datatotal$CWM_logSLA<-(datatotal$CWM_logSLA-mean(datatotal$CWM_logSLA))/sd(datatotal$CWM_logSLA)
datatotal$CWV_logSLA<-(datatotal$CWV_logSLA-mean(datatotal$CWV_logSLA))/sd(datatotal$CWV_logSLA)
datatotal$CWS_logSLA<-(datatotal$CWS_logSLA-mean(datatotal$CWS_logSLA))/sd(datatotal$CWS_logSLA)
datatotal$CWK_logSLA<-(datatotal$CWK_logSLA-mean(datatotal$CWK_logSLA))/sd(datatotal$CWK_logSLA)
#####Zscorring ecosystem functions
datatotal$BGL<-(datatotal$BGL-mean(datatotal$BGL))/sd(datatotal$BGL)
datatotal$FOS<-(datatotal$FOS-mean(datatotal$FOS))/sd(datatotal$FOS)
datatotal$AMP<-(datatotal$AMP-mean(datatotal$AMP))/sd(datatotal$AMP)
datatotal$NTR<-(datatotal$NTR-mean(datatotal$NTR))/sd(datatotal$NTR)
datatotal$I.NDVI<-(datatotal$I.NDVI-mean(datatotal$I.NDVI))/sd(datatotal$I.NDVI)
#####Calculating indices of multifunctionality (M5: 5 functions)
colnames(datatotal)
M5<-rowMeans(datatotal[,c(20,21,22,23,24)])
datatotal<-cbind(datatotal,M5)
#####Log-transfromation of multifunctionality
logM5<-log(datatotal$M5-min(datatotal$M5) 1)
datatotal<-cbind(datatotal,logM5)
library(MuMIn)
mod12<-lm(logM5 ~ LAT SINLONG COSLONG
ARIDITY SLO SAND PH I(PH^2) ELEVATION
CWM_logSLA I(CWM_logSLA^2) CWV_logSLA I(CWV_logSLA^2) CWS_logSLA CWK_logSLA I(CWK_logSLA^2)
CWM_logH I(CWM_logH^2) CWV_logH I(CWV_logH^2) CWS_logH CWK_logH I(CWK_logH^2)
SR
, data=datatotal)
# 这一步要好长时间
dd12<-dredge(mod12, subset = ~ LAT & SINLONG & COSLONG & ARIDITY & SLO & SAND & PH &SR & ELEVATION &
dc(CWM_logSLA,I(CWM_logSLA^2)) & dc(CWV_logSLA,I(CWV_logSLA^2)) & dc(CWK_logSLA,I(CWK_logSLA^2))
& dc(CWM_logH,I(CWM_logH^2)) & dc(CWV_logH,I(CWV_logH^2)) & dc(CWK_logH,I(CWK_logH^2)),
options(na.action = "na.fail"))
subset(dd12,delta<2)
de12<-model.avg(dd12, subset = delta < 2)
summary(de12)
这部分代码和示例数据可以在公众号后台回复20210403获取
接下来是画图
首先是获取画图数据
代码语言：javascript复制load("de12.Rdata")
library(tidyverse)
as.data.frame(de12.summary$coefmat.subset) %>%
mutate(var=rownames(.),
group=sample(LETTERS[1:4],22,replace = T)) %>%
rows_delete(tibble(var="(Intercept)")) -> mydf
这里的分组文件我先随便构造了
添加显著性星号并设置因子水平
代码语言：javascript复制mydf %>%
arrange(group,Estimate) %>%
mutate(signi=case_when(
`Pr(>|z|)` > 0.05 ~ '',
`Pr(>|z|)` < 0.05 & `Pr(>|z|)` >= 0.01 ~ '*',
`Pr(>|z|)` < 0.01 & `Pr(>|z|)` >= 0.001 ~ '**',
`Pr(>|z|)` < 0.001 ~ '***'
)) %>%
mutate(var=fct_relevel(var,var)) -> mydf1
最后是画图代码
代码语言：javascript复制library(ggplot2)
library(ggh4x)
library(see)
ggplot(mydf1,aes(x=Estimate,y=var))
geom_point(aes(color=group),
show.legend = F,
size=5)
xlim(-0.2,0.2)
labs(y=NULL)
theme_minimal()
theme(panel.grid = element_blank(),
axis.text.y = element_blank(),
axis.line.x = element_line(),
axis.ticks.x = element_line())
guides(x=guide_axis_truncated(trunc_lower = -0.2,
trunc_upper = 0.2))
geom_linerange(aes(xmin=Estimate-`Std. Error`,
xmax=Estimate `Std. Error`,
color=group),
show.legend = F)
geom_text(aes(y=1:21,x=0.1,label=var),hjust=0)
geom_text(aes(y=1:21,x=0.19,label=signi))
scale_color_material_d()
最终结果如下
image.png
当然还有一些细节有待调整 我们下期推文一并介绍
画图的示例数据和代码可以在公众号后台回复20210910获取
最后祝关注公众号的老师们 教师节快乐！科研顺利！文章一个劲儿的发！
欢迎大家关注我的公众号
小明的数据分析笔记本
小明的数据分析笔记本 公众号 主要分享：1、R语言和python做数据分析和数据可视化的简单小例子；2、园艺植物相关转录组学、基因组学、群体遗传学文献阅读笔记；3、生物信息学入门学习资料及自己的学习笔记！
