Modified CP model

11 posts / 0 new

Wed, 02/26/2020 - 14:20

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Modified CP model

Attachment	Size
Model	44.06 KB

Dear all,

I am trying to fit a model (I am not sure if it will be possible) like the one in the attached figure.

I am trying to do it by modifying this script (CP model with one latent factor):

# Fit Common Pathway ACE Model
# ------------------------------------------------------------------------------                                           
nl        <- 1
 
# Matrices ac, cc, and ec to store a, c, and e path coefficients for latent phenotype(s)
pathAl    <- mxMatrix( type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels=labLower("al",nl), lbound=.00001, name="al" )
pathCl    <- mxMatrix( type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels=labLower("cl",nl), lbound=.00001, name="cl" )
pathEl    <- mxMatrix( type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels=labLower("el",nl), lbound=.00001, name="el" )
 
# Matrix and Algebra for constraint on variance of latent phenotype
covarLP   <- mxAlgebra( expression=al %*% t(al) + cl %*% t(cl) + el %*% t(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=.2, labels=labFull("fl",nv,nl), name="fl" )
 
# Matrices A, C, and E compute variance components
covA      <- mxAlgebra( expression=fl %&% (al %*% t(al)) + as %*% t(as), name="A" )
covC      <- mxAlgebra( expression=fl %&% (cl %*% t(cl)) + cs %*% t(cs), name="C" )
covE      <- mxAlgebra( expression=fl %&% (el %*% t(el)) + es %*% t(es), name="E" )

And then fixing the loading factor to zero for the first variable. However, I am struggling to figure out how to create the paths for A C and E to the loading factor and the first variable (lines in red).

I would be really grateful if you could help me with this. I have read other threads but I am still confused about how to do it.

Thank you so much in advance.
With all good wishes.

Tue, 03/17/2020 - 16:39

AdminRobK

Offline

Joined: 01/24/2014 - 12:15

Imagine that there's another

Imagine that there's another latent factor that is indicated only by variable 1. Variable 1's loading onto that new latent factor is fixed to 1.0, and the latent factor's variance is entirely due to the left-hand set of A, C, & E. Then, you'd have two latent factors, and their covariance matrix--say, covarLP-- would be a function of the Cholesky-style A, C, & E paths. Then, the model-expected covariance matrix for the 3 traits would be the sum of the common and unique covariance matrices. The common covariance matrix is equal to the matrix of loadings times covarLP times the transpose of the matrix of loadings. The unique covariance matrix is a function of the "residual" A, C, & E. Note that, for model identification, you will not be able to freely estimate both loadings onto the right-hand common factor.

Alternately, you could always specify your model with MxPaths, in a RAM-type or LISREL-type MxModel. You would still need to make some MxAlgebras to calculate all the quantities of interest, though.

Thu, 03/19/2020 - 18:49

(Reply to #2) #3

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Hi Robert,

Thank you so much for your response. After, thinking a lot about the best model I have tried to fit this model (see attached figure). I think this kind of model gives me all the information that I need.
Here you can see the full script:

# Select Continuous Variables
varsc     <- c('SL')                   # list of continuous variables names
nvc       <- 1                         # number of continuous variables
ntvc      <- nvc*2                     # number of total continuous variables
conVars   <- paste(varsc,c(rep(1,nvc),rep(2,nvc)),sep="")
 
# Select Ordinal Variables
nth       <- 1                         # number of thresholds
varso     <- c('NE',"LB")                   # list of ordinal variables names
nvo       <- 2                         # number of ordinal variables
ntvo      <- nvo*2                     # number of total ordinal variables
ordVars   <- paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="")
ordData   <- twinData
 
# Select Variables for Analysis
vars      <- c('SL','NE',"LB")              # list of variables names
nv        <- nvc+nvo                   # number of variables
ntv       <- nv*2                      # number of total variables
oc        <- c(0,1,1)                    # 0 for continuous, 1 for ordinal variables COMPROBAR
selVars   <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")
 
# Select Covariates for Analysis
covVars   <- c('age', "Sex1" , "Sex2")
nc        <- 3                         # number of covariates
 
# Select Data for Analysis
mzData    <- subset(ordData, Zyg==1|Zyg==3, c(selVars, covVars))
dzData    <- subset(ordData, Zyg==2|Zyg==4|Zyg==5, c(selVars, covVars))
mzDataF   <- mzData 
dzDataF   <- dzData
mzDataF$NE1 <- mxFactor(mzDataF$NE1, levels =c(0,1))
mzDataF$NE2 <- mxFactor(mzDataF$NE2, levels =c(0,1))
dzDataF$NE1 <- mxFactor(dzDataF$NE1, levels =c(0,1))
dzDataF$NE2 <- mxFactor(dzDataF$NE2, levels =c(0,1))
mzDataF$LB1 <- mxFactor(mzDataF$LB1, levels =c(0,1))
mzDataF$LB2 <- mxFactor(mzDataF$LB2, levels =c(0,1))
dzDataF$LB1 <- mxFactor(dzDataF$LB1, levels =c(0,1))
dzDataF$LB2 <- mxFactor(dzDataF$LB2, levels =c(0,1))
 
 
# Set Starting Values
frMV      <- c(TRUE,FALSE,FALSE)             # free status for variables
frCvD     <- diag(frMV,ntv,ntv)        # lower bounds for diagonal of covariance matrix
frCvD[lower.tri(frCvD)] <- TRUE        # lower bounds for below diagonal elements
frCvD[upper.tri(frCvD)] <- TRUE        # lower bounds for above diagonal elements
frCv      <- matrix(as.logical(frCvD),4)
svMe      <- c(2,0,0)                   # start value for means
svPa      <- 1                        # start value for path coefficient
svPaD     <- vech(diag(svPa,nv,nv))    # start values for diagonal of covariance matrix
svPe      <- 1                        # start value for path coefficient for e
svPeD     <- vech(diag(svPe,nv,nv))    # start values for diagonal of covariance matrix
lbPa      <- .000000000001                     # start value for lower bounds
lbPaD     <- diag(lbPa,nv,nv)          # lower bounds for diagonal of covariance matrix
lbPaD[lower.tri(lbPaD)] <- -10         # lower bounds for below diagonal elements
lbPaD[upper.tri(lbPaD)] <- NA          # lower bounds for above diagonal elements
svLTh     <- 0                     # start value for first threshold
svITh     <- 0                         # start value for increments
svTh      <- 0    # start value for thresholds
lbTh      <- 0     # lower bounds for thresholds
 
# Create Labels
labMe     <- paste("mean",vars,sep="_")
labTh     <- paste(paste("t",1:nth,sep=""),rep(varso,each=nth),sep="_")
labBe     <- labFull("beta",nc,1)
 
# ------------------------------------------------------------------------------
# PREPARE MODEL
 
# ACE Model
# Create Matrices for Covariates and linear Regression Coefficients
defAge    <- mxMatrix( type="Full", nrow=1, ncol=1, free=FALSE, labels=c("data.age"), name="Age" )
pathB1     <- mxMatrix( type="Full", nrow=nc, ncol=1, free=TRUE, values=.01, label=c("b11","b12", "b13"), name="b1" )
defSex    <- mxMatrix( type="Full", nrow=1, ncol=6, free=FALSE, labels=c("data.Sex1", "data.Sex1", "data.Sex1","data.Sex2", "data.Sex2", "data.Sex2"), name="Sex" )
pathB2     <- mxMatrix( type="Full", nrow=1, ncol=6, free=TRUE, values=.01, label=c("b21", "b22","b23", "b21", "b22", "b23"), name="b2" )
 
# Create Algebra for expected Mean Matrices
meanG     <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMe, name="meanG" )
expMean   <- mxAlgebra( expression= meanG + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex , name="expMean" )
thinG     <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labTh, name="thinG" )
inc       <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" )
threG     <- mxAlgebra( expression= inc %*% thinG, name="threG" )
 
