Attachment | Size |
---|---|
2-IP Model Script.doc | 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 )
**************************************************************************************************
(model 3) Run the Independent Pathway ACE Model
**************************************************************************************************
IPACEFit2 <- mxRun(IPACE2Model)
IPACE2Summ <- summary(IPACEFit2)
IPACESumm
I think what you're after could be accomplished like 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.
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.
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:
That should resolve this error.
Wowww!!! Thanks so much!!! That was exactly the problem. Now it's running ;-)