R/confintAsymp.R

In quantPsych/RMediation: Mediation Analysis Confidence Intervals

confintAsymp <- function(mu, Sigma, quant = NULL, alpha = 0.05, type = "asymp", plot = FALSE, plotCI = FALSE) {
  # Sigma <- vech.reverse(Sigma)
  q1 <- quant
  fx <- deriv(quant, names(mu), func = TRUE)
  muList <- as.list(mu)
  grad1 <- do.call("fx", muList)
  grad <- as.vector(attr(grad1, "gradient"))
  quantSE <- sqrt(t(grad) %*% Sigma %*% grad)
  quant <- parse(text = sub("~", "", quant))
  quantMean <- eval(quant, muList)
  q.a <- qnorm(alpha / 2, 0, 1)
  CI <- c(quantMean + q.a * quantSE, quantMean - q.a * quantSE)
  res <- list(CI, quantMean, quantSE)
  names(res) <- c(paste((1 - alpha / 2) * 100, "% ", "CI", sep = ""), "Estimate", "SE")
  attr(res, "quant") <- q1
  outer <- FALSE # outer position
  mcex <- .8

  if (plot & quantSE > 40 * .Machine$double.eps) {
    if (type == "all") {
      f <- Vectorize(function(x) dnorm(x, quantMean, quantSE), "x")
      curve(f, lwd = 2, col = "blue", lty = 2, add = TRUE)
      if (plotCI) {
        max1 <- par("usr")[2]
        min1 <- par("usr")[1]
        yci <- par("usr")[3] - .1 * diff(par("usr")[3:4]) / 25
        arrows(CI[1], yci, CI[2], yci, length = 0, angle = 90, code = 3, cex = 1.5, lwd = 2, col = "blue", lty = 2)
        points(quantMean, yci, pch = 19, cex = 1.5, col = "blue")
        # New
        mtext(paste("LL=", round(CI[1], 3)), side = 3, line = 1, outer = outer, at = max1 - 1 * (max1 - min1) / 9, col = "blue", cex = mcex, adj = 0)
        mtext(paste("UL=", round(CI[2], 3)), side = 3, line = 2, outer = outer, at = max1 - 1 * (max1 - min1) / 9, col = "blue", cex = mcex, adj = 0)
        mtext("Asymptotic Normal", side = 3, line = 3, outer = outer, at = max1 - 1 * (max1 - min1) / 9, col = "blue", cex = mcex, adj = 0)
      }
    }

    if (type %in% c("asymp", "Asymp")) {
      f <- Vectorize(function(x) dnorm(x, quantMean, quantSE), "x")
      curve(f, from = quantMean - 2.58 * quantSE, to = quantMean + 2.58 * quantSE, lwd = 2, xlab = "", ylab = "", axes = FALSE)
      range1 <- c(quantMean - 2.58 * quantSE, quantMean + 2.58 * quantSE)
      max1 <- range1[2]
      min1 <- range1[1]
      xrange <- seq(min1, max1, length = 7)
      xrange <- pretty(xrange, n = 9)
      mtext(quant, 1, 4.5)
      axis(1, xrange, xrange, line = 2.5)
      axis(2, line = 1.1)
      mtext(paste("LL=", round(CI[1], 3)), side = 3, line = 1, outer = outer, at = max1 - 4 * (max1 - min1) / 9, col = "black", cex = mcex)
      mtext(paste("UL=", round(CI[2], 3)), side = 3, line = 2, outer = outer, at = max1 - 4 * (max1 - min1) / 9, col = "black", cex = mcex)
      mtext("Asymptotic", side = 3, line = 3, outer = outer, at = max1 - 4 * (max1 - min1) / 9, col = "black", cex = mcex)

      if (plotCI) {
        yci <- par("usr")[3] - 1.2 * diff(par("usr")[3:4]) / 25
        arrows(CI[1], yci, CI[2], yci, length = 0, angle = 90, code = 3, cex = 1.5, lwd = 2)
        points(quantMean, yci, pch = 19, cex = 1.5)
      } # if
    }
  } # if

  if (plot & quantSE <= 40 * .Machine$double.eps) {
    if (type %in% c("asymp", "Asymp")) {
      #     range1 <- c(quantMean-2.58*quantSE, quantMean+2.58*quantSE)
      #     max1 <- range1[2]
      #     min1 <- range1[1]
      #     xrange<-round(seq(min1,max1,length=7),1)
      #     axis(1,c(xrange, round(c(quantMean,CI),2)));axis(2)
      #     smidge <- par("cin")*abs(par("tcl"))
      #     text(max1-0*(max1-min1)/9,(par("usr")[4]),adj=0, bquote(mu== .(round(quantMean,3)) ))
      #     text(max1-0*(max1-min1)/9,(par("usr")[4]-1*par("cxy")[2]),adj=0, bquote(sigma== .(round(quantSE,3)) ))
      cat("
        \n *******************************************
Note: Because the SE =", quantSE, ", the plot of the sampling distribution using aymptotic normal is not available.
        \n *******************************************\n")
    }

    if (type == "all") {
      max1 <- par("usr")[2]
      min1 <- par("usr")[1]
      yci <- par("usr")[3] + 2.5 * diff(par("usr")[3:4]) / 25
      # arrows(CI[1],yci,CI[2],yci,length=0,angle=90,code=3,cex=1.5,lwd=2, col="blue",lty=2)
      # points(quantMean,yci,pch=19,cex=1.5, col="blue")
      #       text(max1+0*(max1-min1)/9,(par("usr")[4]-2.5*par("cxy")[2]),adj=0, paste("LL=",round(CI[1],3)), col="blue")
      #       text(max1+0*(max1-min1)/9,(par("usr")[4]-3.5*par("cxy")[2]),adj=0, paste("UL=",round(CI[2],3)), col="blue")
      mtext(paste("LL=", round(CI[1], 3)), side = 3, line = 2, outer = outer, at = max1 - 1 * (max1 - min1) / 9, col = "blue", cex = mcex)
      mtext(paste("UL=", round(CI[2], 3)), side = 3, line = 3, outer = outer, at = max1 - 1 * (max1 - min1) / 9, col = "blue", cex = mcex)
      mtext("Asymptotic", side = 3, line = 4, outer = outer, at = max1 - 1 * (max1 - min1) / 9, col = "blue", cex = mcex)
      cat("Because the SE,", quantSE, ", is near zero, the plot of the sampling distribution using aymptotic normal is not available.\n")
    }
  }

  return(res)
}