R/MCres_comb.R

In koohyun-kwon/rdadapt:

#' Specification Data Frame Generation
#'
#' Generates all possible combinations of specifications.
#'
#' @param spec a vector of strings, containing specification information.
#'
#' @return \code{data.frame} object containing binary indicators for all possible
#' combinations of specifications.
#' @export
#'
#' @examples spec.grid(c("unif", "homo"))
spec.grid <- function(spec){

  spec.list <- list()
  spec.len <- length(spec)

  for(i in 1:spec.len){

    spec.list[[spec[i]]] <- c(1, 0)
  }

  res <- expand.grid(spec.list)
  return(res)
}

#' Regression Function Evaluation
#'
#' Calculates the regression function value given a specification.
#'
#' When \code{which.fun} is 5 or 6, \code{x} should be a matrix with two columns.
#'
#' @param which.fun a positive integer; curretly supports 1:6.
#' @inheritParams f.1
#'
#' @return vector of evaluated values with length \code{length(x)}
#' @export
#'
#' @examples f.eval(2, runif(100,-1,1), 1, 1)
f.eval <- function(which.fun, x, C, th){

  if(which.fun == 1){

    res <- f.1(x, C, th)
  }else if(which.fun == 2){

    res <- f.2(x, C, th)
  }else if(which.fun == 3){

    res <- f.3(x, C, th)
  }else if(which.fun == 4){

    res <- f.4(x, C, th)
  }else if(which.fun == 5){

    res <- f.5(x, C, th)
  }else if(which.fun == 6){

    res <- f.6(x, C, th)
  }

  return(res)
}

#' Data Generation for a Given Specification
#'
#' Generates running variables, and corresponding standard deviations and
#' regression function values given a specification.
#'
#' @param spec.ind \code{data.frame} with a single rows, with \code{names(spec.ind)} containing
#' \code{c("unif", "homo", "trueCsmall")}.
#' @param which.fun which regression function to evaluate; supports \code{1:6}.
#' @param n sample size.
#' @param th true RD parameter.
#' @param C.big the value of C corresponding to the large C.
#' @param C.small the value of C corresponding to the small C.
#' @param sig.scale scaling parameter for the standard deviation.
#' @param seed the random number seed used to generate regressors; the default is
#' \code{seed = 1}.
#'
#' @return a list with components \code{x} (running variables),
#' \code{sig} (standard deviations), and \code{fx} (regression function values).
#' @export
#'
#' @examples spec.ind <- data.frame(unif = 1, homo = 0, trueCsmall = 0)
#' spec.data(spec.ind, 1, 100, 1, 3, 1, 1/2)
#' spec.data(spec.ind, 6, 100, 1, 3, 1, 1/2)
spec.data <- function(spec.ind, which.fun, n, th, C.big, C.small, sig.scale, seed = 1){

  d <- ifelse(which.fun <= 4, 1, 2) # reg. functions 1 ~ 4 correspond to 1-dim
  set.seed(1)

  if(d == 1){
    if(spec.ind$unif == 1){
      x <- stats::runif(n, -1, 1)
    }else{
      x <- 2 * (stats::rbeta(n, 2, 2)) - 1
    }
  }else if(d == 2){
    if(spec.ind$unif == 1){
      x <- cbind(stats::runif(n, -1, 1), stats::runif(n, -1, 1))
    }else{
      x <- cbind(2 * (stats::rbeta(n, 2, 2)) - 1, 2 * (stats::rbeta(n, 2, 2)) - 1)
    }
  }

  if(spec.ind$homo == 1){
    sig <- sig.scale * rep(1, n)
  }else{
    if(d == 1){
      sig <- sig.scale * stats::dnorm(x) / stats::dnorm(0)
    }else if(d == 2){
      sig <- sig.scale * mvtnorm::dmvnorm(x) / mvtnorm::dmvnorm(c(0, 0))
    }
  }

  if(spec.ind$trueCsmall == 1){
    trueC <- C.small
  }else{
    trueC <- C.big
  }

  fx <- f.eval(which.fun, x, trueC, th)

  res <- list(x = x, sig = sig, fx = fx)
  return(res)
}

