New 2.3.0 version- Multigroup. Error in `[.data.frame`(data, , var) : undefined columns selected

[BFadrique]

I am using the exact same code to perform a multigroup analysis before and after the version update. Now I get the "Error in [.data.frame(data, , var) : undefined columns selected" I have the data=dataframe in the model but it does not solve the problem.

model<- psem(lm(y ~ x+q+p, dataframe), lm(z ~ y+u+x+q+p,dataframe), data=dataframe)

summary(model, .progressBar = T) ##works multigroup(model, group = "Group") ##doesn't work. Group is on dataframe.

Before the update the code worked without problems

[ShanKothari]

I just wanted to note that I encountered the same issue in piecewiseSEM 2.3.0. It runs fine in piecewiseSEM 2.1.2 (but alas, then I run into the object 'ret' not found error with random effects).

Sometime when I'm less in a rush, I'll try to put together a good example, but for now perhaps it helps to know that this is not just a single user's issue. Thanks to Dr. Lefcheck and others for all the work you put into this package!

[KatSym]

I am running into the same problem. SEM model runs fine and I get the summary, but when I add a group with multilevel I get "Error in [.data.frame(data, , var) : undefined columns selected"

Thanks!

[nehamohanbabu]

I am running into a similar issue while using mixed models within psem(). Just wondering if anyone has found a solution to this error message?

[I have already tried to include "data=data_frame" within psem() and have tried different versions of piecewiseSEM] Thanks in advance!

[liaojiahui-r]

I ran into this similar problem! this is how I solved it:

In short, the issues is due to the multigroup function remove the data from modelList (i.e., modelList <- piecewiseSEM:::removeData(modelList, formulas = 1)); the basisSet function, however, require modelList with data (i.e., modelList$data).

Example:

library(piecewiseSEM) library(lme4) library(glmmTMB) library(nlme)

# Generate data set.seed(1) dat <- data.frame( x1 = rnorm(100), y1 = rnorm(100), y2 = rnorm(100), category = LETTERS[1:2], random = rep(letters[1:4], each = 25) )

# nlme: nlme.model <- psem( lme(y1 ~ x1, random = ~ 1 | random, dat), lme(y2 ~ y1, random = ~ 1 | random, dat), data = dat ) # NO! multigroup(nlme.model, group = "category") # with the modified multigroup2 function... then, YES! multigroup2(nlme.model, group = "category")

# glmmTMB: glmmTMB.model <- psem( glmmTMB(y1 ~ x1 + (1 | random), dat, REML = FALSE), glmmTMB(y2 ~ y1 + (1 | random), dat, REML = FALSE), data = dat ) # A simple trick: glmmTMB.model$data$category <- dat$category # ... then, NO! multigroup(glmmTMB.model, group = "category") # with the modified multigroup2 function... then, YES! multigroup2(glmmTMB.model, group = "category") # similar as nlme!

Below, are two modified functions: baseSet2 and multigroup2, respectively:

