R/mtc.network.R

Defines functions plot.mtc.network summary.mtc.network mtc.study.treatment.matrix print.mtc.network mtc.network.graph graph.create edges.create mtc.comparisons mtc.nr.comparisons mtc.treatment.pairs has.indirect.evidence coerce.factor mtc.study.design mtc.studies.list mtc.merge.data as.treatment.factor mtc.network.validate mtc.validate.data.re mtc.validate.data.ab fix.network mtc.network check.duplicated.treatments remove.onearm standardize.data treatment.id.to.description standardize.treatments

Documented in mtc.network plot.mtc.network

standardize.treatments <- function(treatments) {
  treatments[['id']] <- factor(treatments[['id']], levels=treatments[['id']])
  treatments[['description']] <- as.character(treatments[['description']])
  rownames(treatments) <- as.character(treatments[['id']])
  treatments[order(treatments[['id']]), , drop=FALSE]
}

treatment.id.to.description <- function(network, ids) {
  all.ts <- network[['treatments']][['description']]
  names(all.ts) <- as.character(network[['treatments']][['id']])
  all.ts[ids]
}

standardize.data <- function(data, treatment.levels, re.order=FALSE) {
  data[['study']] <- as.factor(data[['study']])
  data[['treatment']] <- factor(as.character(data[['treatment']]), levels=treatment.levels)
  if (re.order) {
    na <- sapply(data[['study']], function(study) { sum(data[['study']] == study) })
    bl <- !is.na(data[['diff']])
    data <- data[order(na, data[['study']], bl, data[['treatment']]), , drop=FALSE]
  } else {
    data <- data[order(data[['study']], data[['treatment']]), , drop=FALSE]
  }
  if (nrow(data) > 0) {
    rownames(data) <- seq(1:nrow(data))
  }
  data
}

remove.onearm <- function(data, warn=FALSE) {
  # Remove 1-arm studies
  sel <- as.logical(sapply(data[,'study'], function(study) {
    sum(data[,'study'] == study) > 1
  }))

  if (warn && !all(sel)) {
    warning(paste("Removing", sum(!sel), "one-arm studies:",
                  paste(data[!sel,'study'], collapse=", ")))
  }
  data[sel, , drop=FALSE]
}

check.duplicated.treatments <- function(network, warn=FALSE) {
  duplicates <- FALSE
  for (study in mtc.studies.list(network)[['values']]) {
    design <- rle(sort(as.vector(mtc.study.design(network, study))))
    dup.v <- design[['values']][design[['lengths']] > 1]
    dup.l <- design[['lengths']][design[['lengths']] > 1]
    if (any(design[['lengths']] > 1)) {
      duplicates <- TRUE
      if (warn) warning(paste(dup.v, "occurs", dup.l, "times in", study))
    }
  }
  duplicates
}

mtc.network <- function(
  data.ab=data, treatments=NULL, description="Network",
  data.re=NULL, studies=NULL, data=NULL
) {
  if (is.null(data.ab) && is.null(data.re)) {
    stop("Either `data.ab' or `data.re' (or both) must be specified")
  }
  # standardize the data
  if (!is.null(data.ab)) {
    data.ab <- remove.onearm(as.data.frame(data.ab), warn=TRUE)
    mtc.validate.data.ab(data.ab)
  }
  if (!is.null(data.re)) {
    data.re <- remove.onearm(as.data.frame(data.re), warn=TRUE)
    mtc.validate.data.re(data.re)
  }

  # standardize the treatments
  if (is.null(treatments)) {
    data.treatments <- vector(mode="character")
    if (!is.null(data.ab)) {
      data.treatments <- c(data.treatments, as.character(data.ab[['treatment']]))
    }
    if (!is.null(data.re)) {
      data.treatments <- c(data.treatments, as.character(data.re[['treatment']]))
    }
    treatments <- sort(unique(data.treatments))
  }
  if (is.list(treatments) && !is.data.frame(treatments)) {
    treatments <- as.data.frame(do.call(rbind, treatments))
  }
  if (is.character(treatments) || is.factor(treatments)) {
    treatments <- data.frame(id=treatments, description=treatments)
  }
  treatments <- standardize.treatments(treatments)

  # standardize the study data
  if (!is.null(studies)) {
    if (!is.data.frame(studies)) {
      stop("studies must be a data frame")
    }
  }

  network <- list(
    description=description,
    treatments=treatments)

  if (!is.null(data.ab) && nrow(data.ab) > 0) {
    network <- c(network, list(data.ab=standardize.data(data.ab, levels(treatments[['id']]))))
  }
  if (!is.null(data.re) && nrow(data.re) > 0) {
    network <- c(network, list(data.re=standardize.data(data.re, levels(treatments[['id']]), re.order=TRUE)))
  }

  if (!is.null(studies) && nrow(studies) > 0) {
    network <- c(network, list(studies=studies[order(studies[['study']]), ,drop=FALSE]))
  }

  mtc.network.validate(network)

  class(network) <- "mtc.network"
  network
}

