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]]
______________________________________________
R-help@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