# Create Matrices for Path Coefficients
pathA     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("a",nv),  name="a" ) 
pathC     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("c",nv),  name="c" )
pathE     <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("e",nv),  name="e" )
 
# Create Algebra for Variance Comptwonts
covA      <- mxAlgebra( expression=a %*% t(a), name="A" )
covC      <- mxAlgebra( expression=c %*% t(c), name="C" ) 
covE <- mxMatrix(
  type="Symm",nrow=nv,
  free=c(T,T,T,F,T,F),
  values=c(0.8,0,0,1,0,1),
  labels=c("e_1_1", "e_2_1", "e_3_1", "e_2_2", "e_3_2", "e_3_3"),
  lbound=c(0.0001,rep(NA,5)),
  name="E")
 
# Create Algebra for expected Variance/Covariance Matrices in MZ & DZ twins
covP      <- mxAlgebra( expression= A+C+E, name="V" )
covMZ     <- mxAlgebra( expression= A+C, name="cMZ" )
covDZ     <- mxAlgebra( expression= 0.5%x%A+ C, 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 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=solve(sqrt(I*A))%&%A, name ="rA" ) #cov2cor()
corC      <- mxAlgebra( expression=solve(sqrt(I*C))%&%C, name ="rC" )
corE      <- mxAlgebra( expression=solve(sqrt(I*E))%&%E, name ="rE" )
 
