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