The data of the paper titled “Minimally Invasive Percutaneous Nephrolithotomy versus Retrograde Intrarenal Surgery for Upper Urinary Stones: A Systematic Review and Meta-Analysis” by Jiang H et al (2017) were re-analyzed to reproduce the MH estimates presented in the paper. In addition, RR, RD,OR and Peto OR methods were also presented.

stone <- read.table("stone.txt", sep = ",", header = TRUE)
stone
##                 Study event_mini n_mini event_rirs n_rirs
## 1    Ozgor et al 2016         45     56         43     56
## 2       Gu et al 2013         30     30         26     29
## 3       Hu et al 2016         87    104         53     80
## 4    Kirac et al 2013         33     37         32     36
## 5    Knoll et al 2011         25     25         18     21
## 6    Kruck et al 2013        137    172         84    108
## 7    Kumar et al 2015         39     41         37     43
## 8      Lee et al 2015         30     35         32     33
## 9      Pan et al 2013         57     59         40     56
## 10  Sabnis et al 2012         32     32         31     32
## 11 Wilhelm et al 2015         23     25         24     25
## 12    Zeng et al 2015         38     53         23     53
## 13   Zhang et al 2014         30     32         37     44
# The events
xE = stone$event_mini
# the total number
nE = stone$n_mini
# the risk
pE = xE/nE
# print the risk
pE
##  [1] 0.8035714 1.0000000 0.8365385 0.8918919 1.0000000 0.7965116 0.9512195
##  [8] 0.8571429 0.9661017 1.0000000 0.9200000 0.7169811 0.9375000
xC = stone$event_rirs
nC = stone$n_rirs
pC = xC/nC
pC
##  [1] 0.7678571 0.8965517 0.6625000 0.8888889 0.8571429 0.7777778 0.8604651
##  [8] 0.9696970 0.7142857 0.9687500 0.9600000 0.4339623 0.8409091
# The risk-ratios
ES = pE/pC
# print the RR
ES
##  [1] 1.0465116 1.1153846 1.2626996 1.0033784 1.1666667 1.0240864 1.1054713
##  [8] 0.8839286 1.3525424 1.0322581 0.9583333 1.6521739 1.1148649
#confidence intervel approach for statistical significance
# calculate the log risk-ratio
lnES = log(ES)
# print the log risk-ratio
lnES
##  [1]  0.045462374  0.109199292  0.233251941  0.003372684  0.154150680
##  [6]  0.023800877  0.100271782 -0.123379021  0.301986061  0.031748698
## [11] -0.042559614  0.502091944  0.108733200
# calculate the variance
VlnES = 1/stone$event_mini - 1/stone$n_mini + 1/stone$event_rirs - 1/stone$n_rirs
# print the variance
VlnES
##  [1] 0.009763750 0.003978780 0.008246793 0.006748225 0.007936508
##  [6] 0.004130819 0.005021995 0.005708874 0.007737564 0.001008065
## [11] 0.005144928 0.032058201 0.006383088
# upper CI
ciup = lnES+1.96*sqrt(VlnES)
print(ciup)
##  [1] 0.23913329 0.23283133 0.41124317 0.16438195 0.32876136 0.14977292
##  [7] 0.23916921 0.02471288 0.47439437 0.09397876 0.09802756 0.85302611
## [13] 0.26532589
# lower CI
cilow = lnES-1.96*sqrt(VlnES)
print(cilow)
##  [1] -0.14820854 -0.01443275  0.05526071 -0.15763658 -0.02046000
##  [6] -0.10217116 -0.03862564 -0.27147093  0.12957775 -0.03048136
## [11] -0.18314679  0.15115778 -0.04785949
# then transform back to the original scale
cat("The low CI is:", exp(cilow),"\n\n")
## The low CI is: 0.8622513 0.9856709 1.056816 0.8541601 0.9797479 0.902875 0.9621108 0.7622574 1.138348 0.9699785 0.8326459 1.16318 0.9532677
cat("The upper CI is:", exp(ciup),"\n\n")
## The upper CI is: 1.270148 1.262169 1.508692 1.178664 1.389246 1.16157 1.270193 1.025021 1.607041 1.098536 1.102993 2.346738 1.303856
#The p value apparoach
# The z-value
z = lnES/sqrt(VlnES)
# Print the z-values
z
##  [1]  0.46009104  1.73119053  2.56851866  0.04105641  1.73033709
##  [6]  0.37031804  1.41494841 -1.63292438  3.43308663  0.99995800
## [11] -0.59334605  2.80423030  1.36096440
pval = 2*(1-pnorm(abs(z)))
cat("p-values = ", pval, sep=" ", "\n\n")
## p-values =  0.6454509 0.08341779 0.01021342 0.9672509 0.08357007 0.7111455 0.1570836 0.1024849 0.0005967514 0.3173308 0.5529496 0.005043686 0.1735249
#The postpower approach
# load the pwr library
library(pwr)
# calculate the power and print it
pow.study = pwr.2p2n.test(ES.h(pE,pC),n1=stone$n_mini,
n2=stone$n_rirs,sig.level=0.05)
pow.study
## 
##      difference of proportion power calculation for binomial distribution (arcsine transformation) 
## 
##               h = 0.087106467, 0.654908852, 0.407349632, 0.009612818, 0.775193373, 0.045777435, 0.320254057, 0.425255282, 0.757542651, 0.355421202, 0.170797263, 0.581348612, 0.315190573
##              n1 = 56, 30, 104, 37, 25, 172, 41, 35, 59, 32, 25, 53, 32
##              n2 = 56, 29, 80, 36, 21, 108, 43, 33, 56, 32, 25, 53, 44
##       sig.level = 0.05
##           power = 0.07467296, 0.71052225, 0.78207415, 0.05019317, 0.74501817, 0.06607345, 0.31138493, 0.41797081, 0.98215860, 0.29555242, 0.09270927, 0.84913082, 0.27360542
##     alternative = two.sided
## 
## NOTE: different sample sizes
#the risk ratio meta analysis step by step
VlnES
##  [1] 0.009763750 0.003978780 0.008246793 0.006748225 0.007936508
##  [6] 0.004130819 0.005021995 0.005708874 0.007737564 0.001008065
## [11] 0.005144928 0.032058201 0.006383088
# Inverse weighting
w = 1/VlnES
w
##  [1] 102.41966 251.33333 121.25926 148.18710 126.00000 242.08273 199.12406
##  [8] 175.16588 129.23964 992.00000 194.36620  31.19327 156.66399
# the inverse of variance for each study
fwi = 1/VlnES
fwi
##  [1] 102.41966 251.33333 121.25926 148.18710 126.00000 242.08273 199.12406
##  [8] 175.16588 129.23964 992.00000 194.36620  31.19327 156.66399
# the total weight
fw = sum(fwi)
fw
## [1] 2869.035
# the relative weight for each study
rw = fwi/fw
rw
##  [1] 0.03569829 0.08760204 0.04226482 0.05165050 0.04391720 0.08437775
##  [7] 0.06940454 0.06105393 0.04504638 0.34576084 0.06774619 0.01087239
## [13] 0.05460512
# the weighted mean
flnES = sum(lnES*rw)
flnES
## [1] 0.06252015
# the variance for the weighted mean
var= 1/fw;
var
## [1] 0.0003485492
# the RR
exp(flnES)
## [1] 1.064516
# the lower and upper CI bounds
exp(flnES-1.96*sqrt(var))
## [1] 1.026267
exp(flnES+1.96*sqrt(var))
## [1] 1.10419
# estimate heterogeneity statistic Q
Q = sum(fwi*lnES^2)-(sum(fwi*lnES))^2/fw
Q
## [1] 29.247
# Get the degrees of freedom
K = dim(stone)[1]; df = K -1
df
## [1] 12
# calculate the statistical significance: p-value
pval4Q = pchisq(Q, df, lower.tail=F)
pval4Q
## [1] 0.003620415
#risk ratio meta analysis with R package meta
# call "metabin" for RR meta-analysis
library(meta)
## Loading 'meta' package (version 4.8-2).
## Type 'help(meta)' for a brief overview.
RR.stone = metabin(event_mini,n_mini,event_rirs,n_rirs,studlab=Study,
data=stone, method="Inverse", sm="RR")
# print the analysis summary
summary(RR.stone)
## Number of studies combined: k = 13
## 
##                          RR           95%-CI    z  p-value
## Fixed effect model   1.0640 [1.0261; 1.1033] 3.35   0.0008
## Random effects model 1.0858 [1.0199; 1.1559] 2.58   0.0099
## 
## Quantifying heterogeneity:
##  tau^2 = 0.0071; H = 1.56 [1.15; 2.12]; I^2 = 59.0% [24.3%; 77.8%]
## 
## Test of heterogeneity:
##      Q d.f.  p-value
##  29.24   12   0.0036
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
## - Continuity correction of 0.5 in studies with zero cell frequencies
forest(RR.stone)