# Create Data Objects for Multiple Groups
dataMZ    <- mxData( observed=mzDataF, type="raw" )
dataDZ    <- mxData( observed=dzDataF, type="raw" )
 
# Create Expectation Objects for Multiple Groups
expMZ     <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="threG", threshnames=ordVars )
expDZ     <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="threG", threshnames=ordVars )
funML     <- mxFitFunctionML()
 
# Create Model Objects for Multiple Groups
pars      <- list(pathB1,pathB2, meanG, thinG, inc, threG, matI, invSD,
                  pathA, pathC, covA, covC, covE, covP, corA, corC, corE)
defs      <- list(defAge, defSex)
modelMZ   <- mxModel( name="MZ", pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML )
modelDZ   <- mxModel( name="DZ", pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML )
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
 
# Create Algebra for Variance Components
rowVC     <- rep('VC',nv)
colVC     <- rep(c('A','C','E','SA','SC','SE'),each=nv)
estVC     <- mxAlgebra( expression=cbind(A,C,E,A/V,C/V,E/V), name="VC", dimnames=list(rowVC,colVC))
 
# Create Confidence Interval Objects
ciACE     <- mxCI( c ("VC[1,10]", "VC[1,16]","MZ.rA", "MZ.rC", "MZ.rE","MZ.A","MZ.C","MZ.E","VC[2,11]","VC[2,17]", "VC[3,12]", "VC[3,18]" ))
 
# Build Model with Confidence Intervals
modelACE  <- mxModel( "twoACEja", pars, modelMZ, modelDZ, multi, estVC, ciACE )
 
# ------------------------------------------------------------------------------
# RUN MODEL
 
# Run ACE Model
set.seed(180316) 
fitACE    <- mxTryHardOrdinal( modelACE, intervals=F, extraTries=11 )
sumACE    <- summary( fitACE )
 
# Compare with Saturated Model
mxCompare( fit, fitACE )
 
# Print Goodness-of-fit Statistics & Parameter Estimates
fitGofs(fitACE)
fitEsts(fitACE)
 
# ------------------------------------------------------------------------------
# RUN SUBMODELS
 
# Run AE model
modelAE   <- mxModel( fitACE, name="twoAEja" )
modelAE   <- omxSetParameters( modelAE, labels=labLower("c",nv), free=FALSE, values=0 )
set.seed(180316)
fitAE     <- mxTryHardOrdinal( modelAE, intervals=F, extraTries=11 )
mxCompare( fitACE, fitAE )
fitGofs(fitAE)
fitEsts(fitAE)
 
# Run CE model
modelCE   <- mxModel( fitACE, name="twoCEja" )
modelCE   <- omxSetParameters( modelCE, labels=labLower("a",nv), free=FALSE, values=0 )
fitCE     <- mxTryHardOrdinal( modelCE, intervals=F, extraTries=31 )
mxCompare( fitACE, fitCE )
fitGofs(fitCE)
fitEsts(fitCE)
 
