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.