1

I am using the rdd (regression discontinuity design) package's DCdensity function and would like to change how the plot looks. I went to the help file for this function and found the following, under plot: "The user may wrap this function in additional graphical options to modify the plot." How can I go about doing this in practice?

I note that my question is the same as R plot options in rdd::DCdensity. However, the one answer to this old question is unsatisfactory to me because I do not want to change global plot options, but want to change them locally/specifically for each application (for example, changing the vertical line).

Here is a MWE (directly from the help file):

library(rdd)
x<-runif(1000,-1,1)
x<-x+2*(runif(1000,-1,1)>0&x<0)
DCdensity(x,0,plot=TRUE)

Here are two options I would like to implement:

Change the x-axis:

xlim=c(-0.5,0.5)

and add a vertical line at the cutoff: 

abline(v=0)
bill999
  • 2,147
  • 8
  • 51
  • 103
  • 2
    Copy anything needed from an older question. Explain _why_ it is unsatisfactory since it has an accepted answer. Read how to [edit] and then include a [MCVE]. – IRTFM Sep 14 '18 at 15:12

1 Answers1

2

I think with 'wrapping' they meant something like the following:

library(rdd)

myDCdensity <- function(runvar, cutpoint, my_abline = 0, my_title = "Default"){

  # get the default plot
  myplot <- DCdensity(runvar, cutpoint)

  # 'additional graphical options to modify the plot'
  abline(v = my_abline)
  title(main = my_title)

  # return
  return(myplot)
}

# run to verify
x <- runif(1000, -1, 1)
x <- x + 2 * (runif(1000, -1, 1) > 0 & x < 0)
myDCdensity(x, 0)
myDCdensity(x, 0, my_abline = 0.4, my_title = "Some other title")

You mentioned that you'd also like to set xlim. Here the situation is more complicated as this option is passed to the plot function and it seems you cannot modify it once the plot is drawn, so 'wrapping' won't help you. If you need to control xlim as well, the easiest is to write your own function based on the original code:

1) get the original function code with rdd::DCdensity (no parentheses)

2) modify setting xlim and save as your own function:

DCdensity2 <- function (runvar, cutpoint, bin = NULL, bw = NULL, verbose = FALSE, 
          plot = TRUE, ext.out = FALSE, htest = FALSE, my_xlim = c(-0.5,0.5)) # my_xlim param added
{
  runvar <- runvar[complete.cases(runvar)]
  rn <- length(runvar)
  rsd <- sd(runvar)
  rmin <- min(runvar)
  rmax <- max(runvar)
  if (missing(cutpoint)) {
  if (verbose) 
  cat("Assuming cutpoint of zero.\n")
  cutpoint <- 0
  }
  if (cutpoint <= rmin | cutpoint >= rmax) {
  stop("Cutpoint must lie within range of runvar")
  }
  if (is.null(bin)) {
  bin <- 2 * rsd * rn^(-1/2)
  if (verbose) 
  cat("Using calculated bin size: ", sprintf("%.3f", 
  bin), "\n")
  }
  l <- floor((rmin - cutpoint)/bin) * bin + bin/2 + cutpoint
  r <- floor((rmax - cutpoint)/bin) * bin + bin/2 + cutpoint
  lc <- cutpoint - (bin/2)
  rc <- cutpoint + (bin/2)
  j <- floor((rmax - rmin)/bin) + 2
  binnum <- round((((floor((runvar - cutpoint)/bin) * bin + 
  bin/2 + cutpoint) - l)/bin) + 1)
  cellval <- rep(0, j)
  for (i in seq(1, rn)) {
  cnum <- binnum[i]
  cellval[cnum] <- cellval[cnum] + 1
  }
  cellval <- (cellval/rn)/bin
  cellmp <- seq(from = 1, to = j, by = 1)
  cellmp <- floor(((l + (cellmp - 1) * bin) - cutpoint)/bin) * 
  bin + bin/2 + cutpoint
  if (is.null(bw)) {
  leftofc <- round((((floor((lc - cutpoint)/bin) * bin + 
  bin/2 + cutpoint) - l)/bin) + 1)
  rightofc <- round((((floor((rc - cutpoint)/bin) * bin + 
  bin/2 + cutpoint) - l)/bin) + 1)
  if (rightofc - leftofc != 1) {
  stop("Error occurred in bandwidth calculation")
  }
  cellmpleft <- cellmp[1:leftofc]
  cellmpright <- cellmp[rightofc:j]
  P.lm <- lm(cellval ~ poly(cellmp, degree = 4, raw = T), 
  subset = cellmp < cutpoint)
  mse4 <- summary(P.lm)$sigma^2
  lcoef <- coef(P.lm)
  fppleft <- 2 * lcoef[3] + 6 * lcoef[4] * cellmpleft + 
  12 * lcoef[5] * cellmpleft * cellmpleft
  hleft <- 3.348 * (mse4 * (cutpoint - l)/sum(fppleft * 
  fppleft))^(1/5)
  P.lm <- lm(cellval ~ poly(cellmp, degree = 4, raw = T), 
  subset = cellmp >= cutpoint)
  mse4 <- summary(P.lm)$sigma^2
  rcoef <- coef(P.lm)
  fppright <- 2 * rcoef[3] + 6 * rcoef[4] * cellmpright + 
  12 * rcoef[5] * cellmpright * cellmpright
  hright <- 3.348 * (mse4 * (r - cutpoint)/sum(fppright * 
  fppright))^(1/5)
  bw = 0.5 * (hleft + hright)
  if (verbose) 
  cat("Using calculated bandwidth: ", sprintf("%.3f", 
  bw), "\n")
  }
  if (sum(runvar > cutpoint - bw & runvar < cutpoint) == 0 | 
  sum(runvar < cutpoint + bw & runvar >= cutpoint) == 0) 
  stop("Insufficient data within the bandwidth.")
  if (plot) {
  d.l <- data.frame(cellmp = cellmp[cellmp < cutpoint], 
  cellval = cellval[cellmp < cutpoint], dist = NA, 
  est = NA, lwr = NA, upr = NA)
  pmin <- cutpoint - 2 * rsd
  pmax <- cutpoint + 2 * rsd
  for (i in 1:nrow(d.l)) {
  d.l$dist <- d.l$cellmp - d.l[i, "cellmp"]
  w <- kernelwts(d.l$dist, 0, bw, kernel = "triangular")
  newd <- data.frame(dist = 0)
  pred <- predict(lm(cellval ~ dist, weights = w, data = d.l), 
  interval = "confidence", newdata = newd)
  d.l$est[i] <- pred[1]
  d.l$lwr[i] <- pred[2]
  d.l$upr[i] <- pred[3]
  }
  d.r <- data.frame(cellmp = cellmp[cellmp >= cutpoint], 
  cellval = cellval[cellmp >= cutpoint], dist = NA, 
  est = NA, lwr = NA, upr = NA)
  for (i in 1:nrow(d.r)) {
  d.r$dist <- d.r$cellmp - d.r[i, "cellmp"]
  w <- kernelwts(d.r$dist, 0, bw, kernel = "triangular")
  newd <- data.frame(dist = 0)
  pred <- predict(lm(cellval ~ dist, weights = w, data = d.r), 
  interval = "confidence", newdata = newd)
  d.r$est[i] <- pred[1]
  d.r$lwr[i] <- pred[2]
  d.r$upr[i] <- pred[3]
  }
  plot(d.l$cellmp, d.l$est, lty = 1, lwd = 2, col = "black", # xlim set here based on the parameter
  type = "l", xlim = my_xlim, ylim = c(min(cellval[cellmp <= 
  pmax & cellmp >= pmin]), max(cellval[cellmp <= 
  pmax & cellmp >= pmin])), xlab = NA, ylab = NA, 
  main = NA)
  lines(d.l$cellmp, d.l$lwr, lty = 2, lwd = 1, col = "black", 
  type = "l")
  lines(d.l$cellmp, d.l$upr, lty = 2, lwd = 1, col = "black", 
  type = "l")
  lines(d.r$cellmp, d.r$est, lty = 1, lwd = 2, col = "black", 
  type = "l")
  lines(d.r$cellmp, d.r$lwr, lty = 2, lwd = 1, col = "black", 
  type = "l")
  lines(d.r$cellmp, d.r$upr, lty = 2, lwd = 1, col = "black", 
  type = "l")
  points(cellmp, cellval, type = "p", pch = 20)
  }
  cmp <- cellmp
  cval <- cellval
  padzeros <- ceiling(bw/bin)
  jp <- j + 2 * padzeros
  if (padzeros >= 1) {
  cval <- c(rep(0, padzeros), cellval, rep(0, padzeros))
  cmp <- c(seq(l - padzeros * bin, l - bin, bin), cellmp, 
  seq(r + bin, r + padzeros * bin, bin))
  }
  dist <- cmp - cutpoint
  w <- 1 - abs(dist/bw)
  w <- ifelse(w > 0, w * (cmp < cutpoint), 0)
  w <- (w/sum(w)) * jp
  fhatl <- predict(lm(cval ~ dist, weights = w), newdata = data.frame(dist = 0))[[1]]
  w <- 1 - abs(dist/bw)
  w <- ifelse(w > 0, w * (cmp >= cutpoint), 0)
  w <- (w/sum(w)) * jp
  fhatr <- predict(lm(cval ~ dist, weights = w), newdata = data.frame(dist = 0))[[1]]
  thetahat <- log(fhatr) - log(fhatl)
  sethetahat <- sqrt((1/(rn * bw)) * (24/5) * ((1/fhatr) + 
  (1/fhatl)))
  z <- thetahat/sethetahat
  p <- 2 * pnorm(abs(z), lower.tail = FALSE)
  if (verbose) {
  cat("Log difference in heights is ", sprintf("%.3f", 
  thetahat), " with SE ", sprintf("%.3f", sethetahat), 
  "\n")
  cat("  this gives a z-stat of ", sprintf("%.3f", z), 
  "\n")
  cat("  and a p value of ", sprintf("%.3f", p), "\n")
  }
  if (ext.out) 
  return(list(theta = thetahat, se = sethetahat, z = z, 
  p = p, binsize = bin, bw = bw, cutpoint = cutpoint, 
  data = data.frame(cellmp, cellval)))
  else if (htest) {
  structure(list(statistic = c(z = z), p.value = p, method = "McCrary (2008) sorting test", 
  parameter = c(binwidth = bin, bandwidth = bw, cutpoint = cutpoint), 
  alternative = "no apparent sorting"), class = "htest")
  }
  else return(p)
}

3) run to verify:

DCdensity2(x, 0)
DCdensity2(x, 0, my_xlim = c(-5, 5))
pieca
  • 2,463
  • 1
  • 16
  • 34