#' Data Generation for All Specifications
#'
#' Generates running variables, and corresponding standard deviations and
#' regression function values for all specifications.
#'
#' @param spec.all \code{data.frame} object created by \code{spec.grid} function.
#' @inheritParams spec.data
#'
#' @return a list with \code{nrow(spec.all)} components, where each of them contains
#' data corresponding to each specification.
#' @export
#'
#' @examples spec.all <- spec.grid(c("unif", "homo", "trueCsmall"))
#' spec.data.all(spec.all, 4, 100, 1, 3, 1, 1/2)
#' spec.data.all(spec.all, 5, 100, 1, 3, 1, 1/2)
spec.data.all <- function(spec.all, which.fun, n, th, C.big, C.small, sig.scale,
                          seed = 1){

  res.list <- list()
  s.len <- nrow(spec.all)

  for(i in 1:s.len){

    spec.ind <- spec.all[i, ]
    res.list[[i]] <- spec.data(spec.ind, which.fun, n, th, C.big, C.small, sig.scale,
                               seed)
  }

  return(list(data = res.list, spec.all = spec.all))

}


#' CI Length and Coverage Result Generation
#'
#' Generates CI length and coverage results for each DGP specification.
#'
#' This function is used for MC simulations.
#'
#' @param specdata.list a list generated by \code{spec.data.all} function.
#' @param eps N(0,1) random noises.
#' @param CI.method method to be used to generate the CI.
#' @param th true RD parameter value.
#' @param spec.etc list of other specification parameters.
#'
#' @return a list with two components; \code{len} contains the length result,
#' and \code{cov} contains the coverage indicator.
#' @export
#'
#' @examples spec.etc <- list(alpha = 0.05, mon_ind = 1,
#' C.small = 1, C.large = 3, se.method.RDH = "nn", se.initial.RDH = "nn", M.RDH = 1,
#' se.method = "nn.test", se.init = "S.test", t.dir = "left", C.l = 1/2, C.u = 2, C = 3)
#' spec.all <- spec.grid(c("unif", "homo", "trueCsmall"))
#' specdata.list <- spec.data.all(spec.all, 4, 100, 1, 3, 1, 1/2)
#' eps <- rnorm(100)
#' spec.res.gen(specdata.list, eps, "RDH", 1, spec.etc)
#' spec.res.gen(specdata.list, eps, "RDR", 1, spec.etc)
#' spec.res.gen(specdata.list, eps, "adpt", 1, spec.etc)
#' spec.res.gen(specdata.list, eps, "mm.largeC", 1, spec.etc)
#' spec.res.gen(specdata.list, eps, "RDR.L", 1, spec.etc)
#' spec.etc$mon_ind <- c(1, 2)
#' specdata.list <- spec.data.all(spec.all, 5, 100, 1, 3, 1, 1/2)
#' spec.res.gen(specdata.list, eps, "mm.largeC", 1, spec.etc)
spec.res.gen <- function(specdata.list, eps, CI.method, th, spec.etc){

  s.len <- nrow(specdata.list$spec.all) # Number of simuation designs
  res.cov <- numeric(s.len)
  res.len <- numeric(s.len)
  res.sd <- numeric(s.len)

  for(i in 1:s.len){

    data.i <- specdata.list$data[[i]]

    x <- data.i$x
    sig <- data.i$sig

    if(is.vector(x)){
      d <- 1
      x <- matrix(x, ncol = 1)
    }else if(is.matrix(x)){
      d <- ncol(x)
    }

    if(d == 1){
      t.ind <- x[, 1] < 0
    }else if(d == 2){
      t.ind <- x[, 1] < 0 & x[, 2] < 0
    }

    xt <- x[t.ind, ]
    xc <- x[!t.ind, ]
    sig.t <- sig[t.ind]
    sig.c <- sig[!t.ind]

    fx <- data.i$fx
    y <- fx + eps * sig
    yt <- y[t.ind]
    yc <- y[!t.ind]

    ci.res <- CI.gen(CI.method, xt, xc, yt, yc, sig.t, sig.c, spec.etc)
    ci.l <- ci.res$ci[1]
    ci.u <- ci.res$ci[2]

    if(ci.u == Inf){
      res.len[i] <- th - ci.l # excess length of a one-sided confidence interval
    }else{
      res.len[i] <- ci.u - ci.l
    }

    res.cov[i] <- (th < ci.u) * (th > ci.l)

    res.sd[i] <- ci.res$sd
  }

  res <- list(len = res.len, cov = res.cov, sd = res.sd)
  return(res)
}

