# Fix mean latent variable

5 posts / 0 new
Offline
Joined: 07/20/2016 - 13:13
Fix mean latent variable

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

# 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")

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.

Offline
Joined: 03/01/2013 - 14:09
I think it's a transformation

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!

Offline
Joined: 07/20/2016 - 13:13
Thank you!

Dear Mike,

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!

Offline
Joined: 03/01/2013 - 14:09
Yes, these are equivalent models

Juan

Looks good! But beware of over interpretation of factor means.

Cheers
Mike

Offline
Joined: 07/20/2016 - 13:13
Thank you so much, Mike!

Thank you so much, Mike!