Dear R users, I am trying to simulate some multitrait-multimethod models using the packages sem and psych but whatever I do to deal with models which do not converge I always get stuck and get error messages such as these: "Error in summary.sem(M1) : coefficient covariances cannot be computed" "Error in solve.default(res$hessian) : System ist f?r den Rechner singul?r: reziproke Konditionszahl = 8.61916e-17" I am aware that these models could not be computed but it is ok, itis part of what I am trying to show with the simulation, that under specifications x the models converge more easily than under specifications y. When models do not converge I just let R write a string with -1s andhave expected the simulation to go on. But it doesn't happen! Instead of that the computations using mapply to fill matrix rows with the results of the single simulations break up and the whole simulation stops. How could I bring R just to spring the undesired solutions and go on up to the end? Best, Gui # Simulation MTMM myMaxN <- 1 # myRep <- 1 # number of replications in each experimental cell traitLoads <- c(0.3, 0.5, 0.7) # loads of observed variables on trait factors traitCorrs <- c(0.0, 0.4, 0.7) # correlations between traits methodLoads <- c(0.2, 0.3, 0.4) # loads of observed variables on method factors methdCorrs <- c(0.0, 0.2, 0.4) # correlations between methods SampleSize <- 500 # Sample size myMaxIter <- 500 # Maximal number of interactions in every model estimation nCond <- length(traitLoads)* length(traitCorrs)* length(methodLoads)* length(methdCorrs) myRes <- as.numeric(gl(nCond, 1, myRep*nCond)) myloadTrait <- as.numeric(gl(length(traitLoads), 1, length(myRes))) mycorrTrait <- as.numeric(gl(length(traitCorrs), length(traitLoads), length(myRes))) myloadMethd <- as.numeric(gl(length(methodLoads), length(traitLoads) * length(traitCorrs), length(myRes))) mycorrMethd <- as.numeric(gl(length(methdCorrs), length(traitLoads) * length(traitCorrs) * length(methodLoads), length(myRes))) theTotalReplications <- myRes ##### ######## BIG FUNCTION ####### ##### sizeControlGroup <- function(theTotalReplications) { # Big function for running the whole simulation traitLoad <- traitLoads[myloadTrait[theTotalReplications]] traitCorr <- traitCorrs[mycorrTrait[theTotalReplications]] methodLoad <- methodLoads[myloadMethd[theTotalReplications]] methdCorr <- methdCorrs[mycorrMethd[theTotalReplications]] fx = matrix(c( rep(traitLoad, 4), rep(0, 16), rep(traitLoad, 4), rep(0, 16), rep(traitLoad, 4), rep(0, 16), rep(traitLoad, 4), rep(c(methodLoad, 0, 0, 0), 4), rep(c(0, methodLoad, 0, 0), 4), rep(c(0, 0, methodLoad, 0), 4), rep(c(0, 0, 0, methodLoad), 4)), ncol = 8) Phi = matrix(c(1, traitCorr, traitCorr, traitCorr, rep(0,4), traitCorr, 1, traitCorr, traitCorr, rep(0,4), traitCorr, traitCorr, 1, traitCorr, rep(0,4), traitCorr, traitCorr, traitCorr, 1, rep(0,4), rep(0,4),1, methdCorr, methdCorr, methdCorr, methdCorr, rep(0,4),1, methdCorr, methdCorr, methdCorr, methdCorr, rep(0,4),1, methdCorr, methdCorr, methdCorr, methdCorr, rep(0,4),1), ncol = 8) mmtm <- sim.structure(fx, Phi, n = SampleSize, raw=T) correMatrix <- mmtm$r colnames(correMatrix) <- c(paste("item", seq(1:16), sep = "")) rownames(correMatrix) <- c(paste("item", seq(1:16), sep = "")) corForModel = correMatrix # establishes the correlation matrix corForModel for model estimation M1 <- try(sem(CTM1, corForModel, SampleSize, maxiter = myMaxIter), silent = FALSE) # SEM model estimation) # tries to estimate the CT(M-1) model M2 <- try(sem(CTCM, corForModel, SampleSize, maxiter = myMaxIter), silent = FALSE) # SEM model estimation) # tries to estimate the CTCM model if(M1$convergence > 2){ convergenceM1 <- c(0) resultsM1 <- as.numeric(c(rep(-1, 12))) } else { # else needs to be in the same line as the last command myModlChiM1 <- try(summary(M1)) convergenceM1 <- as.numeric(M1$convergence) chiM1 <- as.numeric(myModlChiM1$chisq) dfM1 <- as.numeric(myModlChiM1$df) chiM0 <- as.numeric(myModlChiM1$chisqNull) dfM0 <- as.numeric(myModlChiM1$dfNull) GFIM1 <- as.numeric(myModlChiM1$GFI) AGFIM1 <- as.numeric(myModlChiM1$AGFI) RMSEAM1 <- as.numeric(myModlChiM1$RMSEA[1]) CFIM1 <- as.numeric(myModlChiM1$CFI) BICM1 <- as.numeric(myModlChiM1$BIC) SRMRM1 <- as.numeric(myModlChiM1$SRMR) iterM1 <- as.numeric(myModlChiM1$iterations) resultsM1 <- as.numeric(c(convergenceM1, chiM1, dfM1, chiM0, dfM0, GFIM1, AGFIM1, RMSEAM1, CFIM1, BICM1, SRMRM1, iterM1)) } if(M2$convergence > 2){ convergenceM2 <- c(0) resultsM2 <- as.numeric(c(rep(-1, 12))) } else { # else needs to be in the same line as the last command myModlChiM2 <- try(summary(M2)) convergenceM2 <- as.numeric(M2$convergence) chiM2 <- as.numeric(myModlChiM2$chisq) dfM2 <- as.numeric(myModlChiM2$df) chiM0 <- as.numeric(myModlChiM2$chisqNull) dfM0 <- as.numeric(myModlChiM2$dfNull) GFIM2 <- as.numeric(myModlChiM2$GFI) AGFIM2 <- as.numeric(myModlChiM2$AGFI) RMSEAM2 <- as.numeric(myModlChiM2$RMSEA[1]) CFIM2 <- as.numeric(myModlChiM2$CFI) BICM2 <- as.numeric(myModlChiM2$BIC) SRMRM2 <- as.numeric(myModlChiM2$SRMR) iterM2 <- as.numeric(myModlChiM2$iterations) resultsM2 <- as.numeric(c(convergenceM2, chiM2, dfM2, chiM0, dfM0, GFIM2, AGFIM2, RMSEAM2, CFIM2, BICM2, SRMRM2, iterM2)) } designparameters <- c(traitLoad, traitCorr, methodLoad, methdCorr) myResults <- c(designparameters, SampleSize, resultsM1, resultsM2) #, convergenceM3, resultsM3) # return(myResults) } # End of function sizeControlGroup ############ Loop for replications ########################################################## totalRepeats = 100 for(myRepeats in 1:totalRepeats){ myTests <- mapply(sizeControlGroup, myRes, SIMPLIFY = F) # effectSize myTests <- matrix(unlist(myTests), nc=length(myTests[[1]]), byrow=T) colnames(myTests) <- c("trL", "trCorr", "mthL", "methCorr", "n", "convM1", "ChiM1", "dfM1", "Chim0", "dfm0", "GFIM1", "AGFIM1", "RMSEAM1", "CFIM1", "BICM1", "SRMRM1", "iterM1","ConvM2", "ChiM2", "dfM2", "Chim0", "dfm0", "GFIM2", "AGFIM2", "RMSEAM2", "CFIM2", "BICM2", "SRMRM2", "iterM2") write.table(myTests, paste("C:\\Dokumente und Einstellungen\\simulations\\rep", myRepeats, sep=""), sep = " ", row.names F) } -- View this message in context: r.789695.n4.nabble.com/sem-psych-tp2321572p2321572.html Sent from the R help mailing list archive at Nabble.com.
Dear R users, I am trying to simulate some multitrait-multimethod models using the packages sem and psych but whatever I do to deal with models which do not converge I always get stuck and get error messages such as these: "Error in summary.sem(M1) : coefficient covariances cannot be computed" "Error in solve.default(res$hessian) : System ist für den Rechner singulär: reziproke Konditionszahl = 8.61916e-17" I am aware that these models could not be computed but it is ok, itis part of what I am trying to show with the simulation, that under specifications x the models converge more easily than under specifications y. When models do not converge I just let R write a string with -1s andhave expected the simulation to go on. But it doesn't happen! Instead of that the computations using mapply to fill matrix rows with the results of the single simulations break up and the whole simulation stops. How could I bring R just to spring the undesired solutions and go on up to the end? Best, Gui # Simulation MTMM myMaxN <- 1 # myRep <- 1 # number of replications in each experimental cell traitLoads <- c(0.3, 0.5, 0.7) # loads of observed variables on trait factors traitCorrs <- c(0.0, 0.4, 0.7) # correlations between traits methodLoads <- c(0.2, 0.3, 0.4) # loads of observed variables on method factors methdCorrs <- c(0.0, 0.2, 0.4) # correlations between methods SampleSize <- 500 # Sample size myMaxIter <- 500 # Maximal number of interactions in every model estimation nCond <- length(traitLoads)* length(traitCorrs)* length(methodLoads)* length(methdCorrs) myRes <- as.numeric(gl(nCond, 1, myRep*nCond)) myloadTrait <- as.numeric(gl(length(traitLoads), 1, length(myRes))) mycorrTrait <- as.numeric(gl(length(traitCorrs), length(traitLoads), length(myRes))) myloadMethd <- as.numeric(gl(length(methodLoads), length(traitLoads) * length(traitCorrs), length(myRes))) mycorrMethd <- as.numeric(gl(length(methdCorrs), length(traitLoads) * length(traitCorrs) * length(methodLoads), length(myRes))) theTotalReplications <- myRes ##### ######## BIG FUNCTION ####### ##### sizeControlGroup <- function(theTotalReplications) { # Big function for running the whole simulation traitLoad <- traitLoads[myloadTrait[theTotalReplications]] traitCorr <- traitCorrs[mycorrTrait[theTotalReplications]] methodLoad <- methodLoads[myloadMethd[theTotalReplications]] methdCorr <- methdCorrs[mycorrMethd[theTotalReplications]] fx = matrix(c( rep(traitLoad, 4), rep(0, 16), rep(traitLoad, 4), rep(0, 16), rep(traitLoad, 4), rep(0, 16), rep(traitLoad, 4), rep(c(methodLoad, 0, 0, 0), 4), rep(c(0, methodLoad, 0, 0), 4), rep(c(0, 0, methodLoad, 0), 4), rep(c(0, 0, 0, methodLoad), 4)), ncol = 8) Phi = matrix(c(1, traitCorr, traitCorr, traitCorr, rep(0,4), traitCorr, 1, traitCorr, traitCorr, rep(0,4), traitCorr, traitCorr, 1, traitCorr, rep(0,4), traitCorr, traitCorr, traitCorr, 1, rep(0,4), rep(0,4),1, methdCorr, methdCorr, methdCorr, methdCorr, rep(0,4),1, methdCorr, methdCorr, methdCorr, methdCorr, rep(0,4),1, methdCorr, methdCorr, methdCorr, methdCorr, rep(0,4),1), ncol = 8) mmtm <- sim.structure(fx, Phi, n = SampleSize, raw=T) correMatrix <- mmtm$r colnames(correMatrix) <- c(paste("item", seq(1:16), sep = "")) rownames(correMatrix) <- c(paste("item", seq(1:16), sep = "")) corForModel = correMatrix # establishes the correlation matrix corForModel for model estimation M1 <- try(sem(CTM1, corForModel, SampleSize, maxiter = myMaxIter), silent = FALSE) # SEM model estimation) # tries to estimate the CT(M-1) model M2 <- try(sem(CTCM, corForModel, SampleSize, maxiter = myMaxIter), silent = FALSE) # SEM model estimation) # tries to estimate the CTCM model if(M1$convergence > 2){ convergenceM1 <- c(0) resultsM1 <- as.numeric(c(rep(-1, 12))) } else { # else needs to be in the same line as the last command myModlChiM1 <- try(summary(M1)) convergenceM1 <- as.numeric(M1$convergence) chiM1 <- as.numeric(myModlChiM1$chisq) dfM1 <- as.numeric(myModlChiM1$df) chiM0 <- as.numeric(myModlChiM1$chisqNull) dfM0 <- as.numeric(myModlChiM1$dfNull) GFIM1 <- as.numeric(myModlChiM1$GFI) AGFIM1 <- as.numeric(myModlChiM1$AGFI) RMSEAM1 <- as.numeric(myModlChiM1$RMSEA[1]) CFIM1 <- as.numeric(myModlChiM1$CFI) BICM1 <- as.numeric(myModlChiM1$BIC) SRMRM1 <- as.numeric(myModlChiM1$SRMR) iterM1 <- as.numeric(myModlChiM1$iterations) resultsM1 <- as.numeric(c(convergenceM1, chiM1, dfM1, chiM0, dfM0, GFIM1, AGFIM1, RMSEAM1, CFIM1, BICM1, SRMRM1, iterM1)) } if(M2$convergence > 2){ convergenceM2 <- c(0) resultsM2 <- as.numeric(c(rep(-1, 12))) } else { # else needs to be in the same line as the last command myModlChiM2 <- try(summary(M2)) convergenceM2 <- as.numeric(M2$convergence) chiM2 <- as.numeric(myModlChiM2$chisq) dfM2 <- as.numeric(myModlChiM2$df) chiM0 <- as.numeric(myModlChiM2$chisqNull) dfM0 <- as.numeric(myModlChiM2$dfNull) GFIM2 <- as.numeric(myModlChiM2$GFI) AGFIM2 <- as.numeric(myModlChiM2$AGFI) RMSEAM2 <- as.numeric(myModlChiM2$RMSEA[1]) CFIM2 <- as.numeric(myModlChiM2$CFI) BICM2 <- as.numeric(myModlChiM2$BIC) SRMRM2 <- as.numeric(myModlChiM2$SRMR) iterM2 <- as.numeric(myModlChiM2$iterations) resultsM2 <- as.numeric(c(convergenceM2, chiM2, dfM2, chiM0, dfM0, GFIM2, AGFIM2, RMSEAM2, CFIM2, BICM2, SRMRM2, iterM2)) } designparameters <- c(traitLoad, traitCorr, methodLoad, methdCorr) myResults <- c(designparameters, SampleSize, resultsM1, resultsM2) #, convergenceM3, resultsM3) # return(myResults) } # End of function sizeControlGroup ############ Loop for replications ########################################################## totalRepeats = 100 for(myRepeats in 1:totalRepeats){ myTests <- mapply(sizeControlGroup, myRes, SIMPLIFY = F) # effectSize myTests <- matrix(unlist(myTests), nc=length(myTests[[1]]), byrow=T) colnames(myTests) <- c("trL", "trCorr", "mthL", "methCorr", "n", "convM1", "ChiM1", "dfM1", "Chim0", "dfm0", "GFIM1", "AGFIM1", "RMSEAM1", "CFIM1", "BICM1", "SRMRM1", "iterM1","ConvM2", "ChiM2", "dfM2", "Chim0", "dfm0", "GFIM2", "AGFIM2", "RMSEAM2", "CFIM2", "BICM2", "SRMRM2", "iterM2") write.table(myTests, paste("C:\\Dokumente und Einstellungen\\wood\\Eigene Dateien\\Wood\\papers\\simulations\\rep", myRepeats, sep=""), sep = " ", row.names = F) } [[alternative HTML version deleted]]