#' Length Generation for One-sided CIs
#'
#' Calculate excess lengths of one-sided CIs for different values of true Lipschitz coefficients.
#'
#' Modification for multi-dim hasn't been made yet.
#'
#' @inheritParams spec.res.gen
#' @inheritParams spec.data
#' @param which.spec which design is used; the default is \code{which.spec = 1}.
#'
#' @return vector of excess lengths and coverage indicators corresponding to
#' true Lipschitz coefficient values contained in \code{spec.etc$C.true.vec}.
#' @export
#'
#' @examples spec.etc <- list(alpha = 0.05, mon_ind = 1,
#' C.small = 1, C.large = 3, se.method.RDH = "nn", se.initial.RDH = "nn", M.RDH = 1,
#' se.method = "nn", se.init = "Silverman", t.dir = "left", C.l = 1/2, C.u = 2, C = 3,
#' C.true.vec = seq(1/2, 3, length.out = 2), adpt.orc = TRUE)
#' n <- 100
#' th <- 1
#' eps <- rnorm(n)
#' which.fun <- 1
#' sig.scale <- 1/2
#' CI.method <- "adpt.one"
#' adpt.res.gen(eps, CI.method, th, spec.etc, which.fun, n, sig.scale)
#' spec.etc$adpt.orc <- FALSE
#' adpt.res.gen(eps, CI.method, th, spec.etc, which.fun, n, sig.scale)
adpt.res.gen <- function(eps, CI.method, th, spec.etc, which.fun, n,
                         sig.scale, which.spec = 1){

  spec.all <- spec.grid(c("unif", "homo", "trueCsmall"))
  spec <- spec.all[which.spec, ]

  adpt.orc <- spec.etc$adpt.orc

  C.true.vec <- spec.etc$C.true.vec
  c.len <- length(C.true.vec)
  res.len <- numeric(c.len)
  res.cov <- numeric(c.len)

  for(i in 1:c.len){

    C.true <- C.true.vec[i]
    obs.data <- spec.data(spec, which.fun, n, th, C.true, C.true, sig.scale)

    x <- obs.data$x
    xt <- matrix(x[x < 0], ncol = 1) # only allow matrix input
    xc <- matrix(x[x >= 0], ncol = 1)

    sig <- obs.data$sig
    sig.t <- sig[x < 0]
    sig.c <- sig[x >= 0]

    fx <- obs.data$fx
    y <- fx + eps * sig
    yt <- y[x < 0]
    yc <- y[x >= 0]

    if(adpt.orc == T){
      C.one <- C.true
    }else{
      C.one <- C.true.vec[length(C.true.vec)]
    }

    ci.res <- CI.gen(CI.method, xt, xc, yt, yc, sig.t, sig.c, spec.etc, C.one)
    ci.l <- ci.res$ci[1]

    res.len[i] <- th - ci.l
    res.cov[i] <- (th > ci.l)
  }

  return(list(len = res.len, cov = res.cov))
}