fix.network <- function(network) {
  # move data into data.ab for legacy networks
  if (is.null(network[['data.ab']]) && !is.null(network[['data']])) {
    network[['data.ab']] <- network[['data']]
    network[['data']] <- NULL
  }
  network
}

mtc.validate.data.ab <- function(data) {
  # No longer validating arm-based data at this stage
  # This is now defined at the likelihood/link level
}

mtc.validate.data.re <- function(data) {
  # These checks are still essential, because they check more then just the presence / absence of
  # columns.

  columns <- colnames(data)
  reColumns <- c('diff', 'std.err')
  if (!all(reColumns %in% columns)) {
    stop(paste('data.re must contain columns: ', paste(reColumns, collapse=', ')))
  }

  baselineCount <- sapply(unique(data[['study']]), function(study) { sum(is.na(data[['diff']][data[['study']] == study])) })
  if (!all(baselineCount == 1)) {
    stop('Each study in data.re must have a unique baseline arm (diff=NA)')
  }
  if (!all(!is.na(data[['std.err']][!is.na(data[['diff']])]))) {
    stop('All non-baseline arms in data.re must have std.err specified')
  }

  studies <- unique(data[['study']])
  ma <- sapply(studies, function(study) { sum(data[['study']] == study) > 2 })
  ma <- studies[ma]
  ma.studies <- data[['study']] %in% ma
  nobaseline <- is.na(data[['std.err']][ma.studies])
  if (!all(!nobaseline)) {
    stop(paste('All multi-arm trials (> 2 arms) must have the std.err of the baseline specified.',
               'Constraint violated by:',
               paste(data[['study']][ma.studies][nobaseline], collapse=", ")))
  }

  se <- data[['std.err']][!is.na(data[['std.err']])]
  if (!all(se > 0.0)) {
    stop(paste('In data.re, std.err must be either positive or NA'))
  }

  if (!all(sapply(ma, function(study) {
    ref <- data[['study']] == study & is.na(data[['diff']])
    rel <- data[['study']] == study & !is.na(data[['diff']])
    all(data[['std.err']][ref] < data[['std.err']][rel])
  }))) {
    stop("In data.re, std.err of the reference arm must < std.err of the relative effects")
  }
    
}