#Analysis with risk difference
# the risk difference
ESRD = pE-pC
ESRD
##  [1]  0.035714286  0.103448276  0.174038462  0.003003003  0.142857143
##  [6]  0.018733850  0.090754396 -0.112554113  0.251815981  0.031250000
## [11] -0.040000000  0.283018868  0.096590909
# calculate the variance
VarRD = pE*(1-pE)/nE + pC*(1-pC)/nC
VarRD
##  [1] 0.0060017310 0.0031981631 0.0041097475 0.0053494503 0.0058309038
##  [6] 0.0025426963 0.0039239375 0.0043889894 0.0041993862 0.0009460449
## [11] 0.0044800000 0.0084633624 0.0048715318
# calculate the standard error
SERD = sqrt(VarRD)
SERD
##  [1] 0.07747084 0.05655230 0.06410731 0.07313994 0.07636035 0.05042516
##  [7] 0.06264134 0.06624945 0.06480267 0.03075784 0.06693280 0.09199653
## [13] 0.06979636
# call metabin for RD meta-analysis
RD.stone = metabin(event_mini,n_mini,event_rirs,n_rirs,studlab=Study,
data=stone, method="Inverse", sm="RD")
# print the summary
summary(RD.stone)
## Number of studies combined: k = 13
## 
##                          RD           95%-CI    z  p-value
## Fixed effect model   0.0691 [0.0346; 0.1037] 3.93 < 0.0001
## Random effects model 0.0771 [0.0208; 0.1334] 2.68   0.0073
## 
## Quantifying heterogeneity:
##  tau^2 = 0.0063; H = 1.59 [1.17; 2.16]; I^2 = 60.5% [27.6%; 78.5%]
## 
## Test of heterogeneity:
##      Q d.f.  p-value
##  30.40   12   0.0024
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
## - Continuity correction of 0.5 in studies with zero cell frequencies
forest(RD.stone)

