[R] sem & psych

[chegadesalzburg]

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: http://r.789695.n4.nabble.com/sem-psych-tp2321572p2321572.html
Sent from the R help mailing list archive at Nabble.com.