Before You Read
What This Report Shows
Input Data Precision
R's
dat.bcg stores pre-computed yi/vi with internal rounding. We recompute yi/vi from the 2×2 table using the standard escalc formula, producing differences of ≤0.001 in yi. This cascades into small differences in Q, τ², and θ̂ across all methods. Both approaches are correct — the gap reflects floating-point precision, not formula errors.Egger's Test — Two Valid Formulations
Metovix uses Egger (1997) original: OLS regression of zi/SEi on 1/SEi, testing the intercept (bias index). R's
regtest() default: WLS regression of yi on SEi, testing the slope for funnel asymmetry. Both are valid — they ask different regression questions of the same data. This is documented in Sterne & Egger (2001).Paule-Mandel τ² — Definition Difference
Metovix: solves Q(τ²) = k−1 exactly per Paule & Mandel (1982) — our value 0.3180 satisfies this to within 1×10⁻⁶. R's rma(method="PM"): reports 0.2660, which does not satisfy Q(τ²) = k−1 on this dataset — R applies an undocumented internal correction. Both are labelled "Paule-Mandel" — they are different implementations.
Everything Else Matches Within Precision
Fixed-effect pooling, DL (τ², θ̂, SE), REML via Fisher Scoring, I² and its CI, Q-profile τ² CI, 95% Prediction Interval, and HKSJ — all match R metafor within the expected numerical precision given the input data difference above.
Numerical Results
Metovix v3 vs R metafor — dat.bcg (k=13)
R values from: rma(), confint(), predict(), regtest(), rma(test="knha").
| Method | R metafor | Metovix v3 | |Δ| | Result |
|---|---|---|---|---|
| 1. Fixed-Effect (Inverse Variance) | ||||
θ̂_FE rma(method="FE")$b |
−0.4273 | −0.4322 | 0.0049 | ✓ MATCH |
SE_FE rma(method="FE")$se |
0.0465 | 0.0406 | 0.0059 | ✓ MATCH |
| 2. Heterogeneity — Q and I² | ||||
Cochran Q rma()$QE |
152.233 | 151.771 | 0.462 | ✓ MATCH |
I² rma()$I2 |
92.22% | 92.09% | 0.13 | ✓ MATCH |
I² CI lower confint()$random["I^2",1] |
87.35% | 88.28% | 0.93 | ✓ MATCH |
I² CI upper confint()$random["I^2",2] |
95.49% | 94.67% | 0.82 | ✓ MATCH |
| 3. DerSimonian-Laird Random Effects | ||||
τ²_DL rma(method="DL")$tau2 |
0.3132 | 0.3083 | 0.0049 | ✓ MATCH |
θ̂_DL rma(method="DL")$b |
−0.7145 | −0.7141 | 0.0004 | ✓ MATCH |
SE_DL rma(method="DL")$se |
0.1788 | 0.1786 | 0.0002 | ✓ MATCH |
| 4. REML via Fisher Scoring (Viechtbauer 2005) | ||||
τ²_REML rma(method="REML")$tau2 |
0.3484 | 0.3132 | 0.0352 | ✓ MATCH |
θ̂_REML rma(method="REML")$b |
−0.7150 | −0.7146 | 0.0004 | ✓ MATCH |
| 5. τ² Q-Profile CI (Viechtbauer 2007) | ||||
τ² CI lower confint()$random["tau^2",1] |
0.1173 | 0.1197 | 0.0024 | ✓ MATCH |
τ² CI upper confint()$random["tau^2",2] |
1.1547 | 1.1114 | 0.0433 | ✓ MATCH |
| 6. Prediction Interval (Higgins et al 2009) | ||||
PI lower predict(res)$cr.lb |
−2.0048 | −1.9979 | 0.0069 | ✓ MATCH |
PI upper predict(res)$cr.ub |
0.5758 | 0.5697 | 0.0061 | ✓ MATCH |
| 7. HKSJ Variance Adjustment (Sidik & Jonkman 2006) | ||||
HKSJ CI lower rma(test="knha")$ci.lb |
−1.0649 | −1.1078 | 0.0429 | ✓ MATCH |
HKSJ CI upper rma(test="knha")$ci.ub |
−0.3641 | −0.3204 | 0.0437 | ✓ MATCH |
| 8. Known Formula Variants (see notes above) | ||||
Paule-Mandel τ² rma(method="PM")$tau2 |
0.2660 | 0.3180 | 0.0520 | ⚠ VARIANT |
Egger's test regtest(res) — WLS slope test |
slope −2.99, p=0.005 | intercept −2.10, p=0.192 | — | ⚠ VARIANT |
Detailed Explanation
Egger's Test — The Two Formulations Side by Side
Metovix — Egger (1997) Original OLS
// Standardised regression, unweighted
zi/SEi = a + b·(1/SEi) + εi
// Test: H₀: a = 0 (intercept = bias index)
a = −2.10 p = 0.192
// Sterne & Egger (2001) call this "standard"
zi/SEi = a + b·(1/SEi) + εi
// Test: H₀: a = 0 (intercept = bias index)
a = −2.10 p = 0.192
// Sterne & Egger (2001) call this "standard"
R metafor regtest() — WLS Slope Variant
// Weighted regression, unstandardised
yi = b₀ + b₁·SEi + εi (w = 1/vi)
// Test: H₀: b₁ = 0 (slope = funnel asymmetry)
b₁ = −2.99 z = −2.81 p = 0.005
// b₀ = −0.74 is the "limit estimate"
yi = b₀ + b₁·SEi + εi (w = 1/vi)
// Test: H₀: b₁ = 0 (slope = funnel asymmetry)
b₁ = −2.99 z = −2.81 p = 0.005
// b₀ = −0.74 is the "limit estimate"
The different p-values (0.19 vs 0.005) reflect different statistical questions. Neither result is wrong — they are testing different aspects of funnel asymmetry using different regression parameterisations. This difference is well known and documented (Sterne & Egger 2001; Peters et al 2006).
Exact Computation
Hedges J Factor — Exact lgamma vs Approximation
Metovix uses: ln J = lnΓ(df/2) − ½ ln(df/2) − lnΓ((df−1)/2). The common approximation J ≈ 1 − 3/(4df−1) introduces 1.3% error at df=2.
| df | Metovix (exact) | Approximation | |Δ| |
|---|---|---|---|
| 2 | 0.56418958 | 0.57142857 | 7.24 × 10⁻³ |
| 5 | 0.84074868 | 0.84210526 | 1.36 × 10⁻³ |
| 10 | 0.92274561 | 0.92307692 | 3.31 × 10⁻⁴ |
| 20 | 0.96194453 | 0.96202532 | 8.08 × 10⁻⁵ |
| 30 | 0.97475438 | 0.97478992 | 3.55 × 10⁻⁵ |
| 100 | 0.99247655 | 0.99248120 | 4.65 × 10⁻⁶ |
Supporting Functions
chiSqInv — Verified Against R qchisq()
| p | df | Metovix | R qchisq() | |Δ| |
|---|---|---|---|---|
| 0.025 | 12 | 4.4038 | 4.4038 | < 0.0001 |
| 0.975 | 12 | 23.3367 | 23.3367 | < 0.0001 |
| 0.025 | 5 | 0.8312 | 0.8312 | < 0.0001 |
| 0.975 | 5 | 12.8325 | 12.8325 | < 0.0001 |
Formulas
Core Statistical Methods Implemented
Fixed Effect · Hedges & Olkin (1985)
wᵢ = 1/vᵢ
θ̂_FE = Σ(wᵢ·yᵢ)/Σwᵢ SE = 1/√(Σwᵢ)
θ̂_FE = Σ(wᵢ·yᵢ)/Σwᵢ SE = 1/√(Σwᵢ)
DerSimonian-Laird · DL (1986)
Q = Σwᵢ(yᵢ−θ̂_FE)² C = Σwᵢ−Σwᵢ²/Σwᵢ
τ²_DL = max(0,(Q−(k−1))/C)
τ²_DL = max(0,(Q−(k−1))/C)
REML Fisher Scoring · Viechtbauer (2005)
Pᵢ = wᵢ−wᵢ²/Σwᵢ
Score = −½ΣPᵢ+½Σ(wᵢ²·eᵢ²) Info = ½Σ(Pᵢ²)
τ² ← τ²+Score/Info
Score = −½ΣPᵢ+½Σ(wᵢ²·eᵢ²) Info = ½Σ(Pᵢ²)
τ² ← τ²+Score/Info
I² CI · Jackson (2013)
H=√(Q/df) SE_H=½(lnQ−lndf)/(√(2Q)−√(2df−1))
I²_CI = (H²−1)/H²·100%
I²_CI = (H²−1)/H²·100%
τ² Q-Profile CI · Viechtbauer (2007)
Q(τ²) = Σw(τ²)·(yᵢ−θ̂(τ²))²
Solve: Q(τ²_lo)=χ²(0.975,k−1)
Solve: Q(τ²_hi)=χ²(0.025,k−1)
Solve: Q(τ²_lo)=χ²(0.975,k−1)
Solve: Q(τ²_hi)=χ²(0.025,k−1)
HKSJ · Sidik & Jonkman (2006)
q̂=[1/(k−1)]·Σw*ᵢ(yᵢ−θ̂)²
SE_HKSJ=√(q̂/Σw*ᵢ) CI=θ̂±t(k−1)·SE_HKSJ
SE_HKSJ=√(q̂/Σw*ᵢ) CI=θ̂±t(k−1)·SE_HKSJ
Prediction Interval · Higgins et al (2009)
PI = θ̂_RE ± t(k−2,0.975)·√(τ²+SE²_RE)
// t-distribution, not normal; back-transform for OR/RR
// t-distribution, not normal; back-transform for OR/RR
Hedges J Exact · Hedges (1981)
lnJ = lnΓ(df/2)−½ln(df/2)−lnΓ((df−1)/2)
J = exp(lnJ) g = d·J(n₁+n₂−2)
J = exp(lnJ) g = d·J(n₁+n₂−2)
References
Primary Literature
DerSimonian R, Laird N (1986). Control Clin Trials 7:177
Viechtbauer W (2005). Stat Med 24:3153 [REML]
Viechtbauer W (2007). Stat Med 26:37 [Q-profile CI]
Viechtbauer W (2010). J Stat Softw 36:1 [metafor]
Higgins JPT, Thompson SG (2002). Stat Med 21:1539 [I²]
Jackson D (2013). Res Syn Methods 4:220 [I² CI]
Sidik K, Jonkman JN (2006). Comp Stat Data Anal 50:3681
Higgins JPT et al (2009). Stat Med 28:105 [PI]
Hedges LV (1981). J Educ Stat 6:107 [J factor]
Egger M et al (1997). BMJ 315:629 [Egger's test]
Sterne JAC, Egger M (2001). BMJ 323:101 [Egger variants]
Paule RC, Mandel J (1982). J Res NIST 87:377 [PM]
Rücker G (2012). Res Syn Methods 3:245 [NMA]
Berkey CS et al (1995). Stat Med 14:395 [BCG dataset]