mtc.network.validate <- function(network) {
  # Check that there is some data
  stopifnot(nrow(network[['treatments']]) > 0)
  stopifnot(nrow(network[['data.ab']]) > 0 || nrow(network[['data.re']]) > 0)

  # Check that the treatments are correctly cross-referenced and have valid names
  all.treatments <- forcats::fct_c(
    if (is.null(network[['data.ab']])) factor() else network[['data.ab']][['treatment']],
    if (is.null(network[['data.re']])) factor() else network[['data.re']][['treatment']])
  stopifnot(all(all.treatments %in% network[['treatments']][['id']]))
  # stopifnot(all(network[['treatments']][['id']] %in% all.treatments)) -- disabled for node-splitting
  invalidId <- grep('^[a-zA-Z0-9_]*$', network[['treatments']][['id']], invert=TRUE)
  if (length(invalidId) > 0) {
    stop(paste0("\n", paste0(' Treatment name "',
      network[['treatments']][['id']][invalidId], '" invalid.', collapse="\n"),
      '\nTreatment names may only contain letters, digits, and underscore (_).'))
  }

  # Check that studies are not duplicated between [['data.ab']] and [['data.re']]
  if (!is.null(network[['data.ab']]) && !is.null(network[['data.re']])) {
    dup.study <- intersect(unique(network[['data.ab']][['study']]), unique(network[['data.re']][['study']]))
    if (length(dup.study) > 0) {
      stop(paste('Studies', paste(dup.study, collapse=", "), 'occur in both data and data.re'))
    }
  }

  # Check that the data frame has a sensible combination of columns
  if (!is.null(network[['data.ab']])) {
    mtc.validate.data.ab(network[['data.ab']])
  }

  # Check data.re is well formed
  if (!is.null(network[['data.re']])) {
    mtc.validate.data.re(network[['data.re']])
  }

  # If studies are given, they must match the treatments in the data tables
  if (!is.null(network[['studies']])) {
    data.studies <- c(as.character(network[['data.ab']][['study']]), as.character(network[['data.re']][['study']]))
    if (!setequal(data.studies, as.character(network[['studies']][['study']]))) {
      stop(paste('The studies data frame must match the studies in data.ab and data.re'))
    }
  }

  check.duplicated.treatments(network)
}

as.treatment.factor <- function(x, network) {
  v <- network[['treatments']][['id']]
  if (is.numeric(x)) {
    factor(x, levels=1:nlevels(v), labels=levels(v))
  } else if (is.factor(x) || is.character(x)) {
    x <- as.character(x)
    factor(x, levels=levels(v))
  }
}

mtc.merge.data <- function(network) {
  data.frame(
    study=c(
      as.character(network[['data.ab']][['study']]),
      as.character(network[['data.re']][['study']])),
    treatment=as.treatment.factor(c(
      network[['data.ab']][['treatment']],
      network[['data.re']][['treatment']]), network),
    stringsAsFactors=FALSE)
}

mtc.studies.list <- function(network) {
  rle(as.character(mtc.merge.data(network)[['study']]))
}

mtc.study.design <- function(network, study) {
  data <- mtc.merge.data(network)
  data[['treatment']][data[['study']] == study]
}

coerce.factor <- function(x, prototype) {
  factor(x, levels=1:nlevels(prototype), labels=levels(prototype))
}

# See nodesplit-auto draft for definition
has.indirect.evidence <- function(network, t1, t2) {
  has.both <- sapply(mtc.studies.list(network)[['values']], function(study) {
    all(c(t1, t2) %in% mtc.study.design(network, study))
  })

  data <- mtc.merge.data(network)
  not.both <- !has.both[data[['study']]]
  data <- data[not.both, , drop=FALSE]

  if (nrow(data) > 0) {
    n <- mtc.network(data)
    g <- mtc.network.graph(n)
    all(c(t1, t2) %in% V(g)$name) && as.logical(is.finite(shortest.paths(as.undirected(g), t1, t2)))
  } else {
    FALSE
  }
}

mtc.treatment.pairs <- function(treatments) {
  n <- length(treatments)
  t1 <- do.call(forcats::fct_c, lapply(1:(n-1), function(i) { rep(treatments[i], n - i) }))
  t2 <- do.call(forcats::fct_c, lapply(1:(n-1), function(i) { treatments[(i+1):n] }))
  data.frame(t1=t1, t2=t2)
}

## Get all direct comparisons with the number of studies
## measuring the direct comparison. This is like mtc.comparisons
## but in addition the numbers are appended as an additional column.
mtc.nr.comparisons <- function(network) {
  m <- mtc.study.treatment.matrix(network)
  comp <- mtc.comparisons(network)
  cm <- as.matrix(comp)
  nr <- aaply(cm, 1, function(co) {
    sum(rowSums(m[,co, drop=FALSE]) == 2)
  })
  cbind(comp, nr)
}

