Count data analysis

Introduction

This document describes how to use the main functions of NMA to run a single network meta-analysis with non-survival data.

Example

First load the required packages.

library(NMA)
library(dplyr)
library(purrr)

Settings

Define the BUGS parameters for MCMC. This is not necessary, but recommended, because there are default values for these.

bugs_params <-
  list(
    PROG = "openBugs",  # which version of BUGS to use to run the MCMC
    N.BURNIN = 10,#00,  # number of steps to throw away
    N.SIMS = 150,#0,    # total number of simulations
    N.CHAINS = 2,       # number of chains
    N.THIN = 1,         # thinning rate
    PAUSE = TRUE)

Define the scenario we will use for the analysis.

RANDOM <- FALSE             # is this a random effects model?
REFTX <- "X"                # reference treatment
data_type <- "count_data"   # which type of data to use
label_name <- "label_name"

Read in datasets

Lets read in the each data set. In another article we will show how to do this in one function call by including a Reference file in the data folder which contains the meta data of how to read in the study data.

file_name <- here::here(file.path("inst", "extdata"))

countData <-
  read.csv(paste0(file_name, "/count_data_test.csv"),
           header = TRUE,
           as.is = TRUE)

The format of the data should be fairly self-explanatory and looks like the following. For the hazard ratio data

study	treatment	E	r
Y	B	68364.23	143
X	A	48992.76	159
X	B	45774.40	121
Z	B	68364.23	143
Z	C	68364.23	143

More information about these data is available in the help documentation which can be accessed with e.g. help(countData).

Build model

Now we can create the NMA object to use in the modelling. The workflow is to first create this separately to actually doing the fitting. This then means that we can perform modified fits but we don’t have to redo any of the preparatory work.

nma_model <-
  new_NMA(countData = countData,
          bugs_params = bugs_params,
          is_random = RANDOM,
          data_type = data_type,
          refTx = REFTX ,
          label = "",
          endpoint = "")

nma_model
#> $dat
#> $dat$inits
#> function() {
#>     list(
#>       beta = c(NA, rnorm(nTx - 1, 0, 2)),
#>       sd = 0.1,
#>       alpha = rnorm(nStudies),
#>       d = c(NA, rnorm(nTx - 1, 0, 2)),  ##TODO: can we remove duplication?
#>       mu = rnorm(nStudies),
#>       baseLod = 0) %>% 
#>       .[param_names]  # filter redundant
#>   }
#> <bytecode: 0x00000000135948a8>
#> <environment: 0x00000000135a10f8>
#> 
#> $dat$countData
#>   study treatment        E   r
#> 1     Y         B 68364.23 143
#> 2     X         A 48992.76 159
#> 3     X         B 45774.40 121
#> 4     Z         B 68364.23 143
#> 5     Z         C 68364.23 143
#> 
#> $dat$bugsData
#> $dat$bugsData$mu_beta
#> [1] 0
#> 
#> $dat$bugsData$prec_beta
#> [1] 1e-06
#> 
#> $dat$bugsData$mu_alpha
#> [1] 0
#> 
#> $dat$bugsData$prec_alpha
#> [1] 1e-06
#> 
#> $dat$bugsData$t
#>      [,1] [,2]
#> [1,]    2   NA
#> [2,]    1    2
#> [3,]    2    3
#> 
#> $dat$bugsData$r
#>      [,1] [,2]
#> [1,]  143   NA
#> [2,]  159  121
#> [3,]  143  143
#> 
#> $dat$bugsData$E
#>          [,1]     [,2]
#> [1,] 68364.23       NA
#> [2,] 48992.76 45774.40
#> [3,] 68364.23 68364.23
#> 
#> $dat$bugsData$nt
#> [1] 3
#> 
#> $dat$bugsData$na
#> [1] 1 2 2
#> 
#> $dat$bugsData$ns
#> [1] 3
#> 
#> 
#> $dat$txList
#> [1] "A" "B" "C"
#> 
#> 
#> $data_type
#> [1] "count_data"
#> 
#> $bugs_params
#> $bugs_params$PROG
#> [1] "openBugs"
#> 
#> $bugs_params$N.BURNIN
#> [1] 10
#> 
#> $bugs_params$N.SIMS
#> [1] 150
#> 
#> $bugs_params$N.CHAINS
#> [1] 2
#> 
#> $bugs_params$N.THIN
#> [1] 1
#> 
#> $bugs_params$PAUSE
#> [1] TRUE
#> 
#> $bugs_params$run_bugs
#> [1] TRUE
#> 
#> 
#> $bugs_fn
#> function(...)
#>         R2OpenBUGS::bugs(...)
#> <bytecode: 0x0000000016b49030>
#> <environment: 0x0000000016b543f0>
#> 
#> $is_random
#> [1] FALSE
#> 
#> $refTx
#> [1] "X"
#> 
#> $effectParam
#> [1] NA
#> 
#> $modelParams
#>  count_data 
#> "totresdev" 
#> 
#> $label
#> [1] ""
#> 
#> $endpoint
#> [1] ""
#> 
#> attr(,"class")
#> [1] "nma"
#> attr(,"CALL")
#> attr(,"CALL")$countData
#> countData
#> 
#> attr(,"CALL")$bugs_params
#> bugs_params
#> 
#> attr(,"CALL")$is_random
#> RANDOM
#> 
#> attr(,"CALL")$data_type
#> data_type
#> 
#> attr(,"CALL")$refTx
#> REFTX
#> 
#> attr(,"CALL")$label
#> [1] ""
#> 
#> attr(,"CALL")$endpoint
#> [1] ""

We can view the network graph.

library(sna)
plotNetwork(nma_model)

Run MCMC

The NMA MCMC function calls the appropriate BUGS model.

nma_res <- NMA_run(nma_model, save = FALSE)
#> ====== RUNNING BUGS MODEL

nma_res
#> Inference for Bugs model at "C:/Users/n8tha/Documents/R/NMA/inst/FE_count.txt", 
#>  2 chains, each with 160 iterations (first 10 discarded)
#>  n.sims = 300 iterations saved
#>              mean      sd 2.5%  25%  50%     75%  97.5% Rhat n.eff
#> totresdev 55302.0 87204.9 4584 8617 8704 37409.5 239200  1.6     7
#> deviance  55336.2 87204.4 4618 8651 8738 37439.5 239300  1.6     7
#> 
#> 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 = Dbar-Dhat)
#> pD = 50530.0 and DIC = 105900.0
#> DIC is an estimate of expected predictive error (lower deviance is better).

Plot and tables

BUGS plots are available for diagnosing the performance.

diagnostics(nma_res, save = TRUE)

Different NMA tables can be created. They can provide a record of the analysis.

write_data_to_file(nma_model)
write_results_table(nma_model, nma_res)

pairwiseTable(nma_model, nma_res)

Currently available NMA plots are a treatment effect forest plot of posterior samples and a rank probability grid.

txEffectPlot(nma_model, nma_res)
rankProbPlot(nma_model, nma_res)

It’s also possible to create all of the BUGS and output plots and table functions together and write to an analysis folder.

nma_outputs(nma_res2, save = TRUE)