Dear all,
I am working on a CP model (using a reference script from Boulder). I would like to know if it is possible to fix (or freely estimate) the mean of the latent variables (here the latent factor) using a model specification like this one (see below). Would it be possible by fixing one of the manifest means? I know it is possible (and must be done) to fix the variance to one so that the model is identified but I am also curious if the mean of latent variables can be fixed (or freely estimated). If so how could I do it?
# Select Variables for Analysis vars <- c('family','happy','life','anxdep','somatic','social') nv <- 6 # number of variables ntv <- nv*2 # number of total variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Random Subset to reduce time to Fit Examples testData <- head(nl,n=500) # Select Data for Analysis mzData <- subset(testData, zyg2==1, selVars) dzData <- subset(testData, zyg2==2, selVars) # Generate Descriptive Statistics round(colMeans(mzData,na.rm=TRUE),4) round(colMeans(dzData,na.rm=TRUE),4) round(cov(mzData,use="complete"),4) round(cov(dzData,use="complete"),4) # Set Starting Values svMe <- round(colMeans(mzData,na.rm=TRUE)[1:nv],1) # start value for means svPa <- .25 # start value for path coefficient svPe <- .5 # start value for path coefficient for e lbPa <- .0001 # lower bound for path coefficient lbPaD <- valDiagLU(lbPa,-10,NA,nv) # lower bound for diagonal, lower & upper elements of covariance matrix # ------------------------------------------------------------------------------ # PREPARE MODEL # ACE Model # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labVars("mean",vars), name="meanG" ) # Create Matrices for Variance Components covA <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), label=labLower("VA",nv), name="VA" ) covC <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), label=labLower("VC",nv), name="VC" ) covE <- mxMatrix( type="Symm", nrow=nv, ncol=nv, free=TRUE, values=valDiag(svPa,nv), label=labLower("VE",nv), name="VE" ) # Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins covP <- mxAlgebra( expression= VA+VC+VE, name="V" ) covMZ <- mxAlgebra( expression= VA+VC, name="cMZ" ) covDZ <- mxAlgebra( expression= 0.5%x%VA+VC, name="cDZ" ) expCovMZ <- mxAlgebra( expression= rbind( cbind(V, cMZ), cbind(t(cMZ), V)), name="expCovMZ" ) expCovDZ <- mxAlgebra( expression= rbind( cbind(V, cDZ), cbind(t(cDZ), V)), name="expCovDZ" ) # Create Data Objects for Multiple Groups dataMZ <- mxData( observed=mzData, type="raw" ) dataDZ <- mxData( observed=dzData, type="raw" ) # Create Expectation Objects for Multiple Groups expMZ <- mxExpectationNormal( covariance="expCovMZ", means="meanG", dimnames=selVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list(meanG, covA, covC, covE, covP ) modelMZ <- mxModel( name="MZ", pars, covMZ, expCovMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( name="DZ", pars, covDZ, expCovDZ, dataDZ, expDZ, funML ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Create Algebra for Standardization matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") # Calculate genetic and environmental correlations corA <- mxAlgebra( expression=cov2cor(VA), name ="rA" ) corC <- mxAlgebra( expression=cov2cor(VC), name ="rC" ) corE <- mxAlgebra( expression=cov2cor(VE), name ="rE" ) # Build Model with Confidence Intervals calc <- list( matI, invSD, corA, corC, corE ) modelACE <- mxModel( "mulACEvc", pars, modelMZ, modelDZ, multi, calc ) # ------------------------------------------------------------------------------ # RUN MODEL # Run ACE Model fitACE <- mxRun( modelACE, intervals=F ) sumACE <- summary( fitACE ) parameterSpecifications(fitACE) # ACE Covariance Matrices & Proportions of Variance Matrices matACEcov <- c("VA","VC","VE","V","VA/V","VC/V","VE/V") labACEcov <- c("covA","covC","covE","Var","stCovA","stCovC","stCovE") formatOutputMatrices(fitACE, matACEcov, labACEcov, vars,4) # ACE Correlation Matrices matACEcor <- c("solve(sqrt(I*VA)) %&% VA","solve(sqrt(I*VC)) %&% VC","solve(sqrt(I*VE)) %&% VE") labACEcor <- c("corA","corC","corE") formatOutputMatrices(fitACE, matACEcor, labACEcor, vars, 4) # Fit Common Pathway ACE Model # ------------------------------------------------------------------------------ nl <- 1 svFl <- .6 svFs <- .4 svFe <- .7 svfl <- .2 # Matrices Al, Cl, and El for variance components for latent phenotype(s) covAl <- mxMatrix( type="Symm", nrow=nl, ncol=nl, free=TRUE, values=svFl, label=labLower("Al",nl), name="Al" ) covCl <- mxMatrix( type="Symm", nrow=nl, ncol=nl, free=TRUE, values=svFl, label=labLower("Cl",nl), name="Cl" ) covEl <- mxMatrix( type="Symm", nrow=nl, ncol=nl, free=TRUE, values=svFl, label=labLower("El",nl), name="El" ) # Matrix and Algebra for constraint on variance of latent phenotype covarLP <- mxAlgebra( expression= Al+Cl+El, name="CovarLP" ) varLP <- mxAlgebra( expression= diag2vec(CovarLP), name="VarLP" ) unit <- mxMatrix( type="Unit", nrow=nl, ncol=1, name="Unit") varLP1 <- mxConstraint( expression=VarLP == Unit, name="varLP1") # Matrix f for factor loadings on latent phenotype pathFl <- mxMatrix( type="Full", nrow=nv, ncol=nl, free=TRUE, values=svfl, labels=labFull("fl",nv,nl), name="fl" ) # Matrices As, Cs, and Es for variance components for specific factors covAs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=svFs, labels=labDiag("As",nv), name="As" ) covCs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=svFs, labels=labDiag("Cs",nv), name="Cs" ) covEs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=svFe, labels=labDiag("Es",nv), name="Es" ) # Matrices A, C, and E compute variance components covA <- mxAlgebra( expression=fl %&% Al + As, name="VA" ) covC <- mxAlgebra( expression=fl %&% Cl + Cs, name="VC" ) covE <- mxAlgebra( expression=fl %&% El + Es, name="VE" ) # Create Model Objects for Multiple Groups pars <- list(meanG, matI, invSD, covAl, covCl, covEl, covarLP, varLP, unit, pathFl,covAs, covCs, covEs, covA, covC, covE, covP) modelMZ <- mxModel( name="MZ", pars, covMZ, expCovMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( name="DZ", pars, covDZ, expCovDZ, dataDZ, expDZ, funML ) multi <- mxFitFunctionMultigroup( c("MZ","DZ") ) # Build & Run Model modelCP <- mxModel( "mulCPc", pars, varLP1, modelMZ, modelDZ, multi ) fitCP <- mxRun(modelCP, intervals=F ) fitCP <- mxTryHard(fitCP) sumCP <- summary( fitCP ) mxCompare( fitACE, fitCP ) fitGofs(fitCP) parameterSpecifications(fitCP) # Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices matCPpaths <- c("Al","Cl","El","iSD %*% fl","iSD %*% As","iSD %*% Cs","iSD %*% Es") labCPpaths <- c("stPathAl","stPathCl","stPathEl","stPathFl","stPathAs","stPathCs","stPathEs") formatOutputMatrices(fitCP, matCPpaths, labCPpaths, vars, 4)
Thank you so much in advance.
With all good wishes.
Printer-friendly version
Hi Juan! I think you can fit a model with a latent mean estimated and fixing one of the observed variables' means. There may be a way to simply transform the results of one model into the other, as they should fit the same. I worry that doing so could create difficulties with interpretation. Essentially, the latent factor's mean is being scaled to be the mean-fixed item's mean when freely estimated, divided by the common pathway factor loading to that item. Possibly, comparison between groups on this latent mean would be informative. However, that should be considered in light of which item's mean is fixed as the results could vary as a function of which item is selected to be fixed. And we certainly don't want capitalizing on chance by trying them all to find the one that when fixed supports the user's hypothesis the best!
Dear Mike,
Thank you so much for your response. This is really helpful.
I have read your paper (Estabrook & Neale, 2013; A Comparison of Factor Score Estimation Methods in the Presence of Missing Data: Reliability and an Application to Nicotine Dependence) and adapted the script in the Appendix.
This would be the Factor Model with mean 0 and variance 1:
and here the model with one fixed mean and free mean for the factor mean.
The fit is the same. Would that be correct? Mu would be the factor mean?
Thank you so much!
Juan
Looks good! But beware of over interpretation of factor means.
Cheers
Mike
Thank you so much, Mike!