basisSet2 <- function (modelList.with.data, direction = NULL, interactions = FALSE) { amat <- getDAG(modelList.with.data) b <- lapply(1:nrow(amat), function(i) { lapply(i:ncol(amat), function(j) { if (amat[i, j] != 0 | i == j) NULL else { cvars <- c(rownames(amat)[amat[, rownames(amat)[i], drop = FALSE] == 1], rownames(amat)[amat[, rownames(amat)[j], drop = FALSE] == 1]) cvars <- cvars[!duplicated(cvars)] c(rownames(amat)[i], rownames(amat)[j], cvars) } }) }) b <- unlist(b, recursive = FALSE) b <- b[!sapply(b, is.null)] if (length(b) > 0) { b <- piecewiseSEM:::filterExogenous(b, modelList.with.data, amat) b <- piecewiseSEM:::filterSmoothed(b, modelList.with.data) b <- piecewiseSEM:::filterExisting(b, modelList.with.data) b <- piecewiseSEM:::filterInteractions(b, interactions) b <- piecewiseSEM:::removeCerror(b, modelList.with.data) b <- piecewiseSEM:::reverseAddVars(b, modelList.with.data, amat) b <- piecewiseSEM:::reverseNonLin(b, modelList.with.data, amat) b <- piecewiseSEM:::fixCatDir(b, modelList.with.data) # Here is the cause of the problem! if (!is.null(direction)) b <- piecewiseSEM:::specifyDir(b, direction) } class(b) <- "basisSet" return(b) }

and,

multigroup2 <- function( modelList, group, standardize = "scale", standardize.type = "latent.linear", test.type = "III") { name <- deparse(match.call()$modelList) data <- modelList$data

modelList.with.data <- modelList # My solution modelList <- piecewiseSEM:::removeData(modelList, formulas = 1) # piecewiseSEM:::fixCatDir(b, modelList) caused the issue

intModelList <- lapply(modelList, function(i) { rhs2 <- paste(paste(piecewiseSEM:::all.vars_trans(i)[-1], "", group), collapse = " + " ) i <- update(i, formula(paste(". ~ ", rhs2))) return(i) }) newModelList <- lapply(unique(data[, group]), function(i) { update(as.psem(modelList), data = data[data[, group] == i, ] ) }) names(newModelList) <- unique(data[, group]) coefsList <- lapply( newModelList, coefs, standardize, standardize.type, test.type ) names(coefsList) <- unique(data[, group]) coefTable <- coefs( modelList, standardize, standardize.type, test.type ) anovaTable <- anova(as.psem(intModelList))[[1]] anovaInts <- anovaTable[grepl(":", anovaTable$Predictor), ] global <- anovaInts[anovaInts$P.Value >= 0.05, c( "Response", "Predictor" )] global$Predictor <- sub(":", "\1", sub(group, "\1", global$Predictor)) if (nrow(global) == nrow(anovaInts)) { newCoefsList <- list(global = coefTable) } else { newCoefsList <- lapply(names(coefsList), function(i) { ct <- as.matrix(coefsList[[i]]) idx <- which( apply(ct[, 1:2], 1, paste, collapse = "") %in% apply(global[, 1:2], 1, paste, collapse = "") ) ct[idx, ] <- as.matrix(coefTable[idx, ]) ct <- cbind(ct, ifelse(1:nrow(ct) %in% idx, "c", "" )) for (j in 1:nrow(ct)) { if (ct[j, ncol(ct)] == "c") { model <- modelList[[which(sapply( piecewiseSEM:::listFormula(modelList), function(x) { piecewiseSEM:::all.vars.merMod(x)[1] == ct[ j, "Response" ] } ))]] data. <- data[data[, group] == i, ] sd.x <- piecewiseSEM:::GetSDx(model, modelList, data., standardize) sd.x <- sd.x[which(names(sd.x) == ct[j, "Predictor"])] sd.y <- piecewiseSEM:::GetSDy( model, data., standardize, standardize.type ) new.coef <- as.numeric(ct[j, "Estimate"]) (sd.x / sd.y) ct[j, "Std.Estimate"] <- ifelse(length(new.coef) > 0, round(as.numeric(new.coef), 4), "-") } } ct <- as.data.frame(ct) ct[is.na(ct)] <- "-" names(ct)[(ncol(ct) - 1):ncol(ct)] <- "" return(ct) }) names(newCoefsList) <- names(coefsList) }

if (nrow(global) == nrow(anovaInts)) { gof <- piecewiseSEM::fisherC(modelList) } else {

b <- piecewiseSEM:::basisSet(modelList) # Here is the cause of the problem!

b <- basisSet2(modelList.with.data)
cf <- coefTable[coefTable$Response %in% global$Response & coefTable$Predictor %in% global$Predictor, ]
b <- lapply(
  X = b,
  FUN = function(i) {
    for (j in 3:length(i)) {
      value <- cf[cf$Response == i[2] & cf$Predictor == i[j], "Estimate"]
      if (length(value) != 0) {
        i[j] <- paste0("offset(", value, "*", i[j], ")")
      }
    }
    return(i)
  }
)
if (length(b) == 0) {
  b <- NULL
}
gof <- fisherC(modelList, basis.set = b)

}

ret <- list( name = name, group = group, Cstat = gof, global = global, anovaInts = anovaInts, group.coefs = newCoefsList ) class(ret) <- "multigroup.psem" return(ret) }

Hope it helps!

[misternickel]

@liaojiahui-r , can you explain, please, where I should introduce code, suggested by you, to make multigroup option viable?

[HectorHer-Alo]

Thanks a lot for your new functions, it worked for me using nlme::lme

I ran into this similar problem! this is how I solved it:

In short, the issues is due to the multigroup function remove the data from modelList (i.e., modelList <- piecewiseSEM:::removeData(modelList, formulas = 1)); the basisSet function, however, require modelList with data (i.e., modelList$data).

Example:

library(piecewiseSEM) library(lme4) library(glmmTMB) library(nlme)

# Generate data set.seed(1) dat <- data.frame( x1 = rnorm(100), y1 = rnorm(100), y2 = rnorm(100), category = LETTERS[1:2], random = rep(letters[1:4], each = 25) )

# nlme: nlme.model <- psem( lme(y1 ~ x1, random = ~ 1 | random, dat), lme(y2 ~ y1, random = ~ 1 | random, dat), data = dat ) # NO! multigroup(nlme.model, group = "category") # with the modified multigroup2 function... then, YES! multigroup2(nlme.model, group = "category")

# glmmTMB: glmmTMB.model <- psem( glmmTMB(y1 ~ x1 + (1 | random), dat, REML = FALSE), glmmTMB(y2 ~ y1 + (1 | random), dat, REML = FALSE), data = dat ) # A simple trick: glmmTMB.model$data$category <- dat$category # ... then, NO! multigroup(glmmTMB.model, group = "category") # with the modified multigroup2 function... then, YES! multigroup2(glmmTMB.model, group = "category") # similar as nlme!

Below, are two modified functions: baseSet2 and multigroup2, respectively:

basisSet2 <- function (modelList.with.data, direction = NULL, interactions = FALSE) { amat <- getDAG(modelList.with.data) b <- lapply(1:nrow(amat), function(i) { lapply(i:ncol(amat), function(j) { if (amat[i, j] != 0 | i == j) NULL else { cvars <- c(rownames(amat)[amat[, rownames(amat)[i], drop = FALSE] == 1], rownames(amat)[amat[, rownames(amat)[j], drop = FALSE] == 1]) cvars <- cvars[!duplicated(cvars)] c(rownames(amat)[i], rownames(amat)[j], cvars) } }) }) b <- unlist(b, recursive = FALSE) b <- b[!sapply(b, is.null)] if (length(b) > 0) { b <- piecewiseSEM:::filterExogenous(b, modelList.with.data, amat) b <- piecewiseSEM:::filterSmoothed(b, modelList.with.data) b <- piecewiseSEM:::filterExisting(b, modelList.with.data) b <- piecewiseSEM:::filterInteractions(b, interactions) b <- piecewiseSEM:::removeCerror(b, modelList.with.data) b <- piecewiseSEM:::reverseAddVars(b, modelList.with.data, amat) b <- piecewiseSEM:::reverseNonLin(b, modelList.with.data, amat) b <- piecewiseSEM:::fixCatDir(b, modelList.with.data) # Here is the cause of the problem! if (!is.null(direction)) b <- piecewiseSEM:::specifyDir(b, direction) } class(b) <- "basisSet" return(b) }

and,

multigroup2 <- function( modelList, group, standardize = "scale", standardize.type = "latent.linear", test.type = "III") { name <- deparse(match.call()$modelList) data <- modelList$data

modelList.with.data <- modelList # My solution modelList <- piecewiseSEM:::removeData(modelList, formulas = 1) # piecewiseSEM:::fixCatDir(b, modelList) caused the issue

intModelList <- lapply(modelList, function(i) { rhs2 <- paste(paste(piecewiseSEM:::all.varstrans(i)[-1], "*", group), collapse = " + " ) i <- update(i, formula(paste(". ~ ", rhs2))) return(i) }) newModelList <- lapply(unique(data[, group]), function(i) { update(as.psem(modelList), data = data[data[, group] == i, ] ) }) names(newModelList) <- unique(data[, group]) coefsList <- lapply( newModelList, coefs, standardize, standardize.type, test.type ) names(coefsList) <- unique(data[, group]) coefTable <- coefs( modelList, standardize, standardize.type, test.type ) anovaTable <- anova(as.psem(intModelList))[[1]] anovaInts <- anovaTable[grepl(":", anovaTable$Predictor), ] global <- anovaInts[anovaInts$P.Value >= 0.05, c( "Response", "Predictor" )] global$Predictor <- sub(":", "\1", sub(group, "\1", global$Predictor)) if (nrow(global) == nrow(anovaInts)) { newCoefsList <- list(global = coefTable) } else { newCoefsList <- lapply(names(coefsList), function(i) { ct <- as.matrix(coefsList[[i]]) idx <- which( apply(ct[, 1:2], 1, paste, collapse = "") %in% apply(global[, 1:2], 1, paste, collapse = "_") ) ct[idx, ] <- as.matrix(coefTable[idx, ]) ct <- cbind(ct, ifelse(1:nrow(ct) %in% idx, "c", "" )) for (j in 1:nrow(ct)) { if (ct[j, ncol(ct)] == "c") { model <- modelList[[which(sapply( piecewiseSEM:::listFormula(modelList), function(x) { piecewiseSEM:::all.vars.merMod(x)[1] == ct[ j, "Response" ] } ))]] data. <- data[data[, group] == i, ] sd.x <- piecewiseSEM:::GetSDx(model, modelList, data., standardize) sd.x <- sd.x[which(names(sd.x) == ct[j, "Predictor"])] sd.y <- piecewiseSEM:::GetSDy( model, data., standardize, standardize.type ) new.coef <- as.numeric(ct[j, "Estimate"]) * (sd.x / sd.y) ct[j, "Std.Estimate"] <- ifelse(length(new.coef) > 0, round(as.numeric(new.coef), 4), "-") } } ct <- as.data.frame(ct) ct[is.na(ct)] <- "-" names(ct)[(ncol(ct) - 1):ncol(ct)] <- "" return(ct) }) names(newCoefsList) <- names(coefsList) }

if (nrow(global) == nrow(anovaInts)) { gof <- piecewiseSEM::fisherC(modelList) } else { # b <- piecewiseSEM:::basisSet(modelList) # Here is the cause of the problem! b <- basisSet2(modelList.with.data) cf <- coefTable[coefTable$Response %in% global$Response & coefTable$Predictor %in% global$Predictor, ] b <- lapply( X = b, FUN = function(i) { for (j in 3:length(i)) { value <- cf[cf$Response == i[2] & cf$Predictor == i[j], "Estimate"] if (length(value) != 0) { i[j] <- paste0("offset(", value, "*", i[j], ")") } } return(i) } ) if (length(b) == 0) { b <- NULL } gof <- fisherC(modelList, basis.set = b) }

ret <- list( name = name, group = group, Cstat = gof, global = global, anovaInts = anovaInts, group.coefs = newCoefsList ) class(ret) <- "multigroup.psem" return(ret) }

Hope it helps!

[daniela-pazpez]

I have tried the code using lm and it still does not work- is there a way of using it with lm?