# Run E model
modelE    <- mxModel( fitAE, name="twoEja" )
modelE    <- omxSetParameters( modelE, labels=labLower("a",nv), free=FALSE, values=0 )
fitE      <- mxTryHardOrdinal( modelE, intervals=F )
mxCompare( fitAE, fitE )
fitGofs(fitE)
fitEsts(fitE)
 
#-----------------------------------------------
#COMMON AND INDEPENDT PATHWAYMODELS   
# Fit Independent Pathway ADE Model
# ------------------------------------------------------------------------------                                           
nf        <- 1       # number of factors
 
# Matrices ac, cc, and ec to store a, d, and e path coefficients for common factors
pathAc    <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.6,lbound=.00001, labels=labFull("ac",nv,nf), name="ac" )
pathCc    <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.6,lbound=.00001, labels=labFull("cc",nv,nf), name="cc" )
pathEc    <- mxMatrix( type="Full", nrow=nv, ncol=nf, free=TRUE, values=.6,lbound=.00001, labels=labFull("ec",nv,nf), name="ec" )
 
# Matrices as, ds, and es to store a, d, and e path coefficients for specific factors
pathAs    <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=4, labels=labDiag("as",nv), lbound=.00001, name="as" )
pathCs    <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=4, labels=labDiag("cs",nv), lbound=.00001, name="cs" )
pathEs    <- mxMatrix( type="Diag", nrow=nv, ncol=nv, free=TRUE, values=5, labels=labDiag("es",nv), lbound=.00001, name="es" )
 
# Matrices A, D, and E compute variance components
covA      <- mxAlgebra( expression=ac %*% t(ac) + as %*% t(as), name="A" )
covC      <- mxAlgebra( expression=cc %*% t(cc) + cs %*% t(cs), name="C" )
covE      <- mxAlgebra( expression=ec %*% t(ec) + es %*% t(es), name="E" )
 
# Create Model Objects for Multiple Groups
pars      <- list(pathB1, pathB2,meanG, matI,thinG, invSD,inc, threG, matI,
                  pathAc, pathCc, pathEc, pathAs, pathCs, pathEs, covA, covC, covE, covP, corA, corC, corE)
modelMZ   <- mxModel( name="MZ", pars, defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML )
modelDZ   <- mxModel( name="DZ", pars, defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML )
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
 
# Build & Run Model 
modelIP   <- mxModel( "mulIPc", pars, modelMZ, modelDZ, multi )
fitIP     <- mxRun( modelIP, intervals=F )
sumIP     <- summary( fitIP )
mxCompare( fitIP,fitACE )
fitGofs(fitIP)
fitEsts(fitIP)
parameterSpecifications(fitIP)
 
# Fit Common Pathway ACE Model
# ------------------------------------------------------------------------------                                           
nl        <- 2
 
# Matrices ac, cc, and ec to store a, c, and e path coefficients for latent phenotype(s)
pathAl    <- mxMatrix( type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels=labLower("al",nl), lbound=.00001, name="al" )
pathCl    <- mxMatrix( type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels=labLower("cl",nl), lbound=.00001, name="cl" )
pathEl    <- mxMatrix( type="Lower", nrow=nl, ncol=nl, free=TRUE, values=.6, labels=labLower("el",nl), lbound=.00001, name="el" )
 