# Get all comparisons with direct evidence from the data set.
# Returns a (sorted) data frame with two columns (t1 and t2).
mtc.comparisons <- function(network) {
  ## Generate all pair-wise comparisons from each "design"
  data <- mtc.merge.data(network)
  ## Identify the unique "designs" (treatment combinations)
  design <- function(study) { mtc.study.design(network, study) }
  designs <- unique(lapply(unique(data[['study']]), design))
  comparisons <- do.call(rbind, lapply(designs, mtc.treatment.pairs))
  ## Make sure we include each comparison in only one direction
  swp <- as.character(comparisons[['t1']]) > as.character(comparisons[['t2']])
  tmp <- comparisons[['t1']]
  comparisons[['t1']][swp] <- comparisons[['t2']][swp]
  comparisons[['t2']][swp] <- tmp[swp]

  ## Ensure the output comparisons are always in the same order
  comparisons <- comparisons[order(comparisons[['t1']], comparisons[['t2']]), ,drop=FALSE]
  comparisons[['t1']] <- as.treatment.factor(comparisons[['t1']], network)
  comparisons[['t2']] <- as.treatment.factor(comparisons[['t2']], network)
  ## Ensure the output comparisons are unique
  comparisons <- unique(comparisons)
  row.names(comparisons) <- NULL
  comparisons
}

edges.create <- function(e, ...) {
  ed <- t(matrix(c(e[['t1']], e[['t2']]), ncol=2))
  edges(as.vector(ed), weight=e[['nr']], ...)
}

graph.create <- function(v, e, ...) {
  g <- graph.empty()
  g <- g + vertex(levels(v))
  g <- g + edges.create(e, ...)
  g
}

mtc.network.graph <- function(network, include.nr.comparisons=FALSE) {
    comp <- if (include.nr.comparisons) mtc.nr.comparisons else mtc.comparisons
    comparisons <- comp(network)
    treatments <-   network[['treatments']][['id']]
    graph.create(treatments, comparisons, arrow.mode=0)
}

## mtc.network class methods
print.mtc.network <- function(x, ...) {
  x <- fix.network(x)
  cat("MTC dataset: ", x[['description']], "\n", sep="")
  if (!is.null(x[['data.ab']])) {
    cat('Arm-level data: \n')
    print(x[['data.ab']])
  }
  if (!is.null(x[['data.re']])) {
    cat('Relative effect data: \n')
    print(x[['data.re']])
  }
  if (!is.null(x[['studies']])) {
    cat('Study-level data: \n')
    print(x[['studies']])
  }
}

mtc.study.treatment.matrix <- function(network) {
  object <- fix.network(network)
  data <- mtc.merge.data(object)
  studies <- unique(data[['study']])
  m <- laply(object[['treatments']][['id']], function(treatment) {
    laply(studies, function(study) {
      any(data[['study']] == study & data[['treatment']] == treatment)
    }, .drop=TRUE)
  }, .drop=FALSE)
  m <- t(m)
  colnames(m) <- object[['treatments']][['id']]
  m
}

summary.mtc.network <- function(object, ...) {
  object <- fix.network(object)
  m <- mtc.study.treatment.matrix(object)
  x <- as.factor(apply(m, 1, sum))
  levels(x) <- sapply(levels(x), function(y) { paste(y, "arm", sep="-") })
  list("Description"=paste("MTC dataset: ", object[['description']], sep=""),
       "Studies per treatment"=apply(m, 2, sum),
       "Number of n-arm studies"=summary(x),
       "Studies per treatment comparison"=mtc.nr.comparisons(object)
       )
}

plot.mtc.network <- function(x, layout=igraph::layout.circle, dynamic.edge.width=TRUE, use.description=FALSE, ...) {
  x <- fix.network(x)
  g <- mtc.network.graph(x, TRUE)
  labels <- if (use.description)
                treatment.id.to.description(x, x[['treatments']][['id']])
            else
                as.character(x[['treatments']][['id']])

  if (dynamic.edge.width) {
    igraph::plot.igraph(g, layout=layout, edge.width=E(g)$weight, vertex.label=labels, ...)
  } else {
    igraph::plot.igraph(g, layout=layout, vertex.label=labels, ...)
  }
}

Try the gemtc package in your browser

Any scripts or data that you put into this service are public.

gemtc documentation built on May 15, 2021, 1:07 a.m.