R语言_效应量effect size

本文主要介绍通过effect size来展示实验处理或控制因子对一些实验指标的影响。这种effect size，不管是实际差值（raw mean difference）或者其变型（ln RR）、还是标准化的差值（standardized mean difference），相比于传统的方差分析表格和柱状图等都更为清晰、明了，如下面的NM和GCB文章结果：
https://www.nature.com/articles/s41564-022-01147-3

Wu et al. Nat Microbiol

https://onlinelibrary.wiley.com/doi/10.1111/gcb.16651

Petro et al. GCB

示例数据

#   导入数据
> str(da)
'data.frame':   20 obs. of  34 variables:
 $ W       : chr  "C" "C" "C" "C" ...
 $ ANPP    : num  222 277 197 164 141 ...
 $ BNPP    : num  1043 776 866 1351 529 ...
 $ NPP     : num  1265 754 1063 1515 670 ...
 $ Richness: num  15 17 14 17 14 15 17 14 17 14 ...
 $ Grass   : num  44 41 40.7 47.9 44.5 ...
 $ others  : num  56 59 59.3 52.1 55.5 ...
> table(da$W)   # 处理组与对照组

   C W2.5 
  10   10

下图是原始数据的均值差异，红色圆点示 (Treatment - Control)>0，黑色反之。

Raw_difference.jpg

1. ln RR + "metafor"显著性检验

借用meta分析的效应值方法，展示处理效应对所测指标的影响

计算 log RR (yi)

> library(metafor)
> da7 = escalc(measure = "ROM", 
+              m1i = meanW, sd1i = sdW,  n1i = nW,
+              m2i = meanC, sd2i = sdC,  n2i = nC,  data = da6)
> da7$order = factor(rownames(da7), levels = rev(rownames(da7)))
> da8 = da7[order(da7$order),]
> head(da8)

             meanC     meanW         sdC        sdW nC nW      yi     vi    order 
NO3      17.967609  4.924877 11.08805328 2.32914353 10 10 -1.2943 0.0604      NO3 
NH4      30.818814 27.123295 10.51116292 8.28958950 10 10 -0.1277 0.0210      NH4 
pltR     17.200000 15.400000  1.27633261 1.08337258 10 10 -0.1105 0.0010     pltR 
ShannonT  6.824683  6.677340  0.07435578 0.09345605 10 10 -0.0218 0.0000 ShannonT 
CaC       3.013501  2.825361  0.27431139 0.83272419 10 10 -0.0645 0.0095      CaC 
SOC      25.691747 27.265316  2.13963451 3.25220321 10 10  0.0594 0.0021      SOC

提取置信区间、显著性

# rma()计算随机效应模型
> res8 = list()
> for (i in 1:33) {
+   res8[[i]] <-  rma(yi, vi, data = da8[i,])
+ }
> res9= data.frame(
+   da8,
+   ci.lb = list.stack(list.select(res8, ci.lb)),
+   ci.ub = list.stack(list.select(res8, ci.ub)),
+   pvale = list.stack(list.select(res8, pval))
+     )%>%
+   mutate(gp = ifelse(pval<0.05, "solid", "blank"))

作图

> res9$order = factor(res9$order, levels = rev(res9$order))
> ggplot(res9,aes(x=order, y=yi, color= gp)) +  
+   geom_point(size = 4, shape = 19)+
+   scale_color_manual(values = c( "grey","navyblue"))+
+   geom_errorbar(aes(x=order,ymin=ci.lb, ymax=ci.ub), width=.15, linewidth =0.5)+
+   geom_text(aes(x= order, y = 1.2,label = ifelse(pval<0.05,"*", "")), size = 10)+
+   geom_hline(aes(yintercept=0))+
+   xlab("")+ylab("log RR") +
+   theme_bw()+
+   theme(legend.position = "none")+
+   coord_flip()

# metafor内置的forest函数成图也附在最后
> res7 <- rma(yi, vi, data = da8)
> forest(res7, slab = rownames(da8), 
             cex=1, header=TRUE)

实际均值-lnRR-ggplot.jpg

实际均值-lnRR-metafor包.jpg

2. regression coefficients + LMM显著性

而在上述NM文章中，effect size采用LMM模型的回归系数，以线性混合模型结果作为处理效应的显著性。详见原文:
"Linear mixed-effects models (LMMs) for determining the sources of variations in hierarchical biological data were first employed to test the effects of treatments and their interactions on soil biogeochemistry and plant communities. In these models, the regression coefficients represent the directions and magnitudes of the treatment effects, namely effect sizes (β)."

首先，以Wu et al.文章数据展示

# data and code from Wu et al. 2022 NM article 
# https://github.com/Linwei-Wu/warming_soil_biodiversity

