Calculate the DIC for a model fitted using the nma() function.
Usage
dic(x, penalty = c("pD", "pV"), ...)Arguments
- x
A fitted model object, inheriting class stan_nma
- penalty
The method for estimating the effective number of parameters, used to penalise model fit in the DIC. Either
"pD"(the default), or"pV". For survival likelihoods only"pV"is currently available.- ...
Other arguments (not used)
Value
A nma_dic object.
See also
print.nma_dic() for printing details, plot.nma_dic() for
producing plots of residual deviance contributions.
Examples
## Smoking cessation
# \donttest{
# Run smoking FE NMA example if not already available
if (!exists("smk_fit_FE")) example("example_smk_fe", run.donttest = TRUE)
#>
#> exmp__> # Set up network of smoking cessation data
#> exmp__> head(smoking)
#> studyn trtn trtc r n
#> 1 1 1 No intervention 9 140
#> 2 1 3 Individual counselling 23 140
#> 3 1 4 Group counselling 10 138
#> 4 2 2 Self-help 11 78
#> 5 2 3 Individual counselling 12 85
#> 6 2 4 Group counselling 29 170
#>
#> exmp__> smk_net <- set_agd_arm(smoking,
#> exmp__+ study = studyn,
#> exmp__+ trt = trtc,
#> exmp__+ r = r,
#> exmp__+ n = n,
#> exmp__+ trt_ref = "No intervention")
#>
#> exmp__> # Print details
#> exmp__> smk_net
#> A network with 24 AgD studies (arm-based).
#>
#> ------------------------------------------------------- AgD studies (arm-based) ----
#> Study Treatment arms
#> 1 3: No intervention | Group counselling | Individual counselling
#> 2 3: Group counselling | Individual counselling | Self-help
#> 3 2: No intervention | Individual counselling
#> 4 2: No intervention | Individual counselling
#> 5 2: No intervention | Individual counselling
#> 6 2: No intervention | Individual counselling
#> 7 2: No intervention | Individual counselling
#> 8 2: No intervention | Individual counselling
#> 9 2: No intervention | Individual counselling
#> 10 2: No intervention | Self-help
#> ... plus 14 more studies
#>
#> Outcome type: count
#> ------------------------------------------------------------------------------------
#> Total number of treatments: 4
#> Total number of studies: 24
#> Reference treatment is: No intervention
#> Network is connected
#>
#> exmp__> ## No test:
#> exmp__> # Fitting a fixed effect model
#> exmp__> smk_fit_FE <- nma(smk_net, ## Don't show:
#> exmp__+ refresh = if (interactive()) 200 else 0,
#> exmp__+ ## End(Don't show)
#> exmp__+ trt_effects = "fixed",
#> exmp__+ prior_intercept = normal(scale = 100),
#> exmp__+ prior_trt = normal(scale = 100))
#>
#> exmp__> smk_fit_FE
#> A fixed effects NMA with a binomial likelihood (logit link).
#> Inference for Stan model: binomial_1par.
#> 4 chains, each with iter=2000; warmup=1000; thin=1;
#> post-warmup draws per chain=1000, total post-warmup draws=4000.
#>
#> mean se_mean sd 2.5% 25% 50%
#> d[Group counselling] 0.84 0.00 0.17 0.50 0.72 0.84
#> d[Individual counselling] 0.77 0.00 0.06 0.65 0.73 0.76
#> d[Self-help] 0.23 0.00 0.13 -0.02 0.14 0.23
#> lp__ -5859.37 0.09 3.75 -5867.83 -5861.63 -5859.06
#> 75% 97.5% n_eff Rhat
#> d[Group counselling] 0.95 1.17 2321 1
#> d[Individual counselling] 0.80 0.88 1744 1
#> d[Self-help] 0.31 0.47 2622 1
#> lp__ -5856.74 -5853.06 1894 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Sep 24 08:57:45 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at
#> convergence, Rhat=1).
#>
#> exmp__> ## End(No test)
#> exmp__>
#> exmp__> ## Don't show:
#> exmp__> if (requireNamespace("pkgdown", quietly = TRUE) && pkgdown::in_pkgdown()) {
#> exmp__+ assign("smk_net", smk_net, .GlobalEnv)
#> exmp__+ assign("smk_fit_FE", smk_fit_FE, .GlobalEnv)
#> exmp__+ }
#>
#> exmp__> ## End(Don't show)
#> exmp__>
#> exmp__>
#> exmp__>
# }
# \donttest{
# Run smoking RE NMA example if not already available
if (!exists("smk_fit_RE")) example("example_smk_re", run.donttest = TRUE)
#>
#> exmp__> # Set up network of smoking cessation data
#> exmp__> head(smoking)
#> studyn trtn trtc r n
#> 1 1 1 No intervention 9 140
#> 2 1 3 Individual counselling 23 140
#> 3 1 4 Group counselling 10 138
#> 4 2 2 Self-help 11 78
#> 5 2 3 Individual counselling 12 85
#> 6 2 4 Group counselling 29 170
#>
#> exmp__> smk_net <- set_agd_arm(smoking,
#> exmp__+ study = studyn,
#> exmp__+ trt = trtc,
#> exmp__+ r = r,
#> exmp__+ n = n,
#> exmp__+ trt_ref = "No intervention")
#>
#> exmp__> # Print details
#> exmp__> smk_net
#> A network with 24 AgD studies (arm-based).
#>
#> ------------------------------------------------------- AgD studies (arm-based) ----
#> Study Treatment arms
#> 1 3: No intervention | Group counselling | Individual counselling
#> 2 3: Group counselling | Individual counselling | Self-help
#> 3 2: No intervention | Individual counselling
#> 4 2: No intervention | Individual counselling
#> 5 2: No intervention | Individual counselling
#> 6 2: No intervention | Individual counselling
#> 7 2: No intervention | Individual counselling
#> 8 2: No intervention | Individual counselling
#> 9 2: No intervention | Individual counselling
#> 10 2: No intervention | Self-help
#> ... plus 14 more studies
#>
#> Outcome type: count
#> ------------------------------------------------------------------------------------
#> Total number of treatments: 4
#> Total number of studies: 24
#> Reference treatment is: No intervention
#> Network is connected
#>
#> exmp__> ## No test:
#> exmp__> # Fitting a random effects model
#> exmp__> smk_fit_RE <- nma(smk_net, ## Don't show:
#> exmp__+ refresh = if (interactive()) 200 else 0,
#> exmp__+ ## End(Don't show)
#> exmp__+ trt_effects = "random",
#> exmp__+ prior_intercept = normal(scale = 100),
#> exmp__+ prior_trt = normal(scale = 100),
#> exmp__+ prior_het = normal(scale = 5))
#>
#> exmp__> smk_fit_RE
#> A random effects NMA with a binomial likelihood (logit link).
#> Inference for Stan model: binomial_1par.
#> 4 chains, each with iter=2000; warmup=1000; thin=1;
#> post-warmup draws per chain=1000, total post-warmup draws=4000.
#>
#> mean se_mean sd 2.5% 25% 50%
#> d[Group counselling] 1.11 0.01 0.44 0.28 0.81 1.10
#> d[Individual counselling] 0.85 0.01 0.24 0.38 0.69 0.84
#> d[Self-help] 0.48 0.01 0.40 -0.29 0.23 0.48
#> lp__ -5767.81 0.19 6.30 -5781.15 -5771.84 -5767.54
#> tau 0.84 0.00 0.18 0.55 0.71 0.82
#> 75% 97.5% n_eff Rhat
#> d[Group counselling] 1.39 1.99 2259 1
#> d[Individual counselling] 1.00 1.34 1413 1
#> d[Self-help] 0.74 1.24 1985 1
#> lp__ -5763.33 -5756.53 1114 1
#> tau 0.94 1.26 1343 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Sep 24 08:57:50 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at
#> convergence, Rhat=1).
#>
#> exmp__> ## End(No test)
#> exmp__>
#> exmp__> ## Don't show:
#> exmp__> if (requireNamespace("pkgdown", quietly = TRUE) && pkgdown::in_pkgdown()) {
#> exmp__+ assign("smk_net", smk_net, .GlobalEnv)
#> exmp__+ assign("smk_fit_RE", smk_fit_RE, .GlobalEnv)
#> exmp__+ }
#>
#> exmp__> ## End(Don't show)
#> exmp__>
#> exmp__>
#> exmp__>
# }
# \donttest{
# Compare DIC of FE and RE models
(smk_dic_FE <- dic(smk_fit_FE))
#> Residual deviance: 267.2 (on 50 data points)
#> pD: 27.1
#> DIC: 294.2
(smk_dic_RE <- dic(smk_fit_RE)) # substantially better fit
#> Residual deviance: 53.8 (on 50 data points)
#> pD: 43.7
#> DIC: 97.5
# Plot residual deviance contributions under RE model
plot(smk_dic_RE)
# Check for inconsistency using UME model
# }
# \donttest{
# Run smoking UME NMA example if not already available
if (!exists("smk_fit_RE_UME")) example("example_smk_ume", run.donttest = TRUE)
#>
#> exmp__> # Set up network of smoking cessation data
#> exmp__> head(smoking)
#> studyn trtn trtc r n
#> 1 1 1 No intervention 9 140
#> 2 1 3 Individual counselling 23 140
#> 3 1 4 Group counselling 10 138
#> 4 2 2 Self-help 11 78
#> 5 2 3 Individual counselling 12 85
#> 6 2 4 Group counselling 29 170
#>
#> exmp__> smk_net <- set_agd_arm(smoking,
#> exmp__+ study = studyn,
#> exmp__+ trt = trtc,
#> exmp__+ r = r,
#> exmp__+ n = n,
#> exmp__+ trt_ref = "No intervention")
#>
#> exmp__> # Print details
#> exmp__> smk_net
#> A network with 24 AgD studies (arm-based).
#>
#> ------------------------------------------------------- AgD studies (arm-based) ----
#> Study Treatment arms
#> 1 3: No intervention | Group counselling | Individual counselling
#> 2 3: Group counselling | Individual counselling | Self-help
#> 3 2: No intervention | Individual counselling
#> 4 2: No intervention | Individual counselling
#> 5 2: No intervention | Individual counselling
#> 6 2: No intervention | Individual counselling
#> 7 2: No intervention | Individual counselling
#> 8 2: No intervention | Individual counselling
#> 9 2: No intervention | Individual counselling
#> 10 2: No intervention | Self-help
#> ... plus 14 more studies
#>
#> Outcome type: count
#> ------------------------------------------------------------------------------------
#> Total number of treatments: 4
#> Total number of studies: 24
#> Reference treatment is: No intervention
#> Network is connected
#>
#> exmp__> ## No test:
#> exmp__> # Fitting an unrelated mean effects (inconsistency) model
#> exmp__> smk_fit_RE_UME <- nma(smk_net, ## Don't show:
#> exmp__+ refresh = if (interactive()) 200 else 0,
#> exmp__+ ## End(Don't show)
#> exmp__+ consistency = "ume",
#> exmp__+ trt_effects = "random",
#> exmp__+ prior_intercept = normal(scale = 100),
#> exmp__+ prior_trt = normal(scale = 100),
#> exmp__+ prior_het = normal(scale = 5))
#>
#> exmp__> smk_fit_RE_UME
#> A random effects NMA with a binomial likelihood (logit link).
#> An inconsistency model ('ume') was fitted.
#> Inference for Stan model: binomial_1par.
#> 4 chains, each with iter=2000; warmup=1000; thin=1;
#> post-warmup draws per chain=1000, total post-warmup draws=4000.
#>
#> mean se_mean sd 2.5%
#> d[Group counselling vs. No intervention] 1.12 0.02 0.80 -0.39
#> d[Individual counselling vs. No intervention] 0.91 0.01 0.27 0.39
#> d[Self-help vs. No intervention] 0.33 0.01 0.61 -0.87
#> d[Individual counselling vs. Group counselling] -0.31 0.01 0.60 -1.47
#> d[Self-help vs. Group counselling] -0.63 0.02 0.72 -2.06
#> d[Self-help vs. Individual counselling] 0.16 0.02 1.08 -2.00
#> lp__ -5765.22 0.22 6.48 -5779.03
#> tau 0.94 0.01 0.24 0.57
#> 25% 50% 75%
#> d[Group counselling vs. No intervention] 0.59 1.07 1.61
#> d[Individual counselling vs. No intervention] 0.73 0.89 1.08
#> d[Self-help vs. No intervention] -0.05 0.33 0.70
#> d[Individual counselling vs. Group counselling] -0.69 -0.31 0.07
#> d[Self-help vs. Group counselling] -1.09 -0.63 -0.18
#> d[Self-help vs. Individual counselling] -0.52 0.17 0.85
#> lp__ -5769.23 -5764.85 -5760.70
#> tau 0.78 0.90 1.06
#> 97.5% n_eff Rhat
#> d[Group counselling vs. No intervention] 2.79 2108 1
#> d[Individual counselling vs. No intervention] 1.47 1069 1
#> d[Self-help vs. No intervention] 1.54 1858 1
#> d[Individual counselling vs. Group counselling] 0.88 2104 1
#> d[Self-help vs. Group counselling] 0.83 2226 1
#> d[Self-help vs. Individual counselling] 2.28 2896 1
#> lp__ -5753.82 906 1
#> tau 1.47 949 1
#>
#> Samples were drawn using NUTS(diag_e) at Tue Sep 24 08:57:56 2024.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at
#> convergence, Rhat=1).
#>
#> exmp__> ## End(No test)
#> exmp__>
#> exmp__> ## Don't show:
#> exmp__> if (requireNamespace("pkgdown", quietly = TRUE) && pkgdown::in_pkgdown()) {
#> exmp__+ assign("smk_net", smk_net, .GlobalEnv)
#> exmp__+ assign("smk_fit_RE_UME", smk_fit_RE_UME, .GlobalEnv)
#> exmp__+ }
#>
#> exmp__> ## End(Don't show)
#> exmp__>
#> exmp__>
#> exmp__>
# }
# \donttest{
# Compare DIC
smk_dic_RE
#> Residual deviance: 53.8 (on 50 data points)
#> pD: 43.7
#> DIC: 97.5
(smk_dic_RE_UME <- dic(smk_fit_RE_UME)) # no difference in fit
#> Residual deviance: 53.6 (on 50 data points)
#> pD: 44.8
#> DIC: 98.4
# Compare residual deviance contributions
plot(smk_dic_RE, smk_dic_RE_UME, show_uncertainty = FALSE)
# }