R Markdown
library(metamisc)
## Registered S3 methods overwritten by 'ggplot2':
## method from
## [.quosures rlang
## c.quosures rlang
## print.quosures rlang
######### Validation of prediction models with a binary outcome #########
data(EuroSCORE)
# Calculate the c-statistic and its standard error
est1 <- ccalc(cstat = c.index, cstat.se = se.c.index, cstat.cilb = c.index.95CIl,
cstat.ciub = c.index.95CIu, N = n, O = n.events, data = EuroSCORE, slab = Study)
est1
## theta theta.se theta.cilb theta.ciub theta.source
## Nashef 0.8095 0.01377576 0.7820000 0.8360000 c-statistic
## Biancari 0.8670 0.03520473 0.7980000 0.9360000 c-statistic
## Di Dedda 0.8100 0.03571494 0.7400000 0.8800000 c-statistic
## Chalmers 0.7900 0.01000000 0.7704004 0.8095996 c-statistic
## Grant 0.8080 0.00800000 0.7923203 0.8236797 c-statistic
## Carneo 0.8500 0.01000000 0.8304004 0.8695996 c-statistic
## Kunt 0.7200 0.05100000 0.6200418 0.8199582 c-statistic
## Kirmani 0.8180 0.00700000 0.8042803 0.8317197 c-statistic
## Howell 0.6700 0.02972796 0.6117343 0.7282657 c-statistic
## Wang 0.7200 0.01500000 0.6906005 0.7493995 c-statistic
## Borde 0.7200 0.09044574 0.5427296 0.8972704 c-statistic
## Qadir 0.8400 0.02338189 0.7941723 0.8858277 c-statistic
## Spiliopoulos 0.7700 0.06700000 0.6386824 0.9013176 c-statistic
## Wendt 0.7200 0.03400000 0.6533612 0.7866388 c-statistic
## Laurent 0.7700 0.06100000 0.6504422 0.8895578 c-statistic
## Wang.1 0.6420 0.07100000 0.5028426 0.7811574 c-statistic
## Nishida 0.7697 0.04247895 0.6864428 0.8529572 c-statistic
## Barilli 0.8000 0.01500000 0.7706005 0.8293995 c-statistic
## Barilli.1 0.8200 0.02000000 0.7808007 0.8591993 c-statistic
## Paparella 0.8300 0.01200000 0.8064804 0.8535196 c-statistic
## Carosella 0.7600 0.05600000 0.6502420 0.8697580 c-statistic
## Borracci 0.8560 0.03300000 0.7913212 0.9206788 c-statistic
## Osnabrugge 0.7700 0.01000000 0.7504004 0.7895996 c-statistic
## theta.se.source
## Nashef Confidence Interval
## Biancari Confidence Interval
## Di Dedda Confidence Interval
## Chalmers Standard Error
## Grant Standard Error
## Carneo Standard Error
## Kunt Standard Error
## Kirmani Standard Error
## Howell Newcombe (Method 4)
## Wang Standard Error
## Borde Newcombe (Method 4)
## Qadir Newcombe (Method 4)
## Spiliopoulos Standard Error
## Wendt Standard Error
## Laurent Standard Error
## Wang.1 Standard Error
## Nishida Newcombe (Method 4)
## Barilli Standard Error
## Barilli.1 Standard Error
## Paparella Standard Error
## Carosella Standard Error
## Borracci Standard Error
## Osnabrugge Standard Error
# Calculate the logit c-statistic and its standard error
est2 <- ccalc(cstat = c.index, cstat.se = se.c.index, cstat.cilb = c.index.95CIl,
cstat.ciub = c.index.95CIu, N = n, O = n.events, data = EuroSCORE, slab = Study,
g = "log(cstat/(1-cstat))")
est2
## theta theta.se theta.cilb theta.ciub theta.source
## Nashef 1.4467646 0.08964514 1.27735968 1.6287622 c-statistic
## Biancari 1.8746898 0.33390703 1.37384090 2.6827324 c-statistic
## Di Dedda 1.4500102 0.24144872 1.04596856 1.9924302 c-statistic
## Chalmers 1.3249254 0.06027728 1.20678413 1.4430667 c-statistic
## Grant 1.4370667 0.05156766 1.33599594 1.5381374 c-statistic
## Carneo 1.7346011 0.07843137 1.58087839 1.8883237 c-statistic
## Kunt 0.9444616 0.25297619 0.44863739 1.4402858 c-statistic
## Kirmani 1.5028556 0.04701900 1.41070011 1.5950112 c-statistic
## Howell 0.7081851 0.13445481 0.44465847 0.9717116 c-statistic
## Wang 0.9444616 0.07440476 0.79863096 1.0902923 c-statistic
## Borde 0.9444616 0.44863958 0.06514418 1.8237790 c-statistic
## Qadir 1.6582281 0.17397241 1.31724841 1.9992077 c-statistic
## Spiliopoulos 1.2083112 0.37831733 0.46682285 1.9497996 c-statistic
## Wendt 0.9444616 0.16865079 0.61391213 1.2750111 c-statistic
## Laurent 1.2083112 0.34443817 0.53322480 1.8833976 c-statistic
## Wang.1 0.5840553 0.30891592 -0.02140877 1.1895194 c-statistic
## Nishida 1.2066180 0.23963947 0.73693330 1.6763027 c-statistic
## Barilli 1.3862944 0.09375000 1.20254774 1.5700410 c-statistic
## Barilli.1 1.5163475 0.13550136 1.25076971 1.7819253 c-statistic
## Paparella 1.5856273 0.08504607 1.41894004 1.7523145 c-statistic
## Carosella 1.1526795 0.30701754 0.55093618 1.7544228 c-statistic
## Borracci 1.7824571 0.26771807 1.25773930 2.3071748 c-statistic
## Osnabrugge 1.2083112 0.05646527 1.09764130 1.3189811 c-statistic
## theta.se.source
## Nashef Confidence Interval
## Biancari Confidence Interval
## Di Dedda Confidence Interval
## Chalmers Standard Error
## Grant Standard Error
## Carneo Standard Error
## Kunt Standard Error
## Kirmani Standard Error
## Howell Newcombe (Method 4)
## Wang Standard Error
## Borde Newcombe (Method 4)
## Qadir Newcombe (Method 4)
## Spiliopoulos Standard Error
## Wendt Standard Error
## Laurent Standard Error
## Wang.1 Standard Error
## Nishida Newcombe (Method 4)
## Barilli Standard Error
## Barilli.1 Standard Error
## Paparella Standard Error
## Carosella Standard Error
## Borracci Standard Error
## Osnabrugge Standard Error
# Display the results of all studies in a forest plot
plot(est1)
data(Collins)
data("Daniels")
data(DVTipd)
str(DVTipd)
## 'data.frame': 500 obs. of 16 variables:
## $ sex : num 0 1 0 1 0 0 1 0 1 0 ...
## $ malign : num 0 0 0 0 0 0 0 0 0 0 ...
## $ par : num 0 0 1 0 0 0 0 0 0 0 ...
## $ surg : num 0 0 0 0 0 0 0 0 1 0 ...
## $ tend : num 1 1 0 1 1 0 0 1 1 1 ...
## $ oachst : num 0 0 0 0 0 0 0 0 0 0 ...
## $ leg : num 1 0 0 0 0 1 1 0 0 0 ...
## $ notraum : num 1 1 1 1 1 0 0 1 0 1 ...
## $ calfdif3: num 0 0 0 0 0 0 0 0 0 0 ...
## $ pit : num 0 0 0 0 0 1 0 1 1 1 ...
## $ vein : num 0 0 0 0 1 0 0 0 0 1 ...
## $ altdiagn: num 1 0 1 1 1 0 1 1 1 1 ...
## $ histdvt : num 0 1 0 0 0 0 1 0 0 0 ...
## $ ddimdich: num 1 0 0 0 0 1 1 0 1 1 ...
## $ dvt : num 0 0 0 0 0 0 0 0 0 0 ...
## $ study : Factor w/ 4 levels "a","b","c","d": 1 4 1 4 1 4 4 4 4 2 ...
summary(apply(DVTipd,2,as.factor))
## sex malign par surg tend oachst leg notraum calfdif3
## 0:299 0:460 0:446 0:446 0:192 0:472 0:313 0:103 0:315
## 1:201 1: 40 1: 54 1: 54 1:308 1: 28 1:187 1:397 1:185
##
##
## pit vein altdiagn histdvt ddimdich dvt study
## 0:185 0:411 0:226 0:417 0:190 0:418 a:146
## 1:315 1: 89 1:274 1: 83 1:310 1: 82 b:102
## c:116
## d:136
## Develop a prediction model to predict presence of DVT
model.dvt <- glm("dvt~sex+oachst+malign+surg+notraum+vein+calfdif3+ddimdich", family=binomial, data=DVTipd)
summary(model.dvt)
##
## Call:
## glm(formula = "dvt~sex+oachst+malign+surg+notraum+vein+calfdif3+ddimdich",
## family = binomial, data = DVTipd)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.5424 -0.5687 -0.2874 -0.1260 2.7104
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.1664 0.6365 -8.117 4.76e-16 ***
## sex 0.8146 0.2825 2.883 0.00393 **
## oachst 0.4324 0.6227 0.694 0.48739
## malign 0.5679 0.4025 1.411 0.15826
## surg 0.1002 0.4111 0.244 0.80734
## notraum 0.3351 0.3700 0.906 0.36513
## vein 0.4831 0.3186 1.516 0.12939
## calfdif3 1.1841 0.2819 4.200 2.67e-05 ***
## ddimdich 2.6081 0.5310 4.911 9.04e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 446.24 on 499 degrees of freedom
## Residual deviance: 345.98 on 491 degrees of freedom
## AIC: 363.98
##
## Number of Fisher Scoring iterations: 6
data(DVTmodels)
data("EuroSCORE")
data(Fibrinogen)
b <- log(Fibrinogen$HR)
b.se <- ((log(Fibrinogen$HR.975) - log(Fibrinogen$HR.025))/(2*qnorm(0.975)))
n.total <- Fibrinogen$N.total
d.total <- Fibrinogen$N.events
fat(b=b, b.se=b.se)
## Call: fat(b = b, b.se = b.se)
##
## Fixed effect summary estimate: 0.4186
##
## test for funnel plot asymmetry: t =1.9021, df = 29, p = 0.0671
fat(b=b, b.se=b.se, d.total=d.total, method="D-FIV")
## Call: fat(b = b, b.se = b.se, d.total = d.total, method = "D-FIV")
##
## Fixed effect summary estimate: 0.4186
##
## test for funnel plot asymmetry: t =1.6847, df = 29, p = 0.1028
# Note that many tests are also available via metafor
require(metafor)
## Loading required package: metafor
## Loading required package: Matrix
## Loading 'metafor' package (version 2.1-0). For an overview
## and introduction to the package please type: help(metafor).
##
## Attaching package: 'metafor'
## The following object is masked from 'package:metamisc':
##
## forest
fat(b=b, b.se=b.se, n.total=n.total, method="M-FIV")
## Call: fat(b = b, b.se = b.se, n.total = n.total, method = "M-FIV")
##
## Fixed effect summary estimate: 0.4186
##
## test for funnel plot asymmetry: t =-1.4275, df = 29, p = 0.1641
regtest(x=b, sei=b.se, ni=n.total, model="lm", predictor="ni")
##
## Regression Test for Funnel Plot Asymmetry
##
## model: weighted regression with multiplicative dispersion
## predictor: sample size
##
## test for funnel plot asymmetry: t = -1.4275, df = 29, p = 0.1641
str(Fibrinogen)
## 'data.frame': 31 obs. of 5 variables:
## $ N.total : num 3134 615 838 1074 7271 ...
## $ N.events: num 17 23 25 25 25 29 31 38 61 100 ...
## $ HR : num 2.62 2.19 1.38 1.82 1.49 1.62 1.92 3.04 1.71 1.26 ...
## $ HR.025 : num 1.39 1.36 0.78 1.03 0.88 0.87 1.17 1.14 1.22 0.95 ...
## $ HR.975 : num 4.93 3.52 2.44 3.21 2.54 3 3.16 8.1 2.4 1.67 ...
plot(Fibrinogen)
data(DVTipd)
# Internal-external cross-validation of a pre-specified model 'f'
f <- dvt ~ histdvt + ddimdich + sex + notraum
fit <- metapred(DVTipd, strata = "study", formula = f, scope = f, family = binomial)
fit
## Call: metapred(data = DVTipd, strata = "study", formula = f, scope = f,
## family = binomial)
##
## Started with model:
## dvt ~ histdvt + ddimdich + sex + notraum
##
## Generalizability:
## unchanged
## 1 0.1493671
##
## Cross-validation stopped after 0 steps, as Final model:
## Meta-analytic model of prediction models estimated in 4 strata. Coefficients:
## (Intercept) histdvt ddimdich sex notraum
## -4.1180636 0.6174010 1.4400776 0.9647970 0.3658235
data("Framingham")
data("Kertai")
data(Daniels)
fit <- riley(Daniels,control=list(maxit=10000))
logLik(fit)
## 'log Lik.' -48.85119 (df=5)
data(DVTipd)
# Scope
f <- dvt ~ histdvt + ddimdich + sex + notraum
# Internal-external cross-validation of a pre-specified model 'f'
fit <- metapred(DVTipd, strata = "study", formula = f, scope = f, family = binomial)
fit
## Call: metapred(data = DVTipd, strata = "study", formula = f, scope = f,
## family = binomial)
##
## Started with model:
## dvt ~ histdvt + ddimdich + sex + notraum
##
## Generalizability:
## unchanged
## 1 0.1493671
##
## Cross-validation stopped after 0 steps, as Final model:
## Meta-analytic model of prediction models estimated in 4 strata. Coefficients:
## (Intercept) histdvt ddimdich sex notraum
## -4.1180636 0.6174010 1.4400776 0.9647970 0.3658235
# Let's try to simplify model 'f' in order to improve its external validity
metapred(DVTipd, strata = "study", formula = f, family = binomial)
## Call: metapred(data = DVTipd, strata = "study", formula = f, family = binomial)
##
## Started with model:
## dvt ~ histdvt + ddimdich + sex + notraum
##
## Generalizability:
## unchanged
## 1 0.1493671
##
## Generalizability:
## ddimdich histdvt notraum sex
## 1 0.1407603 0.1450895 0.1363067 0.1496112
##
## Continued with model:
## dvt ~ histdvt + ddimdich + sex
##
## Generalizability:
## ddimdich histdvt sex
## 1 0.1415329 0.1342381 0.1374618
##
## Continued with model:
## dvt ~ ddimdich + sex
##
## Generalizability:
## ddimdich sex
## 1 0.1402638 0.1357657
##
## Cross-validation stopped after 3 steps, as no improvement was possible. Final model:
## Meta-analytic model of prediction models estimated in 4 strata. Coefficients:
## (Intercept) ddimdich sex
## -3.6187987 1.5175573 0.8784071
# We can also try to build a generalizable model from scratch
# Some additional examples:
metapred(DVTipd, strata = "study", formula = dvt ~ 1, scope = f, family = binomial) # Forwards
## Warning in metafor::rma(coefficients[, col], variances[, col], method =
## method, : Ratio of largest to smallest sampling variance extremely large.
## May not be able to obtain stable results.
## Warning in metafor::rma(coefficients[, col], variances[, col], method =
## method, : Ratio of largest to smallest sampling variance extremely large.
## May not be able to obtain stable results.
## Call: metapred(data = DVTipd, strata = "study", formula = dvt ~ 1,
## scope = f, family = binomial)
##
## Started with model:
## dvt ~ 1
##
## Generalizability:
## unchanged
## 1 0.1435791
##
## Generalizability:
## ddimdich histdvt notraum sex
## 1 0.1357657 0.1443994 0.1425857 0.1402638
##
## Continued with model:
## dvt ~ ddimdich
##
## Generalizability:
## histdvt notraum sex
## 1 0.1374618 0.1454819 0.1342381
##
## Continued with model:
## dvt ~ ddimdich + sex
##
## Generalizability:
## histdvt notraum
## 1 0.1363067 0.1450895
##
## Cross-validation stopped after 3 steps, as no improvement was possible. Final model:
## Meta-analytic model of prediction models estimated in 4 strata. Coefficients:
## (Intercept) ddimdich sex
## -3.6187987 1.5175573 0.8784071
metapred(DVTipd, strata = "study", formula = f, scope = f, family = binomial) # no selection
## Call: metapred(data = DVTipd, strata = "study", formula = f, scope = f,
## family = binomial)
##
## Started with model:
## dvt ~ histdvt + ddimdich + sex + notraum
##
## Generalizability:
## unchanged
## 1 0.1493671
##
## Cross-validation stopped after 0 steps, as Final model:
## Meta-analytic model of prediction models estimated in 4 strata. Coefficients:
## (Intercept) histdvt ddimdich sex notraum
## -4.1180636 0.6174010 1.4400776 0.9647970 0.3658235
metapred(DVTipd, strata = "study", formula = f, max.steps = 0, family = binomial) # no selection
## Call: metapred(data = DVTipd, strata = "study", formula = f, max.steps = 0,
## family = binomial)
##
## Started with model:
## dvt ~ histdvt + ddimdich + sex + notraum
##
## Generalizability:
## unchanged
## 1 0.1493671
##
## Cross-validation stopped after 0 steps, as max.steps was reached. Final model:
## Meta-analytic model of prediction models estimated in 4 strata. Coefficients:
## (Intercept) histdvt ddimdich sex notraum
## -4.1180636 0.6174010 1.4400776 0.9647970 0.3658235
metapred(DVTipd, strata = "study", formula = f, recal.int = TRUE, family = binomial)
## Call: metapred(data = DVTipd, strata = "study", formula = f, recal.int = TRUE,
## family = binomial)
##
## Started with model:
## dvt ~ histdvt + ddimdich + sex + notraum
##
## Generalizability:
## unchanged
## 1 0.1260719
##
## Generalizability:
## ddimdich histdvt notraum sex
## 1 0.1370616 0.1251329 0.1264966 0.1304285
##
## Continued with model:
## dvt ~ ddimdich + sex + notraum
##
## Generalizability:
## ddimdich notraum sex
## 1 0.1362763 0.1254852 0.1296413
##
## Cross-validation stopped after 2 steps, as no improvement was possible. Final model:
## Meta-analytic model of prediction models estimated in 4 strata. Coefficients:
## (Intercept) ddimdich sex notraum
## -4.2145310 1.4298534 0.9124567 0.4564684
metapred(DVTipd, strata = "study", formula = f, meta.method = "REML", family = binomial)
## Call: metapred(data = DVTipd, strata = "study", formula = f, meta.method = "REML",
## family = binomial)
##
## Started with model:
## dvt ~ histdvt + ddimdich + sex + notraum
##
## Generalizability:
## unchanged
## 1 0.1482474
##
## Generalizability:
## ddimdich histdvt notraum sex
## 1 0.1408751 0.1432147 0.1346666 0.1483425
##
## Continued with model:
## dvt ~ histdvt + ddimdich + sex
##
## Generalizability:
## ddimdich histdvt sex
## 1 0.1414789 0.1318535 0.1360003
##
## Continued with model:
## dvt ~ ddimdich + sex
##
## Generalizability:
## ddimdich sex
## 1 0.1402132 0.1337662
##
## Cross-validation stopped after 3 steps, as no improvement was possible. Final model:
## Meta-analytic model of prediction models estimated in 4 strata. Coefficients:
## (Intercept) ddimdich sex
## -3.6187987 1.7130967 0.8784071
# By default, metapred assumes the first column is the outcome.
newdat <- data.frame(dvt=0, histdvt=0, ddimdich=0, sex=1, notraum=0)
fitted <- predict(fit, newdata = newdat)
######### Validation of prediction models with a binary outcome #########
data(EuroSCORE)
# Calculate the total O:E ratio and its standard error
est1 <- oecalc(O = n.events, E = e.events, N = n, data = EuroSCORE, slab = Study)
est1
## theta theta.se theta.cilb theta.ciub theta.source O
## Nashef 1.0450450 0.06716203 0.9134099 1.1766802 O, E and N 232
## Biancari 0.6086957 0.11345371 0.3863305 0.8310608 O, E and N 28
## Di Dedda 1.2058824 0.18475130 0.8437765 1.5679882 O, E and N 41
## Chalmers 0.7318008 0.05203645 0.6298112 0.8337903 O, E and N 191
## Grant 0.9214541 0.03320253 0.8563783 0.9865298 O, E and N 746
## Carneo 1.2573099 0.08328543 1.0940735 1.4205464 O, E and N 215
## Kunt 4.8571429 0.79922240 3.2906957 6.4235900 O, E and N 34
## Kirmani 1.4134367 0.05935803 1.2970971 1.5297763 O, E and N 547
## Howell 0.8571429 0.08588255 0.6888162 1.0254696 O, E and N 90
## Wang 0.7793103 0.05131185 0.6787410 0.8798797 O, E and N 226
## Borde 0.8000000 0.28056169 0.2501092 1.3498908 O, E and N 8
## Qadir 1.0270270 0.11555260 0.8005481 1.2535060 O, E and N 76
## Spiliopoulos 1.5555556 0.40204098 0.7675697 2.3435414 O, E and N 14
## Wendt 1.3235294 0.19309081 0.9450784 1.7019804 O, E and N 45
## Laurent 2.5714286 0.58846311 1.4180621 3.7247951 O, E and N 18
## Wang.1 0.6190476 0.17032315 0.2852204 0.9528749 O, E and N 13
## Nishida 0.9705882 0.16279815 0.6515097 1.2896668 O, E and N 33
## Barilli 0.6885246 0.04710205 0.5962063 0.7808429 O, E and N 210
## Barilli.1 1.2019231 0.10340170 0.9992595 1.4045867 O, E and N 125
## Paparella 1.1029412 0.06211634 0.9811954 1.2246870 O, E and N 300
## Carosella 2.2500000 0.73637626 0.8067290 3.6932710 O, E and N 9
## Borracci 1.3125000 0.28036848 0.7629879 1.8620121 O, E and N 21
## Osnabrugge 0.6830357 0.02064915 0.6425641 0.7235073 O, E and N 1071
## E N
## Nashef 222.00 5553
## Biancari 46.00 1027
## Di Dedda 34.00 1090
## Chalmers 261.00 5576
## Grant 809.59 23740
## Carneo 171.00 3798
## Kunt 7.00 428
## Kirmani 387.00 15497
## Howell 105.00 933
## Wang 290.00 11170
## Borde 10.00 498
## Qadir 74.00 2004
## Spiliopoulos 9.00 216
## Wendt 34.00 1066
## Laurent 7.00 314
## Wang.1 21.00 818
## Nishida 34.00 461
## Barilli 305.00 12201
## Barilli.1 104.00 1670
## Paparella 272.00 6191
## Carosella 4.00 250
## Borracci 16.00 503
## Osnabrugge 1568.00 50588
# Calculate the log of the total O:E ratio and its standard error
est2 <- oecalc(O = n.events, E = e.events, N = n, data = EuroSCORE, slab = Study, g = "log(OE)")
est2
## theta theta.se theta.cilb theta.ciub theta.source
## Nashef 0.04405999 0.06426711 -0.081901240 0.17002122 O, E and N
## Biancari -0.49643689 0.18638824 -0.861751123 -0.13112265 O, E and N
## Di Dedda 0.18721154 0.15320840 -0.113071397 0.48749448 O, E and N
## Chalmers -0.31224698 0.07110740 -0.451614919 -0.17287904 O, E and N
## Grant -0.08180235 0.03603276 -0.152425252 -0.01117944 O, E and N
## Carneo 0.22897447 0.06624097 0.099144553 0.35880439 O, E and N
## Kunt 1.58045038 0.16454579 1.257946559 1.90295419 O, E and N
## Kirmani 0.34602411 0.04199553 0.263714374 0.42833385 O, E and N
## Howell -0.15415068 0.10019631 -0.350531831 0.04223047 O, E and N
## Wang -0.24934592 0.06584264 -0.378395127 -0.12029672 O, E and N
## Borde -0.22314355 0.35070211 -0.910507050 0.46421995 O, E and N
## Qadir 0.02666825 0.11251174 -0.193850721 0.24718721 O, E and N
## Spiliopoulos 0.44183275 0.25845491 -0.064729568 0.94839507 O, E and N
## Wendt 0.28030197 0.14589084 -0.005638818 0.56624275 O, E and N
## Laurent 0.94446161 0.22884677 0.495930190 1.39299303 O, E and N
## Wang.1 -0.47957308 0.27513739 -1.018832454 0.05968629 O, E and N
## Nishida -0.02985296 0.16773143 -0.358600527 0.29889460 O, E and N
## Barilli -0.37320425 0.06841012 -0.507285614 -0.23912288 O, E and N
## Barilli.1 0.18392284 0.08603021 0.015306718 0.35253896 O, E and N
## Paparella 0.09798041 0.05631881 -0.012402434 0.20836325 O, E and N
## Carosella 0.81093022 0.32727834 0.169476459 1.45238397 O, E and N
## Borracci 0.27193372 0.21361408 -0.146742192 0.69060962 O, E and N
## Osnabrugge -0.38120813 0.03023143 -0.440460642 -0.32195562 O, E and N
## O E N
## Nashef 232 222.00 5553
## Biancari 28 46.00 1027
## Di Dedda 41 34.00 1090
## Chalmers 191 261.00 5576
## Grant 746 809.59 23740
## Carneo 215 171.00 3798
## Kunt 34 7.00 428
## Kirmani 547 387.00 15497
## Howell 90 105.00 933
## Wang 226 290.00 11170
## Borde 8 10.00 498
## Qadir 76 74.00 2004
## Spiliopoulos 14 9.00 216
## Wendt 45 34.00 1066
## Laurent 18 7.00 314
## Wang.1 13 21.00 818
## Nishida 33 34.00 461
## Barilli 210 305.00 12201
## Barilli.1 125 104.00 1670
## Paparella 300 272.00 6191
## Carosella 9 4.00 250
## Borracci 21 16.00 503
## Osnabrugge 1071 1568.00 50588
# Display the results of all studies in a forest plot
plot(est1)
data(Fibrinogen)
b <- log(Fibrinogen$HR)
b.se <- ((log(Fibrinogen$HR.975) - log(Fibrinogen$HR.025))/(2*qnorm(0.975)))
n.total <- Fibrinogen$N.total
# A very simple funnel plot
plot(fat(b=b, b.se=b.se))
# Add custom tickmarks for the X-axis
plot(fat(b=b, b.se=b.se, n.total=n.total, method="M-FIV"), xlab="Hazard ratio", xaxt="n")
axis(1, at=c(log(0.5), log(1), log(1.5), log(2), log(3)), labels=c(0.5, 1, 1.5, 2,3))
data(EuroSCORE)
# Calculate the c-statistic and its standard error
est1 <- ccalc(cstat = c.index, cstat.se = se.c.index, cstat.cilb = c.index.95CIl,
cstat.ciub = c.index.95CIu, N = n, O = n.events, data = EuroSCORE, slab = Study)
plot(est1)
# Calculate the total O:E ratio and its standard error
est2 <- oecalc(O = n.events, E = e.events, N = n, data = EuroSCORE, slab = Study)
plot(est2)
data(Scheidler)
#Perform the analysis
fit <- riley(Scheidler[which(Scheidler$modality==1),])
summary(fit)
## $call
## riley.default(X = Scheidler[which(Scheidler$modality == 1), ])
##
## $confints
## Estimate SE 2.5 % 97.5 %
## beta1 -0.01731291 0.23398704 -0.475919069 0.4412933
## beta2 -2.32166611 0.17377316 -2.662255249 -1.9810770
## psi1 0.71181410 0.21932344 0.281948058 1.1416801
## psi2 0.38103153 0.19690972 -0.004904432 0.7669675
## rho 0.70119871 0.17024870 0.405271507 0.8641157
## Sens 0.49567188 0.05849238 0.383216240 0.6085671
## FPR 0.08934441 0.01413853 0.065237670 0.1212041
##
## attr(,"class")
## [1] "summary.riley"
plot(fit)
require(ggplot2)
## Loading required package: ggplot2
plot(fit, sort="none", theme=theme_gray())
data(Roberts)
# Frequentist random-effects meta-analysis
fit <- with(Roberts, uvmeta(r=SDM, r.se=SE, labels=rownames(Roberts)))
plot(fit)
data(EuroSCORE)
fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl,
cstat.ciub=c.index.95CIu, N=n, O=n.events, data=EuroSCORE)
fit
## Summary c-statistic with 95% confidence and (approximate) 95% prediction interval:
##
## Estimate CIl CIu PIl PIu
## 0.7888603 0.7648784 0.8110005 0.6792568 0.8682736
##
## Number of studies included: 23
plot(fit)
library(ggplot2)
plot(fit, theme=theme_grey())
data(DVTipd)
DVTipd$cluster <- 1:4 # Add a fictional clustering to the data set.
# Suppose we estimated the model in three studies:
DVTipd123 <- DVTipd[DVTipd$cluster <= 3, ]
mp <- metamisc:::metapred(DVTipd123, strata = "cluster", f = dvt ~ vein + malign,
family = binomial)
# and now want to recalibrate it for the fourth:
DVTipd4 <- DVTipd[DVTipd$cluster == 4, ]
metamisc:::recalibrate(mp, newdata = DVTipd4)
## Call: metamisc:::recalibrate(object = mp, newdata = DVTipd4)
##
## Started with model:
## dvt ~ vein + malign
##
## Generalizability:
## unchanged
## 1 0.1333288
##
## Generalizability:
## malign vein
## 1 0.1344872 0.1318589
##
## Continued with model:
## dvt ~ malign
##
## Generalizability:
## malign
## 1 0.1344533
##
## Cross-validation stopped after 2 steps, as no improvement was possible. Final model:
## Meta-analytic model of prediction models estimated in 3 strata. Coefficients:
## (Intercept) malign
## -1.775619 1.145096
data(Scheidler)
data(Daniels)
data(Kertai)
# Meta-analysis of potential surrogate markers data
# The results obtained by Riley (2008) were as follows:
# beta1 = -0.042 (SE = 0.063),
# beta2 = 14.072 (SE = 4.871)
# rho = -0.759
fit1 <- riley(Daniels) #maxit reached, try again with more iterations
## Warning in rileyES(X, optimization = optimization, control = control, pars
## = pars.default, : Iteration limit had been reached.
fit1 <- riley(Daniels, control=list(maxit=10000))
summary(fit1)
## $call
## riley.default(X = Daniels, control = list(maxit = 10000))
##
## $confints
## Estimate SE 2.5 % 97.5 %
## beta1 0.005298983 0.06479973 -0.12170615 0.1323041
## beta2 13.505678310 4.99256707 3.72042666 23.2909300
## psi1 0.134785102 0.09190903 -0.04535329 0.3149235
## psi2 18.076027226 4.00992776 10.21671324 25.9353412
## rho -0.748689375 0.15266774 -0.89430728 -0.4596724
##
## attr(,"class")
## [1] "summary.riley"
# Meta-analysis of prognostic test studies
fit2 <- riley(Kertai)
fit2
## Call:
## riley.default(X = Kertai)
##
## Coefficients
## beta1 beta2 psi1 psi2 rho
## 0.8164679 -0.9715821 0.3499043 0.7692122 0.1537878
##
## Degrees of Freedom: 9 Residual
# Meta-analysis of computed tomography data
ds <- Scheidler[which(Scheidler$modality==1),]
fit3 <- riley(ds)
fit3
## Call:
## riley.default(X = ds)
##
## Coefficients
## beta1 beta2 psi1 psi2 rho
## -0.01731291 -2.32166611 0.71181410 0.38103153 0.70119871
##
## Degrees of Freedom: 29 Residual
data(DVTipd)
DVTipd$cluster <- letters[1:4] # Add a fictional clustering to the data.
mp <- metapred(DVTipd, strata = "cluster", formula = dvt ~ histdvt + ddimdich, family = binomial)
subset(mp) # best cross-validated model
## Prediction models estimated in 4 strata. Coefficients:
## (Intercept) ddimdich
## a -3.198673 2.135779
## b -3.555348 2.251292
## c -19.566069 18.540216
## d -3.891820 2.947359
##
## Meta-analytic models, estimated in 4 fold combinations. Coefficients:
## (Intercept) ddimdich
## b, c, d -3.724268 2.600774
## a, c, d -3.432711 2.421694
## a, b, d -3.463478 2.377572
## a, b, c -3.318460 2.176257
##
## Cross-validation at stratum level yields the following performance:
## val.strata estimate se var ci.lb ci.ub
## b, c, d a 0.1285223 0.01986149 0.0003944788 0.08921088 0.1678338
## a, c, d b 0.1293544 0.01729276 0.0002990394 0.09512723 0.1635816
## a, b, d c 0.1123561 0.01717570 0.0002950048 0.07836058 0.1463516
## a, b, c d 0.1297595 0.01948263 0.0003795729 0.09119796 0.1683211
## measure
## b, c, d mse
## a, c, d mse
## a, b, d mse
## a, b, c mse
##
## Generalizability:
## 1
## 0.1249981
subset(mp, select = "global") # Final model fitted on all strata.
## Meta-analytic model of prediction models estimated in 4 strata. Coefficients:
## (Intercept) ddimdich
## -3.463480 2.377574
subset(mp, step = 1) # The best model of step 1
## Prediction models estimated in 4 strata. Coefficients:
## (Intercept) ddimdich
## a -3.198673 2.135779
## b -3.555348 2.251292
## c -19.566069 18.540216
## d -3.891820 2.947359
##
## Meta-analytic models, estimated in 4 fold combinations. Coefficients:
## (Intercept) ddimdich
## b, c, d -3.724268 2.600774
## a, c, d -3.432711 2.421694
## a, b, d -3.463478 2.377572
## a, b, c -3.318460 2.176257
##
## Cross-validation at stratum level yields the following performance:
## val.strata estimate se var ci.lb ci.ub
## b, c, d a 0.1285223 0.01986149 0.0003944788 0.08921088 0.1678338
## a, c, d b 0.1293544 0.01729276 0.0002990394 0.09512723 0.1635816
## a, b, d c 0.1123561 0.01717570 0.0002950048 0.07836058 0.1463516
## a, b, c d 0.1297595 0.01948263 0.0003795729 0.09119796 0.1683211
## measure
## b, c, d mse
## a, c, d mse
## a, b, d mse
## a, b, c mse
##
## Generalizability:
## 1
## 0.1249981
subset(mp, step = 1, model = "histdvt") # The model in which histdvt was removed, in step 1.
## Prediction models estimated in 4 strata. Coefficients:
## (Intercept) ddimdich
## a -3.198673 2.135779
## b -3.555348 2.251292
## c -19.566069 18.540216
## d -3.891820 2.947359
##
## Meta-analytic models, estimated in 4 fold combinations. Coefficients:
## (Intercept) ddimdich
## b, c, d -3.724268 2.600774
## a, c, d -3.432711 2.421694
## a, b, d -3.463478 2.377572
## a, b, c -3.318460 2.176257
##
## Cross-validation at stratum level yields the following performance:
## val.strata estimate se var ci.lb ci.ub
## b, c, d a 0.1285223 0.01986149 0.0003944788 0.08921088 0.1678338
## a, c, d b 0.1293544 0.01729276 0.0002990394 0.09512723 0.1635816
## a, b, d c 0.1123561 0.01717570 0.0002950048 0.07836058 0.1463516
## a, b, c d 0.1297595 0.01948263 0.0003795729 0.09119796 0.1683211
## measure
## b, c, d mse
## a, c, d mse
## a, b, d mse
## a, b, c mse
##
## Generalizability:
## 1
## 0.1249981
data(Roberts)
# Frequentist random-effects meta-analysis
fit1 <- with(Roberts, uvmeta(r=SDM, r.se=SE, labels=rownames(Roberts)))
summary(fit1)
## Call:
## uvmeta.default(r = SDM, r.se = SE, labels = rownames(Roberts))
##
## Random effects summary: 0.36195 (SE: 0.0859)
## Tau squared: 0.01322 (SE: 0.03431)
plot(fit1) #show a forest plot
fit1
## Summary estimate with 95% confidence and (approximate) 95% prediction interval:
##
## Estimate CIl CIu PIl PIu
## 0.36194806 0.17636909 0.54752703 0.04923206 0.67466405
data(Collins)
#Extract effect size and error variance
r <- Collins$logOR
vars <- Collins$SE**2
#Frequentist random-effects meta-analysis
fit1 <- uvmeta(r,vars)
fit1
## Summary estimate with 95% confidence and (approximate) 95% prediction interval:
##
## Estimate CIl CIu PIl PIu
## -0.5391532 -1.1186728 0.0403663 -2.2361189 1.1578124
summary(fit1)
## Call:
## uvmeta.default(r = r, r.se = vars)
##
## Random effects summary: -0.53915 (SE: 0.25131)
## Tau squared: 0.45186 (SE: 0.24831)
######### Validation of prediction models with a binary outcome #########
data(EuroSCORE)
# Meta-analysis of the c-statistic (random effects)
fit <- valmeta(cstat=c.index, cstat.se=se.c.index, cstat.cilb=c.index.95CIl,
cstat.ciub=c.index.95CIu, cstat.cilv=0.95, N=n, O=n.events,
slab=Study, data=EuroSCORE)
plot(fit)
# Nearly identical results when we need to estimate the SE
valmeta(cstat=c.index, N=n, O=n.events, slab=Study, data=EuroSCORE)
## Summary c-statistic with 95% confidence and (approximate) 95% prediction interval:
##
## Estimate CIl CIu PIl PIu
## 0.7889020 0.7650864 0.8109000 0.6818676 0.8669518
##
## Number of studies included: 23
# Two-stage meta-analysis of the total O:E ratio (random effects)
valmeta(measure="OE", O=n.events, E=e.events, N=n, data=EuroSCORE)
## Summary Total O:E ratio with 95% confidence and (approximate) 95% prediction interval:
##
## Estimate CIl CIu PIl PIu
## 1.1075973 0.8998973 1.3632352 0.4295250 2.8561122
##
## Number of studies included: 23
valmeta(measure="OE", O=n.events, E=e.events, data=EuroSCORE)
## Summary Total O:E ratio with 95% confidence and (approximate) 95% prediction interval:
##
## Estimate CIl CIu PIl PIu
## 1.1059784 0.8990028 1.3606056 0.4316383 2.8338269
##
## Number of studies included: 23
valmeta(measure="OE", Po=Po, Pe=Pe, N=n, data=EuroSCORE)
## Summary Total O:E ratio with 95% confidence and (approximate) 95% prediction interval:
##
## Estimate CIl CIu PIl PIu
## 1.1230955 0.9212978 1.3690944 0.4549877 2.7722586
##
## Number of studies included: 23
# One-stage meta-analysis of the total O:E ratio (random effects)
valmeta(measure="OE", O=n.events, E=e.events, data=EuroSCORE, method="ML", pars=list(model.oe="poisson/log"))
## Warning in valmeta(measure = "OE", O = n.events, E = e.events, data =
## EuroSCORE, : The Sidik-Jonkman-Hartung-Knapp correction cannot be applied
## Summary with 95% confidence and (approximate) 95% prediction interval:
##
## Estimate CIl CIu PIl PIu
## 1.0890504 0.9002780 1.3253873 0.4449741 2.6653930
##
## Number of studies included: 23
data(Framingham)
# Meta-analysis of total O:E ratio after 10 years of follow-up
valmeta(measure="OE", Po=Po, Pe=Pe, N=n, data=Framingham)
## Warning in rma(yi = ds$theta, sei = ds$theta.se, data = ds, method =
## method, : Studies with NAs omitted from model fitting.
## Summary Total O:E ratio with 95% confidence and (approximate) 95% prediction interval:
##
## Estimate CIl CIu PIl PIu
## 0.5781061 0.4400900 0.7594053 0.1935434 1.7267794
##
## Number of studies included: 16
library("mada")
## Loading required package: mvtnorm
## Loading required package: ellipse
##
## Attaching package: 'ellipse'
## The following object is masked from 'package:graphics':
##
## pairs
## Loading required package: mvmeta
## This is mvmeta 0.4.11. For an overview type: help('mvmeta-package').
##
## Attaching package: 'mvmeta'
## The following object is masked from 'package:metafor':
##
## blup
##
## Attaching package: 'mada'
## The following object is masked from 'package:metafor':
##
## forest
## The following object is masked from 'package:metamisc':
##
## forest
AuditC6 <- data.frame(TP = c(47, 126, 19, 36, 130, 84),
FN = c(9, 51, 10, 3, 19, 2),
FP = c(101, 272, 12, 78, 211, 68),
TN = c(738, 1543, 192, 276, 959, 89))
AuditC6
## TP FN FP TN
## 1 47 9 101 738
## 2 126 51 272 1543
## 3 19 10 12 192
## 4 36 3 78 276
## 5 130 19 211 959
## 6 84 2 68 89
AuditC6$names <- c("Study 1", "Study 2", "Study 4", "Study 4", "Study 5", "Study 6")
data("AuditC")
tail(AuditC)
## TP FN FP TN
## 9 59 5 55 136
## 10 142 50 571 2788
## 11 137 24 107 358
## 12 57 3 103 437
## 13 34 1 21 56
## 14 152 51 88 264
madad(AuditC)
## Descriptive summary of AuditC with 14 primary studies.
## Confidence level for all calculations set to 95 %
## Using a continuity correction of 0.5 if applicable
##
## Diagnostic accuracies
## sens 2.5% 97.5% spec 2.5% 97.5%
## [1,] 0.833 0.716 0.908 0.879 0.855 0.899
## [2,] 0.711 0.640 0.772 0.850 0.833 0.866
## [3,] 0.650 0.471 0.795 0.939 0.898 0.964
## [4,] 0.912 0.785 0.967 0.779 0.733 0.819
## [5,] 0.870 0.807 0.915 0.819 0.796 0.840
## [6,] 0.971 0.912 0.991 0.566 0.489 0.641
## [7,] 0.993 0.934 0.999 0.790 0.754 0.822
## [8,] 0.999 0.994 1.000 0.480 0.468 0.492
## [9,] 0.915 0.822 0.962 0.711 0.643 0.770
## [10,] 0.738 0.672 0.795 0.830 0.817 0.842
## [11,] 0.849 0.786 0.896 0.769 0.729 0.805
## [12,] 0.943 0.854 0.979 0.809 0.773 0.840
## [13,] 0.958 0.838 0.990 0.724 0.616 0.811
## [14,] 0.748 0.684 0.802 0.749 0.702 0.792
##
## Test for equality of sensitivities:
## X-squared = 272.3603, df = 13, p-value = <2e-16
## Test for equality of specificities:
## X-squared = 2204.8, df = 13, p-value = <2e-16
##
##
## Diagnostic OR and likelihood ratios
## DOR 2.5% 97.5% posLR 2.5% 97.5% negLR 2.5% 97.5%
## [1,] 36.379 17.587 75.251 6.897 5.556 8.561 0.190 0.106 0.339
## [2,] 13.913 9.818 19.717 4.736 4.100 5.470 0.340 0.270 0.429
## [3,] 28.600 11.133 73.469 10.660 5.862 19.384 0.373 0.229 0.608
## [4,] 36.732 11.925 113.144 4.127 3.320 5.129 0.112 0.041 0.306
## [5,] 30.361 18.440 49.986 4.817 4.201 5.523 0.159 0.105 0.240
## [6,] 44.162 12.077 161.486 2.240 1.868 2.687 0.051 0.015 0.173
## [7,] 515.729 31.687 8393.755 4.730 4.009 5.581 0.009 0.001 0.145
## [8,] 1388.854 86.752 22234.711 1.922 1.876 1.968 0.001 0.000 0.022
## [9,] 26.607 10.523 67.271 3.167 2.506 4.001 0.119 0.053 0.266
## [10,] 13.768 9.865 19.216 4.341 3.879 4.857 0.315 0.249 0.400
## [11,] 18.716 11.572 30.271 3.679 3.079 4.397 0.197 0.136 0.284
## [12,] 69.444 23.113 208.648 4.927 4.099 5.922 0.071 0.026 0.196
## [13,] 60.442 10.948 333.681 3.477 2.411 5.014 0.058 0.012 0.277
## [14,] 8.850 5.949 13.165 2.982 2.448 3.632 0.337 0.264 0.430
##
## Correlation of sensitivities and false positive rates:
## rho 2.5 % 97.5 %
## 0.677 0.228 0.888
madad(AuditC, level = 0.80)
## Descriptive summary of AuditC with 14 primary studies.
## Confidence level for all calculations set to 80 %
## Using a continuity correction of 0.5 if applicable
##
## Diagnostic accuracies
## sens 10% 90% spec 10% 90%
## [1,] 0.833 0.761 0.887 0.879 0.864 0.893
## [2,] 0.711 0.665 0.752 0.850 0.839 0.860
## [3,] 0.650 0.533 0.751 0.939 0.914 0.957
## [4,] 0.912 0.838 0.955 0.779 0.749 0.806
## [5,] 0.870 0.831 0.901 0.819 0.805 0.833
## [6,] 0.971 0.938 0.987 0.566 0.516 0.616
## [7,] 0.993 0.964 0.999 0.790 0.767 0.812
## [8,] 0.999 0.997 1.000 0.480 0.472 0.488
## [9,] 0.915 0.860 0.950 0.711 0.667 0.751
## [10,] 0.738 0.696 0.777 0.830 0.821 0.838
## [11,] 0.849 0.809 0.881 0.769 0.743 0.793
## [12,] 0.943 0.892 0.970 0.809 0.786 0.829
## [13,] 0.958 0.892 0.985 0.724 0.655 0.784
## [14,] 0.748 0.707 0.784 0.749 0.719 0.778
##
## Test for equality of sensitivities:
## X-squared = 272.3603, df = 13, p-value = <2e-16
## Test for equality of specificities:
## X-squared = 2204.8, df = 13, p-value = <2e-16
##
##
## Diagnostic OR and likelihood ratios
## DOR 10% 90% posLR 10% 90% negLR 10% 90%
## [1,] 36.379 22.618 58.513 6.897 5.987 7.944 0.190 0.130 0.277
## [2,] 13.913 11.077 17.475 4.736 4.310 5.204 0.340 0.293 0.396
## [3,] 28.600 15.433 53.001 10.660 7.210 15.760 0.373 0.271 0.513
## [4,] 36.732 17.603 76.651 4.127 3.579 4.757 0.112 0.058 0.216
## [5,] 30.361 21.914 42.063 4.817 4.405 5.268 0.159 0.121 0.208
## [6,] 44.162 18.918 103.093 2.240 1.989 2.523 0.051 0.023 0.113
## [7,] 515.729 83.223 3195.961 4.730 4.245 5.270 0.009 0.002 0.056
## [8,] 1388.854 226.547 8514.419 1.922 1.892 1.952 0.001 0.000 0.008
## [9,] 26.607 14.508 48.797 3.167 2.718 3.690 0.119 0.070 0.201
## [10,] 13.768 11.072 17.121 4.341 4.033 4.672 0.315 0.270 0.368
## [11,] 18.716 13.667 25.630 3.679 3.275 4.134 0.197 0.155 0.250
## [12,] 69.444 33.825 142.573 4.927 4.369 5.557 0.071 0.036 0.138
## [13,] 60.442 19.778 184.715 3.477 2.737 4.417 0.058 0.021 0.161
## [14,] 8.850 6.826 11.474 2.982 2.621 3.392 0.337 0.287 0.395
##
## Correlation of sensitivities and false positive rates:
## rho 2.5 % 97.5 %
## 0.677 0.228 0.888
AuditC.d <- madad(AuditC)
AuditC.d$fpr
## $fpr
## [1] 0.12083333 0.15005507 0.06097561 0.22112676 0.18061486 0.43354430
## [7] 0.20988806 0.52006770 0.28906250 0.17008929 0.23068670 0.19131238
## [13] 0.27564103 0.25070822
##
## $fpr.ci
## 2.5% 97.5%
## [1,] 0.10050071 0.1446182
## [2,] 0.13436927 0.1672182
## [3,] 0.03560834 0.1024938
## [4,] 0.18106960 0.2671547
## [5,] 0.15963660 0.2036817
## [6,] 0.35875648 0.5114869
## [7,] 0.17753840 0.2463665
## [8,] 0.50762738 0.5324832
## [9,] 0.22957510 0.3568250
## [10,] 0.15776394 0.1831681
## [11,] 0.19473290 0.2710443
## [12,] 0.16038974 0.2265879
## [13,] 0.18879306 0.3835508
## [14,] 0.20834216 0.2984416
forest(madad(AuditC), type = "sens")
forest(madad(AuditC), type = "spec")
rs <- rowSums(AuditC)
weights <- 4 * rs / max(rs)
crosshair(AuditC, xlim = c(0,0.6), ylim = c(0.4,1), col = 1:14, lwd = weights)
ROCellipse(AuditC, pch = "")
points(fpr(AuditC), sens(AuditC))
(fit.DOR.DSL <- madauni(AuditC))
## Call:
## madauni(x = AuditC)
##
## DOR tau^2
## 26.337 0.311
(fit.DOR.MH <- madauni(AuditC, method = "MH"))
## Call:
## madauni(x = AuditC, method = "MH")
##
## DOR
## 17.93335
summary(fit.DOR.DSL)
## Call:
## madauni(x = AuditC)
##
## Estimates:
## DSL estimate 2.5 % 97.5 %
## DOR 26.337 17.971 38.596
## lnDOR 3.271 2.889 3.653
## tau^2 0.311 0.000 3.787
## tau 0.557 0.000 1.946
##
## Cochran's Q: 19.683 (13 df, p = 0.103)
## Higgins' I^2: 33.955%
forest(fit.DOR.DSL)
(fit.phm.homo <- phm(AuditC, hetero = FALSE))
## Call:
## phm.default(data = AuditC, hetero = FALSE)
##
## Coefficients:
## theta
## 0.004586893
(fit.phm.het <- phm(AuditC))
## Call:
## phm.default(data = AuditC)
##
## Coefficients:
## theta taus_sq
## 0.084631351 0.003706143
summary(fit.phm.homo)
## Call:
## phm.default(data = AuditC, hetero = FALSE)
##
## Estimate 2.5 % 97.5 %
## theta 0.004586893 0.003508507 0.00566528
##
## Log-likelihood: -61.499 on 1 degrees of freedom
## AIC: 125
## BIC: 125.6
##
## Chi-square goodness of fit test (Adjusted Profile Maximum
## Likelihood under homogeneity)
##
## data: x
## Chi-square = 222.47, df = 1, p-value < 2.2e-16
##
##
## AUC 2.5 % 97.5 % pAUC 2.5 % 97.5 %
## 0.995 0.997 0.994 0.994 0.995 0.992
summary(fit.phm.het)
## Call:
## phm.default(data = AuditC)
##
## Estimate 2.5 % 97.5 %
## theta 0.084631351 0.047449859 0.121812844
## taus_sq 0.003706143 -0.001277798 0.008690085
##
## Log-likelihood: 31.121 on 2 degrees of freedom
## AIC: -58.2
## BIC: -57
##
## Chi-square goodness of fit test (Adjusted Profile Maximum
## Likelihood under heterogeneity)
##
## data: x
## Chi-square = 13.726, df = 2, p-value = 0.3185
##
##
## AUC 2.5 % 97.5 % pAUC 2.5 % 97.5 %
## 0.922 0.955 0.891 0.891 0.937 0.848
plot(fit.phm.het, xlim = c(0,0.6), ylim = c(0.4,1))
ROCellipse(AuditC, add = TRUE)
(fit.reitsma <- reitsma(AuditC))
## Call: reitsma.default(data = AuditC)
##
## Fixed-effects coefficients:
## tsens tfpr
## (Intercept) 2.0997 -1.2637
##
## 14 studies, 2 fixed and 3 random-effects parameters
## logLik AIC BIC
## 31.5640 -53.1279 -46.4669
summary(fit.reitsma)
## Call: reitsma.default(data = AuditC)
##
## Bivariate diagnostic random-effects meta-analysis
## Estimation method: REML
##
## Fixed-effects coefficients
## Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
## tsens.(Intercept) 2.100 0.338 6.215 0.000 1.438 2.762
## tfpr.(Intercept) -1.264 0.174 -7.249 0.000 -1.605 -0.922
## sensitivity 0.891 - - - 0.808 0.941
## false pos. rate 0.220 - - - 0.167 0.285
##
## tsens.(Intercept) ***
## tfpr.(Intercept) ***
## sensitivity
## false pos. rate
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Variance components: between-studies Std. Dev and correlation matrix
## Std. Dev tsens tfpr
## tsens 1.175 1.000 .
## tfpr 0.638 0.854 1.000
##
## logLik AIC BIC
## 31.564 -53.128 -46.467
##
## AUC: 0.887
## Partial AUC (restricted to observed FPRs and normalized): 0.861
##
## HSROC parameters
## Theta Lambda beta sigma2theta sigma2alpha
## -0.083 3.262 -0.610 0.695 0.218
plot(fit.reitsma, sroclwd = 2, main = "SROC curve (bivariate model) for AUDIT-C data")
points(fpr(AuditC), sens(AuditC), pch = 2)
legend("bottomright", c("data", "summary estimate"), pch = c(2,1))
legend("bottomleft", c("SROC", "conf. region"), lwd = c(2,1))
data("IAQ")
data("SAQ")
# both datasets contain more than one 2x2-table per study
# reduce (somewhat arbitrarily) to one row per study by using the first coded table only:
IAQ1 <- subset(IAQ, IAQ$result_id == 1)
SAQ1 <- subset(SAQ, SAQ$result_id == 1)
fit.IAQ <- reitsma(IAQ1)
fit.SAQ <- reitsma(SAQ1)
plot(fit.IAQ, xlim = c(0,.5), ylim = c(.5,1), main = "Comparison of IAQ and SAQ")
lines(sroc(fit.SAQ), lty = 2)
ROCellipse(fit.SAQ, lty = 2, pch = 2, add = TRUE)
points(fpr(IAQ1), sens(IAQ1), cex = .5)
points(fpr(SAQ1), sens(SAQ1), pch = 2, cex = 0.5)
legend("bottomright", c("IAQ", "SAQ"), pch = 1:2, lty = 1:2)
summary(fit.IAQ)
## Call: reitsma.default(data = IAQ1)
##
## Bivariate diagnostic random-effects meta-analysis
## Estimation method: REML
##
## Fixed-effects coefficients
## Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
## tsens.(Intercept) 2.791 0.341 8.184 0.000 2.123 3.460
## tfpr.(Intercept) -3.341 0.293 -11.413 0.000 -3.915 -2.768
## sensitivity 0.942 - - - 0.893 0.970
## false pos. rate 0.034 - - - 0.020 0.059
##
## tsens.(Intercept) ***
## tfpr.(Intercept) ***
## sensitivity
## false pos. rate
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Variance components: between-studies Std. Dev and correlation matrix
## Std. Dev tsens tfpr
## tsens 1.022 1.000 .
## tfpr 0.827 0.736 1.000
##
## logLik AIC BIC
## 40.176 -70.351 -65.372
##
## AUC: 0.983
## Partial AUC (restricted to observed FPRs and normalized): 0.953
##
## HSROC parameters
## Theta Lambda beta sigma2theta sigma2alpha
## -0.601 6.225 -0.211 0.734 0.447
summary(fit.SAQ)
## Call: reitsma.default(data = SAQ1)
##
## Bivariate diagnostic random-effects meta-analysis
## Estimation method: REML
##
## Fixed-effects coefficients
## Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
## tsens.(Intercept) 1.682 0.472 3.564 0.000 0.757 2.607
## tfpr.(Intercept) -2.460 0.240 -10.267 0.000 -2.930 -1.990
## sensitivity 0.843 - - - 0.681 0.931
## false pos. rate 0.079 - - - 0.051 0.120
##
## tsens.(Intercept) ***
## tfpr.(Intercept) ***
## sensitivity
## false pos. rate
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Variance components: between-studies Std. Dev and correlation matrix
## Std. Dev tsens tfpr
## tsens 1.797 1.000 .
## tfpr 0.899 0.499 1.000
##
## logLik AIC BIC
## 32.819 -55.637 -48.309
##
## AUC: 0.949
## Partial AUC (restricted to observed FPRs and normalized): 0.869
##
## HSROC parameters
## Theta Lambda beta sigma2theta sigma2alpha
## -1.143 4.667 -0.692 1.212 1.618
data("smoking")
# again reduce to one result per study:
smoking1 <- subset(smoking, smoking$result_id == 1)
summary(smoking1$type)
## IAQ SAQ
## 10 16
fit.smoking.type <- reitsma(smoking1, formula = cbind(tsens, tfpr) ~ type)
summary(fit.smoking.type)
## Call: reitsma.default(data = smoking1, formula = cbind(tsens, tfpr) ~
## type)
##
## Bivariate diagnostic random-effects meta-analysis
## Estimation method: REML
##
## Fixed-effects coefficients
## Estimate Std. Error z Pr(>|z|) 95%ci.lb 95%ci.ub
## tsens.(Intercept) 2.813 0.491 5.735 0.000 1.852 3.775
## tsens.typeSAQ -1.166 0.634 -1.838 0.066 -2.409 0.077
## tfpr.(Intercept) -3.337 0.311 -10.733 0.000 -3.946 -2.727
## tfpr.typeSAQ 0.882 0.389 2.269 0.023 0.120 1.645
##
## tsens.(Intercept) ***
## tsens.typeSAQ .
## tfpr.(Intercept) ***
## tfpr.typeSAQ *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Variance components: between-studies Std. Dev and correlation matrix
## Std. Dev tsens tfpr
## tsens 1.508 1.000 .
## tfpr 0.875 0.551 1.000
##
## logLik AIC BIC
## 70.721 -127.441 -113.783
fit.smoking.ml.type <- reitsma(smoking1, formula = cbind(tsens, tfpr) ~ type, method = "ml")
fit.smoking.ml.intercept <- reitsma(smoking1, formula = cbind(tsens, tfpr) ~ 1, method = "ml")
anova(fit.smoking.ml.type, fit.smoking.ml.intercept)
## Likelihood-ratio test
## Model 1: cbind(tsens, tfpr) ~ type
## Model 2: cbind(tsens, tfpr) ~ 1
##
## ChiSquared Df Pr(>ChiSquared)
## 13.25 2 0.00133 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit.smoking1 <- reitsma(smoking1, method = "ml")
fit.smoking2 <- reitsma(smoking1, alphasens = 0, alphafpr = 2, method = "ml")
AIC(fit.smoking1)
## [1] -120.0473
AIC(fit.smoking2)
## [1] -120.1002
set.seed(2017)
library("bamdit")
data("glas")
head(glas)
## tp n1 fp n2 Author cutoff(U/ml) marker
## 1 1 2 15 52 Kirollos <NA> BTA
## 2 17 60 9 70 Johnston <NA> BTA
## 3 8 28 7 34 Murphy <NA> BTA
## 4 19 47 8 30 Landman <NA> BTA
## 5 33 41 27 304 Leyh <NA> BTA
## 6 8 12 12 35 Chong <NA> BTA
plotdata(glas, group = glas$marker, max.size = 20)
glas.t <- glas[glas$marker == "Telomerase", 1:4]
plotdata(glas.t)
glas.m1 <- metadiag(glas.t, re = "normal", re.model = "DS", link = "logit", sd.Fisher.rho = 1.7, nr.burnin = 1000, nr.iterations = 10000, nr.chains = 4, r2jags = TRUE)
## module glm loaded
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 20
## Unobserved stochastic nodes: 28
## Total graph size: 208
##
## Initializing model
summary(glas.m1, digits = 3)
## Inference for Bugs model at "4", fit using jags,
## 4 chains, each with 10000 iterations (first 1000 discarded)
## n.sims = 36000 iterations saved
## mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
## deviance 80.088 5.458 71.291 76.185 79.430 83.314 92.484 1.002 3300
## fp.new 11.947 13.119 0.000 2.000 7.000 18.000 46.000 1.001 36000
## mu.D 3.099 0.511 2.018 2.792 3.117 3.430 4.068 1.003 1400
## mu.S -0.572 0.718 -1.933 -1.045 -0.589 -0.127 0.944 1.002 2000
## rho -0.868 0.139 -0.991 -0.958 -0.913 -0.828 -0.480 1.001 22000
## se.new 0.758 0.126 0.445 0.695 0.778 0.845 0.942 1.001 28000
## se.pool 0.777 0.041 0.691 0.752 0.779 0.804 0.852 1.001 36000
## sigma.D 1.531 0.547 0.774 1.158 1.427 1.793 2.857 1.003 1100
## sigma.S 2.446 0.737 1.399 1.934 2.313 2.806 4.218 1.002 1900
## sp.new 0.761 0.258 0.083 0.640 0.866 0.959 0.998 1.001 36000
## sp.pool 0.849 0.076 0.652 0.815 0.864 0.901 0.950 1.004 1700
## tp.new 37.862 6.951 21.000 34.000 39.000 43.000 48.000 1.001 34000
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 14.9 and DIC = 95.0
## DIC is an estimate of expected predictive error (lower deviance is better).
library("R2jags")
## Loading required package: rjags
## Loading required package: coda
## Linked to JAGS 4.3.0
## Loaded modules: basemod,bugs,dic,glm
##
## Attaching package: 'R2jags'
## The following object is masked from 'package:coda':
##
## traceplot
attach.jags(glas.m1)
cor(se.pool, sp.pool)
## [,1]
## [1,] -0.424367
plot(glas.m1, level = c(0.5, 0.75, 0.95), parametric.smooth = TRUE)
plotsesp(glas.m1)
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 2 rows containing missing values (geom_bar).
bsroc(glas.m1, level = c(0.025, 0.5, 0.975), plot.post.bauc = TRUE, fpr.x = seq(0.01, 0.75, 0.01), lower.auc = 0.01, upper.auc = 0.75, partial.AUC = FALSE)
##
## ------------------------------------------------------------
## These results are based on the following random effects model:
## ------------------------------------------------------------
## Link function: logit
## Random Effects distribution: Bivariate Normal
## Parametrization: Differences and Sums
## Splitting study weights: No
## ------------------------------------------------------------
##
## ------------------------------------------------------------
## Posteriors for the parameters of the Bayesian SROC Curve
## ------------------------------------------------------------
## Mean SD 2.5% 25% 50% 75% 97.5%
## A 3.102 0.507 2.039 2.795 3.119 3.431 4.068
## B -0.547 0.151 -0.839 -0.646 -0.550 -0.453 -0.236
##
## ------------------------------------------------------------
## Summary results for the Bayesian Area Under the Curve (BAUC)
## ------------------------------------------------------------
## Mean SD 2.5% 25% 50% 75% 97.5%
## BAUC 0.628 0.032 0.551 0.611 0.632 0.65 0.678
## ------------------------------------------------------------
library("ggplot2")
library("GGally")
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
library("R2jags")
attach.jags(glas.m1)
hyper.post <- data.frame(mu.D, mu.S, sigma.D, sigma.S, rho)
ggpairs(hyper.post, title = "Hyper-Posteriors", lower = list(continuous = "density"))
glas.m2 <- metadiag(glas.t, re = "sm", link = "logit", sd.Fisher.rho = 1.7, df.estimate = TRUE, split.w = TRUE, nr.burnin = 10000, nr.iterations = 100000, nr.chains = 1, r2jags = TRUE)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 20
## Unobserved stochastic nodes: 51
## Total graph size: 333
##
## Initializing model
glas.m2
## Inference for Bugs model at "5", fit using jags,
## 1 chains, each with 1e+05 iterations (first 10000 discarded)
## n.sims = 90000 iterations saved
## mean sd 2.5% 25% 50% 75% 97.5%
## deviance 80.4 5.6 71.3 76.4 79.8 83.8 93.2
## df 7.3 4.1 3.1 4.1 5.9 9.4 18.2
## fp.new 14.5 13.6 0.0 4.0 10.0 22.0 48.0
## mu.D 2.6 0.5 1.5 2.3 2.7 3.0 3.5
## mu.S -0.1 0.6 -1.3 -0.5 -0.1 0.4 1.3
## p.w1[1] 0.6 0.5 0.0 0.0 1.0 1.0 1.0
## p.w1[2] 0.6 0.5 0.0 0.0 1.0 1.0 1.0
## p.w1[3] 0.5 0.5 0.0 0.0 1.0 1.0 1.0
## p.w1[4] 0.5 0.5 0.0 0.0 1.0 1.0 1.0
## p.w1[5] 0.6 0.5 0.0 0.0 1.0 1.0 1.0
## p.w1[6] 0.5 0.5 0.0 0.0 1.0 1.0 1.0
## p.w1[7] 0.6 0.5 0.0 0.0 1.0 1.0 1.0
## p.w1[8] 0.6 0.5 0.0 0.0 1.0 1.0 1.0
## p.w1[9] 0.6 0.5 0.0 0.0 1.0 1.0 1.0
## p.w1[10] 0.6 0.5 0.0 0.0 1.0 1.0 1.0
## p.w2[1] 0.6 0.5 0.0 0.0 1.0 1.0 1.0
## p.w2[2] 0.5 0.5 0.0 0.0 0.0 1.0 1.0
## p.w2[3] 0.5 0.5 0.0 0.0 1.0 1.0 1.0
## p.w2[4] 0.6 0.5 0.0 0.0 1.0 1.0 1.0
## p.w2[5] 0.7 0.4 0.0 0.0 1.0 1.0 1.0
## p.w2[6] 0.5 0.5 0.0 0.0 0.0 1.0 1.0
## p.w2[7] 0.7 0.4 0.0 0.0 1.0 1.0 1.0
## p.w2[8] 0.5 0.5 0.0 0.0 1.0 1.0 1.0
## p.w2[9] 0.5 0.5 0.0 0.0 1.0 1.0 1.0
## p.w2[10] 0.7 0.5 0.0 0.0 1.0 1.0 1.0
## rho -0.9 0.1 -1.0 -1.0 -0.9 -0.8 -0.5
## se.new 0.8 0.1 0.4 0.7 0.8 0.8 0.9
## se.pool 0.8 0.0 0.7 0.8 0.8 0.8 0.9
## sigma.D 1.5 0.6 0.7 1.1 1.4 1.8 3.0
## sigma.S 2.3 0.8 1.1 1.7 2.1 2.7 4.1
## sp.new 0.7 0.3 0.1 0.6 0.8 0.9 1.0
## sp.pool 0.8 0.1 0.6 0.7 0.8 0.8 0.9
## tp.new 38.0 7.0 20.0 35.0 39.0 43.0 48.0
## w1[1] 1.5 2.2 0.4 0.8 1.1 1.6 5.0
## w1[2] 1.4 1.7 0.4 0.7 1.1 1.6 4.7
## w1[3] 1.4 1.9 0.4 0.7 1.0 1.6 4.5
## w1[4] 1.4 2.0 0.4 0.7 1.1 1.6 4.7
## w1[5] 1.6 2.1 0.4 0.8 1.1 1.7 5.4
## w1[6] 1.4 1.9 0.4 0.7 1.0 1.5 4.5
## w1[7] 1.5 2.1 0.4 0.8 1.1 1.7 5.2
## w1[8] 1.5 2.1 0.4 0.8 1.1 1.6 5.0
## w1[9] 1.5 2.2 0.4 0.8 1.1 1.6 5.1
## w1[10] 1.5 2.1 0.4 0.8 1.1 1.6 5.2
## w2[1] 1.5 1.9 0.4 0.8 1.1 1.6 4.8
## w2[2] 1.3 1.5 0.4 0.7 1.0 1.4 3.9
## w2[3] 1.3 1.7 0.4 0.7 1.0 1.5 4.2
## w2[4] 1.6 2.1 0.5 0.8 1.2 1.8 5.6
## w2[5] 2.0 3.8 0.5 1.0 1.4 2.2 7.3
## w2[6] 1.3 1.8 0.4 0.7 1.0 1.4 3.9
## w2[7] 2.4 4.8 0.5 1.0 1.5 2.4 9.7
## w2[8] 1.3 1.7 0.4 0.7 1.0 1.5 4.0
## w2[9] 1.3 1.5 0.4 0.7 1.0 1.5 4.1
## w2[10] 1.8 2.8 0.5 0.9 1.3 2.0 6.5
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 15.8 and DIC = 96.3
## DIC is an estimate of expected predictive error (lower deviance is better).
plotw(m = glas.m2)
glas.t[c(5, 7, 10), ]
## tp n1 fp n2
## 38 40 57 1 138
## 40 23 42 0 12
## 43 37 44 22 29
dat.hat <- data.frame(tpr = glas.t[ , 1]/glas.t[ , 2], fpr = glas.t[ , 3]/glas.t[, 4], n = glas.t[ , 2] + glas.t[ , 4])
dat.hat[c(5, 7, 10), ]
## tpr fpr n
## 5 0.7017544 0.007246377 195
## 7 0.5476190 0.000000000 54
## 10 0.8409091 0.758620690 73
plotcompare(m1 = glas.m1, m2 = glas.m2, m1.name = "Binomial + Normal", m2.name = "Binomial + Scale mixtures", level = 0.95)
data("ct")
gr <- with(ct, factor(design, labels = c("Retrospective study", "Prospective study")))
plotdata(ct, group = gr, y.lo = 0.75, x.up = 0.75, alpha.p = 0.5, max.size = 20)
ct.m <- metadiag(ct, re = "sm", link = "logit", df.estimate = TRUE, split.w = TRUE, nr.burnin = 1000, nr.iterations = 10000, nr.chains = 4, r2jags = TRUE)
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 102
## Unobserved stochastic nodes: 215
## Total graph size: 1399
##
## Initializing model
plotw(m = ct.m, group = gr)
m1.ct <- metadiag(ct[ct$design == 1, 1:4])
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 44
## Unobserved stochastic nodes: 52
## Total graph size: 388
##
## Initializing model
m2.ct <- metadiag(ct[ct$design == 2, 1:4])
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 58
## Unobserved stochastic nodes: 66
## Total graph size: 493
##
## Initializing model
plotcompare(m1.ct, m2.ct, m1.name = "Retrospective design", m2.name = "Prospective design", group = gr, limits.x = c(0, 0.75), limits.y = c(0.65, 1))
summary(m1.ct)
## Inference for Bugs model at "7", fit using jags,
## 2 chains, each with 10000 iterations (first 1000 discarded)
## n.sims = 18000 iterations saved
## mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
## deviance 161.515 9.376 144.899 154.977 160.844 167.439 181.622 1.001 18000
## fp.new 5.403 6.948 0.000 1.000 3.000 7.000 26.000 1.002 6800
## mu.D 5.940 0.373 5.210 5.696 5.938 6.184 6.676 1.002 1200
## mu.S 0.447 0.438 -0.425 0.158 0.447 0.728 1.320 1.001 18000
## rho -0.340 0.220 -0.711 -0.502 -0.358 -0.195 0.129 1.001 18000
## se.new 0.939 0.068 0.757 0.925 0.960 0.978 0.995 1.006 18000
## se.pool 0.959 0.010 0.937 0.954 0.960 0.966 0.976 1.002 2100
## sigma.D 1.403 0.358 0.825 1.152 1.359 1.606 2.227 1.001 18000
## sigma.S 1.870 0.412 1.204 1.583 1.823 2.100 2.826 1.001 14000
## sp.new 0.892 0.133 0.488 0.865 0.940 0.975 0.996 1.011 2500
## sp.pool 0.937 0.019 0.893 0.927 0.940 0.950 0.967 1.001 6300
## tp.new 46.947 3.795 37.000 46.000 48.000 49.000 50.000 1.005 18000
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 44.0 and DIC = 205.5
## DIC is an estimate of expected predictive error (lower deviance is better).
summary(m2.ct)
## Inference for Bugs model at "8", fit using jags,
## 2 chains, each with 10000 iterations (first 1000 discarded)
## n.sims = 18000 iterations saved
## mean sd 2.5% 25% 50% 75% 97.5% Rhat n.eff
## deviance 189.820 8.900 174.021 183.559 189.331 195.595 208.399 1.001 18000
## fp.new 3.195 3.500 0.000 1.000 2.000 4.000 13.000 1.005 1700
## mu.D 6.091 0.317 5.500 5.880 6.081 6.297 6.741 1.002 2000
## mu.S 0.001 0.263 -0.520 -0.169 -0.001 0.176 0.518 1.002 1500
## rho -0.293 0.230 -0.699 -0.462 -0.309 -0.139 0.180 1.001 4400
## se.new 0.944 0.044 0.832 0.930 0.955 0.971 0.989 1.001 18000
## se.pool 0.954 0.008 0.936 0.949 0.954 0.960 0.969 1.001 18000
## sigma.D 1.266 0.284 0.801 1.064 1.235 1.431 1.914 1.001 8700
## sigma.S 0.978 0.223 0.608 0.821 0.957 1.111 1.475 1.001 18000
## sp.new 0.936 0.062 0.767 0.919 0.954 0.975 0.993 1.009 2500
## sp.pool 0.954 0.010 0.933 0.948 0.954 0.960 0.971 1.002 1400
## tp.new 47.177 2.704 40.000 46.000 48.000 49.000 50.000 1.001 18000
##
## For each parameter, n.eff is a crude measure of effective sample size,
## and Rhat is the potential scale reduction factor (at convergence, Rhat=1).
##
## DIC info (using the rule, pD = var(deviance)/2)
## pD = 39.6 and DIC = 229.4
## DIC is an estimate of expected predictive error (lower deviance is better).