R/ci.R

Defines functions ci

Documented in ci

#' 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, 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 satisft 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
#' @export
#' @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{[email protected]@gatech.edu} and David P. MacKinnon \email{[email protected]@asu.edu}
#' @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:10.3758/s13428-011-0076-x
#' @seealso \code{\link{medci}} \code{\link{RMediation-package}}

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(class(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 margines 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)
}

Try the RMediation package in your browser

Any scripts or data that you put into this service are public.

RMediation documentation built on May 2, 2019, 9:41 a.m.