#meta analysis with Odds ratio
oddsE = pE/(1-pE)
oddsE
##  [1]  4.090909       Inf  5.117647  8.250000       Inf  3.914286 19.500000
##  [8]  6.000000 28.500000       Inf 11.500000  2.533333 15.000000
oddsC = pC/(1-pC)
oddsC
##  [1]  3.3076923  8.6666667  1.9629630  8.0000000  6.0000000  3.5000000
##  [7]  6.1666667 32.0000000  2.5000000 31.0000000 24.0000000  0.7666667
## [13]  5.2857143
OR = oddsE/oddsC
OR
##  [1]  1.2367865        Inf  2.6071032  1.0312500        Inf  1.1183673
##  [7]  3.1621622  0.1875000 11.4000000        Inf  0.4791667  3.3043478
## [13]  2.8378378
LogOR= log(OR)
LogOR
##  [1]  0.21251646         Inf  0.95823973  0.03077166         Inf
##  [6]  0.11186990  1.15125602 -1.67397643  2.43361336         Inf
## [11] -0.73570679  1.19523912  1.04304244
VLogOR = 1/334+ 1/(4995-334)+1/418+1/(5006-418)
VLogOR
## [1] 0.005818863
SE.LogOR=sqrt(VLogOR)
SE.LogOR
## [1] 0.07628147
OR.stone = metabin(event_mini,n_mini,event_rirs,n_rirs,studlab=Study,
data=stone, method="Inverse", sm="OR")
summary(OR.stone)
## Number of studies combined: k = 13
## 
##                          OR           95%-CI    z  p-value
## Fixed effect model   1.8879 [1.3769; 2.5886] 3.95 < 0.0001
## Random effects model 2.0215 [1.2433; 3.2868] 2.84   0.0045
## 
## Quantifying heterogeneity:
##  tau^2 = 0.2719; H = 1.31 [1.00; 1.82]; I^2 = 42.0% [0.0%; 69.8%]
## 
## Test of heterogeneity:
##      Q d.f.  p-value
##  20.68   12   0.0552
## 
## Details on meta-analytical method:
## - Inverse variance method
## - DerSimonian-Laird estimator for tau^2
## - Continuity correction of 0.5 in studies with zero cell frequencies
forest(OR.stone)

