Получение R в квадрате из многоуровневой модели со смешанными эффектами в metafor

Я провожу мета-анализ в 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 , но есть теперь способы сделать это. Мне было интересно, знает ли кто-нибудь хороший способ сделать что-то подобное с помощью метаанализа? У меня есть рецензенты, которые просят об этом, и я не уверен, стоит ли мне просто сказать им, что это невозможно или нет.

Заранее спасибо за помощь!


person Phil_Martin    schedule 12.03.2014    source источник
comment
Этот вопрос, вероятно, лучше задать здесь: stats.stackexchange.com.   -  person Andrew Cassidy    schedule 12.03.2014
comment
Этот вопрос о CrossValidated может быть полезен. Другой пакет, но проблемы те же.   -  person Matt Parker    schedule 12.03.2014
comment
В этом сообщении блога говорится, что статья Накагавы, на которую вы ссылаетесь, содержит сценарий R. для расчета их статистики - вы это проверили?   -  person Matt Parker    schedule 12.03.2014
comment
Ха! Это твой блог. Я сейчас просто молчу.   -  person Matt Parker    schedule 12.03.2014


Ответы (1)


Во-первых, вы не совсем используете правильный синтаксис для функции 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)
person Wolfgang    schedule 12.03.2014
comment
Спасибо за это Вольфгангу, очень помог. Кажется, он отлично работает с большинством моих моделей, но для нескольких явно хорошо поддерживаемых моделей (относительно низкий AICc) пропорциональное уменьшение общей дисперсии оказывается отрицательным. Любые идеи? Я проведу расследование и опубликую здесь, если я выясню, в чем проблема. - person Phil_Martin; 13.03.2014
comment
Вы имеете в виду, что модель с модераторами / предикторами имеет заметно более низкий AICc, чем модель только с перехватом, и все же пропорциональное уменьшение общей дисперсии отрицательно? Я полагаю, что это могло произойти - в зависимости от размера вашего набора данных эти компоненты дисперсии могут быть не очень точно оценены, и это может привести к таким нелогичным результатам. Вы также можете подумать о вычислении одного из значений псевдоквадрата R на основе логарифмических вероятностей. См., Например: ats.ucla.edu/stat/mult_pkg /faq/general/Psuedo_RSquareds.htm - person Wolfgang; 13.03.2014