
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]]

