From 75f2135c9f2f140fada4198bf78a9fd09637486e Mon Sep 17 00:00:00 2001 From: Mike Cheung Date: Wed, 10 Jul 2024 09:33:50 +0800 Subject: [PATCH] Fix Mmatrix defaults in lavaan2RAM --- DESCRIPTION | 4 +- NEWS | 9 ++ R/create.mxModel.R | 27 +++- R/lavaan2RAM.R | 280 +++++++++++++++++++++-------------------- man/metaSEM-package.Rd | 4 +- 5 files changed, 180 insertions(+), 144 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 013e07c..7556632 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,8 +1,8 @@ Package: metaSEM Type: Package Title: Meta-Analysis using Structural Equation Modeling -Version: 1.4.0 -Date: 2024-05-15 +Version: 1.4.2 +Date: 2024-07-10 Depends: R (>= 3.4.0), OpenMx Imports: Matrix, MASS, ellipse, graphics, stats, utils, mvtnorm, numDeriv, lavaan Suggests: metafor, semPlot, R.rsp, testthat, matrixcalc diff --git a/NEWS b/NEWS index 55ee577..37efe4c 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,12 @@ +Release 1.4.2 (Jul 10, 2024) +==================================== +* Fix the means of the latent variables at 0 in lavaan2RAM(). + +Release 1.4.1 (Jun 17, 2024) +==================================== +* Add "<~" operator in lavaan2RAM(). +* Fix startvalues not assigned in create.mxModel(). + Release 1.4.0 (May 15, 2024) ==================================== * Release to CRAN. diff --git a/R/create.mxModel.R b/R/create.mxModel.R index 9f90bd2..18b09a0 100644 --- a/R/create.mxModel.R +++ b/R/create.mxModel.R @@ -21,6 +21,7 @@ create.mxModel <- function(model.name="mxModel", RAM=NULL, data=NULL, ## Extract all observed and latent variable names var.names <- colnames(Fmatrix$values) + ## Without raw data if (is.null(data)) { ## Without means @@ -43,7 +44,8 @@ create.mxModel <- function(model.name="mxModel", RAM=NULL, data=NULL, ## Create an incomplete model, which will be used to store other mx objects. mx.model <- mxModel(model.name, mx.data, expFun, mxFitFunctionML()) - + + ## Note. startvalues may overwrite the starting values in RAM ## Collate the starting values from RAM and add them to startvalues para.labels <- c(Amatrix$labels[Amatrix$free], Smatrix$labels[Smatrix$free], Mmatrix$labels[Mmatrix$free]) @@ -52,9 +54,22 @@ create.mxModel <- function(model.name="mxModel", RAM=NULL, data=NULL, ## Name the starting values with names, which is consistent with the startvalues names(para.values) <- para.labels para.values <- as.list(para.values) - ## Remove starting values from para.values if they are overlapped with startvalues - para.values[names(para.values) %in% names(startvalues)] <- NULL - startvalues <- c(startvalues, para.values) + + ## Prepare startvalues + if (is.null(startvalues)) {## Note. startvalues may overwrite the starting values in RAM + startvalues <- para.values + } else { + ## Remove starting values from para.values if they are overlapped with startvalues + para.values[names(para.values) %in% names(startvalues)] <- NULL + startvalues <- c(startvalues, para.values) + + ## Replace the startvalues in Amatrix, Smatrix, and Mmatrix + for (i in seq_along(startvalues)) { + Amatrix$values[Amatrix$labels==names(startvalues)[i]] <- startvalues[[i]] + Smatrix$values[Smatrix$labels==names(startvalues)[i]] <- startvalues[[i]] + Mmatrix$values[Mmatrix$labels==names(startvalues)[i]] <- startvalues[[i]] + } + } ## Extract a local copy for ease of reference ## Remove starting values for ease of matching @@ -174,7 +189,6 @@ create.mxModel <- function(model.name="mxModel", RAM=NULL, data=NULL, } if (run) { - ## Default is z mx.fit <- tryCatch(mxRun(mx.model, intervals=(intervals.type=="LB"), suppressWarnings=TRUE, silent=TRUE, ...), @@ -362,6 +376,9 @@ plot.mxRAMmodel <- function(x, manNames=NULL, latNames=NULL, S <- x$mx.fit@matrices$Smatrix$values F <- x$mx.fit@matrices$Fmatrix$values M <- x$mx.fit@matrices$Mmatrix$values + ## Fixed when M is NULL, i.e., no mean structure + if (is.null(M)) M <- matrix(0, nrow=1, ncol=ncol(A)) + RAM <- x$RAM ## When there are definition variables, data in the first role are used in diff --git a/R/lavaan2RAM.R b/R/lavaan2RAM.R index 19b4ad5..8b77f1c 100644 --- a/R/lavaan2RAM.R +++ b/R/lavaan2RAM.R @@ -2,167 +2,177 @@ lavaan2RAM <- function(model, obs.variables = NULL, A.notation="ON", S.notation= M.notation="mean", A.start=0.1, S.start=0.5, M.start=0, auto.var = TRUE, std.lv = TRUE, ngroups=1, ...) { - ## if (!requireNamespace("lavaan", quietly=TRUE)) - ## stop("\"lavaan\" package is required for this function.") + ## if (!requireNamespace("lavaan", quietly=TRUE)) + ## stop("\"lavaan\" package is required for this function.") - ## Default: fix the latent independent variables at 1 - my.model <- lavaan::lavaanify(model, fixed.x = FALSE, auto.var=auto.var, - std.lv=std.lv, ngroups=ngroups, ...) + ## Default: fix the latent independent variables at 1 + my.model <- lavaan::lavaanify(model, fixed.x = FALSE, auto.var=auto.var, + std.lv=std.lv, ngroups=ngroups, ...) - ## Maximum no. of groups - max.gp <- max(my.model$group) + ## Maximum no. of groups + max.gp <- max(my.model$group) - ## Empty list to store the matrices per group - out <- list() + ## Empty list to store the matrices per group + out <- list() - for (gp in seq_len(max.gp)) { - ## Select the ith group - mod <- my.model[my.model$group==gp, ] + for (gp in seq_len(max.gp)) { + ## Select the ith group + mod <- my.model[my.model$group==gp, ] - ## if (any...) is required to avoid error when there is no element for the assignment - ## set the starting values in A if NA - if (any((mod$op=="=~"|mod$op=="~")&is.na(mod$ustart))) { - mod[(mod$op=="=~"|mod$op=="~")&is.na(mod$ustart), ]$ustart <- A.start - } + ## if (any...) is required to avoid error when there is no element for the assignment + ## set the starting values in A if NA + if (any((mod$op=="=~"|mod$op=="~")&is.na(mod$ustart))) { + mod[(mod$op=="=~"|mod$op=="~")&is.na(mod$ustart), ]$ustart <- A.start + } - ## set the starting values in M if NA - if (any(mod$op=="~1"&is.na(mod$ustart))) { - mod[mod$op=="~1"&is.na(mod$ustart), ]$ustart <- M.start - } + ## set the starting values in M if NA + if (any(mod$op=="~1"&is.na(mod$ustart))) { + mod[mod$op=="~1"&is.na(mod$ustart), ]$ustart <- M.start + } - ## set the starting values in S and free parameters if NA (variances) - if (any(mod$op=="~~"&is.na(mod$ustart)&(mod$lhs==mod$rhs))) { - mod[mod$op=="~~"&is.na(mod$ustart)&(mod$lhs==mod$rhs), ]$ustart <- S.start - } + ## set the starting values in S and free parameters if NA (variances) + if (any(mod$op=="~~"&is.na(mod$ustart)&(mod$lhs==mod$rhs))) { + mod[mod$op=="~~"&is.na(mod$ustart)&(mod$lhs==mod$rhs), ]$ustart <- S.start + } - ## Set the starting values in S and free parameters if NA (covariances) - if (any(mod$op=="~~"&is.na(mod$ustart)&(mod$lhs!=mod$rhs))) { - mod[mod$op=="~~"&is.na(mod$ustart)&(mod$lhs!=mod$rhs), ]$ustart <- 0 - } + ## Set the starting values in S and free parameters if NA (covariances) + if (any(mod$op=="~~"&is.na(mod$ustart)&(mod$lhs!=mod$rhs))) { + mod[mod$op=="~~"&is.na(mod$ustart)&(mod$lhs!=mod$rhs), ]$ustart <- 0 + } - ## all variables - ## Removed sort(); otherwise, the variables will be arranged as x1, x10, x2, x3... - all.var <- unique(c(mod$lhs, mod$rhs)) + ## all variables + ## Removed sort(); otherwise, the variables will be arranged as x1, x10, x2, x3... + all.var <- unique(c(mod$lhs, mod$rhs)) - ## latent variables: (with indicators) - latent <- unique(mod[mod$op== "=~", ]$lhs) - ## observed variables: not latent - observed <- all.var[!(all.var %in% latent)] - ## remove empty string "" when there are mean structure - observed <- observed[observed !=""] + ## latent variables: (with indicators) + latent <- unique(mod[mod$op== "=~", ]$lhs) + ## observed variables: not latent + observed <- all.var[!(all.var %in% latent)] + ## remove empty string "" when there are mean structure + observed <- observed[observed !=""] - ## check whether observed in model = observed in argument - if (!is.null(obs.variables)) { - if (!identical(sort(observed), sort(obs.variables))) { - stop("Names in \"obs.variables\" do not agree with those in model.\n") - } else { - ## arrange the observed variables according to obs.var argument - observed <- obs.variables - } - } + ## check whether observed in model = observed in argument + if (!is.null(obs.variables)) { + if (!identical(sort(observed), sort(obs.variables))) { + stop("Names in \"obs.variables\" do not agree with those in model.\n") + } else { + ## arrange the observed variables according to obs.var argument + observed <- obs.variables + } + } - ## if there are latent variables - if (length(latent)>0) { - ## arrange variable list as observed + latent - all.var <- c(observed, latent) - } else { - all.var <- observed - } + ## if there are latent variables + if (length(latent)>0) { + ## arrange variable list as observed + latent + all.var <- c(observed, latent) + } else { + all.var <- observed + } + + no.lat <- length(latent) + no.obs <- length(observed) + no.all <- no.lat+no.obs - no.lat <- length(latent) - no.obs <- length(observed) - no.all <- no.lat+no.obs - - Amatrix <- matrix(0, ncol=no.all, nrow=no.all, dimnames=list(all.var, all.var)) - Smatrix <- matrix(0, ncol=no.all, nrow=no.all, dimnames=list(all.var, all.var)) - Mmatrix <- matrix(0, nrow=1, ncol=no.all, dimnames=list(1, all.var)) + Amatrix <- matrix(0, ncol=no.all, nrow=no.all, dimnames=list(all.var, all.var)) + Smatrix <- matrix(0, ncol=no.all, nrow=no.all, dimnames=list(all.var, all.var)) + ## Fixed Mmatrix by setting the default to estimate the means; + ## otherwise, the mean structure is misspecified. + ## Defaults: observed variables are free, whereas latent variables are fixed at 0 + Mmatrix <- matrix(c(paste0("0*", observed, "mean"), rep("0", no.lat)), + nrow=1, ncol=no.all, dimnames=list(1, all.var)) - ## Prepare the labels - for (i in seq_len(nrow(mod))) { - ## if there is no label - if (mod[i, ]$label=="") { - switch(mod[i, ]$op, - "=~" = mod[i, ]$label <- paste0(mod[i, ]$rhs, A.notation, mod[i, ]$lhs), - "~" = mod[i, ]$label <- paste0(mod[i, ]$lhs, A.notation, mod[i, ]$rhs), - "~~" = mod[i, ]$label <- paste0(mod[i, ]$lhs, S.notation, mod[i, ]$rhs), - "~1" = mod[i, ]$label <- paste0(mod[i, ]$lhs, M.notation)) + ## Prepare the labels + for (i in seq_len(nrow(mod))) { + ## if there is no label + if (mod[i, ]$label=="") { + switch(mod[i, ]$op, + "=~" = mod[i, ]$label <- paste0(mod[i, ]$rhs, A.notation, mod[i, ]$lhs), + "<~" = mod[i, ]$label <- paste0(mod[i, ]$lhs, A.notation, mod[i, ]$rhs), + "~" = mod[i, ]$label <- paste0(mod[i, ]$lhs, A.notation, mod[i, ]$rhs), + "~~" = mod[i, ]$label <- paste0(mod[i, ]$lhs, S.notation, mod[i, ]$rhs), + "~1" = mod[i, ]$label <- paste0(mod[i, ]$lhs, M.notation)) } - } + } - ## Replace NA to 0 in ustart if there are still NA - ## mod$ustart[is.na(mod$ustart)] <- 0 + ## Replace NA to 0 in ustart if there are still NA + ## mod$ustart[is.na(mod$ustart)] <- 0 - ## keys in as.mxMatrix format - key <- with(mod, ifelse(free==0, yes=ustart, no=paste(ustart, label, sep="*"))) + ## keys in as.mxMatrix format + key <- with(mod, ifelse(free==0, yes=ustart, no=paste(ustart, label, sep="*"))) - for (i in seq_len(nrow(mod))) { - my.line <- mod[i, ] - switch(my.line$op, - ## lhs: IV; rhs: DV - "=~" = Amatrix[my.line$rhs, my.line$lhs] <- key[i], - ## lhs: DV; rhs: IV - "~" = Amatrix[my.line$lhs, my.line$rhs] <- key[i], - "~~" = Smatrix[my.line$lhs, my.line$rhs] <- - Smatrix[my.line$rhs, my.line$lhs] <- key[i], - ## means - "~1" = Mmatrix[1, my.line$lhs] <- key[i] - ) ## from switch - } ## from for loop - - Fmatrix <- create.Fmatrix(c(rep(1, no.obs), rep(0, no.lat)), as.mxMatrix=FALSE) - dimnames(Fmatrix) <- list(observed, all.var) + for (i in seq_len(nrow(mod))) { + my.line <- mod[i, ] + switch(my.line$op, + ## lhs: IV; rhs: DV + "=~" = Amatrix[my.line$rhs, my.line$lhs] <- key[i], + ## lhs: DV; rhs: IV + "<~" = Amatrix[my.line$lhs, my.line$rhs] <- key[i], + "~" = Amatrix[my.line$lhs, my.line$rhs] <- key[i], + "~~" = Smatrix[my.line$lhs, my.line$rhs] <- + Smatrix[my.line$rhs, my.line$lhs] <- key[i], + ## means + "~1" = Mmatrix[1, my.line$lhs] <- key[i] + ) ## from switch + } ## from for loop - out[[gp]] <- list(A=Amatrix, S=Smatrix, F=Fmatrix, M=Mmatrix) - } + Fmatrix <- create.Fmatrix(c(rep(1, no.obs), rep(0, no.lat)), as.mxMatrix=FALSE) + dimnames(Fmatrix) <- list(observed, all.var) + + out[[gp]] <- list(A=Amatrix, S=Smatrix, F=Fmatrix, M=Mmatrix) + } - ## Add group names, 1, 2, 3... to the list - names(out) <- seq_along(out) + ## Add group names, 1, 2, 3... to the list + names(out) <- seq_along(out) - ## If there are constraints such as .p1.==.p2.; remove them first - ## otherwise, .p1.==.p2. will create an empty list in mxalgebra - if (length(grep("^\\.", my.model$lhs)) >0 ) { - my.model <- my.model[-grep("^\\.", my.model$lhs), ] - } + ## If there are constraints such as .p1.==.p2.; remove them first + ## otherwise, .p1.==.p2. will create an empty list in mxalgebra + if (length(grep("^\\.", my.model$lhs)) >0 ) { + my.model <- my.model[-grep("^\\.", my.model$lhs), ] + } - ## Check if there are constraints or algebras - if (any(my.model$group==0)) { + ## Check if there are constraints or algebras + if (any(my.model$group==0)) { - ## An empty list to store mxAlgebra and mxConstraint - mxalgebra <- list() - con_index <- 1 # A counter for no. of constraints, started from 1 + ## An empty list to store mxAlgebra and mxConstraint + mxalgebra <- list() + con_index <- 1 # A counter for no. of constraints, started from 1 - ## Constraints and algebras only - y <- my.model[my.model$group==0, , drop=FALSE] - - for (i in seq_len(nrow(y))) { - switch(y[i, 'op'], - ## mxAlgebra - ":=" = { eval(parse(text=paste0(y[i,'lhs'], - "<- mxAlgebra(", y[i,'rhs'], ", name=\"", y[i,'lhs'], "\")"))) - eval(parse(text=paste0("mxalgebra <- c(mxalgebra, ", y[i, 'lhs'], "=", y[i, 'lhs'], ")"))) }, - ## Default condition to test if there are mxConstraints - if (y[i, 'op'] %in% c("==", ">", "<")) { - eval(parse(text=paste0("constraint", con_index, " <- mxConstraint(", y[i, 'lhs'], - y[i, 'op'], y[i, 'rhs'], ", name=\"constraint", con_index, "\")"))) - eval(parse(text=paste0("mxalgebra <- c(mxalgebra, constraint", - con_index, "=constraint", con_index, ")"))) - con_index <- con_index + 1 - } - ) ## from switch - } ## from for loop + ## Constraints and algebras only + y <- my.model[my.model$group==0, , drop=FALSE] + + for (i in seq_len(nrow(y))) { + switch(y[i, 'op'], + ## mxAlgebra + ":=" = { eval(parse(text=paste0(y[i,'lhs'], + "<- mxAlgebra(", y[i,'rhs'], + ", name=\"", y[i,'lhs'], "\")"))) + eval(parse(text=paste0("mxalgebra <- + c(mxalgebra, ", y[i, 'lhs'], + "=", y[i, 'lhs'], ")"))) }, + ## Default condition to test if there are mxConstraints + if (y[i, 'op'] %in% c("==", ">", "<")) { + eval(parse(text=paste0("constraint", con_index, " <- mxConstraint(", y[i, 'lhs'], + y[i, 'op'], y[i, 'rhs'], ", name=\"constraint", con_index, "\")"))) + eval(parse(text=paste0("mxalgebra <- c(mxalgebra, constraint", + con_index, "=constraint", con_index, ")"))) + con_index <- con_index + 1 + } + ) ## from switch + } ## from for loop - ## Append mxalgebra to out[[1]] - out[[1]] <- list(A=out[[1]]$A, S=out[[1]]$S, F=out[[1]]$F, M=out[[1]]$M, mxalgebras=mxalgebra) - } ## End if - ## else { - ## out[[1]] <- list(A=out[[1]]$A, S=out[[1]]$S, F=out[[1]]$F, M=out[[1]]$M) - ## } + ## Append mxalgebra to out[[1]] + out[[1]] <- list(A=out[[1]]$A, S=out[[1]]$S, F=out[[1]]$F, M=out[[1]]$M, + mxalgebras=mxalgebra) + } ## End if + ## else { + ## out[[1]] <- list(A=out[[1]]$A, S=out[[1]]$S, F=out[[1]]$F, M=out[[1]]$M) + ## } - ## Output the first list instead of a list of one item when there is only 1 group - if (max.gp==1) { - out <- out[[1]] - } + ## Output the first list instead of a list of one item when there is only 1 group + if (max.gp==1) { + out <- out[[1]] + } - out + out } diff --git a/man/metaSEM-package.Rd b/man/metaSEM-package.Rd index eee3cbc..200a51c 100644 --- a/man/metaSEM-package.Rd +++ b/man/metaSEM-package.Rd @@ -14,8 +14,8 @@ \tabular{ll}{ Package: \tab metaSEM\cr Type: \tab Package\cr -Version: \tab 1.4.0\cr -Date: \tab 2024-05-15\cr +Version: \tab 1.4.2\cr +Date: \tab 2024-07-10\cr License: \tab GPL (>=2)\cr LazyLoad: \tab yes\cr }