R/sim_intervals.R

In FrankLef/fciR: Data and functions used with Fundamentals of Causal Inference by Babette Brumback

#' Simulate a sampling distribution
#'
#' Simulate a sampling distribution.
#'
#' Simulate a sampling distribution using base R. Section 2.4.
#'
#' @param nsim Nb of replicates.
#' @param n Nb of observations.
#' @param seed Integer, the seed used for random numbers.
#'
#' @importFrom stats rnorm var
#'
#' @return List with percentage of bad and good.
#'
#' @examples
#' \dontrun{
#' sim_intervals()
#' }
#' @export
sim_intervals <- function(nsim = 1e4, n = 1e4, seed = 1009) {
  # get the statistics from the 2 samplings
  set.seed(seed)
  out <- replicate(n = nsim, expr = {
    y <- rnorm(n)
    x <- rnorm(n)
    muhaty <- mean(y)  # mean of the y sample
    sehaty <- sqrt(var(y) / n)  # standard error of the y sample
    lci <- muhaty - 1.96 * sehaty  # lower bound of the y sample
    uci <- muhaty + 1.96 * sehaty  # upper bound of the y sample
    muhatx <- mean(x)  # mean of the x sample
    # is mean(x) within the interval from y?
    bad <- (muhatx >= lci) & (muhatx <= uci)
    # is the true mean of 0 within the interval from y?
    good <- (0 >= lci) & (0 <= uci)
    # return the results
    c("bad" = bad, "good" = good)
  })
  # convert to data.frame
  out <- data.frame(t(out))
  pbad <- sum(out$bad) / nrow(out)  # proportion in bad interval
  pgood <- sum(out$good) / nrow(out)  # proportion in good interval
  # return results in a list
  list("bad" = pbad, "good" = pgood)
}

#' @rdname sim_intervals
#' @examples
#' \dontrun{
#' sim()
#' }
#' @export
sim <- sim_intervals