joelkuiper/convert.R

## convert.R
library(gemtc)

convert <- function(re) {
    options(stringsAsFactors = FALSE)  # Because R
    results <- data.frame()
    studies <- unique(re$study)

    entry <- function(study, treatment, diff, std.err) {
        list(study = study, treatment = treatment, diff = diff, std.err = std.err)
    }

    as.entry <- function(row) {
        entry(row$study, row$treatment, row$diff, row$std.err)
    }

    for(studyId in studies) {
        data <- subset(re, study == studyId)
        if(nrow(data) == 1) {
            ## Two-arm study
            base <- entry(studyId, data$base, NA, NA)
            treatment <- as.entry(data[1, ])
            results <- do.call("rbind", list(results, base, treatment))

        } else {
            ## Multi-arm study
            ## These entries are similar, but require standard error of the mean for the baseline
            ## We try to compute this if possible, if not we approximate it

            ## Guess at most common base, which is the one most referred to
            base.table <- table(data$base)
            base.trt <- names(base.table)[[which.max(base.table)]]

            treatments <- subset(data, base == base.trt)

            calc.base.se <- function(d.ab.se, d.ac.se, d.bc.se) {
                ## Let $Var(D_{BC}) = Var(D_{AB}) + Var(D_{AC}) - 2Var(Y_A)$, thus
                ## $se_B = \sqrt{(se^2_{AB} + se^2_{CB} - se^2_{AC}) / 2}$
                sqrt((d.ab.se^2 + d.ac.se^2 - d.bc.se^2) / 2)
            }

            approx.active.comparison <- function(d.ab, d.ac) {
                ## Try to use approximate Var(D_bc) using number of participants, if available (formula 9, Brooks et.al.)
                ## assuming the standard error is proportional to 1/sqrt(n)

                stopifnot(d.ab$base.n == d.ac$base.n) ## A is base should have same number of participants
                n.a <- d.ab$base.n
                n.b <- d.ab$treatment.n
                n.c <- d.ac$treatment.n

                if(all(!is.na(c(n.a, n.b, n.c)))) {
                    sqrt(((d.ab$std.err^2 + d.ac$std.err^2) * (1/n.b + 1/n.c)) / (1/n.b + 1/n.c + 2/n.a))
                } else {
                    NA
                }
            }

            base.se <- function(b, trt1, trt2) {
                d.ab <- subset(data, base == b & treatment == trt1)
                d.ac <- subset(data, base == b & treatment == trt2)
                d.bc <- subset(data, base == trt2 & treatment == trt1)

                se <- NA
                if(nrow(d.bc) != 0) {
                    ##  We have enough data to compute baseline se
                    se <- calc.base.se(d.ab$std.err, d.ac$std.err, d.bc$std.err)
                } else {
                    ## We need to approximate D_bc
                    d.bc.std.err <- approx.active.comparison(d.ab, d.ac)
                    if(!is.na(d.bc.std.err)) {
                        se <- calc.base.se(d.ab$std.err, d.ac$std.err, d.bc.std.err)
                    }
                }
                se
            }

            ## Find the combinations of treatments that we could use (e.g. BC, BD)
            combinations <- t(combn(as.character(treatments$treatment), 2))

            ## Try methods by Brooks et.al. to computer or approximate baseline se,
            ## We take the median of all these values, NAs ommited
            se <- median(apply(combinations, 1, function(row) { base.se(base.trt, row[[1]], row[[2]]) }), na.rm=T)


            ## Append the treatments associated with the base
            if(is.finite(se)) { # not NaN, NA or infinite
                results <- rbind(results, entry(studyId, base.trt, NA, se))
                for(entry in by(treatments, 1:nrow(treatments), as.entry)) {
                    results <- rbind(results, entry)
                }
            } else {
                warning(paste("could not calculate baseline std.err from data for study", studyId, "omitting"))
            }
        }
    }

    options(stringsAsFactors = TRUE) # Because R, but lets not break compatibility
    results
}

