写在前面最近觉得有些疲惫，可能是由于长时间没有产出的原因。反省了一下，今年得出几篇看得过去的成果，要不然无颜面对江东父老呀。前几天师兄带我们集体交流学习，提到了使用效应值的形式展示处理对关键指标的影响。主要看了王岭老师团队发表的那篇NC和一篇GCB。总体而言，相对于以前单纯的展示单因素的结果，使用效应值能够更加定性的说明处理对指标的影响，更利于我们表述结果和讨论原因。今天就一起看看怎么在R里计算效应值吧。 
效应值的计算公式效应值的计算通常是使用ln RR来计算，其计算公式如下，小编直接截图GCB中的计算公式，先看原文的可以根据文末的参考文献查找 
数据准备练习数据，小编又编了点儿，值得注意的是，如果我们需要使用自己的实测数据求效应值，最少得9个重复，因为我们表格中的数据分别是处理的均值、标准差、和重复数以及对照的均值、标准差、重复数。文末下载~ 
开始实践我们先读取数据 library(metafor)library(ggplot2)library(glmulti)d1 <- read.csv(file.choose(),header = T)variable.names
读进来的数据是这样的 > d1 TR TM TSD TN CM SCD CN1 A 202.07733 42.130348 3 255.5513 51.77846 32 A 139.57867 39.573730 3 194.7520 18.52936 33 A 160.88000 15.805961 3 238.2107 37.50586 34 B 146.81867 84.565882 3 255.5513 51.77846 35 B 248.60800 71.051697 3 194.7520 18.52936 36 B 89.57867 26.335072 3 238.2107 37.50586 37 C 95.24133 34.520747 3 255.5513 51.77846 38 C 161.84547 3.379314 3 194.7520 18.52936 39 C 110.33467 21.170659 3 238.2107 37.50586 310 D 127.23867 32.082678 3 255.5513 51.77846 311 D 179.39733 14.401514 3 194.7520 18.52936 312 D 160.92533 19.359273 3 238.2107 37.50586 313 E 157.14400 38.824887 3 255.5513 51.77846 314 E 110.25467 1.636035 3 194.7520 18.52936 315 E 141.53333 43.818992 3 238.2107 37.50586 3
然后我们再计算yi和vi d2<-escalc(measure="ROM",data=d1,m1i=TM,sd1i=TSD,n1i=TN,m2i=CM,sd2i=SCD,n2i=CN)head(d2)
计算结果如下，我们的效应值计算就是根据yi和vi算的 > head(d2)
TR TM TSD TN CM SCD CN yi vi 1 A 202.07733 42.13035 3 255.5513 51.77846 3 -0.2348 0.0282 2 A 139.57867 39.57373 3 194.7520 18.52936 3 -0.3331 0.0298 3 A 160.88000 15.80596 3 238.2107 37.50586 3 -0.3925 0.0115 4 B 146.81867 84.56588 3 255.5513 51.77846 3 -0.5542 0.1243 5 B 248.60800 71.05170 3 194.7520 18.52936 3 0.2442 0.0302 6 B 89.57867 26.33507 3 238.2107 37.50586 3 -0.9780 0.0371
然后就是，先整体计算所有处理的效应值 #全组计算r1<-rma(yi,vi, data=d2, method="REML")summary(r1)
计算结果如下，我们需要的绘图结果就是第18行和第19行 > summary(r1)
Random-Effects Model (k = 15; tau^2 estimator: REML)
logLik deviance AIC BIC AICc -4.1848 8.3697 12.3697 13.6478 13.4606
tau^2 (estimated amount of total heterogeneity): 0.0759 (SE = 0.0379)tau (square root of estimated tau^2 value): 0.2755I^2 (total heterogeneity / total variability): 85.37%H^2 (total variability / sampling variability): 6.84
Test for Heterogeneity:Q(df = 14) = 81.5383, p-val < .0001
Model Results:
estimate se zval pval ci.lb ci.ub -0.4385 0.0824 -5.3207 <.0001 -0.6000 -0.2770 ***
---Signif. codes: 0 '***’ 0.001 '**’ 0.01 '*’ 0.05 '.’ 0.1 ' ’ 1
随后是基于每个处理的亚组计算，计算结果同上，就不再赘述 #亚组计算r2<-rma(yi,vi, data=subset(d2,TR=="A"), method="REML")r3<-rma(yi,vi, data=subset(d2,TR=="B"), method="REML")r4<-rma(yi,vi, data=subset(d2,TR=="C"), method="REML")r5<-rma(yi,vi, data=subset(d2,TR=="D"), method="REML")r6<-rma(yi,vi, data=subset(d2,TR=="E"), method="REML")#查看结果summary(r2)summary(r3)summary(r4)summary(r5)summary(r6)
开始实践随后们就可以根据计算结果和显著性绘制一个森林图了，绘制前根据结果吧数据整理成如下形式 
然后我们开始绘图吧，颜色和显著性懒得改，显著性在ppt里自己加就好，颜色可以使用代码修改一下。颜色和显著性的代码可以参考这篇推文。这里边有详细的修改细节。 #森林图的绘制#调用R包library(ggthemes)library(ggplot2)#读取数据dataset <- read.csv(file.choose(),header = T)#森林图绘制a<-ggplot(dataset, aes(Estimate, index))+ geom_point(size=4,color = "black")+ geom_errorbarh(aes(xmax = Estimate +`SE`, xmin = Estimate -`SE`),size= 0.5,height = 0.2, colour = "black") + scale_x_continuous(limits= c(-1, 1))+ geom_vline(aes(xintercept = 0),color="gray",linetype="dashed", size = 0.6) + xlab('效应量 Effect size ')+ ylab(' ')+ theme_few()+ theme(axis.text.x = element_text(size = 24, color = "black"))+ theme(axis.text.y = element_text(size = 24, color = "black"))+ theme(title=element_text(size=24))+ theme_bw()+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")) )#查看图aa#保存至pptlibrary(eoffice)topptx(filename="森林图.pptx")
完整版代码如下library(metafor)library(ggplot2)library(glmulti)d1 <- read.csv(file.choose(),header = T)variable.namesd2<-escalc(measure="ROM",data=d1,m1i=TM,sd1i=TSD,n1i=TN,m2i=CM,sd2i=SCD,n2i=CN)head(d2)#全组计算r1<-rma(yi,vi, data=d2, method="REML")summary(r1)#亚组计算r2<-rma(yi,vi, data=subset(d2,TR=="A"), method="REML")r3<-rma(yi,vi, data=subset(d2,TR=="B"), method="REML")r4<-rma(yi,vi, data=subset(d2,TR=="C"), method="REML")r5<-rma(yi,vi, data=subset(d2,TR=="D"), method="REML")r6<-rma(yi,vi, data=subset(d2,TR=="E"), method="REML")#查看结果summary(r2)summary(r3)summary(r4)summary(r5)summary(r6)
#森林图的绘制#调用R包library(ggthemes)library(ggplot2)#读取数据dataset <- read.csv(file.choose(),header = T)#森林图绘制a<-ggplot(dataset, aes(Estimate, index))+ geom_point(size=4,color = "black")+ geom_errorbarh(aes(xmax = Estimate +`SE`, xmin = Estimate -`SE`),size= 0.5,height = 0.2, colour = "black") + scale_x_continuous(limits= c(-1, 1))+ geom_vline(aes(xintercept = 0),color="gray",linetype="dashed", size = 0.6) + xlab('效应量 Effect size ')+ ylab(' ')+ theme_few()+ theme(axis.text.x = element_text(size = 24, color = "black"))+ theme(axis.text.y = element_text(size = 24, color = "black"))+ theme(title=element_text(size=24))+ theme_bw()+ theme(axis.ticks.length=unit(-0.25, "cm"), axis.text.x = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")), axis.text.y = element_text(margin=unit(c(0.5,0.5,0.5,0.5), "cm")) )#查看图aa#保存至pptlibrary(eoffice)topptx(filename="森林图.pptx")
参考文献： [1] Zhang, M., Delgado-Baquerizo, M., Li, G. et al. Experimental impacts of grazing on grassland biodiversity and function are explained by aridity. Nat Commun 14, 5040 (2023). https:///10.1038/s41467-023-40809-6 [2] Shi, Y., Wang, J., Ao, Y., Han, J., Guo, Z., Liu, X., Zhang, J., Mu, C., & Le Roux, X. (2021). Responses of soil N2O emissions and their abiotic and biotic drivers to altered rainfall regimes and co-occurring wet N deposition in a semi-arid grassland. Global Change Biology, 27, 4894–4908. https:///10.1111/gcb.15792 
|