# Matrix and Algebra for constraint on variance of latent phenotype
covarLP   <- mxAlgebra( expression=al %*% t(al) + cl %*% t(cl) + el %*% t(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=.2, labels=labFull("fl",nv,nl), name="fl" )
 
# Matrices A, C, and E compute variance components
covA      <- mxAlgebra( expression=fl %&% (al %*% t(al)) + as %*% t(as), name="A" )
covC      <- mxAlgebra( expression=fl %&% (cl %*% t(cl)) + cs %*% t(cs), name="C" )
covE      <- mxAlgebra( expression=fl %&% (el %*% t(el)) + es %*% t(es), name="E" )
 
#     Standardized Variance components: TOTAL
stTA    <- mxAlgebra( expression=A/V, name="Th2" )
stTC    <- mxAlgebra( expression=C/V, name="Tc2" )
stTE    <- mxAlgebra( expression=E/V, name="Te2" )
 
#     Standardized Variance components: Latent Factor
stLA    <- mxAlgebra( expression=(al%*%t(al)), name="Lh2" )
stLC    <- mxAlgebra( expression=(cl%*%t(cl)), name="Lc2" )
stLE    <- mxAlgebra( expression=(el%*%t(el)), name="Le2" )
 
#     Standardized Variance components: Common ACE
stCA    <- mxAlgebra( expression=fl%&%(al%*%t(al))/V, name="Ch2" )
stCC    <- mxAlgebra( expression=fl%&%(cl%*%t(cl))/V, name="Cc2" )
stCE    <- mxAlgebra( expression=fl%&%(el%*%t(el))/V, name="Ce2" )
 
#     Standardized Variance components: Specific ACE
stSA    <- mxAlgebra( expression=as%*%t(as)/V, name="Sh2" )
stSC    <- mxAlgebra( expression=cs%*%t(cs)/V, name="Sc2" )
stSE    <- mxAlgebra( expression=es%*%t(es)/V, name="Se2" )
 
#     Standardized Path Estimates: Factor Loadings
 
#stFL <- mxAlgebra( expression=iSD %*% fl, name="stl" )
 
#     Standardized Path Estimates: Specific ACE
stPSA   <- mxAlgebra( expression=as*iSD, name="SpA" )
stPSC   <- mxAlgebra( expression=cs*iSD, name="SpC" )
stPSE   <- mxAlgebra( expression=es*iSD, name="SpE" )
 
# Constrain Variance of Binary Variables
matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )
var1 <- mxConstraint(expression=rbind(V[2,2],V[3,3]) == Unv1, name="Var1")
 
# Create Model Objects for Multiple Groups
pars      <- list(pathB1, pathB2,meanG, matI,thinG, invSD,inc, threG, matI,
                  pathAl, pathCl, pathEl, covarLP, varLP, unit, pathFl, pathAs, pathCs, pathEs, covA, covC, covE, covP, matUnv)
pars2       <- list( stTA, stTC, stTE, stLA, stLC, stLE, stCA, stCC, stCE, stSA, stSC, stSE, stPSA, stPSC, stPSE   )# stl
modelMZ   <- mxModel( name="MZ", pars,pars2,defs, expMean, covMZ, expCovMZ, dataMZ, expMZ, funML )
modelDZ   <- mxModel( name="DZ", pars,pars2,defs, expMean, covDZ, expCovDZ, dataDZ, expDZ, funML )
multi     <- mxFitFunctionMultigroup( c("MZ","DZ") )
 
ciCP     <- mxCI( c ("Th2", "Tc2","Te2", "Lh2", "Lc2","Le2","Ch2","Cc2","Ce2","Sh2", "Sc2", "Se2","SpA","SpC","SpE"  ))
 
# Build & Run Model 
modelCP   <- mxModel( "mulCPc", pars,pars2,  varLP1, modelMZ, modelDZ, multi, ciCP,var1 )
fitCP     <- mxTryHardOrdinal(modelCP, intervals=F, extraTries=100 )
sumCP     <- summary( fitCP )
fitGofs(fitCP)
fitEsts(fitCP)
mxCompare(fitCP,fitACE )
parameterSpecifications(fitCP)
 
modelCP2  <- mxModel( fitCP, name="mulCPc2" )
modelCP2  <- omxSetParameters( modelCP2, labels=c("cl_1_1","cl_2_1","cl_2_2","cs_1_1","cs_2_2","cs_3_3" ), free=FALSE, values=0 )
fitCP2    <- mxTryHardOrdinal( modelCP2, extraTries=11)
mxCompare( fitCP, fitCP2 )
fitGofs(fitCP2)
fitEsts(fitCP2)
 
