# How to add a second latent common A factor in an independent pathway model

5 posts / 0 new
Offline
Joined: 12/23/2012 - 05:28
How to add a second latent common A factor in an independent pathway model
AttachmentSize
34.5 KB

Dear all,
I have run without problems a multivariate IP model for 6 variables.
The thing is that now I would like to add another latent common A factor that loads only on 3 of the 6 variables to see if this "new" 2-IP factor model (two Ac factors) fits better the data than the 1-IP factor model (just one Ac factor for all the observed variables).
But I think I’m doing something wrong because the model doesn't understand that the second common A factor is loading only in the first 3 observed variables. Can anybody help me? Many many thanks!

In order to be easy to read I've attached a doc with the added code marked in bold.

# SCRIPT---------------------------------------

nf <- 1 # number of A, C and E factors

# Create start values

Stmean <-colMeans(mzData[,1:nv],na.rm=TRUE)

# Create Labels for Column and Diagonal Matrices

AcLabs <- paste("ac",1:nv, sep="")
CcLabs <- paste("cc",1:nv, sep="")
EcLabs <- paste("ec",1:nv, sep="")
AsLabs <- paste("as",1:nv, sep="")
CsLabs <- paste("cs",1:nv, sep="")
EsLabs <- paste("es",1:nv, sep="")
mLabs <- paste("m",1:nv, sep="")

# Matrices to store a, c, and e path coefficients for common ACE factors

pathAc <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.6, labels=AcLabs, name="ac" )
pathCc <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.6, labels=CcLabs, name="cc" )
pathEc <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.6, labels=EcLabs, name="ec" )

# BULDING A COMMON GENETIC FACTOR JUST FOR the first three variables

pathAc2 <- mxMatrix( type="Full", nrow=3, ncol=nf, free=TRUE, values=.6, labels=c("ocd_a2", "hd_a2", "bdd_a2"), name="ac2" )

# Matrices as, cs, and es to store a, c, and e path coefficients for specific factors

pathAs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=AsLabs, name="as" )
pathCs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=CsLabs, name="cs" )
pathEs <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=.6, labels=EsLabs, name="es" )

# Matrices A, C, and E compute variance components

covA <- mxAlgebra( expression=ac %% t(ac) + as %% t(as) + ac2 %% t(ac2), name="A" )
covC <- mxAlgebra( expression=cc %
% t(cc) + cs %% t(cs), name="C" )
covE <- mxAlgebra( expression=ec %
% t(ec) + es %*% t(es), name="E" )

# Algebra to compute total variances and standard deviations (diagonal only)

covP <- mxAlgebra( expression=A+C+E, name="V" )
matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I")
invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD")

# Unstandardized & Standardized ACE Variance components (TOTAL)

rowVars <- Vars
colVars <- rep(c('A','C','E','h2', "h22", 'c2','e2'),each=nv)
estVars1<- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="est", dimnames=list(rowVars,colVars))

# Standardized Variance components: Common ACE

colVars2 <- rep(c('h2c',"h22c", 'c2c','e2c'),each=nv)
estVars2 <- mxAlgebra( expression=cbind((ac%%t(ac))/V, (ac2%%t(ac2))/V, (cc%%t(cc))/V, (ec%%t(ec))/V), name="estCom", dimnames=list(rowVars,colVars2))

# Standardized Variance components: Specific ACE

colVars3 <- rep(c('h2s','c2s','e2s'),each=nv)
estVars3 <- mxAlgebra( expression=cbind( (as%%t(as))/V,(cs%%t(cs))/V,(es%*%t(es))/V), name="estSp", dimnames=list(rowVars,colVars3))

# Standardized Path estimates: Common ACE

colVars4 <- rep(c('sdAc',"sdAc2", 'sdCc', 'sdEc'),each=nv)
estVars4 <- mxAlgebra( expression=cbind( vec2diag(ac)iSD, vec2diag(ac2)iSD, vec2diag(cc)iSD, vec2diag(ec)iSD), name="stCom", dimnames=list(rowVars,colVars4))

# Standardized Path estimates: Specific ACE

colVars5 <- rep(c('sdAs','sdCs','sdEs'),each=nv)
estVars5 <- mxAlgebra( expression=cbind( asiSD, csiSD, es*iSD), name="stSp", dimnames=list(rowVars,colVars5))

# Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins

MeanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=Stmean, labels=mLabs, name="expMean" )
covMZ <- mxAlgebra( expression= rbind( cbind(A+C+E , A+C),
cbind(A+C , A+C+E)), name="expCovMZ" )
covDZ <- mxAlgebra( expression= rbind( cbind(A+C+E , 0.5%x%A+C),
cbind(0.5%x%A+C , A+C+E)), name="expCovDZ" )

