[R] sem & psych
John Fox
jfox at mcmaster.ca
Wed Aug 11 22:03:34 CEST 2010
Dear Gui,
Because you're wrapping the calls to summary() in try(), the errors there
shouldn't stop your simulation. It probably makes sense to test whether the
object returned by try() is of class "try-error", and possibly to specify
the argument silent=TRUE, examining the error messages in your code.
I didn't follow what you're trying to do with mapply(), since there appears
to be only one data argument (myRes).
Regards,
John
--------------------------------
John Fox
Senator William McMaster
Professor of Social Statistics
Department of Sociology
McMaster University
Hamilton, Ontario, Canada
web: socserv.mcmaster.ca/jfox
> -----Original Message-----
> From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org]
On
> Behalf Of Guilherme Wood
> Sent: August-11-10 1:05 PM
> To: r-help at r-project.org
> Subject: [R] sem & psych
>
>
> 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
singuldr:
> 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]]
More information about the R-help
mailing list