R/OneWaySA.R

In feralaes/dampack: Package with useful decision-analytic modeling tools

#' One-way sensitivity analysis using linear regression metamodeling
#'
#' This function displays a one-way sensitivity analysis (OWSA) graph 
#' by estimating a linear regression metamodel of a PSA for a given
#' decision-analytic model
#' @param strategies String vector with the name of the strategies
#' @param x Matrix with the model inputs or parameters
#' @param y Matrix with the model outputs
#' @param parm String with the name of the parameter of interest
#' @param range Range of the parameter of interest. Default: NULL range. If
#' range=NULL, the 2.5th and 9.75th percentile of the parameter
#' of interest will be used as lower and upper bounds of the 
#' range, respectively.
#' @param poly.order Order of polynomial for the linear regression metamodel. 
#' Default: 2
#' @param txtsize Font size for ggplot graph. Default: 12
#' @keywords one-way sensitivity analysis; linear regression metamodel
#' @return owsa A `ggplot` object with the OWSA graph of the `parm` on the 
#' outcome of interest
#'
OneWaySA <- function(strategies, y, x, 
                     parm, range = NULL,
                     poly.order = 2,
                     txtsize = 12){
  # Extract parameter column number in x matrix
  par.col <- which(colnames(x)==parm)
  dep <- length(strategies) #Number of dependent variables, i.e., strategies outcomes
  indep <- ncol(x) #Number of independent variables, i.e., parameters
  Sim <- data.frame(y,x)
  #Determine range of of the parameer to be plotted
  if (is.null(range)){ #If user does not define a range
    #Default range given by the domain of the parameter's sample
    #vector to define 400 samples between the 2.5th and 97.5th percentiles
    percentiles = seq(2.5, 97.5, length = 400) 
    j = round(percentiles*(length(x[,par.col])/100)) #indexing vector;j=round(y*n/100) where n is the size of vector of interest
    vector<-sort(x[j, par.col]) 
  }
  else{ #If user defines a range  
    vector <- seq(range[1],range[2],length.out=400)
  }
  
  #Generate a formula by pasting column names for both dependent and independent variables. Imposes a 1 level interaction
  f <- as.formula(paste('cbind(',paste(colnames(Sim)[1:dep], collapse=','), 
                        ') ~ (','poly(', parm,',', poly.order,',raw=TRUE)+',
                        paste(colnames(x)[-par.col], collapse='+'),')'))
  #Run Multiple Multivariate Regression (MMR) Metamodel
  Oway.mlm = lm(f, data=Sim)
  
  # Create data frame with all combinations between both parameters of interest
  OWSA <- data.frame(parm = vector)
  
  #Generate matrix to use for prediction 
  Sim.fit <- matrix(rep(colMeans(x)),
                    nrow = length(vector), 
                    ncol = ncol(x), 
                    byrow = T)
  Sim.fit[, par.col] <- OWSA[, 1]
  # Transform to data frame, the format required for predict
  Sim.fit <- data.frame(Sim.fit)
  # Name data frame's columns with parameters' names
  colnames(Sim.fit) <- colnames(x) #Name data frame's columns with parameters' names
  
  # Predict Outcomes using MMMR Metamodel fit
  Sim.OW = data.frame(predict(Oway.mlm, newdata = Sim.fit))
  colnames(Sim.OW) <- strategies #Name the predicted outcomes columns with strategies names
  
  #Reshape dataframe for ggplot
  Sim.OW = stack(Sim.OW, select=strategies) #
  Sim.OW = cbind(Sim.fit, Sim.OW) #Append parameter's dataframe to predicted outcomes dataframe
  
  #A simple trick I use to define my variables in my functions environment
  #Borrowed from http://stackoverflow.com/questions/5106782/use-of-ggplot-within-another-function-in-r
  Sim.OW$parm<-Sim.OW[,parm];
  
  owsa <- ggplot(data = Sim.OW, aes(x = parm, y = values, color = ind)) +
    geom_line() +
    ggtitle("One-way sensitivity analysis") + #\n Net Health Benefit
    xlab(parm) +
    ylab("E[Outcome]") +
    scale_colour_hue("Strategy", l=50) +
    scale_x_continuous(breaks=number_ticks(6)) + #Adjust for number of ticks in x axis
    scale_y_continuous(breaks=number_ticks(6)) +
    theme_bw(base_size = 14) +
    theme(legend.position = "bottom")
  return(owsa)
}