1

我正在尝试通过 Chen 等人提供的多个单变量荟萃分析运行示例 R 代码以进行多变量荟萃分析，该文章发表于2016 年医学统计;35:1405-22。

不幸的是，示例代码不起作用（下面详细提供）。

关于如何使其顺利运行的任何建议？

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
4

1 回答 1

2

我得到的错误是"ˆ"角色，与角色不同"^"。这是与“真实”ASCII 插入符号的 pdf 印刷不匹配，这是 R 将作为求幂运算符解析的唯一字符。在多个地方修复了非插入符号求幂运算符后，唯一剩余的错误是在函数之外并且对于它的执行不是致命的被抛出

(beta1+beta2)/2

....我认为这可能是其作者的意思：

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

这是在几行之后使用不同的方法获得的，使用 50-50 的加权平均值作为beta.average.

于 2016-04-20T20:39:15.653 回答