3

I'm having a problem with metafor's rma-function, which inappropriately (I suppose!) gives both tauˆ2 and Iˆ2 zero values. This questions was previously asked on Stack Overflow, and then it was suggested using (the default) REML estimator of tauˆ2. However, despite using REML, I'm facing this issue. Here are the summary results of my data from rma:

Random-Effects Model (k = 44; tau^2 estimator: REML)

tau^2 (estimated amount of total heterogeneity): 0 (SE = 0.0136)
tau (square root of estimated tau^2 value):      0
I^2 (total heterogeneity / total variability):   0.00%
H^2 (total variability / sampling variability):  1.00

Test for Heterogeneity:
Q(df = 43) = 7.7261, p-val = 1.0000

Model Results:

estimate      se    zval    pval    ci.lb   ci.ub 
  0.0589  0.0407  1.4451  0.1484  -0.0210  0.1387 

And I think those zero values truly are errors, as Review Manager 5.4 gives the following, reasonable results from the same data:

Heterogeneity: Tau² = 0,01; Chi² = 123,95, df = 43 (P < 0.00001); I² = 65%

Please see my code below. The error most probably is in it, but I just can't figure it out. Could someone help me on this, thank you so much!

metareg <- read.table(text="
effect_size variance    study   n   SE  subgroup
0.0368  0.0806  A-1 230 0.0187  A
-0.0509 0.09    A-2 22  0.064   A
-0.1031 0.1081  A-3 16  0.0822  A
0.0731  0.0683  A-4 113 0.0246  A
0.1105  0.0448  A-5 90  0.0223  A
0.1193  0.0484  A-6 28  0.0416  A
0.1405  0.0524  A-7 50  0.0324  A
0.0157  0.0622  A-8 13  0.0692  A
-0.0528 0.0864  A-9 64  0.0367  A
-0.011  0.0799  A-10    13  0.0784  A
0.0713  0.0768  A-11    40  0.0438  A
-0.1307 0.0996  A-12    17  0.0766  A
-0.0658 0.106   A-13    18  0.0767  A
-0.012  0.097   A-14    86  0.0336  A
-0.0163 0.0914  A-15    50  0.0428  A
0.0479  0.0721  A-16    27  0.0517  A
0.11    0.0517  A-17    151 0.0185  A
0.2944  0.0303  A-18    19  0.0399  A
0.0138  0.0713  A-19    16  0.0668  A
0.0082  0.0884  A-20    32  0.0526  A
0.0726  0.0868  A-21    33  0.0513  A
0.1554  0.0163  A-22    10  0.0404  A
0.1247  0.0396  A-23    36  0.0332  A
0.0954  0.0685  A-24    47  0.0382  A
0.0299  0.1061  B-1 401 0.0163  B
-0.0516 0.1111  B-2 68  0.0404  B
-0.0857 0.1165  B-3 47  0.0498  B
-0.1722 0.1083  B-4 23  0.0686  B
-0.0324 0.0807  B-5 16  0.071   B
-0.0581 0.1168  B-6 49  0.0488  B
0.049   0.0982  B-7 43  0.0478  B
0.0885  0.1067  B-8 56  0.0437  B
-0.0056 0.1138  B-9 214 0.0231  B
-0.1039 0.1183  B-10    53  0.0472  B
0.0065  0.1172  B-11    42  0.0528  B
0.1047  0.0957  B-12    72  0.0365  B
0.1337  0.0868  B-13    46  0.0434  B
-0.2632 0.0983  B-14    66  0.0386  B
0.1366  0.0862  B-15    51  0.0411  B
-0.0518 0.1034  B-16    12  0.0928  B
-0.1573 0.11    B-17    17  0.0804  B
0.2218  0.0998  B-18    18  0.0745  B
0.2097  0.0932  B-19    13  0.0847  B
0.0877  0.1061  B-20    14  0.087   B
", header=TRUE)

res.rma <- rma(effect_size, variance, data=metareg)
res.rma
forest(res.rma, slab = paste(metareg$study))

I'm not sure if it helps, but here are the values of respective "intermediate variables" and (probably) correct tauˆ2 and Iˆ2, when these Review Manager's statistical calculations are done in Excel:

"intermediate variables" and results in Excel

araappa
  • 31
  • 1

2 Answers2

2

The variance variable does not appear to contain the sampling variances of the effect size estimates. The SE variable appears to be sqrt(variance/n) so that would be appropriate if the effect sizes are just means and variance represents the variance of the raw measurements, but rma(effect_size, sei=SE, data=metareg, method="DL") still differs from the RevMan results. rma(effect_size, variance/(n/2), data=metareg, method="DL") comes close and the difference could be due to the values in shown being rounded, but in the end, I am just guessing. I can guarantee you that metafor and RevMan give the exact same results when they are given the same data (see, for example: https://www.metafor-project.org/doku.php/plots:forest_plot_revman), so just make sure you pass the exact same estimates and corresponding sampling variances or standard errors (the latter via argument sei) to rma() and use the DL estimator and the results must match up.

Wolfgang
  • 15,542
  • 1
  • 47
  • 74
0

Wolfgang, my warmest thanks to you! I apologize, but I had an error in SE calculation in my earlier test data (fixed after comma-to-dot decimal separator conversion), which I posted yesterday. After re-calculating this dot-separated table with the following calculations (in parentheses values of A-1 -study for an example):

  • the original 'effect size' is observed cure rate in each study, thus it's a proportion (and the actual effect_size in the example is a derivative from observed cure rate, sharing the same CI-range). (A-1: cure_rate = 0.7913)

  • 95 % CI for each proportion is estimated with Wilson's method (A-1: CI_low = 0.7342, CI_high = 0.8388 -> CI_range = 0.1046)

  • SE = (CI_high - CI_low)/3.92 (A-1: SE = 0.1046 / 3.92 = 0.02668)

So far so good, with these calculations rma(effect_size, sei=SE, data=metareg, method="DL") works perfectly and corresponds to that of RevMan's! However, although I don't need variance anymore, I'd like to know what's the problem with the following calculation of variance:

  • variance = n * SEˆ2 (A-1: 230 * 0.02668 * 0.02668 = 0.1637)

I suppose that this formula corresponds to that what Wolfgang suggested above ("SE variable appears to be sqrt(variance/n) so that would be appropriate if the effect sizes are just means and variance represents the variance of the raw measurements") and that the stated conditions are also fulfilled in this case, or aren't they? Here is the example data with correct SEs:

metareg <- read.table(text="
effect_size variance  study n SE  subgroup
0.0368  0.1638  A-1 230 0.0267  A
-0.011  0.1498  A-2 13  0.1073  A
0.0713  0.1532  A-3 40  0.0619  A
-0.1307 0.186 A-4 17  0.1046  A
-0.0658 0.1976  A-5 18  0.1048  A
-0.012  0.1946  A-6 86  0.0476  A
-0.0163 0.1817  A-7 50  0.0603  A
0.0479  0.1429  A-8 27  0.0727  A
0.11  0.1059  A-9 151 0.0265  A
-0.0509 0.1729  A-10  22  0.0887  A
-0.1031 0.1991  A-11  16  0.1115  A
0.0731  0.1388  A-12  113 0.035 A
0.2944  0.0695  A-13  19  0.0605  A
0.1105  0.0926  A-14  90  0.0321  A
0.0138  0.1381  A-15  16  0.0929  A
0.0082  0.1734  A-16  32  0.0736  A
0.0726  0.1707  A-17  33  0.0719  A
0.1554  0.0501  A-18  10  0.0708  A
0.1193  0.1005  A-19  28  0.0599  A
0.1247  0.0841  A-20  36  0.0483  A
0.0954  0.1382  A-21  47  0.0542  A
0.1405  0.1076  A-22  50  0.0464  A
0.0157  0.1216  A-23  13  0.0967  A
-0.0528 0.173 A-24  64  0.052 A
0.0299  0.2155  B-1 401 0.0232  B
-0.0516 0.2209  B-2 68  0.057 B
-0.0857 0.2287  B-3 47  0.0698  B
-0.1722 0.2055  B-4 23  0.0945  B
-0.0324 0.1536  B-5 16  0.098 B
-0.0581 0.2296  B-6 49  0.0684  B
0.049 0.1936  B-7 43  0.0671  B
0.0885  0.2113  B-8 56  0.0614  B
-0.0056 0.2302  B-9 214 0.0328  B
-0.1039 0.233 B-10  53  0.0663  B
0.0065  0.2289  B-11  42  0.0738  B
0.1047  0.1915  B-12  72  0.0516  B
0.1337  0.1725  B-13  46  0.0612  B
-0.2632 0.1961  B-14  66  0.0545  B
0.1366  0.1719  B-15  51  0.0581  B
-0.0518 0.1854  B-16  12  0.1243  B
-0.1573 0.2033  B-17  17  0.1094  B
0.2218  0.1872  B-18  18  0.102 B
0.2097  0.1709  B-19  13  0.1147  B
0.0877  0.193 B-20  14  0.1174  B
", header=TRUE)

### with variance this does not work appropriately
#
# res.rma <- rma(effect_size, variance, data=metareg)
# res.rma
# forest(res.rma, slab = paste(metareg$study))

### with sei=SE this works and correspond to that of RevMan's

res.rma <- rma(effect_size, sei=SE, data=metareg, method="DL")
res.rma
forest(res.rma, slab = paste(metareg$study))
araappa
  • 31
  • 1
  • So Wolfgang's answer should have been checkmarked and the duplicate on SO should have been deleted. – DWin Feb 09 '22 at 16:26