R/deltamethod.R

In symSEM: Symbolic Computation for Structural Equation Models

## fn <- "log(a)^b+c"
## fn <- c("x^2 + log(y)*z", "x+y+sin(z)")

#' Compute the Variance-Covariance Matrix of Functions using the first-order
#' Delta Method
#'
#' It computes the variance-covariance matrix of functions using the
#' first-order delta method.
#'
#'
#' @param fn A function in character strings or a vector of functions.
#' @param Covvars Variance-covariance matrix of the variables. If it is not
#' specified, they are automatically generated.
#' @param vars A vector of characters of the random variables. If the random
#' variables are not listed in `vars`, they are treated as constants. If `vars`
#' is missing, all names in `RAM` are treated as random variables.
#' @param Var.name Name of the variances.
#' @param Cov.name Name of the covariances.
#' @param simplify Attempt to simplify the output.
#' @return Variance-covariance matrix of the functions.
#' @author Mike W.-L. Cheung <mikewlcheung@@nus.edu.sg>
#' @export
#' @examples
#' \dontrun{
#' 
#' #### Fisher-z-transformation
#' fn  <- "0.5*log((1+r)/(1-r))"
#'
#' ## Sampling variance of r
#' Covvars <- "(1-r^2)^2/n"
#'
#' deltamethod(fn=fn, Covvars=Covvars, vars="r")
#' ## $fn
#' ##     [,1]
#' ## fn1 "0.5*log((r+1)/(1-r))"
#'
#' ## $Covfn
#' ##     fn1
#' ## fn1 "1/n"
#'
#' ## $vars
#' ## [1] "r"
#'
#' ## $Covvars
#' ##   r
#' ## r "(1-r^2)^2/n"
#'
#' ## $Jmatrix
#' ##     r
#' ## fn1 "(0.5*(1-r+r+1)*(1-r))/((1-r)^2*(r+1))"
#'
#' #### Raw mean difference: y_treatment - y_control
#' fn <- "yt - yc"
#'
#' ## Sampling covariance matrix
#' ## S2p: pooled variance
#' ## nt: n_treatment
#' ## nc: n_control
#' Covvars <- matrix(c("S2p/nt", 0,
#'                     0, "S2p/nc"),
#'                   ncol=2, nrow=2)
#'
#' deltamethod(fn=fn, Covvars=Covvars, vars=c("yt", "yc"))
#' ## $fn
#' ##     [,1]
#' ## fn1 "yt-yc"
#'
#' ## $Covfn
#' ##     fn1
#' ## fn1 "(S2p*nt+S2p*nc)/(nt*nc)"
#'
#' ## $vars
#' ## [1] "yt" "yc"
#'
#' ## $Covvars
#' ##    yt       yc
#' ## yt "S2p/nt" "0"
#' ## yc "0"      "S2p/nc"
#'
#' ## $Jmatrix
#' ##     yt  yc
#' ## fn1 "1" "-1"
#'
#' #### log(odds)
#' fn <- "log(p/(1-p))"
#'
#' ## Sampling variance of p
#' Covvars <- "p*(1-p)/n"
#'
#' ## Though it is correct, the simplification does not work well.
#' deltamethod(fn=fn, Covvars=Covvars, vars="p")
#' ## $fn
#' ##     [,1]
#' ## fn1 "log(p/(1-p))"
#'
#' ## $Covfn
#' ##     fn1
#' ## fn1 "(3*p^2-p^3-3*p+1)/((p^4-4*p^3+6*p^2-4*p+1)*p*n)"
#'
#' ## $vars
#' ## [1] "p"
#'
#' ## $Covvars
#' ##   p
#' ## p "(p*(1-p))/n"
#'
#' ## $Jmatrix
#' ##     p
#' ## fn1 "((1-p+p)*(1-p))/((1-p)^2*p)"
#' }
deltamethod <- function(fn, Covvars, vars, Var.name="V", Cov.name="C",
                        simplify=TRUE) {
    
    ## Univariate or multivariate
    fn.p <- length(fn)

    ## function names
    fn.names <- paste0("fn", seq_len(fn.p))

    ## convert it to a column vector
    fn <- matrix(fn, ncol=1)
    rownames(fn) <- fn.names

    fn.S <- as_sym(fn)

    ## Get the variable names
    varlist <- sort(unique(all.names(parse(text=fn), functions=FALSE)))
    if (missing(vars)) {
        vars <- varlist
    } else {
        if (any(!vars %in% varlist)) {
            stop("Some of \"vars\" do not agree with the names in \"fn\".\n")
        }
    }

    ##
    if (missing(Covvars)) {
        ## Variance covariance matrix of x
        Covvars <- outer(vars, vars, function(y, z) paste0(Cov.name, y, z))
        metaSEM::Diag(Covvars) <- paste0(Var.name, vars)
        ## Make it symmetric
        Covvars <- metaSEM::vec2symMat(OpenMx::vech(Covvars))
        Covvars.S <- caracas::as_sym(Covvars)
    } else {
        Covvars.S <- caracas::as_sym(Covvars)
    }

    ## Jacobian matrix
    Jmatrix <- caracas::jacobian(fn.S, vars)

    Covfn <- Jmatrix %*% Covvars.S %*% t(Jmatrix)

    if (simplify) {Covfn <- caracas::simplify(Covfn)}

    Covfn <- caracas::as_character_matrix(Covfn)
    dimnames(Covfn) <- list(fn.names, fn.names)
    Covvars <- as.matrix(Covvars)
    dimnames(Covvars) <- list(vars, vars)

    Jmatrix <- caracas::as_character_matrix(Jmatrix)
    dimnames(Jmatrix) <- list(fn.names, vars)

    list(fn=fn, Covfn=Covfn, vars=vars, Covvars=Covvars, Jmatrix=Jmatrix)
}