Creates marginal a effect ribbon plot from games (http://cran.r-project.org/web/packages/games/games.pdf) model objects

MarginalGames.R

#' Creates a marginal effect or predicted probability ribbon plot for interaction effects estimated from games model objects
#' 
#' \code{MarginalGames} marginal effect or predicted probability ribbon plot for interaction effects estimated from \code{\link{games}} model objects
#' 
#' @param obj fitted model object created by \code{\link{egame12}}.
#' @param b1 character string of the first constitutive variable's name. 
#' @param b2 character string of the second constitutive variable's name.
#' @param b12 character string of the interaction term's name
#' @param X1 numeric fitted value for \code{b1} to simulate for. Note: can only be set if \code{type = "prob"} and can only have one value.
#' @param X2 numeric vector of fitted values of \code{b2} to simulate for. 
#' @param type character specifying the type of quantity of interest to simulate. If \code{type = "marginal"} the marginal effect of a one unit increase in a variable \code{X1} is found. If \code{type = "prob"} then a predicted probability for \code{b1*X1 + b2*X2 + b12*X1*X2} is found. All other variables are fitted at 0.
#' @param ci the proportion of middle simulations to keep. The default is \code{ci = 0.95}, i.e. keep the middle 95 percent. If \code{spin = TRUE} then \code{ci} is the confidence level of the shortest probability interval. Any value from 0 through 1 may be used.
#' @param spin logical, whether or not to keep only the shortest probability interval rather than the middle simulations.
#' @param nsim the number of simulations to run per value of X2. Default is \code{nsim = 1000}.
#' @param colour character string colour of the simulated points. Default is hexadecimal colour \code{colour = '#A6CEE3'}. 
#' @param alpha point alpha (e.g. transparency). Default is \code{alpha = 0.05}. See \code{\link{ggplot2}}.
#' #' @param lsize size of the smoothing line. Default is 2. See \code{\link{ggplot2}}.
#' @param xlab a label for the plot's x-axis.
#' @param ylab a label of the plot's y-axis. 
#' @param title the plot's main title.
#' 
#' @details Note: MarginalGames assumes that the \code{\link{games}} model used the probit link function.
#' 
#' @source
#' Curtis S. Signorino and Brenton Kenkel (2013). games: Statistical Estimation of Game-Theoretic Models. R package version 1.1-0. \url{http://CRAN.R-project.org/package=games}
#' 
#' Brambor, Thomas, William Roberts Clark, and Matt Golder. 2006. ”Understanding Interaction Models: Improving Empirical Analyses.” Political Analysis 14(1): 63-82.
#' 
#' Christopher Gandrud (2013). simPH: Tools for simulating and plotting quantities of interest estimated from Cox Proportional Hazards models. R package version 0.7.5. \url{http://christophergandrud.github.io/simPH/}


MarginalGames <- function(obj, b1, b2, b12, X1, X2, type = "marginal", ci = 0.95, spin = FALSE, nsim = 1000, colour = "#A6CEE3", alpha = 0.5, lsize = 1, xlab = "", ylab = "", title = ""){
  # Load required packages
  require(MSBVAR)
  require(VGAM)
  require(stringr)
  require(simPH)
  require(plyr)
  require(ggplot2)

  typeOpts <- c("marginal", "prob")
  TesttypeOpts <- type %in% typeOpts
  if (!isTRUE(TesttypeOpts)){
    stop("Invalid type. type must be 'marginal' or 'prob'.")
  }
  
  # Parameter estimates & Variance/Covariance matrix
  Coef <- matrix(obj$coefficients)
  VC <- vcov(obj)
  
  # Draw covariate estimates from the multivariate normal distribution      
  Drawn <- rmultnorm(n = nsim, mu = Coef, vmat = VC)
  DrawnDF <- data.frame(Drawn)
  dfn <- names(DrawnDF)
  DropPatterns <- c("u1.1..", "u1.2..", "u2.3..")
  for (i in DropPatterns){
    dfn <- gsub(pattern = i, replacement = "", x = dfn)
  }
  
  bs <- c(b1, b2, b12)
  bpos <- match(bs, dfn)
  
  # Subset data frame to only include interaction constitutive terms and
  Simb <- data.frame(Drawn[, bpos])	
  
  # Calculate quantity of interest
  if (type == "marginal"){
    message("X1 ignored.")
    X2df <- data.frame(X2)
    names(X2df) <- "X2"
    Simb <- merge(Simb, X2df)
    Simb$QIBefore <- Simb[, 1] + (Simb[, 3] * Simb[, 4])
    Simb$QI <- dnorm(Simb$QIBefore)
  }
  else if (type == "prob"){
    if (is.null(X1) | is.null(X2)){
      stop("Both X1 and X2 must be set")
    }
    else if (length(X1) != 1){
      stop("Xl can only have one fitted value.")
    }
    else{
      message("All variables other than X1 and X2 fitted at 0.")
      X12df <- data.frame(X1, X2)
      names(X12df) <- c("X1", "X2")
      Simb <- merge(Simb, X12df)
      Simb$Comb <- Simb[, 4] * Simb[, 5]
      Simb$QIBefore <- (Simb[, 1]*Simb[, 4]) + (Simb[, 2]*Simb[, 5]) + (Simb[, 3]*Simb[, 6])
      Simb$QI <- probit(Simb$QIBefore, inverse = TRUE)
    }
  }
  # Find middle interval
  SimbPerc <- simPH:::IntervalConstrict(Simb = Simb, SubVar = "X2", 
                                        qi = "Marginal Effect",
                                        QI = QI, spin = spin, ci = ci)	
  
  # Shrink data frame to only necessary variables
  objdf <- data.frame(SimbPerc$X2, SimbPerc$QI)
  names(objdf) <- c("Xj", "QI")
  
  # Find ribbon lines
  obj <- simPH:::MinMaxLines(df = objdf)
  ggplot(obj, aes(Xj, Median)) +
    geom_line(size = lsize, colour = colour) +
    geom_ribbon(aes(ymin = Lower50, ymax = Upper50), alpha = alpha, fill = colour) +
    geom_ribbon(aes(ymin = Min, ymax = Max), alpha = alpha, fill = colour) +
    geom_hline(aes(yintercept = 0), linetype = "dotted") +
    xlab(xlab) + ylab(ylab) +
    ggtitle(title) +
    guides(colour = guide_legend(override.aes = list(alpha = 1))) +
    theme_bw(base_size = 15)
}