来自Chen等人的多个单变量meta分析的R代码，2016年医学统计学；35:1405-22

Question

我正在尝试通过Chen等人提供的多个单变量meta分析来运行多变量meta分析的样本R代码，这篇文章发表在Statistics in Medicine 2016;35:1405-22上，题为“使用多个单变量meta分析推断相关效应大小”。

不幸的是，示例代码不起作用(下面将对其进行扩展)。

有没有关于如何让它顺利运行的建议？

install.packages("xmeta")
library(xmeta)

## working example
ys<-matrix(c(1.139434283, 1.446918983, 1.704748092, 0.470003629,
             0.85566611, 1.440361582, 0.186585956, 1.504077397,
             1.540445041, 1.665007764, 3.218875825, 1.299282984,
             0.661398482, 3.283414346, 4.919980926, 1.386294361,
             3.218875825, 2.197224577, 2.268683541, -1.145132304),
           ncol=2)

vars<-matrix(c(0.164999999811686, 0.308823529160345,
               0.0738636362264944, 0.162499999665475,
               0.0838235293119684, 0.137426900570286,
               0.0938352427206998, 0.203703703772108,
               0.40476190440518, 0.169884169901165,
               1.04000000057403, 0.424242423879069,
               0.0947580646552244, 0.34583333355377,
               1.00729926944171, 0.208333333709766,
               2.07999999946466, 0.555555554810304,
               0.367816092120049, 0.188311688602409),ncol=2)
ys
vars

MMoM <- function(ys,vars){
  ys1 = ys[,1]
  ys2 = ys[,2]
  vars1 = vars[,1]
  vars2 = vars[,2]
  w1 = 1/(vars1)
  w2 = 1/(vars2)
  y1.weight = sum(w1*ys1)/sum(w1)
  y2.weight = sum(w2*ys2)/sum(w2)
  n1 = sum(1-1*(vars1 > 10ˆ4))
  n2 = sum(1-1*(vars2 > 10ˆ4))
  Q1 = sum(w1*(ys1-y1.weight)ˆ2)
  Q2 = sum(w2*(ys2-y2.weight)ˆ2)
  tau1.2.hat = max(0, (Q1-(n1-1))/(sum(w1)-sum(w1ˆ2)/sum(w1)))
  tau2.2.hat = max(0, (Q2-(n2-1))/(sum(w2)-sum(w2ˆ2)/sum(w2)))
  w1.star = 1/(vars1 + tau1.2.hat)
  w2.star = 1/(vars2 + tau2.2.hat)
  beta1.hat = sum(w1.star*ys1)/sum(w1.star)
  beta2.hat = sum(w2.star*ys2)/sum(w2.star)
  var.beta1.hat = 1/sum(w1.star)
  var.beta2.hat = 1/sum(w2.star)
  mycov.beta = sum((w1.star/sum(w1.star))*(w2.star/sum(w2.star))
                   *(ys1 - beta1.hat)*(ys2 - beta2.hat))
  beta.hat = c(beta1.hat,beta2.hat)
  sigma.hat = matrix(c(var.beta1.hat,mycov.beta,mycov.beta,
                       var.beta2.hat),nrow = 2, byrow = T)
  result = list(beta.hat=beta.hat,beta.cov=sigma.hat)
  return(result)
}

## Apply the MMoM method to the working example
MMoM.fit = MMoM(ys,vars)
MMoM.fit

## obtain the estimate and standard error of delta = beta1-beta2
myV = c(1, -1)
delta.hat = (myV%*%(MMoM.fit$beta.hat))[1,1]
delta.se = (sqrt(t(myV)%*%(MMoM.fit$beta.cov)%*%myV))[1,1]
delta.hat
delta.se

## obtain the estimate and standard error of beta.average =
(beta1+beta2)/2
myV2 = c(0.5, 0.5)
beta.average = (myV2%*%(MMoM.fit$beta.hat))[1,1]
beta.average.se = (sqrt(t(myV2)%*%(MMoM.fit$beta.cov)%*%myV2))[1,1]
beta.average
beta.average.se

Answer 1

我得到的错误是与"^"字符不同的"ˆ"字符。这是一个pdf排版与“真正的”ASCII插入符号不匹配，这是唯一一个将被R解析为指数运算符的字符。在多个位置修复了非插入符号求幂运算符之后，仅有的一个在函数外部的错误被抛出，因为它的执行不是致命的

(beta1+beta2)/2

……我认为这可能是它的作者的意思是：

> with( MMoM.fit, (beta.hat[1]+beta.hat[2])/2)
[1] 1.55537

这是在几行之后使用另一种方法获得的，该方法使用50-50的加权平均值作为beta.average。