Skip to content

Commit

Permalink
Fix Mmatrix defaults in lavaan2RAM
Browse files Browse the repository at this point in the history
  • Loading branch information
@mikewlcheung
mikewlcheung committed Jul 10, 2024
1 parent 8507c3d commit 75f2135
Show file tree
Hide file tree
Showing 5 changed files with 180 additions and 144 deletions.
4 changes: 2 additions & 2 deletions DESCRIPTION
View file
Original file line number Diff line number Diff line change
@@ -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
Expand Down
9 changes: 9 additions & 0 deletions NEWS
View file
Original file line number Diff line number Diff line change
@@ -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.
Expand Down
27 changes: 22 additions & 5 deletions R/create.mxModel.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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])
Expand All @@ -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
Expand Down Expand Up @@ -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, ...),
Expand Down Expand Up @@ -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
Expand Down
280 changes: 145 additions & 135 deletions R/lavaan2RAM.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
}
Loading

0 comments on commit 75f2135

Please sign in to comment.