如何“在附加图形选项中包裹 [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屋！