Running MCMC simulation
Once the network model is set up, you can perform a MCMC simulation. The overall process is to set and run an appropriate number of simulations, and then check whether the results converge.
First, an example of the fixed effect model will be explained.
· mcmc_b_bin_fe <- mtc.run(model_b_bin_fe, n.adapt=5000, n.iter=10000, thin=20)
In the “mtc.run” function, enter the name of the fixed effect model that has been set. “n.adapt=5000” means to discard no.1-5,000 of the iterations. This is called burn-in, which is to remove a certain part of the beginning of the created random numbers to exclude the effect of the initial values of the algorithm. “n.iter=10000” means to perform 10,000 simulations and “thin=20” means to extract every 20th value.
To summarize the above explanation, the first to 5,000th data are discarded (to reduce the effect of initial values in simulation), simulations are performed 10,000 times with 5,001st to 15,000th values, and every 20th value is extracted (e.g., 5,020, 5,040….).
In the Bayesian analysis, prior distribution considering multichain is input to determine the posterior distribution. The multichain simulation is performed by setting multiple initial values for the prior parameter of prior distribution, that is, the hyperparameter d (e.g., 4 values of -1, 0, 1, and 2). Therefore, because every 20th data are extracted among 10,000 simulations, 500 data points are extracted in each chain.
You can see a more detailed explanation by using the summary command:
· summary(mcmc_b_bin_fe)
MCMC simulation and convergence status
To verify if the MCMC simulation converged well, you can check the following items in combination.
MCMC error
A smaller MCMC error indicates a higher accuracy, which means a good convergence. Therefore, a sufficient sample size should be achieved by performing many simulations and the burn-in process to remove the effect of initial values, and the data extraction interval “thin” should be adjusted appropriately.
Deviance information criterion
The deviance information criterion (DIC) is expressed as DIC=D¯+pD, where D¯ is the sum of residual deviances and pD is an estimated value of the parameter. Thus, the DIC considers both the fitness and complexity of the model, and the smaller the DIC is, the better the model.
Trace plot and density plot
If the trace plot (a graph visually showing the simulation result) has no specific pattern and the chains are entangled, it is considered that the convergence is good. The density plot is a posterior distribution (posterior density function) and if the shapes are significantly different for the same number of simulations, it means the data did not converge well.
Figure 5A is the case when 100 and 500 were input as the total number of iterations with no burn-in. In the first 100 iterations (
Figure 5A left), the four chains have severe variations and are not even, but at approximately 500th iteration (
Figure 5A right), the graph becomes even, with no specific pattern. Therefore, it is desirable to perform at least 500 simulations for each channel. For burn-in to exclude the effect of initial values, 1-100th data points must be discarded without question because they are too uneven, and the number of discarded values must be at least 500. In this example, the values were discarded for up to 5,000 times to minimize the effect of the algorithm.
In
Figure 5B, when the simulations are performed 10,000 times with an extraction interval “thin” of 20, 500 data points are extracted from each channel. If 10 is input to thin (
Figure 5C), 1,000 values are extracted per channel, which increases the sample size and results in a more even distribution of the trace plot. Furthermore, 1,000 samples of post density function were selected because it looks more similar to the normal distribution.
The finally selected model performed burn-in for 5,000 values and 10,000 simulations, and extracted every 10th data point for 1,000 samples per channel.
· mcmc_b_bin_fe <- mtc.run(model_b_bin_fe, n.adapt=5000, n.iter=10000, thin=10)
Compared to the extraction interval thin of 20, the total sample size increased and the MCMC standard error of the treatment group AB decreased from 0.004 to 0.003, thus increasing the accuracy. However, the DIC decreased from 56.24 to 56.63, which is practically insignificant.
Gelman-Rubin statistics and plot
· gelman.diag(mcmc_b_bin_fe)
· gelman.plot(mcmc_b_bin_fe)
The “gelman.diag” command displays the Gelman-Rubin statistics on the console, and the “gelman.plot” command draws the Gelman-Rubin plot. As the number of simulations increases, it approaches 1, and the variations must be stabilized so that it can be said to have converged well.
Selecting the final model for MCMC simulationFor MCMC simulation, a model that converges best should be selected by adjusting the number of chains appropriate for multichain, the number of data for removal of initial effect (burn-in), the number of iterations, and the extraction interval (thin).
For the fixed effect model of this example, 4 chains, 5,000 burnins, 10,000 iterations, and an interval of 10 were selected, to sufficiently remove the effect of initial values, increase the iterations and extraction interval, and minimize the MCMC error and DIC variation with almost no variations and stability of various plots.
However, you should adjust the iterations appropriately because it can take significant time depending on the computer specifications.
Burn-in 5,000, iteration 10,000, thin 10
· mcmc_b_bin_fe <- mtc.run(model_b_bin_fe, n.adapt=5000, n.iter=10000, thin=10)
· plot(mcmc_b_bin_fe)
· summary(mcmc_b_bin_fe)
· gelman.diag(mcmc_b_bin_fe)
· gelman.plot(mcmc_b_bin_fe)
Consistency test
Consistency test in the assumptions of NMA is a critical tool that determines the applicability of NMA results.
· nodesplit_b_bin_fe <- mtc.nodesplit(network_b_bin, linearModel=‘fixed’, n.adapt=5000, n.iter=10000, thin=10)
· plot(nodesplit_b_bin_fe)
· plot(summary(nodesplit_b_bin_fe))
The fixed effect model “nodesplit_b_bin_fe” is created for consistency test by entering the network set-up data in the “mtc.nodesplit” function. The MCMC simulation is also performed.
The variations between treatments and the consistency test results of all individual treatments can be easily seen. As a result of the consistency test, the p-value of treatments E versus D was 0.043, indicating inconsistency. However, no statistical significance was observed in all the other treatments. Therefore, it is desirable to set up a random effect model for a more robust analysis.
Forest plot
Forest plot allows graphical comparison of the effect sizes by treatment group through NMA.
· forest(relative.effect(mcmc_b_bin_fe, t1=“A”), digits=3)
When you enter the final model through MCMC simulation in the forest function, a forest plot with A reference treatment is created (
Figure 6).
The effect size (OR, blood transfusion rate) of every treatment was lower than the placebo, and the 95% credible intervals did not overlap.
In particular, the blood transfusion rate of the combination treatment method (E) was statistically significantly lower compared to those of all the other treatments including single intravenous injection (IV) injection [B, IV(single)], double IV injection [C, IV(double)], and topical application method [D, topical].
Treatment ranking
One of the most important functions of NMA is that the comparative advantages of treatments can be determined. In other words, the cumulative probability that the top priority to the lowest priority treatments are selected can be calculated.
· ranks_b_bin_fe <- rank.probability(mcmc_b_bin_fe, preferredDirection=-1)
· print(ranks_b_bin_fe)
Enter the final MCMC model in the “rank.probability” function. Set the “preferredDirection” to ‘-1’ or ‘1’, depending on whether a smaller effect size indicates a better treatment, or vice versa. In this example, it is set to ‘-1’ because an effect size smaller than that of the reference treatment means a better treatment.
As can be seen in the probability table, the best treatment is E (combination) at 99.8%, followed by C (IV double) at 68.2%, B (IV single), D (topical), and A (placebo).