R code for multiple univariate meta-analyses from Chen et al, Statistics in Medicine 2016;35:1405-22

Question

I am trying to run the sample R code for multivariate meta-analysis through multiple univariate meta-analyses provided by Chen et al, in their article entitled "Inference for correlated effect sizes using multiple univariate meta-analyses", published in Statistics in Medicine 2016;35:1405-22.

Unfortunately, the sample code does not work (it is provided below in extenso).

Any suggestion on how to make it run smoothly?

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

Do you know that you can debug R code ? https://cran.r-project.org/web/packages/debug/index.html if you use RSudio it is even easier — statquant, Apr 20 '16 at 19:57
#joran: I thought to first ask the SO community, and then if this fails contact the authors. — Giuseppe Biondi-Zoccai, Apr 20 '16 at 20:12
@statquant: I have been using RStudio to try working it out, but didn't help much, that's why I am asking in SO. — Giuseppe Biondi-Zoccai, Apr 20 '16 at 20:12

IRTFM · Answer 1 · 2016-04-21T15:34:03.840

The error I got was with the "ˆ" character, different than the "^" character. This was a pdf typographical mismatch with the "real" ASCII caret which is the only character that will be parsed by R as an exponentiation operator. After fixing the non-caret exponentiation operator in multiple places the only remaining error which was outside the function and not fatal for it's execution was thrown by

(beta1+beta2)/2

.... and I think that was probably meant by its authors to be:

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

Which was obtained a couple of line later using a different method using a 50-50 weighted average as beta.average.

R code for multiple univariate meta-analyses from Chen et al, Statistics in Medicine 2016;35:1405-22

1 Answers1