modelCP3  <- mxModel( fitCP2, name="mulCPc3" )
modelCP3  <- omxSetParameters( modelCP3, labels=c("fl_1_2", "al_2_1", "el_2_1" ), free=FALSE, values=0 )
fitCP3    <- mxTryHardOrdinal( modelCP3, intervals=F, extraTries=101)
mxCompare( fitCP2, fitCP3 )
fitGofs(fitCP3)
fitEsts(fitCP3)
 
 
# Generate List of Parameter Estimates and Derived Quantities using formatOutputMatrices
source("GenEpiHelperFunctions.R")
matCPpaths <- c("al","cl","el","iSD %*% fl","iSD %*% as","iSD %*% cs","iSD %*% es")
labCPpaths <- c("stPathAl","stPathCl","stPathEl","stPathFl","stPathAs","stPathCs","stPathEs")
formatOutputMatrices(fitCP3, matCPpaths, labCPpaths, vars, 3)

However, I am not sure of several things:
1-I fixed the nonshared-environmental variances to 1 in the ACEmodel (you helped me in a different post) but I am not sure how to do it in de CP model (common and specific).
I have used this instead to fix the variance to one
matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )
var1 <- mxConstraint(expression=rbind(V[2,2],V[3,3]) == Unv1, name="Var1")

But as I am having problems to get the CI I think it would be better to change the parameterization

2- I am not sure about this but I think I have to fix either the one threshold, or the regression intercept for the dichotomous variables. Am I Wrong?
I think I know how to fix the one threshold but not how to fix the regression intercept.
I have to do the same in the saturated model, right?
Why do I have to fix the variance and one threshold or the intercept? Is there any paper that I could read to learn more about that?
3-When I run the model, I get this warning:
Warning message:
In mxTryHard(model = model, greenOK = greenOK, checkHess = checkHess, :
The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)

I know that it is common with ordinal variables. I have tried with mxTryhardOrdinal and 101 extra tries and different sets of starting values. The point is that I get different values each time (see the attached file), similar estimates but not exactly the same.

4- I am also not sure about this:
Note that, for model identification, you will not be able to freely estimate both loadings onto the right-hand common factor
Should I change anything?

With all good wishes,
Thank you so much in advance!

File attachments:

Estimates.txt

Model.pdf

Fri, 03/20/2020 - 13:00

tbates

Offline

Joined: 07/31/2009 - 14:25

just umxCP(nFac= 2)

Hi,
Looking at your figure, you're just wanting a 2 common-factor CP model.

umx implements this (see ?umxCP)

m1 = umxCP("new", selDVs = c("gff", "fc", "qol", "hap"), sep = "_T", nFac = 2,
        dzData = dzData, mzData = mzData, tryHard = "yes")

You would then delete the paths you want to delete with umxModify()

m2 = umxModify(m1, c("c_cp_r1c1", "c_cp_r2c2"), name="no_general_c")
m3 = umxModify(m2, c("cp_loadings_r1c2"), name="no_F2_on_var1)
# ... etc.

The figure umx creates, the help and the parameters function in umx can help you find/understand what the path labels are.

parameters(m1, patt="cp")

Tue, 03/24/2020 - 08:31

(Reply to #4) #5

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Thanks

Thank you so much all of you again for your help.

I have used UMX (find attached the figure). This is a really good tool!

However, I am not sure how to get the standardized values.

I have used: cp2 = xmu_standardize_CP(cp1) but they do not seem standardized.

I would also like to know how to fix my previous script.

I have modified these lines:
thinG <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinG" )
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="thinG", threshnames=ordVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="thinG", threshnames=ordVars )
To fix the Threshold and the variance is fixed to one with this:
matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )
var1 <- mxConstraint(expression=rbind(V[2,2],V[3,3]) == Unv1, name="Var1")
And I think the intercept is free:
expMean <- mxAlgebra( expression= meanG + cbind(t(b1%%Age),t(b1%%Age))+ b2*Sex , name="expMean" )

So, I think everything should be OK in the model.

