CODEWITHSUNDEEP
rplot
Andrew Taylor
Andrew Taylor

Reputation: 3488

Changing x axis limits and ticks for mada::forest()

if I'm making forest plots for a diagnostic test accuracy meta-analysis:

dat<-data.frame(TP=sample(c(25:100),20,replace=TRUE),
            FN=sample(c(25:100),20,replace=TRUE),
            FP=sample(c(25:100),20,replace=TRUE),
            TN=sample(c(25:100),20,replace=TRUE))

I want to plot the DOR:

library(mada)
fit.DOR.DSL <- madauni(dat)
forest(fit.DOR.DSL)

But I want to customize the x axis upper and lower limit and tick marks. The forest() plot doesn't have any obvious arguments. It allows for passing on additional graphical parameters, and I've tried throwing some random darts like xlim=c(), at=c(), etc. I feel as if I'm throwing darts in the dark.

Anyone know the magical argument I need to throw in here? I have no experience passing "additional graphical parameters".

As always, thanks!

Here's my updated functions, the maintainer, Philipp Doebler, was very helpful, and includes some other stuff I wanted the function to do.

#instead of using forest() for mada package, use forest.madad()
#and forest.madauni(), from this file, they allow for control over
#the x-axis and the additional of (1) summary sens and spec in forest.madad()
#and a verticle line at the summar statistic in either forest.madad() or 
#forest.madauni()

myforest <- 
  function (x, ci, plotci = TRUE, main = "Forest plot", xlab = NULL,
            SumLine=SumLine, digits = 2L, snames = NULL, subset = NULL, pch     = 15, cex = 1, 
            cipoly = NULL, polycol = NA, plotxaxis = TRUE, ...) 
  {
    stopifnot(length(x) == dim(ci)[1], all(!is.na(c(ci, x))), 
              is.logical(plotci))
    if (!is.null(snames)) {
      stopifnot(length(snames) == length(x))
    }
    if (is.null(snames)) {
      snames <- paste("Study", 1:length(x))
    }
    if (!is.null(subset)) {
      stopifnot(length(subset) > 0, is.integer(subset))
    }
    if (is.null(subset)) {
      subset <- 1:length(x)
    }
    if (!is.null(cipoly)) {
      stopifnot(length(cipoly) == length(x))
      cireg <- which(!cipoly)
      cipoly <- which(cipoly)
    }
    if (is.null(cipoly)) {
      cireg <- 1:length(x)
    }
    x <- x[subset]
    lb <- ci[subset, 1]
    ub <- ci[subset, 2]
    snames <- snames[subset]
    if (is.null(xlab)) {
      xlab <- ""
    }
    N <- length(x)
    plotrange <- max(ub) - min(lb)
    if (plotci) {
      xlim <- c(min(lb) - plotrange * 1.2, max(ub) + plotrange * 
                  1.2)
    }
    else {
      xlim <- c(min(lb) - plotrange * 1.2, max(ub) + plotrange * 
                  0.4)
    }
    ylim <- c(0.5, N + 1)
    plot(NA, NA, xlim = xlim, ylim = ylim, xlab = xlab, ylab = "", 
         yaxt = "n", xaxt = "n", xaxs = "i", bty = "n", main = main, 
         ...)
    abline(h = ylim[2], ...)
    if(!is.null(SumLine)) abline(v = x[length(x)],lty=2, ...)

    for (i in cireg) {
      points(x[i], (N:1)[i], pch = pch)
      arrows(lb[i], (N:1)[i], ub[i], angle = 90, code = 3, 
             length = 0.05, ...)
    }
    for (i in cipoly) {
      polygon(x = c(lb[i], x[i], ub[i], x[i]), y = c((N:1)[i], 
                                                     (N:1)[i] + 0.25, (N:1)    [i], (N:1)[i] - 0.25), col = polycol)
    }
    text(x = xlim[1], N:1, labels = snames[1:N], pos = 4, cex = cex)
    if (plotci) {
      citext <- format(round(cbind(x, ci), digits = digits), 
                       nsmall = digits)
      citext <- paste(citext[, 1], " [", citext[, 2], ", ", 
                      citext[, 3], "]", sep = "")
      text(x = xlim[2], N:1, labels = citext, pos = 2, cex = cex)
    }
    if(plotxaxis){
      axis(1, at = round(seq(from = min(lb), to = max(ub), length.out = 5), 
                         digits), cex = cex)
    }
    return(invisible(NULL))
  }

