如何“在其他图形选项中包装[DCdensity]以修改图"? rdd包 [英] How to "wrap [DCdensity] in additional graphical options to modify the plot"? rdd package

问题描述

我正在使用rdd(回归不连续性设计)包的DCdensity函数，并想更改图形的外观.我转到了该功能的帮助文件，并在plot下找到了以下内容:用户可以将此功能包装在其他图形选项中，以修改图."在实践中我该如何去做?

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?

我注意到我的问题与rdd :: DCdensity中的 R绘图选项相同.但是，这个老问题的一个答案对我来说并不令人满意，因为我不想更改全局绘图选项，而是想针对每个应用程序在本地/特定地更改它们(例如，更改垂直线).

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).

这是MWE(直接来自帮助文件):

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:

更改x轴:

xlim=c(-0.5,0.5) 

并在截止处添加一条垂直线:

and add a vertical line at the cutoff:

abline(v=0)

推荐答案

我认为包装"的含义如下:

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")

您提到您还想设置xlim.由于此选项传递给plot函数和xlim，最简单的方法就是根据原始代码编写自己的函数:

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)用rdd::DCdensity获取原始功能代码(不带括号)

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

2)修改设置xlim并另存为自己的功能:

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)进行验证:

DCdensity2(x, 0)
DCdensity2(x, 0, my_xlim = c(-5, 5))

这篇关于如何“在其他图形选项中包装[DCdensity]以修改图"? rdd包的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！