plotCovariateEffect command only plot the frame and the axes but not the effects values

[fersalme]

Hi, first of all, thank you this amazing package!

The plotCovariateEffect command for meta-regression topic doesn't work only plot the frame and the axes but not the effects values (mean, CI)

Here a examples from certolizumab vignette

library(gemtc)

regressor <- list(coefficient='shared',
                  variable='diseaseDuration',
                  control='Placebo')

hy.prior <- mtc.hy.prior(type="std.dev", distr="dhnorm", 0, 9.77)

model <- mtc.model(certolizumab,
                   type="regression",
                   regressor=regressor,
                   hy.prior=hy.prior)

result <- mtc.run(model)

plotCovariateEffect(result, t1 = "Placebo", t2 = "Tocilizumab")

The problem maybe come from treatment.pairs function. I make a quick and dirty plotCovariateEffect_2 function based on plotCovariateEffect

plotCovariateEffect_2 <- function(result, t1, t2, xlim = NULL, ylim = NULL, plot = TRUE) {

  treatment.pairs <- function(t1, t2, ts) {
    if((is.null(t2) || length(t2) == 0) && length(t1) == 1) {
      t2 <- ts[ts != as.numeric(t1)]
    }
    if(length(t1) > length(t2)) t2 <- rep(t2, length.out=length(t1))
    if(length(t2) > length(t1)) t1 <- rep(t1, length.out=length(t2))
    matrix(c(t1, t2), ncol=2)
  }

  extract.comparisons <- function(parameters) {
    x <- parameters[grep("^d\\.", parameters)]
    matrix(unlist(strsplit(x, "\\.")), ncol=3, byrow=TRUE)[ , -1, drop=FALSE]
  }

  regressor <- result[["model"]][["regressor"]]
  if (is.null(xlim)) {
    studies <- result[["model"]][["network"]][["studies"]]
    observed <- studies[, regressor[["variable"]]]
    ctr <- regressor[["center"]]
    xlim <- c(min(observed, na.rm = TRUE), max(observed, 
                                               na.rm = TRUE))
    xvals <- seq(xlim[1], xlim[2], length.out = 7)
  } else {
    xvals <- seq(xlim[1], xlim[2], length.out = 7)
  }
  pairs <- treatment.pairs(t1, t2, result[["model"]][["network"]][["treatments"]][["id"]])
  res <- lapply(xvals, function(xval) {
    re <- relative.effect(result, t1, t2, preserve.extra = FALSE, 
                          covariate = xval)
    samples <- as.matrix(re[["samples"]])
    stats <- t(apply(samples, 2, quantile, probs = c(0.025, 
                                                     0.5, 0.975)))
    comps <- extract.comparisons(rownames(stats))
    data.frame(t1 = comps[, 1], t2 = comps[, 2], median = stats[, 
                                                                "50%"], lower = stats[, "2.5%"], upper = stats[, 
                                                                                                               "97.5%"])
  })
  if (is.null(ylim)) {
    ylim <- c(min(sapply(res, function(stats) {
      min(stats[["lower"]])
    })), max(sapply(res, function(stats) {
      max(stats[["upper"]])
    })))
  }
  for (pair in split(pairs, seq(nrow(pairs)))) {
    yvals <- sapply(res, function(stats) {
      stats[stats[["t1"]] == pair[1] & stats[["t2"]] == 
              pair[2], c("median", "lower", "upper")]
    })
    if (plot) {
      plot(xvals, yvals["median", ], type = "l", xlim = xlim, 
           ylim = ylim, main = "Treatment effect vs. covariate", 
           xlab = regressor[["variable"]], ylab = paste("d", 
                                                        pair[1], pair[2], sep = "."))
      lines(xvals, yvals["lower", ], lty = 2)
      lines(xvals, yvals["upper", ], lty = 2)
    } else {
      return(cbind(xvals, t(yvals)))
    }
  }
}

plotCovariateEffect_2(result, t1 = "Placebo", t2 = "Tocilizumab")
plotCovariateEffect_2(result, t1 = "Placebo", t2 = "Tocilizumab", plot = FALSE)

I hope this maybe help for fix the problem.

Cheers!

R version 4.0.5 (2021-03-31) Rstudio Version 1.4.1106 gemtc Version: 0.8-8 Built: R 4.0.5; x86_64-pc-linux-gnu; 2021-04-20 13:35:16 UTC; unix

[gertvv]

I think you're right it's due to a change to treatment.pairs - R is becoming more strict about comparing and converting factors to numerics, and that behavior is what I used to use to bridge between gemtc's internal representation of treatments and igraph's numerical IDs. So this function had to return numerics for that, but I hadn't realized that in this particular place the factor levels were being compared to strings instead. Fixed on master, thanks for reporting!