Change backend from Ryacas to caracas and use roxygen2

DESCRIPTION

@@ -1,20 +1,23 @@
 Package: symSEM
 Type: Package
 Title: Symbolic Computation for Structural Equation Models  
-Version: 0.1.1
-Date: 2020-09-26
-Depends: R (>= 3.3.0)
-Imports: OpenMx, metaSEM, Ryacas, mvtnorm
-Suggests: testthat
+Version: 0.2
+Date: 2023-05-27
+Depends: caracas
+Imports: OpenMx, metaSEM
+Suggests: testthat, roxygen2
 Authors@R: person(given = "Mike", family = "Cheung", role = c("aut","cre"), email = "mikewlcheung@nus.edu.sg", comment = c(ORCID = "0000-0003-0113-0758")) 
 Maintainer: Mike Cheung <mikewlcheung@nus.edu.sg>
-Description: A collection of functions for symbolic computation using 'Ryacas'
-  package for structural equation models. This package includes
-  functions to calculate model-implied covariance (and correlation)
-  matrix and sampling covariance matrix of functions of variables using
-  the first-order Taylor approximation. Reference: McArdle and McDonald (1984) <doi:10.1111/j.2044-8317.1984.tb00802.x>.
+Description: A collection of functions for symbolic computation using the
+			 'caracas' package for structural equation models and other
+			 statistical analyses. Among its features is the ability
+			 to calculate the model-implied covariance (and correlation)
+			 matrix and the sampling covariance matrix of variable functions
+			 using the delta method.
 License: GPL (>=2)
+Encoding: UTF-8
 LazyLoad: yes
 LazyData: yes
 ByteCompile: yes
 URL: https://github.com/mikewlcheung/symsem
+RoxygenNote: 7.2.3

NAMESPACE

@@ -1,60 +1,10 @@
-exportPattern("^[^\\.]")
-
-importFrom(Ryacas, "%*%")
-importFrom(Ryacas, "%>%")
-importFrom(Ryacas, "diag<-")
-importFrom(Ryacas, Hessian)
-importFrom(Ryacas, Jacobian)
-importFrom(Ryacas, as_r)
-importFrom(Ryacas, as_y)
-importFrom(Ryacas, diag)
-importFrom(Ryacas, get_output_width)
-importFrom(Ryacas, integrate)
-importFrom(Ryacas, lim)
-importFrom(Ryacas, lower.tri)
-importFrom(Ryacas, set_output_width)
-importFrom(Ryacas, simplify)
-importFrom(Ryacas, tex)
-importFrom(Ryacas, upper.tri)
-importFrom(Ryacas, with_value)
-importFrom(Ryacas, y_fn)
-importFrom(Ryacas, y_print)
-importFrom(Ryacas, y_rmvars)
-importFrom(Ryacas, yac)
-importFrom(Ryacas, yac_assign)
-importFrom(Ryacas, yac_cli)
-importFrom(Ryacas, yac_expr)
-importFrom(Ryacas, yac_silent)
-importFrom(Ryacas, yac_str)
-importFrom(Ryacas, yac_symbol)
-importFrom(Ryacas, ysym)
-importFrom(Ryacas, ysym_ls)
-
-importFrom(OpenMx, vech, vechs, mxAlgebra, mxExpectationNormal, mxFitFunctionML, mxModel, mxData, mxRun)
-
-importFrom(metaSEM, as.mxMatrix, as.mxAlgebra)
-
-importFrom(mvtnorm, rmvnorm)
-
-importFrom(stats, rnorm, coef, cov, cor, var, vcov, sd, deriv)
-
-S3method(as.matrix, yac_symbol)
-
-# importFrom(base, as.matrix)
-
-# required by R3.3
-#importFrom("graphics", "abline", "arrows", "layout", "par", "plot",
-#           "points", "polygon")
-
-#importFrom("stats", "as.formula", "coef", "cov2cor", "deriv",
-#           "na.omit", "pchisq", "rnorm", "pnorm", "printCoefmat", "qchisq",
-#           "qnorm", "reshape", "var", "vcov", 
-#           "cov", "cor", "rWishart", "sd", "quantile")
-#importFrom("utils", "read.table")
-
-
-
-
-
-
-
+# Generated by roxygen2: do not edit by hand
+
+export(deltamethod)
+export(impliedS)
+import(caracas)
+importFrom(OpenMx,vech)
+importFrom(OpenMx,vechs)
+importFrom(metaSEM,Diag)
+importFrom(metaSEM,as.symMatrix)
+importFrom(metaSEM,vec2symMat)