#meta analysis using mantel- Haenzel method
OR
##  [1]  1.2367865        Inf  2.6071032  1.0312500        Inf  1.1183673
##  [7]  3.1621622  0.1875000 11.4000000        Inf  0.4791667  3.3043478
## [13]  2.8378378
w0 = (stone$n_mini-stone$event_mini)*stone$event_rirs/(stone$n_mini+stone$n_rirs)
w0
##  [1]  4.2232143  0.0000000  4.8967391  1.7534247  0.0000000 10.5000000
##  [7]  0.8809524  2.3529412  0.6956522  0.0000000  0.9600000  3.2547170
## [13]  0.9736842
TotWeight = sum(w0)
TotWeight
## [1] 30.49132
relWt = w0/TotWeight
relWt
##  [1] 0.13850544 0.00000000 0.16059450 0.05750569 0.00000000 0.34436024
##  [7] 0.02889190 0.07716756 0.02281476 0.00000000 0.03148436 0.10674239
## [13] 0.03193316
OR.MH = sum(relWt*OR)
OR.MH
## [1] NaN
Var.ORMH1 = 1/TotWeight
Var.ORMH1
## [1] 0.03279621
lowCI = OR.MH-1.96*sqrt(Var.ORMH1)
upCI = OR.MH+1.96*sqrt(Var.ORMH1)
cat("MH estimate=", round(OR.MH,4),sep="", "\n\n")
## MH estimate=NaN
cat("The 95% CI for MH =(", round(lowCI, 4), ",",
round(upCI, 4), ")", sep="", "\n\n")
## The 95% CI for MH =(NaN,NaN)
#Emerson's approximation
A = stone$event_mini; B = stone$n_mini-stone$event_mini
C = stone$event_rirs; D = stone$n_rirs - stone$event_rirs
N = A+B+C+D

T1 = (A+D)/N; T2 = (B+C)/N
T3 = A*D/N; T4 = B*C/N
ST3 = sum(T3); ST4 = sum(T4)