# data prepairation
treatused<-read.csv("treatment info.csv",header = T,row.names = 1) 
divindex<-read.csv("16S_alpha.csv",header = T,row.names = 1)  # 360 samples
treatused$year<-treatused$year-2009  # effect of years 
# check
divindex<-divindex[match(row.names(treatused),row.names(divindex)),]  
sum(row.names(divindex) == row.names(treatused))
divindex<-scale(divindex)
library(lme4)
library(car)
# calculate LMM models
divs1<-sapply(1:ncol(divindex),function(j){
  message("Now j=",j," in ",ncol(divindex),". ",date())
  if (length(unique(divindex[,j]))<3){
    result<-rep(NA,38)
  } else {
    div<-data.frame(divtest=divindex[,j],treatused)
    
    fm1<-lmer(divtest~warm*precip*clip+(1|year)+(1|block),data=div)
    
    presult<-car::Anova(fm1,type=2)
    coefs<-coef(summary(fm1))[ , "Estimate"]  ##four coefs
    names(coefs)<-paste0(names(coefs),".mean")
    
    SEvalues<-coef(summary(fm1))[ , "Std. Error"] ##standard errors
    names(SEvalues)<-paste0(names(SEvalues),".se")
    
    tvalues<-coef(summary(fm1))[ , "t value"] ##t values
    names(tvalues)<-paste0(names(tvalues),".t")
    
    chisqP<-c(presult[,1],presult[,3])
    names(chisqP)<-c(paste0(row.names(presult),".chisq"),paste0(row.names(presult),".P"))
    
    
    
    result<-c(coefs,tvalues,SEvalues,chisqP)}
  result
})
colnames(divs1)<-colnames(divindex)
#  上述代码来自Wu et al.文章

# select the microbial "speceis richness" and "PD" indices
library(reshape2)
divs2 <- divs1[c(2:4, 18:20, 32:34), c(1,6)] %>% t() %>% as.data.frame() 
divs3 <- divs1 %>% as.data.frame() %>% rownames_to_column("rown") %>% 
  separate(col = rown, into = c("eff", "type"), sep = "\\." )

dcof <- divs3[2:8,]; dcof$type = NULL
dcof1 <- melt(dcof, value.name = "mean") 
dse <- divs3[18:24,]; dse$type = NULL
dse1 <- melt(dse, value.name = "se") 
dp <- divs3[32:38,]; dp$type = NULL
dp1 <- melt(dp, value.name = "p")

dsum <- data.frame(dcof1, se = dse1[,3], p = dp1[,3])
dsum1 <- dsum[which(dsum$variable %in% c("richness", "PD")),]
dsum2 <- dsum1[c(1:3,8:10),]
# plotting 
ggplot(dsum2,aes(x = eff, y = mean)) +  
  geom_bar(stat='identity')+
  geom_errorbar(aes(x= eff,ymin=mean-se, ymax=mean+se), width=.15, linewidth =0.5)+
  geom_hline(aes(yintercept=0))+
  theme_bw()+
  ylab("Effect size")+
  facet_grid(.~variable)

这里仅展示细菌多样性的结果

Rplot66.png

下面以为本例中的数据作图，展示以LMM模型的coefficient为效应量，se作为误差棒的effect size图，代码跟上面类似，就不展示了（这里的数据没有标准化，所以某些值就比较大）。最终结果与方法1结果还是存在差异：显著性结果基本一致（除了Richness），但是effect的方向还是与raw-difference存在出入(NO3, Richness等)。这说明实验处理效应的方向（raw-difference）与线性混合模型的处理因子系数（coefficient，方法2）、线性混合模型置信区间均值（cl/2+cu/2，方法3）并非是一致的。

Rplot68.png

3. 置信区间均值+LMM显著性

根据LMM模型结果，提取处理因子效应的置信区间，计算上限、下限的均值，作为effect size，以线性混合模型结果作为处理效应的显著性。其结果如下：

> laa <- colnames(da)[c(2:32, 34:35)]
> mod = list()
> pval = vector()
> res = list()
> library(lme4)
> for (i in 1:33) {
+   mod[[i]] <-lmer(get(laa[i]) ~ W + (1|block), data = da)
+   pval[i] <- car::Anova(mod[[i]],type=2)$`Pr(>Chisq)`
+   res[[i]] = confint(mod[[i]])[4,] 
+ }

> datapl <- as.data.frame(do.call(rbind, res)) %>% mutate(var = laa) %>%
+   rename(cl = "2.5 %", cu = "97.5 %") %>%
+   mutate(mean = (cl+cu)/2,
+          cl1 = ifelse(cl > 0, (cl-min(cl))/(max(cl)-min(cl)), (cl-max(cl))/(max(cl)-min(cl))),
+          cu1 = ifelse(cu > 0, (cu-min(cu))/(max(cu)-min(cu)), (cu-max(cu))/(max(cu)-min(cu))),
+          mean1 = (cl1+cu1)/2,
+          pval = cbind(pval)[,1],
+          shape = ifelse(pval<0.05, "solid", "blank"))

> ggplot(datapl,aes(x=var, y=mean1, color = shape)) +  
+   geom_point(size=3.5)+
+   scale_color_manual(values = c("grey","navyblue"))+
+   geom_errorbar(aes(x=var,ymin=cl1, ymax=cu1), width=.15, linewidth =0.5)+
+   geom_text(aes(x= var, y = 1.2,label = ifelse(pval<0.05,"*", "")), size = 10)+
+   scale_shape_manual(values = ifelse(datapl$pval<0.05, 19, 1))+
+   geom_hline(aes(yintercept=0))+
+   xlab("")+ylab("Effect size") +
+   theme_bw()+
+   theme(legend.position = "none")+
+   coord_flip()

置信区间均值-daCW2.5-Var-W.jpg

该结果与前面方法1 结果存在出入。后面再按照方法2，看看以回归系数为effect sizes的结果与实际差异（raw difference或 ln RR）的出入情况。

（更新：2023.04.02）

填坑完毕！个人感觉，显著性方面可以采用metafor包自带的混合模型计算，也可以采用线性混合模型自己提取P值。但是在effect size的选择问题上，raw-difference及其对数变形可能更符合原始的处理效应。如有谬误，还请指正。大家有兴趣的话也可以用自己的数据试试，欢迎一起讨论。
（更新：2023.06.13）