NEWS

@@ -1,7 +1,11 @@
+Release 0.2 (May 27, 2023)
+====================================
+* Change the backend from Ryacas to caracas.
+
 Release 0.1.1 (Sep 26, 2020)
 ====================================
-* Fixed impliedS() to work with "_" character in RAM.
+* Fix impliedS() to work with "_" character in RAM.
 
 Release 0.1 (Jul 16, 2020)
 ====================================
-* Released to CRAN.
+* Release to CRAN.

R/deltamethod.R

@@ -1,85 +1,164 @@
 ## fn <- "log(a)^b+c"
 ## fn <- c("x^2 + log(y)*z", "x+y+sin(z)")
-## Note. Initiate symbols first does not work as there are problems in creating 
-## Better: replace R -> yac, create character matrices, ysym convert
 
-## TODO: ensure dimensions of Covvars and vars are matched
-deltamethod <- function(fn, Covvars, vars, Var.name="V", Cov.name="C", simplify=TRUE) {
+#' 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
+#'
+#' #### 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 into a column vector if multivariate
-    if (fn.p>1) {
-        fn <- matrix(fn, ncol=1)
-    }
-
+    ## 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")
-        }        
+        }
     }
 
-    ## Need to convert all variables into symbols first;
-    ## otherwise, some functions, e.g., Ln, will not be converted from R to yac.
-    #for (i in seq_along(varlist)) {
-    #    eval(parse(text=paste0(varlist[i], " <- sym('", varlist[i],"')")))
-    #}    
-    #fn <- eval(parse(text=fn))
-    fn <- sym(fn)
-
     ##
     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 <- sym(Covvars)
+        Covvars.S <- caracas::as_sym(Covvars)
     } else {
-        Covvars <- sym(Covvars)
+        Covvars.S <- caracas::as_sym(Covvars)
     }
- 
+
     ## Jacobian matrix
-    ## Univariate
-    if (fn.p==1) {
-        J <- deriv(fn, vars)
-        ## Possibly bugs in Ryacas?
-        ## Need to clean up as too many {}
-        J <- sym(as.matrix(J))
-        J <- t(J)
-    } else {
-    ## Multivariate
-        J <- Ryacas::Jacobian(fn, vars)
-        J <- sym(as.matrix(J))
-        # J <- Ryacas::simplify(J)
-    }
+    Jmatrix <- caracas::jacobian(fn.S, vars)
 
-    fn <- as.matrix(fn)
-    rownames(fn) <- fn.names
-
-    Covfn <- J * Covvars * t(J)
+    Covfn <- Jmatrix %*% Covvars.S %*% t(Jmatrix)
 
-    if (simplify) {
-        Covfn <- Ryacas::simplify(Covfn)
-        Jmatrix <- Ryacas::simplify(J)
-    }
+    if (simplify) {Covfn <- caracas::simplify(Covfn)}
 
-    Covfn <- as.matrix(Covfn)
+    Covfn <- caracas::as_character_matrix(Covfn)
     dimnames(Covfn) <- list(fn.names, fn.names)
-
     Covvars <- as.matrix(Covvars)
     dimnames(Covvars) <- list(vars, vars)
 
-    Jmatrix <- as.matrix(J)
+    Jmatrix <- caracas::as_character_matrix(Jmatrix)
     dimnames(Jmatrix) <- list(fn.names, vars)
-    
+
     list(fn=fn, Covfn=Covfn, vars=vars, Covvars=Covvars, Jmatrix=Jmatrix)
 }
+