Var.lnOddsRatio = 0.5*(
(T1*T3)/ST3^2+(T1*T4+T2*T3)/(ST3*ST4)+T2*T4/ST4^2)
Var.lnOddsRatio
##  [1] 0.0027455277 0.0002999103 0.0042330185 0.0010748530 0.0003097344
##  [6] 0.0063248912 0.0009047728 0.0010729507 0.0017353944 0.0001012411
## [11] 0.0004890577 0.0031906271 0.0009677349
Var.lnMH = sum(Var.lnOddsRatio)
Var.lnMH
## [1] 0.02344971
lowCI.lnMH = log(OR.MH)-1.96*sqrt(Var.lnMH)
upCI.lnMH = log(OR.MH)+1.96*sqrt(Var.lnMH)
cat("The 95% CI for log-MH =(", round(lowCI.lnMH, 4), ",",
round(upCI.lnMH, 4), ")", sep="", "\n\n")
## The 95% CI for log-MH =(NaN,NaN)
lowCI.MH = exp(lowCI.lnMH)
upCI.MH = exp(upCI.lnMH)
cat("The 95% CI for MH =(", round(lowCI.MH, 4), ",",
round(upCI.MH, 4), ")", sep="", "\n\n")
## The 95% CI for MH =(NaN,NaN)
# MH OR meta-analysis
ORMH.stone = metabin(event_mini,n_mini,event_rirs,n_rirs,studlab=Study,
data=stone, method="MH", sm="OR")
# print the summary
summary(ORMH.stone)
## Number of studies combined: k = 13
## 
##                          OR           95%-CI    z  p-value
## Fixed effect model   1.9597 [1.4552; 2.6390] 4.43 < 0.0001
## Random effects model 2.0218 [1.2425; 3.2897] 2.83   0.0046
## 
## Quantifying heterogeneity:
##  tau^2 = 0.2736; H = 1.31 [1.00; 1.82]; I^2 = 42.1% [0.0%; 69.9%]
## 
## Test of heterogeneity:
##      Q d.f.  p-value
##  20.73   12   0.0544
## 
## Details on meta-analytical method:
## - Mantel-Haenszel method
## - DerSimonian-Laird estimator for tau^2
## - Continuity correction of 0.5 in studies with zero cell frequencies
forest(ORMH.stone)

#Peto's meta analysis method
# the needed quantities of observed and expected
Oi = A
Ei = (A+B)*(A+C)/N
Vi = (A+B)*(C+D)/N*(A+C)/N*(B+D)/(N-1)

psii= exp( (Oi-Ei)/Vi)
psii
##  [1] 1.2339585 8.2247085 2.5891637 1.0308188 9.8982582 1.1187806 2.8311704
##  [8] 0.2519599 6.6245551 7.3890561 0.4990568 3.1504364 2.4919073
# lower and upper CI bounds
lowCIi = exp( (Oi-Ei-1.96*sqrt(Vi))/Vi )
upCIi = exp( (Oi-Ei + 1.96*sqrt(Vi))/Vi )
lowCIi
##  [1] 0.50235740 0.82163331 1.30978611 0.23968447 0.96880474 0.62078883
##  [7] 0.66493245 0.04770032 2.43218821 0.14660696 0.04949417 1.46355600
## [13] 0.61447245
upCIi
##  [1]   3.031016  82.330925   5.118216   4.433276 101.130301   2.016257
##  [7]  12.054647   1.330888  18.043312 372.411714   5.032060   6.781599
## [13]  10.105582
# pooled Peto's OR
psi = exp( sum(Oi-Ei)/sum(Vi))
psi
## [1] 1.974851
# lower and upper CI bounds
lowCI = exp( ( sum(Oi-Ei) - 1.96*sqrt(sum(Vi)))/sum(Vi) )
upCI = exp( ( sum(Oi-Ei) + 1.96*sqrt(sum(Vi)))/sum(Vi) )
lowCI
## [1] 1.468895
upCI
## [1] 2.655082
# call "metabin" for Peto's OR
ORPeto.stone =metabin(event_mini,n_mini,event_rirs,
                      n_rirs,studlab=Study, data=stone, method="Peto", sm="OR")
# print the summary
summary(ORPeto.stone)
## Number of studies combined: k = 13
## 
##                          OR           95%-CI    z  p-value
## Fixed effect model   1.9749 [1.4689; 2.6551] 4.51 < 0.0001
## Random effects model 2.0541 [1.2690; 3.3247] 2.93   0.0034
## 
## Quantifying heterogeneity:
##  tau^2 = 0.3311; H = 1.43 [1.04; 1.96]; I^2 = 50.8% [7.0%; 73.9%]
## 
## Test of heterogeneity:
##      Q d.f.  p-value
##  24.37   12   0.0181
## 
## Details on meta-analytical method:
## - Peto method
## - DerSimonian-Laird estimator for tau^2
forest(ORPeto.stone)