forest.madauni <- 
  function (x, log = TRUE, plotxaxis = TRUE,SumLine=NULL, snames=NULL, ...) 
  {
    fit <- x
    stopifnot(class(fit) == "madauni")
    descr <- fit$descr
    level <- fit$descr$level
    summ <- summary(fit, level = level)
    forest.x <- fit$theta
forest.ci <- switch(fit$type, DOR = descr$DOR$DOR.ci, negLR =     descr$negLR$negLR.ci, 
                        posLR = descr$posLR$posLR.ci)
    forest.x <- c(forest.x, summ$CIcoef[1, 1])
    forest.ci <- rbind(forest.ci, summ$CIcoef[1, 2:3])
    if (!exists("snames")) {
      snames <- descr$names
      if (is.null(snames)) {
        snames <- paste("Study", 1:descr$nobs)
      }
      snames <- c(snames, paste("Summary (", fit$method, ")", 
                                sep = ""))
    }
    xlab = switch(fit$type, DOR = "diagnostic odds ratio", negLR = "negative     likelihood ratio", 
                  posLR = "positive likelihood ratio")
    if (log) {
      forest.x <- log(forest.x)
      forest.ci <- log(forest.ci)
      xlab <- paste("log", xlab)
    }
    myforest(x = forest.x, ci = forest.ci, snames = snames, SumLine=SumLine,
             xlab = xlab, cipoly = c(rep(FALSE, descr$nobs), TRUE),
             plotxaxis = plotxaxis, 
             ...)
  }


forest.madad <- 
  function (x, IncludeSummary=NULL,SumLine=NULL,type = "sens", log = FALSE,         plotxaxis = TRUE, ...) 
  {
    dat<-x$data
    fitmod<-summary(reitsma(dat))$coefficients[3:4,c(1,5,6)]
    sensEstimate<-fitmod[1,]
    specEstimate<-1-fitmod[2,]
    d <- x
    stopifnot(class(d) == "madad")
    stopifnot(type %in% c("DOR", "sens", "spec", "posLR", "negLR"))
    if (type == "DOR") {
      ci <- d$DOR$DOR.ci
      x <- d$DOR$DOR
    }
    if (type == "sens") {
      if(!is.null(IncludeSummary)) ci <- rbind(d$sens$sens.ci,sensEstimate[2:3])
      if(!is.null(IncludeSummary)) x <- c(d$sens$sens,sensEstimate[1])
      if(is.null(IncludeSummary)) ci <- d$sens$sens.ci
      if(is.null(IncludeSummary)) x <- d$sens$sens
    }
    if (type == "spec") {
      if(!is.null(IncludeSummary)) ci <-  rbind(d$spec$spec.ci,specEstimate[c(3,2)])
      if(!is.null(IncludeSummary)) x <- c(d$spec$spec,specEstimate[1])
      if(is.null(IncludeSummary)) ci <- d$spec$spec.ci
      if(is.null(IncludeSummary)) x <- d$spec$spec
    }
    if (type == "posLR") {
      ci <- d$posLR$posLR.ci
      x <- d$posLR$posLR
    }
    if (type == "negLR") {
      ci <- d$negLR$negLR.ci
      x <- d$negLR$negLR
    }
    if (log) {
      x <- log(x)
      ci <- log(ci)
    }
    if(is.null(IncludeSummary)) polySum<-rep(FALSE, length(x))
    if(!is.null(IncludeSummary)) polySum<-c(rep(FALSE, length(x)-1),TRUE)
    myforest(x, ci, plotxaxis = plotxaxis, cipoly=polySum,SumLine=SumLine,         ...)
  }

And then I created the forest plot with a call such as:

forest.madad(madad(SingleSmall),
             SumLine=TRUE,
             type="sens",
             plotxaxis = FALSE,
             IncludeSummary=TRUE,
             snames=c(paste(c(rep("A",3),rep("E",16)),SingleCutpoint$Author),"Summary Estimate"),
             main="Sensitivity")
axis(1, at = seq(0 ,1, by = .25))

Where I indicate I don't want the function to plot the x axis, and then I call axis afterwards. Note I'm also requesting a summary statistic and a line at the summary estimate. You can ignore those arguments if you'd like.

Upvotes: 2

Views: 1546

Answers (1)

LyzandeR
LyzandeR

Reputation: 37879

Well apparently you cannot change those (or at least not immediately).

If you have a look at the source code of the forestmada function (this is the function that gets called when you use the forest function) you will see that you cannot. Run getAnywhere(forest.madad) to view the forest function and mada::forestmada to view forestmada function. Or just check ?forest.

So, this is what gets called inside the forestmada function:

   if (plotci) {
        xlim <- c(min(lb) - plotrange * 1.2, max(ub) + plotrange * 
            1.2)
    }
    else {
        xlim <- c(min(lb) - plotrange * 1.2, max(ub) + plotrange * 
            0.4)
    }

So even if you include any xlim or ylim arguments these will end up getting overwritten by the above.

Exactly the same with the at argument:

axis(1, at = round(seq(from = min(lb), to = max(ub), length.out = 5), 
    digits), cex = cex)

The only way would be to copy the source code and try to alter it yourself, but this will need a lot of work.

I tried myself a bit but apparently all of these arguments are used in a specific way later on, so whenever I changed anything I got a lot warnings and errors, and apparently this is why the author decided to have them calculated in a specific way.

Hope this helps.

Upvotes: 1

Related Questions