data.ab <- read.table(textConnection('
study treatment responders  sampleSize
01  Salmeterol  1 229
01  Placebo 1 227
02  Fluticasone 4 374
02  Salmeterol  3 372
02  SFC 2 358
02  Placebo 7 361
03  Salmeterol  1 554
03  Placebo 2 270'), header=T)

re <- read.table(textConnection('
study        base   treatment   diff std.err base.n treatment.n
    4     Placebo Fluticasone -0.267   0.203     NA          NA
    5     Placebo         SFC -0.209   0.098   1524        1533
    5     Placebo  Salmeterol -0.154   0.096   1524        1521
    5     Placebo Fluticasone  0.055   0.092   1524        1534
    5  Salmeterol         SFC -0.056   0.100   1521        1533
    5 Fluticasone         SFC -0.264   0.096   1534        1533
'), header=T)

data.re <- convert(re)
data.re

network <- mtc.network(data.ab=data.ab, data.re=data.re)
model <- mtc.model(network,
                   link="cloglog",
                   likelihood="binom",
                   linearModel="fixed",
                   hy.prior=mtc.hy.prior("std.dev", "dunif", 0, "om.scale"))
mtc.run(model) -> results
forest(relative.effect(results, t1="Placebo"))
	library(gemtc)

	convert <- function(re) {
	options(stringsAsFactors = FALSE) # Because R
	results <- data.frame()
	studies <- unique(re$study)

	entry <- function(study, treatment, diff, std.err) {
	list(study = study, treatment = treatment, diff = diff, std.err = std.err)
	}

	as.entry <- function(row) {
	entry(row$study, row$treatment, row$diff, row$std.err)
	}

	for(studyId in studies) {
	data <- subset(re, study == studyId)
	if(nrow(data) == 1) {
	## Two-arm study
	base <- entry(studyId, data$base, NA, NA)
	treatment <- as.entry(data[1, ])
	results <- do.call("rbind", list(results, base, treatment))

	} else {
	## Multi-arm study
	## These entries are similar, but require standard error of the mean for the baseline
	## We try to compute this if possible, if not we approximate it

	## Guess at most common base, which is the one most referred to
	base.table <- table(data$base)
	base.trt <- names(base.table)[[which.max(base.table)]]

	treatments <- subset(data, base == base.trt)

	calc.base.se <- function(d.ab.se, d.ac.se, d.bc.se) {
	## Let $Var(D_{BC}) = Var(D_{AB}) + Var(D_{AC}) - 2Var(Y_A)$, thus
	## $se_B = \sqrt{(se^2_{AB} + se^2_{CB} - se^2_{AC}) / 2}$
	sqrt((d.ab.se^2 + d.ac.se^2 - d.bc.se^2) / 2)
	}

	approx.active.comparison <- function(d.ab, d.ac) {
	## Try to use approximate Var(D_bc) using number of participants, if available (formula 9, Brooks et.al.)
	## assuming the standard error is proportional to 1/sqrt(n)

	stopifnot(d.ab$base.n == d.ac$base.n) ## A is base should have same number of participants
	n.a <- d.ab$base.n
	n.b <- d.ab$treatment.n
	n.c <- d.ac$treatment.n

	if(all(!is.na(c(n.a, n.b, n.c)))) {
	sqrt(((d.ab$std.err^2 + d.ac$std.err^2) * (1/n.b + 1/n.c)) / (1/n.b + 1/n.c + 2/n.a))
	} else {
	NA
	}
	}

	base.se <- function(b, trt1, trt2) {
	d.ab <- subset(data, base == b & treatment == trt1)
	d.ac <- subset(data, base == b & treatment == trt2)
	d.bc <- subset(data, base == trt2 & treatment == trt1)

	se <- NA
	if(nrow(d.bc) != 0) {
	## We have enough data to compute baseline se
	se <- calc.base.se(d.ab$std.err, d.ac$std.err, d.bc$std.err)
	} else {
	## We need to approximate D_bc
	d.bc.std.err <- approx.active.comparison(d.ab, d.ac)
	if(!is.na(d.bc.std.err)) {
	se <- calc.base.se(d.ab$std.err, d.ac$std.err, d.bc.std.err)
	}
	}
	se
	}

	## Find the combinations of treatments that we could use (e.g. BC, BD)
	combinations <- t(combn(as.character(treatments$treatment), 2))

	## Try methods by Brooks et.al. to computer or approximate baseline se,
	## We take the median of all these values, NAs ommited
	se <- median(apply(combinations, 1, function(row) { base.se(base.trt, row[[1]], row[[2]]) }), na.rm=T)


	## Append the treatments associated with the base
	if(is.finite(se)) { # not NaN, NA or infinite
	results <- rbind(results, entry(studyId, base.trt, NA, se))
	for(entry in by(treatments, 1:nrow(treatments), as.entry)) {
	results <- rbind(results, entry)
	}
	} else {
	warning(paste("could not calculate baseline std.err from data for study", studyId, "omitting"))
	}
	}
	}

	options(stringsAsFactors = TRUE) # Because R, but lets not break compatibility
	results
	}

	data.ab <- read.table(textConnection('
	study treatment responders sampleSize
	01 Salmeterol 1 229
	01 Placebo 1 227
	02 Fluticasone 4 374
	02 Salmeterol 3 372
	02 SFC 2 358
	02 Placebo 7 361
	03 Salmeterol 1 554
	03 Placebo 2 270'), header=T)

	re <- read.table(textConnection('
	study base treatment diff std.err base.n treatment.n
	4 Placebo Fluticasone -0.267 0.203 NA NA
	5 Placebo SFC -0.209 0.098 1524 1533
	5 Placebo Salmeterol -0.154 0.096 1524 1521
	5 Placebo Fluticasone 0.055 0.092 1524 1534
	5 Salmeterol SFC -0.056 0.100 1521 1533
	5 Fluticasone SFC -0.264 0.096 1534 1533
	'), header=T)

	data.re <- convert(re)
	data.re

	network <- mtc.network(data.ab=data.ab, data.re=data.re)
	model <- mtc.model(network,
	link="cloglog",
	likelihood="binom",
	linearModel="fixed",
	hy.prior=mtc.hy.prior("std.dev", "dunif", 0, "om.scale"))
	mtc.run(model) -> results
	forest(relative.effect(results, t1="Placebo"))