R/ci.R

In RMediation: Mediation Analysis Confidence Intervals

#' CI for a nonlinear function of coefficients estimates
#'
#' This function returns a (\eqn{1-\alpha})% confidence interval (CI) for a
#' well--defined nonlinear function of the coefficients in single--level and
#' multilevel structural equation models. The \code{ci} function uses the Monte
#' Carlo (\code{type="MC"}) and the asymptotic normal theory
#' (\code{type="asymp"}) with the multivariate delta standard error
#' (Asymptotic--Delta) method (Sobel, 1982) to compute a CI. In addition, for
#' each of the methods, when a user specifies \code{plot=TRUE} and
#' \code{plotCI=TRUE}, a plot of the sampling distribution of the quantity of
#' interest in the \code{quant} argument and an overlaid plot of the CI will be
#' produced. When \code{type="all"} and \code{plot=TRUE}, two overlaid plots of
#' the sampling distributions corresponding to each method will be produced;
#' when \code{plotCI=TRUE}, then the overlaid plots of the CIs for both methods
#' will be displayed as well.
#'
#' @param mu (1) a \link{vector} of means (e.g., coefficient estimates) for the
#'   normal random variables. A user can assign a name to each mean value, e.g.,
#'   \code{mu=c(b1=.1,b2=3)}; otherwise, the coefficient names are assigned
#'   automatically as follows: \code{b1,b2,...}. Or, (2) a \link{lavaan} object.
#' @param Sigma either a covariance matrix or a \link{vector} that stacks all
#'   the columns of the lower triangle variance--covariance matrix one
#'   underneath the other.
#' @param quant quantity of interest, which is a nonlinear/linear function of
#'   the model parameters. Argument \code{quant} is a \link{formula} that
#'   \strong{must} start with the symbol "tilde" (\code{~}): e.g.,
#'   \code{~b1*b2*b3*b4}. The names of coefficients must conform to the names
#'   provided in the argument \code{mu} or to the default names, i.e.,
#'   \code{b1,b2,...}.
#' @param alpha significance level for the CI. The default value is .05.
#' @param type method used to compute a CI. It takes on the values \code{"MC"}
#'   (default) for Monte Carlo, \code{"asymp"} for Asymptotic--Delta, or
#'   \code{"all"} that produces CIs using both methods.
#' @param plot when \code{TRUE}, plot the approximate sampling distribution of
#'   the quantity of interest using the specified method(s) in the argument
#'   \code{type}. The default value is \code{FALSE}. When \code{type="all"},
#'   superimposed density plots generated by both methods are displayed.
#' @param plotCI when \code{TRUE}, overlays a CI plot with error bars on the
#'   density plot of the sampling distribution of \code{quant}. When
#'   \code{type="all"}, the superimposed CI plots generated by both methods are
#'   added to the density plots. Note that to obtain a CI plot, one must also
#'   specify \code{plot="TRUE"}. The default value is \code{FALSE}.
#' @param n.mc Monte Carlo sample size. The default sample size is 1e+6.
#' @param H0 False. If \code{TRUE}, it will estimate the sampling distribution
#'   of \eqn{H_{0}:f(\bm b)=0}. See the arguments \code{mu0} and \code{Sigma0}.
#' @param mu0 a \link{vector} of means (e.g., coefficient estimates) for the
#'   normal random variables that satisfy the null hypothesis \eqn{H_{0}:f(\bm
#'   b)=0}. If it is not provided, smallest z value of \code{mu} is zet to zero.
#' @param Sigma0 either a covariance matrix or a \link{vector} that stacks all
#'   the columns of the lower triangle variance--covariance matrix one
#'   underneath the other. If it is not provided, then \code{Sigma} is used
#'   instead.
#' @param ... additional arguments.
#' @return When \code{type} is \code{"MC"} or \code{"asymp"}, \code{ci} returns
#'   a \link{list} that contains: \item{(\eqn{1-\alpha})% CI}{a vector of lower
#'   and upper confidence limits,} \item{Estimate}{a point estimate of the
#'   quantity of interest,} \item{SE}{standard error of the quantity of
#'   interest,} \item{MC Error}{When \code{type="MC"}, error of the Monte Carlo
#'   estimate.} When \code{type="all"}, \code{ci} returns a \link{list} of two
#'   objects, each of which a \link{list} that contains the results produced by
#'   each method as described above.
#' @keywords regression distribution
#'
#' @examples
#' ci(mu=c(b1=1,b2=.7,b3=.6, b4= .45), Sigma=c(.05,0,0,0,.05,0,0,.03,0,.03),
#' quant=~b1*b2*b3*b4, type="all", plot=TRUE, plotCI=TRUE)
#' #An Example of Conservative Null Sampling Distribution
#' ci(c(b1=.3,b2=.4,b3=.3), c(.01,0,0,.01,0,.02),
#' quant=~b1*b2*b3, type="mc", plot=TRUE, plotCI=TRUE, H0=TRUE, mu0=c(b1=.3,b2=.4,b3=0)  )
#' #An Example of Less Conservative Null Sampling Distribution
#' ci(c(b1=.3,b2=.4,b3=.3), c(.01,0,0,.01,0,.02),
#' quant=~b1*b2*b3, type="mc", plot=TRUE, plotCI=TRUE, H0=TRUE, mu0=c(b1=0,b2=.4,b3=0.1)  )
#' @author Davood Tofighi \email{dtofighi@@gmail.com}
#' @references  Tofighi, D. and MacKinnon, D. P. (2011). RMediation: An R
#'   package for mediation analysis confidence intervals. \emph{Behavior
#'   Research Methods}, \bold{43}, 692--700. \doi{doi:10.3758/s13428-011-0076-x}
#' @seealso \code{\link{medci}} \code{\link{RMediation-package}}
#' @export
#' @importFrom lavaan lav_matrix_vech_reverse
#' @importFrom MASS mvrnorm
#' @note A shiny web application for  Monte Carlo method of this function is available at \url{https://amplab.shinyapps.io/MEDMC/}