But I still have the Warning message:
In mxTryHard(model = model, greenOK = greenOK, checkHess = checkHess, :
The model does not satisfy the first-order optimality conditions to the required accuracy, and no improved point for the merit function could be found during the final linesearch (Mx status RED)

And I get slightly different values.

I am still confused with this as well:
Note that, for model identification, you will not be able to freely estimate both loadings onto the right-hand common factor

Thank you so much in advance for your help
With all good wishes

File attachments:

UMX-CP.jpeg

Tue, 03/24/2020 - 10:34

tbates

Offline

Joined: 07/31/2009 - 14:25

strip_zero =FALSE, umx_standardize(cp); umxTwinMaker()

Glad you liked umx. The figure you attach is showing standardized values automatically, with leading zeros suppressed. add strip=FALSE to show the decimal place.
You mostly should never need to call an "xmu" function directly (they're all marked "internal, not for end users). So this works for any known model:

cp2 = umx_standardize(cp1)

By default, umxSummary standardizes twin models (and doesn't standardise RAM models). Compare what you get with:

umxSummary(cp1, std=FALSE)

re: "Note that, for model identification, you will not be able to freely estimate both loadings onto the right-hand common factor"

Just that the way you had it drawn, you have two As, Cs, and Es for one latent (and too few manifests). Not identified. (wish images showed up in posts)

Rob mentioned that path-based is a good way to go for custom models. If you're happy to use devtools::install_github("tbates/umx") , I'd be happy to get feedback on umxTwinMaker(), which is designed to make this a whole lot easier. You just sketch one person, and the rules for co-twins and MZ vs. DZ

Tue, 03/24/2020 - 13:11

(Reply to #6) #7

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Thank you so much for your

Thank you so much for your prompt response!

If they are standardized then for my first phenotype:
Ac= 0.77^20.53=0.31
Ec=0.64^20.53=0.22
As= 0.46^2=0.21
Es=0.72^2=0.52
Ptot=0.31+0.22+0.21+0.52=1.26
Atot=0.31+0.21/1.26=0.41
Etot=0.22+0.52/1.26=0.59

(Hope I am not doing something wrong)

These results are slightly different from my univariate analyses (A=32% and E=68%)
However, the ones I got with the full script matches well with my univariate analyses.

Thank you so much for your suggestion I will try with umxTwinMaker()

However, I wonder if there would be a way to get exact results with my script. I have tried with 100 extraTries and the results are very reasonable but I am not 100% confident the script is correct and with values changing from one time to another (small differences)

Thank you so much for explaining me the thing about model specification

With all good wishes.

Tue, 03/24/2020 - 14:03

tbates

Offline

Joined: 07/31/2009 - 14:25

multivariate results need not = univariate results

The results from a multivariate model will often differ from from those of a univariate analysis.

There's more power (the model can see all the phenotypes influenced by a variance component, and use that information in its estimates). Likely, in this case, the univariate model overestimated E because of this.

Tue, 03/24/2020 - 16:44

(Reply to #8) #9

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Thank you so much. I have run

Thank you so much. I have run my multivariate AE model with and without covariates and the discrepancy is due to them,

Can I add covariates to the CP model in UMX?

I would also like to compare the results with those from my script. Do you think if I run the model with a lot of extraTries, like 200 the results would be accurate?

Any thoughts about how to avoid the warning message? I have tried it with all optimizers and tons of starting values.

Thank you so much! I was really stuck with this.

Tue, 03/24/2020 - 18:09

#10

tbates

Offline

Joined: 07/31/2009 - 14:25

umx3.5 supports covariate in twins (forward looking statement)

You'll like umx 3.5, development of which pandemic is accelerating development (home alone)

Plan is by early April, all twin models will support covariates, inc. for thresh vars

Sun, 04/26/2020 - 13:34

(Reply to #10) #11

JuanJMV

Offline

Joined: 07/20/2016 - 13:13

Thank you so much for your

Thank you so much for your response.

Looking forward to trying it! Hope it is going to be available soon.

With all good wishes

Main menu

Navigation

You are here

Modified CP model

Main menu

User login

Navigation

You are here

Modified CP model

Search form