Attachment | Size |
---|---|
2-IP Model Script.doc [6] | 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