ci <- function(mu, Sigma, quant, alpha=0.05, type="MC", plot=FALSE, plotCI=FALSE, n.mc = 1e+06, H0=FALSE, mu0=NULL, Sigma0=NULL, ...){
  if(missing(mu) | is.null(mu) ) stop(paste("argument",sQuote("mu"), "must be specified"))
  if(missing(quant) | is.null(quant)) stop(paste("argument",sQuote("quant"), "must be specified"))
  if(!type%in%c("MC","mc","asymp","Asymp","all")) stop(paste("Please enter a ccorrect", sQuote("type"), "argument"))

  if(!isa(mu,"lavaan")){
  if(is.null(Sigma)| missing(Sigma) ) stop(paste("argument",sQuote("Sigma"), "cannot be a NULL value"))
  if(!is.matrix(Sigma)){
  if(length(mu)!= (sqrt(1 + 8 * length(Sigma)) - 1)/2) stop(paste("Please check the length of", sQuote("Sigma"),"and",sQuote("mu"),". If the length(dimension) of the", sQuote("mu"),"vector (",length(mu),") is correct, the stacked lower triangle matrix", sQuote("Sigma"), "must have ",((2*length(mu)+1)^2-1)/8, "elements, instead of", length(Sigma)) )
  Sigma <- lav_matrix_vech_reverse(Sigma) #converts to a symmetric matrix
  }
  if(is.null(names(mu)) ) names(mu) <- paste("b",1:length(mu), sep="") # if mu names is NULL

  if(!all(all.vars(quant) %in% names(mu))) stop(paste("The parameters names in formula", sQuote("quant"), "must match the parameters names provided in", sQuote("mu"),"."))
  }
  else{
    fm1 <- mu
#    if(!fm1@Options[['do.fit']]) fm1 <- update(mu, do.fit=TRUE) # If it's not fitted, refit.
    pEstM1 <- coef(fm1) #parameter estimate
    name1 <- all.vars(quant)
    if(!all(name1 %in% names(pEstM1))) stop(paste("The parameters names in formula", sQuote("quant"), "must match the parameters names provided in lavaan object."))
    mu <- pEstM1[name1]
    ##Cov
    covM1 <- (vcov(fm1)) #covariance of the coef estimates
    Sigma <- covM1[name1,name1]
  }

  if(length(mu)*n.mc > .Machine$integer.max){
    n.mc <- 1e6
    warning(paste("n.mc is too large. It is reset to", n.mc))
  }
  # This line sets the plot margins and other plot parameters
  if(plot){
    op <- par(mar = c(7, 4, 7, 4) + 0.1, xpd=TRUE, ask=FALSE)
  }

  if(type %in% c("mc", "MC")) res <- confintMC(mu=mu, Sigma=Sigma, quant=quant, alpha=alpha, plot=plot, plotCI=plotCI, n.mc = n.mc, H0=H0, mu0=mu0, Sigma0=Sigma0,... )
  if(type %in% c("asymp", "Asymp")) res <- confintAsymp(mu=mu, Sigma=Sigma, quant=quant, alpha=alpha, plot=plot, plotCI=plotCI )
  if(type %in% "all") {
    res1 <- confintMC(mu=mu, Sigma=Sigma, quant=quant, type=type, alpha=alpha, plot=plot, plotCI=plotCI, n.mc = n.mc)
    res2 <- confintAsymp(mu=mu, Sigma=Sigma, quant=quant, type=type, alpha=alpha, plot=plot, plotCI=plotCI)
    res <- list(MC=res1, Asymptotic=res2)
    }
  if(plot) {
    #par(mar = c(5, 4, 4, 2) + 0.1, xpd=FALSE)
    par(op)
}
  return(res)
}