从元分析的混合效应多层模型中获取R平方

8
我正在使用R进行特定治疗对森林的元分析。为此模型,我需要拟合随机效应以解释方法和地点年龄差异,因为这两个是混淆变量,而我不想明确调查它们引起的变异。但据我所知,包[metfor]不允许在多层模型中计算类似于R平方的统计量。无论如何,为了更清楚地描述我的问题,这里有一个模拟数据集。
Log<-data.frame(Method=rep(c("RIL","Conv"),each=10),
     RU=runif(n=20,min=10,max=50),SDU=runif(n=20,5,20),
     NU=round(runif(n=20,10,20),0))
Log$Study<-rep(1:4,each=5)
Log$Age<-rep(c(0,10,15,10),times=5)
RIL<-(Log$RU-(Log$RU*(abs(rnorm(n=20,mean=.6,sd=0.1)))))+(0.5*Log$Age)
Conv<-(Log$RU-(Log$RU*(abs(rnorm(n=20,mean=.2,sd=0.1)))))+(0.2*Log$Age)
Log$RL<-ifelse(Log$Method=="RIL",RIL,Conv)
Log$SDL<-Log$SDU
Log$NL<-Log$NU

#now we perform a meta-analysis using metafor
require(metafor)
ROM<-escalc(data=Log,measure="ROM",m2i=RU,
sd2i=SDU,n2i=NU,m1i=RL,sd1i=SDL,n1i=NL,append=T)
Model1<-rma.mv(yi,vi,random=~(1|Study)+(1|Age),method="ML",data=ROM)
summary(Model1)
forest(Model1)

上述模型是一个空模型,旨在确定截距是否与零有显著差异。在我们的情况下是有的。然而,我还想看看治疗方法的差异是否描述了您可以在此处看到的森林图中所见效应大小的差异。
因此,我运行了这个模型:
Model2<-rma.mv(yi,vi,mods=~Method,random=~(1|Study)+(1|Age),method="ML",data=ROM)
summary(Model2)

看起来不错。

    Multivariate Meta-Analysis Model (k = 20; method: ML)

  logLik  Deviance       AIC       BIC      AICc  
  0.4725   19.8422    7.0550   11.0380    9.7217  

Variance Components: 

outer factor: Age   (nlvls = 3)
inner factor: Study (nlvls = 4)

            estim    sqrt  fixed
tau^2      0.0184  0.1357     no
rho        1.0000             no

Test for Residual Heterogeneity: 
QE(df = 18) = 23.3217, p-val = 0.1785

Test of Moderators (coefficient(s) 2): 
QM(df = 1) = 19.6388, p-val < .0001

Model Results:

           estimate      se     zval    pval    ci.lb    ci.ub     
intrcpt     -0.1975  0.1007  -1.9622  0.0497  -0.3948  -0.0002    *
MethodRIL   -0.4000  0.0903  -4.4316  <.0001  -0.5768  -0.2231  ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

然而,我希望从这个模型中获得一个类似于R平方的拟合度量。人们在过去使用GLMM时遇到了这些问题,但现在有解决方法。我想知道是否有一种类似的方法可以用于元分析?审稿人要求这样做,我不确定是否应该告诉他们这是不可能的。
谢谢您的帮助!

1
这个问题最好在这里问:http://stats.stackexchange.com/ - Andrew Cassidy
这个在CrossValidated上的问题可能会有用:http://stats.stackexchange.com/questions/44676/r-squared-of-mixed-model-using-lmmfit-package。虽然使用了不同的包,但问题是相同的。 - Matt Parker
这篇博客文章建议查看您链接的中川论文,其中包含计算统计数据的R脚本 - 您是否已经查看了它? - Matt Parker
1
哈!那是你的博客。我现在就安静了。 - Matt Parker
1个回答

6

首先,您对rma.mv()函数的语法使用不太正确。针对这两个模型,我认为您实际上想要使用的是:

Model1 <- rma.mv(yi, vi, random = list(~ 1 | Study, ~ 1 | Age), method="ML", data=ROM)
Model2 <- rma.mv(yi, vi, mods = ~ Method, random = list(~ 1 | Study, ~ 1 | Age), method="ML", data=ROM)

现在,关于R平方,您可以计算方差分量的比例减少,作为一种伪R平方值。这只是通常在常规元回归中执行的逻辑扩展。因此,基于上面的模型:
(Model1$sigma2[1] - Model2$sigma2[1]) / Model1$sigma2[1]
(Model1$sigma2[2] - Model2$sigma2[2]) / Model1$sigma2[2]

如果一个值是负数,通常会将其设为零。

如果您想要一个单一的值,也可以通过以下方式计算总方差的比例减少:

(sum(Model1$sigma2) - sum(Model2$sigma2)) / sum(Model1$sigma2)

谢谢你的帮助,Wolfgang。这对我非常有用。它似乎可以很好地适用于我的大多数模型,但对于一些明显支持良好的模型(相对较低的AICc),总方差的比例减少会变成负数。有什么想法吗?如果我找出问题所在,我会进行调查并在这里发布回复。 - Phil_Martin
1
你的意思是说,具有调节变量/预测变量的模型的AICc明显低于仅拦截模型,但总方差的比例减少为负?我想这可能会发生 - 取决于数据集的大小,这些方差分量可能无法被精确估计,从而导致这种反直觉的结果。您还可以考虑根据对数似然值之一计算伪R平方值。例如,请参见:http://www.ats.ucla.edu/stat/mult_pkg/faq/general/Psuedo_RSquareds.htm - Wolfgang

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接