Fix mean latent variable
JuanJMV
Joined: 07/20/2016
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.
I think it's a transformation
Log in or register to post comments
In reply to I think it's a transformation by AdminNeale
Thank you!
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:
require(OpenMx)
V1 <- c(1,2.4,6,8,9,5,6,7,8)
V2 <- c(1,3,6,8.6,9,6,6,7,8)
V3 <- c(1,3,6,5,9.2,5,6,7,8)
V4 <- c(1,7,6,9,9,5.5,6,6,8)
V5 <- c(1,2,2,8,5.4,5,6,5,8)
exampleData <-as.data.frame(cbind(V1,V2,V3,V4,V5))
factorModel <- mxModel("Initial Model",
mxData(exampleData, "raw"),
mxMatrix(type="Full", nrow=1, ncol=5,
free=TRUE, values=0, name="alpha"),
mxMatrix("Full", 1, 1, FALSE, 0, name="mu"),
mxMatrix("Full", 1, 5, TRUE, 0.8, name="lambda"),
mxMatrix("Symm", 1, 1, FALSE, 1, name="phi"),
mxMatrix("Diag", 5, 5, TRUE, 0.6, name="theta"),
mxAlgebra(t(lambda) %*% phi %*% lambda + theta, name="cov"),
mxAlgebra(alpha + mu %*% lambda, name="mean"),
mxExpectationNormal(covariance = "cov", means = "mean", dimnames=names(exampleData)),
mxFitFunctionML())
factorResults <- mxRun(factorModel)
and here the model with one fixed mean and free mean for the factor mean.
require(OpenMx)
V1 <- c(1,2.4,6,8,9,5,6,7,8)
V2 <- c(1,3,6,8.6,9,6,6,7,8)
V3 <- c(1,3,6,5,9.2,5,6,7,8)
V4 <- c(1,7,6,9,9,5.5,6,6,8)
V5 <- c(1,2,2,8,5.4,5,6,5,8)
exampleData <-as.data.frame(cbind(V1,V2,V3,V4,V5))
factorModel <- mxModel("Initial Model",
mxData(exampleData, "raw"),
mxMatrix(type="Full", nrow=1, ncol=5,
free=c(F,T,T,T,T), values=0, name="alpha"),
mxMatrix("Full", 1, 1, TRUE, 0, name="mu"),
mxMatrix("Full", 1, 5, TRUE, 0.8, name="lambda"),
mxMatrix("Symm", 1, 1, FALSE, 1, name="phi"),
mxMatrix("Diag", 5, 5, TRUE, 0.6, name="theta"),
mxAlgebra(t(lambda) %*% phi %*% lambda + theta, name="cov"),
mxAlgebra(alpha + mu %*% lambda, name="mean"),
mxExpectationNormal(covariance = "cov", means = "mean", dimnames=names(exampleData)),
mxFitFunctionML())
factorResults <- mxRun(factorModel)
The fit is the same. Would that be correct? Mu would be the factor mean?
Thank you so much!
Log in or register to post comments
Yes, these are equivalent models
Looks good! But beware of over interpretation of factor means.
Cheers
Mike
Log in or register to post comments
Thank you so much, Mike!
Log in or register to post comments