# Data objects for Multiple Groups

dataMZ <- mxData( observed=mzData, type="raw" )
dataDZ <- mxData( observed=dzData, type="raw" )

# Objective objects for Multiple Groups

objMZ <- mxFIMLObjective( covariance="expCovMZ", means="expMean", dimnames=selVars )
objDZ <- mxFIMLObjective( covariance="expCovDZ", means="expMean", dimnames=selVars )

# Combine Groups

pars <- list( pathAc, pathCc, pathEc, pathAs, pathCs, pathEs, covA, covC, covE, covP, matI, invSD, estVars1, estVars2, estVars3, estVars4, estVars5 )
modelMZ <- mxModel( pars, MeanG, covMZ, dataMZ, objMZ, name="MZ" )
modelDZ <- mxModel( pars, MeanG, covDZ, dataDZ, objDZ, name="DZ" )
minus2ll <- mxAlgebra( expression=MZ.objective + DZ.objective, name="m2LL" )
obj <- mxAlgebraObjective( "m2LL" )
IPACE2Model <- mxModel( "IPACE2", pars, modelMZ, modelDZ, minus2ll, obj )

# **************************************************************************************************

IPACEFit2 <- mxRun(IPACE2Model)
IPACE2Summ <- summary(IPACEFit2)
IPACESumm

Offline
Joined: 03/01/2013 - 11:03
Try this

I think what you're after could be accomplished like this.

freeAc2 <- c(rep(TRUE, 3), rep(FALSE, 3))
valsAc2 <- c(rep(.6, 3), rep(0, 3))
labsAc2 <- c("ocd_a2", "hd_a2", "bdd_a2", rep(NA, 3))
pathAc2 <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=freeAc2, values=valsAc2, labels=labsAc2, name="ac2" )

algAc <- mxAlgebra(cbind(ac, ac2), name='ac3')
covA    <- mxAlgebra( expression=ac3 %*% t(ac3) + as %*% t(as), name="A" )
# Adjust dependencies subsequent to this.

The idea is to create another Ac matrix that is 6x1 with only the first three values freely estimated. Then combine the two Ac matrices into 'ac3' to make a 6x2 Ac matrix, and use that for computing the A matrix. There may be some dependencies later on in your script that still need to be updated, but I believe this is the right track to start you on.

Offline
Joined: 12/23/2012 - 05:28
Thanks so much for this

Thanks so much for this Hunter.

I adapt the script following your code, but still I have a small problem. I think it's because "ac3" is performed using mxAlgebra and not the mxMatrix as the others (ac, ac2, as, and so on). When I tried to run the model appeared the following error: "In Model "MZ" the following are both named entities and free parameters: "ac2" and "ac3". If you are trying to set a path using an mxAlgebra, then refer to the Algebra with square-bracket notation. i.e., instead of labels=" "ac2", and "ac3" ", use: labels= "ac2" and "ac3" [1,1]"

I don't understand the brakets here because they indicate position, right? but what I need is to estimate different path loadings at a time... so if I put [1,1] it will understand just to estimate one path loading. What I tried was to combine the two matrices "ac" and "ac2" into another big matrix (ac3) using:

ac3 <- list (pathAc, pathAc2)
do.call (rbind, ac3)

but it didn't work, and then I tried

algAc <- mx (cbind(ac, ac2), name="ac3")
ac3 <- as.matrix (algAc)

But, nothing seems to work :-S

Do you have any other ideas to solve this? Anything would help me a lot!!
Thanks so much.

Offline
Joined: 07/31/2009 - 15:26
Re-do labels

No problem! I think you're getting on the wrong track with how you're trying to resolve this new error. The error is complaining about ac2 and ac3 being free parameters and named entities because "ac2" and "ac3" appear in the labels of the "ac" matrix and also as the name of the pathAc2 matrix and the name of the algAc algebra. So, later when you're referring to "ac2" OpenMx doesn't know if you mean the "ac2" free parameter of the "ac" matrix or the "ac2" matrix. They have the same name so OpenMx thinks you're trying to set them equal. In this case, the error message has led you astray. I would recommend changing the AcLabs to something like this:

AcLabs   <- paste("ac",1:nv, "param", sep="")

That should resolve this error.

Offline
Joined: 12/23/2012 - 05:28
Wowww!!! Thanks so much!!!

Wowww!!! Thanks so much!!! That was exactly the problem. Now it's running ;-)