Trivariate Model
JuanJMV
Joined: 07/20/2016
I am trying to fit a trivariate model since I think is the best option instead two separately bivariate models but I am not sure.
I have two errors:
1)
Warning messages:
* thinG <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labTh, name="thinG" )
1: data length 3 is not a sub-multiple or multiple of the number of columns 2 for argument 'values' in mxMatrix(type = "Full", nrow = nth, ncol = ntvo, free = TRUE, values = svTh, lbound = lbTh, labels = labTh, name = "thinG")
2: data length 3 is not a sub-multiple or multiple of the number of columns 2 for argument 'lbound' in mxMatrix(type = "Full", nrow = nth, ncol = ntvo, free = TRUE, values = svTh, lbound = lbTh, labels = labTh, name = "thinG")
2)
*fitACE <- mxRun( modelACE, intervals=T )
Error: In model 'twoACEja' the name 'mean_AG_Ln' is used as a free parameter in 'twoACEja.meanG', 'MZ.meanG', and 'DZ.meanG' and as a fixed parameter in 'twoACEja.meanG', 'MZ.meanG', and 'DZ.meanG'
And here the full script
# Select Continuous Variables
varsc <- c('AG','RU') # list of continuous variables names
nvc <- 2 # 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('SL') # list of ordinal variables names
nvo <- 1 # 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
#for (i in ordVars) { ordData[,i] <- cut(twinData[,i],
#breaks=quantile(ordData[,i],(0:(nth+1))/(nth+1),na.rm=TRUE), labels=c(0:nth)) }
# Select Variables for Analysis
vars <- c('AG','RU','SL') # list of variables names
nv <- nvc+nvo # number of variables
ntv <- nv*2 # number of total variables
oc <- c(0,0,1) # 0 for continuous, 1 for ordinal variables
selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")
# Select Covariates for Analysis
#ordData[,'age'] <- ordData[,'age']/100
#ordData <- ordData[-which(is.na(ordData$age)),]
covVars <- c('age', "Sex1" , "Sex2")
nc <- 2 # number of covariates # number of covariates
# Select Data for Analysis
mzData <- subset(ordData, Zyg==1, c(selVars, covVars))
dzData <- subset(ordData, Zyg==2, c(selVars, covVars))
mzDataF <- mzData
dzDataF <- dzData
mzDataF$SL1 <- mxFactor(mzDataF$SL1, levels =c(0,1))
mzDataF$SL2 <- mxFactor(mzDataF$SL2, levels =c(0,1))
dzDataF$SL1 <- mxFactor(dzDataF$SL1, levels =c(0,1))
dzDataF$SL2 <- mxFactor(dzDataF$SL2, levels =c(0,1))
# Generate Descriptive Statistics
apply(mzData,2,myMean)
apply(dzData,2,myMean)
# Set Starting Values
frMV <- c(TRUE,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(1,1,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 <- .0001 # 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 <- -1 # start value for first threshold
svITh <- 1 # start value for increments
svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds
lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # 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"), 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=frMV, 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), lbound=lbPaD, name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("c",nv), lbound=lbPaD, name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("e",nv), lbound=lbPaD, name="e" )
# Create Algebra for Variance Comptwonts
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )
covE <- mxAlgebra( expression=e %*% t(e), 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" )
# Constrain Variance of Binary Variables
matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )
matOc <- mxMatrix( type="Full", nrow=1, ncol=nv, free=FALSE, values=oc, name="Oc" )
var1 <- mxConstraint( expression=diag2vec(Oc%&%V)==Unv1, name="Var1" )
# 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, matUnv, matOc,
pathA, pathC, pathE, 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( "VC[1,1]" )# "VC[1,seq(1,3*nv,nv),(2,seq(1,3*nv,nv)),(2,seq(2,3*nv,nv)))]" )
# Build Model with Confidence Intervals
modelACE <- mxModel( "twoACEja", pars, var1, modelMZ, modelDZ, multi, estVC, ciACE )
# ------------------------------------------------------------------------------
# RUN MODEL
# Run ACE Model
fitACE <- mxRun( modelACE, intervals=T )
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 )
fitAE <- mxRun( modelAE, intervals=T )
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 <- mxRun( modelCE, intervals=T )
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 <- mxRun( modelE, intervals=T )
mxCompare( fitAE, fitE )
fitGofs(fitE)
fitEsts(fitE)
# Print Comparative Fit Statistics
mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) )
Thank you so much I have tried everything but I can not make work it
I have managed to make it
Moreover, I am not sure if I have added the covariables properly.
Here the full script:
# Select Continuous Variables
varsc <- c('AG','RU') # list of continuous variables names
nvc <- 2 # 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 <- 2 # number of thresholds
varso <- c('SL') # list of ordinal variables names
nvo <- 1 # 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('AG','RU','SL') # list of variables names
nv <- nvc+nvo # number of variables
ntv <- nv*2 # number of total variables
oc <- c(0,0,1) # 0 for continuous, 1 for ordinal variables
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 # number of covariates
# Select Data for Analysis
mzData <- subset(ordData, Zyg==1, c(selVars, covVars))
dzData <- subset(ordData, Zyg==2, c(selVars, covVars))
mzDataF <- mzData
dzDataF <- dzData
mzDataF$SL1 <- mxFactor(mzDataF$SL1, levels =c(0,1))
mzDataF$SL2 <- mxFactor(mzDataF$SL2, levels =c(0,1))
dzDataF$SL1 <- mxFactor(dzDataF$SL1, levels =c(0,1))
dzDataF$SL2 <- mxFactor(dzDataF$SL2, levels =c(0,1))
# Set Starting Values
frMV <- c(TRUE,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(1,1,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 <- .0001 # 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 <- -1 # start value for first threshold
svITh <- 1 # start value for increments
svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds
lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # 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), lbound=lbPaD, name="a" )
pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPaD, label=labLower("c",nv), lbound=lbPaD, name="c" )
pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, values=svPeD, label=labLower("e",nv), lbound=lbPaD, name="e" )
# Create Algebra for Variance Comptwonts
covA <- mxAlgebra( expression=a %*% t(a), name="A" )
covC <- mxAlgebra( expression=c %*% t(c), name="C" )
covE <- mxAlgebra( expression=e %*% t(e), 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" )
# Constrain Variance of Binary Variables
matUnv <- mxMatrix( type="Unit", nrow=nvo, ncol=1, name="Unv1" )
matOc <- mxMatrix( type="Full", nrow=1, ncol=nv, free=FALSE, values=oc, name="Oc" )
var1 <- mxConstraint( expression=diag2vec(Oc%&%V)==Unv1, name="Var1" )
# 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, matUnv, matOc,
pathA, pathC, pathE, 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,1]", "MZ.rA", "MZ.rC", "MZ.rE"))
# Build Model with Confidence Intervals
modelACE <- mxModel( "twoACEja", pars, var1, modelMZ, modelDZ, multi, estVC, ciACE )
# ------------------------------------------------------------------------------
# RUN MODEL
# Run ACE Model
fitACE <- mxTryHard( modelACE, intervals=T, extraTries = 31 )
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 )
fitAE <- mxTryHard( modelAE, intervals=T, extraTries = 31 )
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 <- mxTryHard( modelCE, intervals=T, 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 <- mxRun( modelE, intervals=T )
mxCompare( fitAE, fitE )
fitGofs(fitE)
fitEsts(fitE)
# Print Comparative Fit Statistics
mxCompare( fitACE, nested <- list(fitAE, fitCE, fitE) )
Thank you so much again.
Log in or register to post comments
In reply to I have managed to make it by JuanJMV
Only one threshold needed, check other issues
You may be on the right track, but if the variables are binary then you should need at most one threshold in the threshold matrix.
nth <- 1
I couldn't run your script because the variable names don't agree with those in the dataset. Also labFull function is missing.
HTH
Log in or register to post comments
In reply to Only one threshold needed, check other issues by AdminNeale
labFull() is from Hermine's
labFull()is from Hermine's miFunctions.R . I also assume that Juan is using his own data, and not the 'twinData' object that comes with OpenMx.Log in or register to post comments
In reply to I have managed to make it by JuanJMV
ordinal variable
If the ordinal variable is dichotomous, you should only need 1 threshold, which allows you to simplify your script a bit. You also need to either fix the one threshold, or the regression intercept for the dichotomous variable, in order to identify the model. I'd say the simplest modification to make would be to try changing
thinG <- mxMatrix( type="Full", nrow=nth, ncol=ntvo, free=TRUE, values=svTh, lbound=lbTh, labels=labTh, name="thinG" )
to
thinG <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, lbound=lbTh, labels=labTh, name="thinG" )
, changing
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="threG", threshnames=ordVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="threG", threshnames=ordVars )
to
expMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", dimnames=selVars, thresholds="thinG", threshnames=ordVars )
expDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", dimnames=selVars, thresholds="thinG", threshnames=ordVars )
, and changing
pars <- list(pathB1,pathB2, meanG, thinG, inc, threG, matI, invSD, matUnv, matOc,
pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE)
to
pars <- list(pathB1,pathB2, meanG, thinG, matI, invSD, matUnv, matOc,
pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE)
Log in or register to post comments
In reply to ordinal variable by AdminRobK
Thank you! I have made the
thinG <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, lbound=lbTh, labels=labTh, name="thinG" )
Warning message:
data length 3 is not a sub-multiple or multiple of the number of columns 2 for argument 'lbound' in mxMatrix(type = "Full", nrow = 1, ncol = ntvo, free = FALSE, values = 0, lbound = lbTh, labels = labTh, name = "thinG")
Thank you so much!!
Log in or register to post comments
In reply to Thank you! I have made the by JuanJMV
Sorry I forgot to say that I
Thanks
Log in or register to post comments
In reply to Thank you! I have made the by JuanJMV
warning message
That should be benign, since the threshold is fixed. In fact, because it's fixed, you don't need to give it labels or a lower bound. If you want, you could instead do
thinG <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinG" )
Log in or register to post comments
Thank you so much it works
Here you can see my results
ACE
A C E --------------------rA rC rE-----------Rph
SL 0.34 0.35 0.31
AG 0.50 0.12 0.38
RB 0.48 0.22 0.30
SL&AG ---------------------------0.71(0.01,0.99) 0.84(0.121,1) -0.26(-0.56,0.27)------------- 0.38
SL&RB ---------------------------0.12(-0.45,0.67) 0.93(0.21,1) 0.05(-0.39,0.40)----------------0.32
AG&RB ---------------------------0.78(0.59,0.87) 0.98(0.77,1) 0.35(0.24,0.48)-----------------0.66
AE
A E---------------rA rE-------------Rph
SL 0.79 0.21
AG 0.62 0.38
RB 0.71 0.29
SL&AG-------------------0.61(0.46,0.75) -0.21(-0.27,0.07)-------0.37
SL&RB-------------------0.45(0.33,0.56) -0.08(-0.33,0.18)-------0.31
AG&RB-------------------0.81(0.77,0.86) 0.34(0.27,0.42)--------0.66
This is the comparison 27 11475.559 6073 -670.44092 NA NA NA
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoACEja
2 twoACEja twoAEja 21 11487.116 6079 -670.88425 11.556674 6 0.072621853
I am not sure what model should I chose and the criterion to choose one or the other
Also, I need to know if A, C, E are significant in each model. I mean for example in the ACE model for sl C is 0.35 and for ag 0.12 I would like to know if C is significant for SL and AG or only for one of them. I think E is always significant because contain error. I have thought calculate CI and if it contais 0 is not significative but I am not sure.
Finally, Can I compare the rA correlations in each trait (with statistic parameters)? I mean if the rA between 1-2 is bigger than 1-3…
Thank you so much in advance I really appreciate it.
Log in or register to post comments
In reply to Thank you so much it works by JuanJMV
suggestions
Assuming you only want to report the results from a single "best" model, then I recommend selecting the model with the smallest AIC, and reporting confidence intervals for the quantities of interest. BTW, have you fit a saturated model? Is a CE model plausible for your data?
However, I advocate reporting conclusions based on multiple models; see, e.g., my relevant publication for more information.
So just add the needed character strings to the
referenceargument tomxCI(). For instance,ciACE <- mxCI(c("VC[1,1]", "MZ.rA", "MZ.rC", "MZ.rE","MZ.A","MZ.C","MZ.E"))
will request CIs for the elements of the additive-genetic, shared-environmental, and nonshared-environmental covariance matrices, in addition to that for which your script already requests CIs.
It may or may not be, but inferences about E are usually not of interest.
Correct.
Yes. Create a suitable MxAlgebra and request confidence intervals for it, for instance,
rAcompare <- mxAlgebra(cbind(rA[1,2] - rA[1,3], rA[1,2] - rA[2,3], rA[1,3] - rA[2,3]), name="rAcompare")
, place it inside the MZ and DZ submodels, and add "MZ.rAcompare" to the vector of character strings provided for argument
referencetomxCI()when creating objectciACE. This will give you confidence intervals for the differences in genetic correlations.Log in or register to post comments
Once again thank you so much.
I have been working with the saturated model but I cannot make it works. I have tried all possible starting values but It doesn't find a solution. I don't know if I am doing something wrong here is the script.:
# Select Continuous Variables
varsc <- c('AG','RU') # list of continuous variables names
nvc <- 2 # 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('SL') # list of ordinal variables names
nvo <- 1 # 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('AG','RU','SL') # list of variables names
nv <- nvc+nvo # number of variables
ntv <- nv*2 # number of total variables
oc <- c(0,0,1) # 0 for continuous, 1 for ordinal variables
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 # number of covariates
# Select Data for Analysis
mzData <- subset(ordData, Zyg==1, c(selVars, covVars))
dzData <- subset(ordData, Zyg==2, c(selVars, covVars))
mzDataF <- mzData
dzDataF <- dzData
mzDataF$SL1 <- mxFactor(mzDataF$SL1, levels =c(0,1))
mzDataF$SL2 <- mxFactor(mzDataF$SL2, levels =c(0,1))
dzDataF$SL1 <- mxFactor(dzDataF$SL1, levels =c(0,1))
dzDataF$SL2 <- mxFactor(dzDataF$SL2, levels =c(0,1))
# Generate Descriptive Statistics
apply(mzData,2,myMean)
apply(dzData,2,myMean)
#?apply(mzData,2,myCor)
#?apply(dzData,2,myCor)
# Set Starting Values
frMV <- c(TRUE,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(1,1,0) # start value for means
svVa <- c(0.5,0.5,0.5) # start value for variance
lbVa <- .01 # start value for lower bounds
svVas <- diag(svVa,ntv,ntv) # assign start values to diagonal of matrix
lbVas <- diag(lbVa,ntv,ntv) # assign lower bounds to diagonal of matrix
svLTh <- 0 # start value for first threshold
svITh <- 1 # start value for increments
svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds
lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds
# Create Labels
labMeMZ <- paste("meanMZ",selVars,sep="_")
labMeDZ <- paste("meanDZ",selVars,sep="_")
labMeZ <- paste("meanZ",selVars,sep="_")
labCvMZ <- labLower("covMZ",ntv)
labCvDZ <- labLower("covDZ",ntv)
labCvZ <- labLower("covZ",ntv)
labVaMZ <- labDiag("covMZ",ntv)
labVaDZ <- labDiag("covDZ",ntv)
labVaZ <- labDiag("covZ",ntv)
labThMZ <- paste(paste("t",1:nth,"MZ",sep=""),rep(ordVars,each=nth),sep="_")
labThDZ <- paste(paste("t",1:nth,"DZ",sep=""),rep(ordVars,each=nth),sep="_")
labThZ <- paste(paste("t",1:nth,"Z",sep=""),rep(ordVars,each=nth),sep="_")
labBe <- labFull("beta",nc,1)
# ------------------------------------------------------------------------------
# PREPARE MODEL
# Saturated 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 & Threshold Matrices
meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeMZ, name="meanMZ" )
meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMeDZ, name="meanDZ" )
expMeanMZ <- mxAlgebra( expression= meanMZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanMZ" )
expMeanDZ <- mxAlgebra( expression= meanDZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanDZ" )
thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinMZ" )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinDZ" )
inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" )
threMZ <- mxAlgebra( expression= inc %*% thinMZ, name="threMZ" )
threDZ <- mxAlgebra( expression= inc %*% thinDZ, name="threDZ" )
# Create Algebra for expected Covariance Matrices
covMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=svVas, lbound=lbVas, labels=labCvMZ, name="covMZ" )
covDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=svVas, lbound=lbVas, labels=labCvDZ, name="covDZ" )
# 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="covMZ", means="expMeanMZ", dimnames=selVars, thresholds="thinMZ", threshnames=ordVars )
expDZ <- mxExpectationNormal( covariance="covDZ", means="expMeanDZ", dimnames=selVars, thresholds="thinDZ", threshnames=ordVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
pars <- list(pathB1,pathB2, inc )
defs <- list(defAge, defSex)
modelMZ <- mxModel( name="MZ", pars, defs, meanMZ, expMeanMZ, covMZ, thinMZ, dataMZ, expMZ, funML )
modelDZ <- mxModel( name="DZ", pars, defs, meanDZ, expMeanDZ, covDZ, thinDZ, dataDZ, expDZ, funML )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Confidence Interval Objects
ciCor <- mxCI( c('MZ.covMZ','DZ.covDZ' ))
#ciThre <- mxCI( c('MZ.threMZ','DZ.threDZ' ))
# Build Saturated Model with Confidence Intervals
model <- mxModel( "twoSATja", pars, modelMZ, modelDZ, multi, ciCor )
# ------------------------------------------------------------------------------
# RUN MODEL
# Run Saturated Model
fit <- mxTryHard( model, intervals=F, extraTries = 31 )
sum <- summary( fit )
round(fit$output$estimate,4)
fitGofs(fit)
# ------------------------------------------------------------------------------
# RUN SUBMODELS
# Test Significance of Covariates
modelCov <- mxModel( fit, name="testCov" )
modelCov <- omxSetParameters( modelCov, label=labBe, free=FALSE, values=0 )
fitCov <- mxRun( modelCov )
mxCompare( fit, fitCov )
# Constrain expected Thresholds to be equal across twin order
modelEMTVO <- mxModel( fit, name="eqThresTwin" )
svMe <- c(5,0); svVa <- c(2.5,1); svLTh <- 7; svITh <- 1;
for (i in 1:nv) {
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labMeMZ[nv+i],labMeMZ[i]), free=frMV[i], values=svMe[i], newlabels=labMeMZ[i] )
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labMeDZ[nv+i],labMeDZ[i]), free=frMV[i], values=svMe[i], newlabels=labMeDZ[i] )
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labVaMZ[nv+i],labVaMZ[i]), free=frMV[i], values=svVa[i], newlabels=labVaMZ[i] )
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labVaDZ[nv+i],labVaDZ[i]), free=frMV[i], values=svVa[i], newlabels=labVaDZ[i] ) }
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThMZ[nvo*nth+1],labThMZ[1]), free=TRUE, values=svLTh, newlabels=labThMZ[1] )
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThDZ[nvo*nth+1],labThDZ[1]), free=TRUE, values=svLTh, newlabels=labThDZ[1] )
for (i in 2:nth) {for (j in seq(i,nvo*nth,nth)) {
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThMZ[nvo*nth+j],labThMZ[j]), free=TRUE, values=svITh, newlabels=labThMZ[j] )
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThDZ[nvo*nth+j],labThDZ[j]), free=TRUE, values=svITh, newlabels=labThDZ[j] ) }}
fitEMTVO <- mxRun( modelEMTVO, intervals=F )
mxCompare( fit, subs <- list(fitCov, fitEMTVO) )
round(fitEMTVO$output$estimate,4)
fitGofs(fitEMTVO)
# Constrain expected Thresholds to be equal across twin order and zygosity
modelEMTVZ <- mxModel( fitEMTVO, name="eqThresZyg" )
for (i in 1:nv) {
modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labMeMZ[i],labMeDZ[i]), free=frMV[i], values=svMe[i], newlabels=labMeZ[i] )
modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labVaMZ[i],labVaDZ[i]), free=frMV[i], values=svVa[i], newlabels=labVaZ[i] ) }
modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labThMZ[1],labThDZ[1]), free=TRUE, values=svLTh, newlabels=labThZ[1] )
for (i in 2:nth) {for (j in seq(i,nvo*nth,nth)) {
modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labThMZ[j],labThDZ[j]), free=TRUE, values=svITh, newlabels=labThZ[j] ) }}
fitEMTVZ <- mxRun( modelEMTVZ, intervals=F )
mxCompare( fit, subs <- list(fitCov, fitEMTVO, fitEMTVZ) )
round(fitEMTVZ $output$estimate,4)
fitGofs(fitEMTVZ)
Thanks!!!
Log in or register to post comments
In reply to Once again thank you so much. by JuanJMV
specifically?
So, what happens when you try to run it? Also, what do you get from
mxVersion()?Log in or register to post comments
In reply to specifically? by AdminRobK
This is the output:
Begin fit attempt 31 of at maximum 32 tries
Lowest minimum so far: 12836.907293975
OpenMx status code 6 not in list of acceptable status codes, 0
OpenMx status code 6 not in list of acceptable status codes, 0
Not all eigenvalues of Hessian are greater than 0: 787961.505795225, 241712.973930221, 90801.8260774038, 37461.8509249169, 36108.3570058396, 20452.956895597, 16278.6263111995, 13514.5169971915, 12348.4836755664, 7865.0595383285, 7434.41254730365, 7351.21427400012, 6728.35738212829, 5648.20553748093, 4416.33288252297, 3674.66242643252, 3514.28595340669, 3224.37985610757, 3058.62507833283, 2824.51753405145, 2653.72023767685, 2290.52727017331, 2143.71360085118, 1845.17794007907, 1758.30696016368, 1351.05525266071, 1226.13780443413, 1173.91950744854, 867.931406722967, 807.202050148447, 615.874189695127, 502.241971723273, 418.459657188584, 364.354936896746, 36.5272816693136, 30.8297966297892, 25.5256524187929, 16.5921347221574, 12.4923808408779, 11.2477690283236, 5.48989588395323, 5.08592627903294, 4.88676405380493, 2.5967232615173, 2.08173667306169, 1.82647729161902, 1.52624731830076e-05, -8.33106327891148e-05
Begin fit attempt 32 of at maximum 32 tries
Fit attempt worse than current best: 12836.908910494 vs 12836.907293975
Retry limit reached
Start values from best fit:
0.10631708710754,0.0362139006101167,-0.0613654041393264,0.00569082513418116,-0.185581069609575,-0.344162428087438,1.06767187921709,-8.36632390547581,0.942344815595082,1.07312588323597,0.408909142438524,2.2005873612565,1.30673456966261,41.1843125633866,0.313146221049428,0.167776661842215,0.85467485365846,0.242455084108765,0.292201400080796,0.869137784815067,0.220708232633343,0.521510782740128,0.274754384528732,0.118846729885699,3.69004028161039,0.0849444409204678,0.103993475707172,1.09385843177197,-17.2859999506245,0.966720410717319,1.26641979849299,0.441386464784935,3.94902382148772,1.9668200201659,185.096253954615,0.164428802392192,0.108368824726336,1.08268432674425,0.219192741624143,0.173675315487782,2.13781627870702,0.227229733482567,0.578122306941328,0.197451965247457,0.0792066276546675,3.46936267863339,0.0862225357533247,0.088883979511605
And here my version:
OpenMx version: 2.7.9 [GIT v2.7.9]
R version: R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32
Default optimiser: NPSOL
Thank you so much!!!
Log in or register to post comments
In reply to Once again thank you so much. by JuanJMV
Because your ordinal variable
selVars, look atmeanMZ$free
meanDZ$free
Your script only fixes SL's intercept for twin #2, and you fix the intercepts for twin #1's RU and twin #2's AG, both of which should be free. [Edit: to clarify, all the intercepts should be free if you're going to fix the threshold.] A similar situation holds for the diagonal elements of
covMZandcovDZ: only the SL variances should be fixed (and most people fix them to 1.0, but that's arbitrary).Aside from that, are the start values you're using reasonable? You already have estimates of the regression coefficients, and of the MZ and DZ covariance-matrix elements, from your other models. You could use those as start values. You could also run
cov()on your MZ and DZ data matrices, with argumentuse="pair".Since you have an ordinal variable, consider using
mxTryHardOrdinal(). Also, be advised that parameterizing covariance matrices as a symmetric matrix of free parameters doesn't always work very well with with the TryHard functions, since the matrices are not guaranteed to be positive-definite at the start of an attempt.Log in or register to post comments
In reply to Because your ordinal variable by AdminRobK
Thanks!!!
I have made the changes but I am not sure If I have made the changes properly.
Now I think the saturated model is OK here the output.
Solution found
Running final fit, for Hessian and/or standard errors and/or confidence intervals
Warning messages generated from final fit for final fit for Hessian/SEs/CIs
Start values from best fit:
-0.0691057959162112,-0.0521420158570258,0.0384309413852813,-0.240706583661364,-0.298725394917852,-0.173318222954087,2.73579579403297,1.79440658210804,11.4497520581424,2.69020253759794,1.75369578879318,-0.782094367395943,1.04792129978035,0.40055361067745,2.73025605742978,1.60148107506046,61.4186624152812,0.303395019112009,0.163145240710836,1.15259082182226,0.230520829429149,0.288949826490561,1.08859911323993,0.215686628575476,0.517151654758388,0.332738963798313,0.130010886719768,4.50785881955327,0.123359536402215,0.110825292980583,24.1352589955029,0.478348813364634,2.97100594513919,1.95886486626441,-0.402130467155594,2.9284481924813,1.92691918020341,-0.311361676981213,1.25042143753879,0.436052231591234,0.254634500568326,0.120286524326224,0.663430914290249,0.156510657389065,0.103573383603742,0.0785470263867377,0.225754745967759,0.175168905853076,0.128880192532778,0.230800499171072,0.582082242296643,0.213125789106686,0.0840996230298986,0.191733399855721,0.0926517086525632,0.08013306011578...
But in the Test Significance of Covariates I have this warning:
Running testCov with 58 parameters
Warning message: 58 12559.294 6041 477.29376 NA NA NA
In model 'testCov' Optimizer returned a non-zero status code 6. 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)
> mxCompare( fit, fitCov )
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoSATja
2 twoSATja testCov 58 12559.294 6041 477.29353 -0.00023061521 0 NA
Warning message:
In collectStatistics(nextBaseSummary, nextCompareSummary) :
Model 'twoSATja' has the same degrees of freedom as testCov which means the models are not nested and should not be compared with the likelihood ratio test. Compare these models using the information criteria statistics
>
And in the # Constrain expected Thresholds to be equal across twin order
Error in omxSetParameters(modelEMTVO, label = c(labThMZ[nvo * nth + 1], :
The following labels are not present in the model (use 'strict' = FALSE to ignore): 't1MZ_SLEEP2' and 't1MZ_SLEEP1'
The full script here:
# Select Continuous Variables
varsc <- c('AG','RU') # list of continuous variables names
nvc <- 2 # 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('SL') # list of ordinal variables names
nvo <- 1 # 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('AG','RU','SL') # list of variables names
nv <- nvc+nvo # number of variables
ntv <- nv*2 # number of total variables
oc <- c(0,0,1) # 0 for continuous, 1 for ordinal variables
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 # number of covariates
# Select Data for Analysis
mzData <- subset(ordData, Zyg==1, c(selVars, covVars))
dzData <- subset(ordData, Zyg==2, c(selVars, covVars))
mzDataF <- mzData
dzDataF <- dzData
mzDataF$SL1 <- mxFactor(mzDataF$SL1, levels =c(0,1))
mzDataF$SL2 <- mxFactor(mzDataF$SL2, levels =c(0,1))
dzDataF$SL1 <- mxFactor(dzDataF$SL1, levels =c(0,1))
dzDataF$SL2 <- mxFactor(dzDataF$SL2, levels =c(0,1))
# Generate Descriptive Statistics
apply(mzData,2,myMean)
apply(dzData,2,myMean)
# Set Starting Values
frMV <- c(TRUE,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(1,1,0) # start value for means
svVa <- c(0.5,0.5,0.5) # start value for variance
lbVa <- .01 # start value for lower bounds
svVas <- diag(svVa,ntv,ntv) # assign start values to diagonal of matrix
lbVas <- diag(lbVa,ntv,ntv) # assign lower bounds to diagonal of matrix
svLTh <- 0 # start value for first threshold
svITh <- 1 # start value for increments
svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds
lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds
# Create Labels
labMeMZ <- paste("meanMZ",selVars,sep="_")
labMeDZ <- paste("meanDZ",selVars,sep="_")
labMeZ <- paste("meanZ",selVars,sep="_")
labCvMZ <- labLower("covMZ",ntv)
labCvDZ <- labLower("covDZ",ntv)
labCvZ <- labLower("covZ",ntv)
labVaMZ <- labDiag("covMZ",ntv)
labVaDZ <- labDiag("covDZ",ntv)
labVaZ <- labDiag("covZ",ntv)
labThMZ <- paste(paste("t",1:nth,"MZ",sep=""),rep(ordVars,each=nth),sep="_")
labThDZ <- paste(paste("t",1:nth,"DZ",sep=""),rep(ordVars,each=nth),sep="_")
labThZ <- paste(paste("t",1:nth,"Z",sep=""),rep(ordVars,each=nth),sep="_")
labBe <- labFull("beta",nc,1)
# ------------------------------------------------------------------------------
# PREPARE MODEL
# Saturated 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 & Threshold Matrices
meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeMZ, name="meanMZ" )
meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeDZ, name="meanDZ" )
expMeanMZ <- mxAlgebra( expression= meanMZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanMZ" )
expMeanDZ <- mxAlgebra( expression= meanDZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanDZ" )
thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinMZ" )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinDZ" )
inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" )
threMZ <- mxAlgebra( expression= inc %*% thinMZ, name="threMZ" )
threDZ <- mxAlgebra( expression= inc %*% thinDZ, name="threDZ" )
# Create Algebra for expected Covariance Matrices
covMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=svVas, lbound=lbVas, labels=labCvMZ, name="covMZ" )
covDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=frCv, values=svVas, lbound=lbVas, labels=labCvDZ, name="covDZ" )
# 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="covMZ", means="expMeanMZ", dimnames=selVars, thresholds="thinMZ", threshnames=ordVars )
expDZ <- mxExpectationNormal( covariance="covDZ", means="expMeanDZ", dimnames=selVars, thresholds="thinDZ", threshnames=ordVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
pars <- list(pathB1,pathB2, inc )
defs <- list(defAge, defSex)
modelMZ <- mxModel( name="MZ", pars, defs, meanMZ, expMeanMZ, covMZ, thinMZ, dataMZ, expMZ, funML )
modelDZ <- mxModel( name="DZ", pars, defs, meanDZ, expMeanDZ, covDZ, thinDZ, dataDZ, expDZ, funML )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Confidence Interval Objects
ciCor <- mxCI( c('MZ.covMZ','DZ.covDZ' ))
#ciThre <- mxCI( c('MZ.threMZ','DZ.threDZ' ))
# Build Saturated Model with Confidence Intervals
model <- mxModel( "twoSATja", pars, modelMZ, modelDZ, multi, ciCor )
# ------------------------------------------------------------------------------
# RUN MODEL
# Run Saturated Model
fit <- mxTryHardOrdinal( model, intervals=F, extraTries = 31 )
sum <- summary( fit )
round(fit$output$estimate,4)
fitGofs(fit)
# ------------------------------------------------------------------------------
# RUN SUBMODELS
# Test Significance of Covariates
modelCov <- mxModel( fit, name="testCov" )
modelCov <- omxSetParameters( modelCov, label=labBe, free=FALSE, values=0 )
fitCov <- mxRun( modelCov )
mxCompare( fit, fitCov )
# Constrain expected Thresholds to be equal across twin order
modelEMTVO <- mxModel( fit, name="eqThresTwin" )
svMe <- c(5,0); svVa <- c(2.5,1); svLTh <- 7; svITh <- 1;
for (i in 1:nv) {
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labMeMZ[nv+i],labMeMZ[i]), free=frMV[i], values=svMe[i], newlabels=labMeMZ[i] )
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labMeDZ[nv+i],labMeDZ[i]), free=frMV[i], values=svMe[i], newlabels=labMeDZ[i] )
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labVaMZ[nv+i],labVaMZ[i]), free=frMV[i], values=svVa[i], newlabels=labVaMZ[i] )
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labVaDZ[nv+i],labVaDZ[i]), free=frMV[i], values=svVa[i], newlabels=labVaDZ[i] ) }
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThMZ[nvo*nth+1],labThMZ[1]), free=TRUE, values=svLTh, newlabels=labThMZ[1] )
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThDZ[nvo*nth+1],labThDZ[1]), free=TRUE, values=svLTh, newlabels=labThDZ[1] )
for (i in 2:nth) {for (j in seq(i,nvo*nth,nth)) {
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThMZ[nvo*nth+j],labThMZ[j]), free=TRUE, values=svITh, newlabels=labThMZ[j] )
modelEMTVO <- omxSetParameters( modelEMTVO, label=c(labThDZ[nvo*nth+j],labThDZ[j]), free=TRUE, values=svITh, newlabels=labThDZ[j] ) }}
fitEMTVO <- mxRun( modelEMTVO, intervals=F )
mxCompare( fit, subs <- list(fitCov, fitEMTVO) )
round(fitEMTVO$output$estimate,4)
fitGofs(fitEMTVO)
# Constrain expected Thresholds to be equal across twin order and zygosity
modelEMTVZ <- mxModel( fitEMTVO, name="eqThresZyg" )
for (i in 1:nv) {
modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labMeMZ[i],labMeDZ[i]), free=frMV[i], values=svMe[i], newlabels=labMeZ[i] )
modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labVaMZ[i],labVaDZ[i]), free=frMV[i], values=svVa[i], newlabels=labVaZ[i] ) }
modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labThMZ[1],labThDZ[1]), free=TRUE, values=svLTh, newlabels=labThZ[1] )
for (i in 2:nth) {for (j in seq(i,nvo*nth,nth)) {
modelEMTVZ <- omxSetParameters( modelEMTVZ, label=c(labThMZ[j],labThDZ[j]), free=TRUE, values=svITh, newlabels=labThZ[j] ) }}
fitEMTVZ <- mxRun( modelEMTVZ, intervals=F )
mxCompare( fit, subs <- list(fitCov, fitEMTVO, fitEMTVZ) )
round(fitEMTVZ $output$estimate,4)
fitGofs(fitEMTVZ)
Thank you again and sorry because I am really new and I can not understand some things
Log in or register to post comments
In reply to Thanks!!! by JuanJMV
suggestions
covMZandcovDZdifferently. Try this:covfree <- matrix(TRUE,6,6)
covfree[3,3] <- FALSE
covfree[6,6] <- FALSE
covMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=covfree, values=svVas, lbound=lbVas, labels=labCvMZ, name="covMZ" )
covDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=covfree, values=svVas, lbound=lbVas, labels=labCvDZ, name="covDZ" )
The relevant syntax is:
modelCov <- mxModel( fit, name="testCov" )
modelCov <- omxSetParameters( modelCov, label=labBe, free=FALSE, values=0 )
fitCov <- mxRun( modelCov )
mxCompare( fit, fitCov )
I'm surprised it ran without
omxSetParameters()throwing an error. Your script never actually useslabBein the model. All the regression coefficients' labels start with "b", and not "beta_", as inlabBe. As a result, the call toomxSetParameters()doesn't actually change anything, andfitCovis justfit, re-run from its solution. That's why it has the same degrees of freedom. Change the call toomxSetParameters()to something like this:modelCov <- omxSetParameters( modelCov, labels=c("b11","b12","b13","b21","b22","b23"), free=FALSE, values=0 )
Your script doesn't give any labels to the matrix of thresholds. Try redefining them as:
thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )
To constrain thresholds across twin order (I think this will work):
modelEMTVO <- omxSetParameters(model=fit, labels=c("tauMZ1","tauMZ2","tauDZ1","tauDZ2"), newlabels=c("tauMZ","tauMZ","tauDZ","tauDZ"), name="eqThreshTwin")
modelEMTVO <- omxAssignFirstParameters(modelEMTVO)
To further constrain thresholds across zygosity:
modelEMTVZ <- omxSetParameters(model=fitEMTVO, labels=c("tauMZ","tauDZ"), newlabels=c("tau","tau"), name="eqThresZyg")
modelEMTVZ <- omxAssignFirstParameters(modelEMTVZ)
Log in or register to post comments
Thank you so much!!!
I have made the changes. Also I had to change the optimizer to : mxOption(NULL,"Default optimizer","SLSQP") because with the normal optimizer I always got the same error:
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 don't know if it is ok.
Here you can see the results :
base comparison ep minus2LL df AIC diffLL diffdf p 60 11413.788 6039 -664.21205 NA NA NA
1 twoSATja
2 twoSATja testCov 54 11479.071 6045 -610.92935 6.5282698e+01 6 3.7769626e-12
3 twoSATja eqThreshTwin 58 11413.788 6041 -668.21204 4.0857321e-06 2 9.9999796e-01
4 twoSATja eqThresZyg 58 11413.788 6041 -668.21205 -1.0619260e-08 2 1.0000000e+00
And the full script
# Select Continuous Variables
varsc <- c('AG','RU') # list of continuous variables names
nvc <- 2 # 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('SL') # list of ordinal variables names
nvo <- 1 # 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('AG','RU','SL') # list of variables names
nv <- nvc+nvo # number of variables
ntv <- nv*2 # number of total variables
oc <- c(0,0,1) # 0 for continuous, 1 for ordinal variables
selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="")
#ordData <- ordData[-which(is.na(ordData$age)),]
covVars <- c('age', "Sex1" , "Sex2")
nc <- 3 # number of covariates # number of covariates
# Select Data for Analysis
mzData <- subset(ordData, Zyg==1, c(selVars, covVars))
dzData <- subset(ordData, Zyg==2, c(selVars, covVars))
mzDataF <- mzData
dzDataF <- dzData
mzDataF$SL1 <- mxFactor(mzDataF$SL1, levels =c(0,1))
mzDataF$SL2 <- mxFactor(mzDataF$SL2, levels =c(0,1))
dzDataF$SL1 <- mxFactor(dzDataF$SL1, levels =c(0,1))
dzDataF$SL2 <- mxFactor(dzDataF$SLP2, levels =c(0,1))
# Generate Descriptive Statistics
apply(mzData,2,myMean)
apply(dzData,2,myMean)
# Set Starting Values
frMV <- c(TRUE,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(1,1,0) # start value for means
svVa <- c(0.5,0.5,0.5) # start value for variance
lbVa <- .01 # start value for lower bounds
svVas <- diag(svVa,ntv,ntv) # assign start values to diagonal of matrix
lbVas <- diag(lbVa,ntv,ntv) # assign lower bounds to diagonal of matrix
svLTh <- 0 # start value for first threshold
svITh <- 1 # start value for increments
svTh <- matrix(rep(c(svLTh,(rep(svITh,nth-1)))),nrow=nth,ncol=nv) # start value for thresholds
lbTh <- matrix(rep(c(-3,(rep(0.001,nth-1))),nv),nrow=nth,ncol=nv) # lower bounds for thresholds
# Create Labels
labMeMZ <- paste("meanMZ",selVars,sep="_")
labMeDZ <- paste("meanDZ",selVars,sep="_")
labMeZ <- paste("meanZ",selVars,sep="_")
labCvMZ <- labLower("covMZ",ntv)
labCvDZ <- labLower("covDZ",ntv)
labCvZ <- labLower("covZ",ntv)
labVaMZ <- labDiag("covMZ",ntv)
labVaDZ <- labDiag("covDZ",ntv)
labVaZ <- labDiag("covZ",ntv)
labThMZ <- paste(paste("t",1:nth,"MZ",sep=""),rep(ordVars,each=nth),sep="_")
labThDZ <- paste(paste("t",1:nth,"DZ",sep=""),rep(ordVars,each=nth),sep="_")
labThZ <- paste(paste("t",1:nth,"Z",sep=""),rep(ordVars,each=nth),sep="_")
labBe <- labFull("beta",nc,1)
# ------------------------------------------------------------------------------
# PREPARE MODEL
# Saturated 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 & Threshold Matrices
meanMZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeMZ, name="meanMZ" )
meanDZ <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=TRUE, values=svMe, labels=labMeDZ, name="meanDZ" )
expMeanMZ <- mxAlgebra( expression= meanMZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanMZ" )
expMeanDZ <- mxAlgebra( expression= meanDZ + cbind(t(b1%*%Age),t(b1%*%Age))+ b2*Sex, name="expMeanDZ" )
thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )
inc <- mxMatrix( type="Lower", nrow=nth, ncol=nth, free=FALSE, values=1, name="inc" )
threMZ <- mxAlgebra( expression= inc %*% thinMZ, name="threMZ" )
threDZ <- mxAlgebra( expression= inc %*% thinDZ, name="threDZ" )
# Create Algebra for expected Covariance Matrices
covfree <- matrix(TRUE,6,6)
covfree[3,3] <- FALSE
covfree[6,6] <- FALSE
covMZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=covfree, values=svVas, lbound=lbVas, labels=labCvMZ, name="covMZ" )
covDZ <- mxMatrix( type="Symm", nrow=ntv, ncol=ntv, free=covfree, values=svVas, lbound=lbVas, labels=labCvDZ, name="covDZ" )
# 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="covMZ", means="expMeanMZ", dimnames=selVars, thresholds="thinMZ", threshnames=ordVars )
expDZ <- mxExpectationNormal( covariance="covDZ", means="expMeanDZ", dimnames=selVars, thresholds="thinDZ", threshnames=ordVars )
funML <- mxFitFunctionML()
# Create Model Objects for Multiple Groups
pars <- list(pathB1,pathB2, inc )
defs <- list(defAge, defSex)
modelMZ <- mxModel( name="MZ", pars, defs, meanMZ, expMeanMZ, covMZ, thinMZ, dataMZ, expMZ, funML )
modelDZ <- mxModel( name="DZ", pars, defs, meanDZ, expMeanDZ, covDZ, thinDZ, dataDZ, expDZ, funML )
multi <- mxFitFunctionMultigroup( c("MZ","DZ") )
# Create Confidence Interval Objects
ciCor <- mxCI( c('MZ.covMZ','DZ.covDZ' ))
#ciThre <- mxCI( c('MZ.threMZ','DZ.threDZ' ))
# Build Saturated Model with Confidence Intervals
model <- mxModel( "twoSATja", pars, modelMZ, modelDZ, multi, ciCor )
# ------------------------------------------------------------------------------
# RUN MODEL
# Run Saturated Model
fit <- mxTryHardOrdinal( model, intervals=F, extraTries = 101, finetuneGradient = F,exhaustive = F)
sum <- summary( fit )
round(fit$output$estimate,4)
fitGofs(fit)
# ------------------------------------------------------------------------------
# RUN SUBMODELS
# Test Significance of Covariates
modelCov <- mxModel( fit, name="testCov" )
modelCov <- omxSetParameters( modelCov, labels=c("b11","b12","b13","b21","b22","b23"), free=FALSE, values=0 )
fitCov <- mxTryHardOrdinal( modelCov, extraTries = 31 )
mxCompare( fit, fitCov )
# Constrain expected Thresholds to be equal across twin order
modelEMTVO <- mxModel( fit, name="eqThresTwin" )
modelEMTVO <- omxSetParameters(model=fit, labels=c("tauMZ1","tauMZ2","tauDZ1","tauDZ2"), newlabels=c("tauMZ","tauMZ","tauDZ","tauDZ"), name="eqThreshTwin")
modelEMTVO <- omxAssignFirstParameters(modelEMTVO)
fitEMTVO <- mxTryHardOrdinal( modelEMTVO, intervals=F, extraTries=10 )
mxCompare( fit, subs <- list(fitCov, fitEMTVO) )
round(fitEMTVO$output$estimate,4)
fitGofs(fitEMTVO)
# Constrain expected Thresholds to be equal across twin order and zygosity
modelEMTVZ <- omxSetParameters(model=fitEMTVO, labels=c("tauMZ","tauDZ"), newlabels=c("tau","tau"), name="eqThresZyg")
modelEMTVZ <- omxAssignFirstParameters(modelEMTVZ)
modelEMTVZ <- mxModel( fitEMTVO, name="eqThresZyg" )
fitEMTVZ <- mxTryHardOrdinal( modelEMTVZ, intervals=F, extraTries=10)
mxCompare( fit, subs <- list(fitCov, fitEMTVO, fitEMTVZ) )
round(fitEMTVZ $output$estimate,4)
fitGofs(fitEMTVZ)
I really appreciate it!
Log in or register to post comments
In reply to Thank you so much!!! by JuanJMV
eqThreshTwin vs. eqThresZyg
Something is still not quite right here. The fourth model, "eqThresZyg", should have one less free parameter than the third, "eqThreshTwin". Maybe run
omxGetParameters()on the model named "eqThresZyg" before you run it, to check whether my suggested syntax did what it was supposed to?Log in or register to post comments
In reply to eqThreshTwin vs. eqThresZyg by AdminRobK
You were right I had a
This line was before mxModel
modelEMTVZ <- omxSetParameters(model=fitEMTVO, labels=c("tauMZ","tauDZ"), newlabels=c("tau","tau"), name="eqThresZyg")
Thanks!!
Here you can see the correct results saturated model and saturated vs ADE and ACE
> mxCompare( fit, subs <- list(fitCov, fitEMTVO, fitEMTVZ) ) 60 11413.788 6039 -664.21205 NA NA NA
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoSATja
2 twoSATja testCov 54 11479.071 6045 -610.92935 6.5282698e+01 6 3.7769626e-12
3 twoSATja eqThreshTwin 58 11413.788 6041 -668.21205 1.7364073e-08 2 9.9999999e-01
4 twoSATja eqThresZyg 57 11413.788 6042 -670.21205 4.1691237e-09 3 1.0000000e+00
>
> mxCompare( fit, fitACE ) 60 11413.788 6039 -664.21205 NA NA NA
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoSATja
2 twoSATja twoACEja 27 11475.559 6073 -670.44092 61.771126 34 0.0024790555
> mxCompare( fit, fitADE ) 60 11413.788 6039 -664.21205 NA NA NA
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoSATja
2 twoSATja twoADEja 27 11474.970 6073 -671.02971 61.182339 34 0.0028790296
And here ADE vs AE and ACE vs AE
base comparison ep minus2LL df AIC diffLL diffdf p 27 11474.970 6073 -671.02971 NA NA NA
1 twoADEja
2 twoADEja twoAEja 21 11487.115 6079 -670.88451 12.145205 6 0.05880676
base comparison ep minus2LL df AIC diffLL diffdf p 27 11475.559 6073 -670.44092 NA NA NA
1 twoACEja
2 twoACEja twoAEja 21 11487.116 6079 -670.88425 11.556674 6 0.072621853
So, what model is the best? Should I explain all as you mentioned before?
Also, I would like to know if it is possible:
1-Can I get CI for phenotypic correlations? I have been trying to get but I don't know how since phenotypic correlations are not in VC I mean I have got SA,SD and SC since they are in VC but I don't know how get the CI's phenotypic correlations since to get rPH y use cov2cor.
2-As you told me in order to know if A, C or D are significant I have calculated CI's for A C and D but for example in the ADE model for one variable D(unstandardized) is 0.05(0,0.15) I mean is weird that I don't have any negative value since I think 0.05 should be non-significant. I have a big sample but I don't know if I am doing something wrong.
3- Also, I would like to know where and how should I get the MZ and DZ correlations in the saturated model? with covariables? or in the ACE or ADE model? I mean I need the cross twin and cross trait correlations
Thank you so much. I really need it and I am learning a lot.
Log in or register to post comments
In reply to You were right I had a by JuanJMV
cov2cor() works in MxAlgebras
cov2cor()works in MxAlgebras, too. You can make MxAlgebras for the correlation matrices, e.g.corMZ <- mxAlgebra(cov2cor(covMZ), name="corMZ")
corDZ <- mxAlgebra(cov2cor(covDZ), name="corDZ")
and put them into the MZ and DZ submodels, respectively. Then, request confidence intervals for those two algebras; in your case, I think that would entail adding "MZ.corMZ" and "DZ.corDZ" to the character vector provided for argument
referencetomxCI(). You could also instead request CIs only for specific elements of the correlation matrices (e.g., "MZ.corMZ[1,2]"), since requesting CIs for all the elements is computationally inefficient, since a correlation matrix is symmetric and has a fixed diagonal. This approach should work for all of the models.Does "D" here refer to a path coefficient, or to a variance component? If it's a variance component, then you should not expect to see a negative lower confidence limit, since the Cholesky parameterization prevents variance components from being negative (they are sums of squares of diagonal elements of the lower-triangular Cholesky factor). The fact that OpenMx's profile-likelihood CIs respect constraints on the parameters, and can be asymmetric when necessary, is among their advantages.
Log in or register to post comments
In reply to Thank you so much!!! by JuanJMV
Something I just noticed
meanMZandmeanDZ). One of those two things has to be fixed to identify the model. The "status RED" warning was probably a result of this unidentification. I suggest just fixing the threshold to zero (which would make constraining it equal across zygosity and twin order pointless).Log in or register to post comments
Thank you so much!
1-OK! but using this I have the correlations for MZ1 and MZ2 I mean this:
mxAlgebra 'corDZ'
$formula: cov2cor(covDZ)
$result:
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.00000000 0.67553492 0.35730803 0.34876643 0.41805389 0.35158189
[2,] 0.67553492 1.00000000 0.30300303 0.33465014 0.46008228 0.26706220
[3,] 0.35730803 0.30300303 1.00000000 0.25237671 0.30618880 0.38199174
[4,] 0.34876643 0.33465014 0.25237671 1.00000000 0.64877935 0.32860511
[5,] 0.41805389 0.46008228 0.30618880 0.64877935 1.00000000 0.28502398
[6,] 0.35158189 0.26706220 0.38199174 0.32860511 0.28502398 1.00000000
dimnames: NULL
How can I get the total value for MZ and DZ?
2-I am refering to this value:
A A A D D D E E E SA SA SA SD SD SD
VC 0.7930 0.6586 0.5146 0.2176 -0.0449 0.0040 0.5668 0.1360 -0.0552 0.5027 0.8786 1.1104 0.1380 -0.0600 0.0087
VC 0.6586 0.5482 0.4233 -0.0449 **0.0471** -0.1371 0.1360 0.2326 -0.0004 0.8786 0.6621 1.4809 -0.0600 0.0569 -0.4796
VC 0.5146 0.4233 0.3468 0.0040 -0.1371 0.4905 -0.0552 -0.0004 0.1627 1.1104 1.4809 0.3468 0.0087 -0.4796 0.4905
SE SE SE
VC 0.3593 0.1814 -0.1190
VC 0.1814 0.2809 -0.0013
VC -0.1190 -0.0013 0.1627
and I calculated its CI in order to know if it is significative am I wrong?
3-Also, I would like to know If I can fit a AE model for one variable and ACE for the other 2 variables. is it possible? I should fix C to 0 for one variable? how can I do this only for one variable?
4- for this:
Something I just noticed about the script in post #17: the dichotomous phenotype has free thresholds and free regression intercepts (in meanMZ and meanDZ). One of those two things has to be fixed to identify the model. The "status RED" warning was probably a result of this unidentification. I suggest just fixing the threshold to zero (which would make constraining it equal across zygosity and twin order pointless).
I had thought change this:
thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=TRUE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )
to :
thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )
But I am not sure if it would be correct.
Thank you so much again in adavance, I really appreciate your help.
Log in or register to post comments
In reply to Thank you so much! by JuanJMV
1-OK! but using this I have
Sorry, I don't understand the question?
That sounds good to me. The element you bolded is the dominance-genetic variance component for the second phenotype, so (given the parameterization you're using) its CI should not have a negative lower limit.
Yes, that's possible. The easiest way to do it would be to fix to zero the appropriate path coefficients in the common-environmental Cholesky factor (named "c" in your script), using
omxSetParameters(), in such a way that the phenotype has no common-environmental variance, nor common-environmental covariance with the other traits. I guess you could instead use MxConstraints to force the phenotype's variance and covariance elements to equal zero in the MxAlgebra named "C", but that is a clumsy approach. Be warned that if you drop the effect of C for only one phenotype, the MxAlgebra named "C" will no longer be a valid (positive-definite) covariance matrix, so the results of other MxAlgebras that use it may be rendered uninterpretable.That should work.
Log in or register to post comments
In reply to 1-OK! but using this I have by AdminRobK
Thanks!
1-Sorry I meant the correlations for MZ and DZ I mean the correlations that we would use in the Falconer's formula
2-So, if its CI should not have a negative lower limit how can I know if D (or other parameter) is significant?
3-I have made the change to this :
modelACE <- omxSetParameters(modelACE, labels=c("c_3_1","c_3_2","c_3_3"),free=FALSE, values = 0)
And It works !! :)
Thank you so much!!!
Log in or register to post comments
Hi again,
I have been working in the saturated model, but if I do this:
thinMZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinMZ", labels=c("tauMZ1","tauMZ2") )
thinDZ <- mxMatrix( type="Full", nrow=1, ncol=ntvo, free=FALSE, values=0, name="thinDZ", labels=c("tauDZ1","tauDZ2") )
then Taumz1, Taumz2, TauDZ1 and Taudz2 are not in the model therefore the saturated model and eqThresTwin and eqThresZyg models have the same degrees of freedom so what changes should I do in these models?
Also, I have been thinking in de ACE/AE model about what you told me “I guess you could instead use MxConstraints to force the phenotype's variance and covariance elements to equal zero in the MxAlgebra named "C", but that is a clumsy approach. Be warned that if you drop the effect of C for only one phenotype, the MxAlgebra named "C" will no longer be a valid (positive-definite) covariance matrix, so the results of other MxAlgebras that use it may be rendered uninterpretable.” But, I am not sure how I should do a ACE/AE model I thought
modelACE <- omxSetParameters(modelACE, labels=c("c_3_1","c_3_2","c_3_3"),free=FALSE, values = 0) but I'm not sure if it's both practically and theoretically correct
I made this decision because I think C in one trait is not significant 0.35(0,0.57)
But as I mentioned earlier. I'm not sure if I'm calculating well the intervals to know if it is significant or not. In order to do this I calculated the intervals for VC (VC[1,1], [2,2]…) I mean the unstandardized value. I don’t know if it is correct.
Thank you so much any insight you could give me about this or materials to learn more would be really helpful
Log in or register to post comments
In reply to Hi again, by JuanJMV
1-Sorry I meant the
?? I'm still confused as to which correlations you might want that aren't in the 'corMZ' and 'corDZ' MxAlgebras.
If the CI for a quantity contains zero, then that means that you cannot reject the null hypothesis that the quantity equals zero, i.e. it does not differ significantly from zero.
If you want to drop C from the third phenotype, then that looks reasonable to me.
Are the "eqThresTwin" and "eqThresZyg" models actually of interest? In any event, you could either (1) fix the intercepts of the ordinal phenotype to zero, and free the thresholds, or (2) leave the thresholds fixed, parameterize the ordinal phenotype's intercepts, in the saturated model, in terms of an MZ1 intercept, an MZ2 intercept, a DZ1 intercept, and a DZ2 intercept, and finally fit an analogous pair of additional models that reduce the four free intercept parameters to two (to equate across twin order) and reduce the two free parameters to one (to equate across zygosity).
Log in or register to post comments
Hi!
1- Sorry but my English is bad and my knowledge about Open Mx too ☹ If I use
corMZ <- mxAlgebra(cov2cor(covMZ), name="corMZ")
corDZ <- mxAlgebra(cov2cor(covDZ), name="corDZ")
Then I got this matrix:
dimnames: NULL
mxAlgebra 'expCovDZ'
$formula: rbind(cbind(V, cDZ), cbind(t(cDZ), V))
$result:
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1.60428543 0.76958186 0.48280016 0.59664363 0.41088273 0.40671628
[2,] 0.76958186 0.84383543 0.29644796 0.41088273 0.38962653 0.26024648
[3,] 0.48280016 0.29644796 1.00000003 0.40671628 0.26024648 0.52439786
[4,] 0.59664363 0.41088273 0.40671628 1.60428543 0.76958186 0.48280016
[5,] 0.41088273 0.38962653 0.26024648 0.76958186 0.84383543 0.29644796
[6,] 0.40671628 0.26024648 0.52439786 0.48280016 0.29644796 1.00000003
I would like to know the correlation for MZ and DZ ( rMZ and rDZ) I am not sure if I should get from the saturated model or ACE model and with or without covariables. For example, rMZ could be 0.40 and for DZ 0.20. It should be a heritability about 40%. I don’t know if I am expressing properly
2- So, If I have 0 in one interval the parameter is non-significant, right? I expected for example intervals such as (-0.425, 0.100) I mean not exactly the zero
3- If I do :modelACE <- omxSetParameters(modelACE, labels=c("c_3_1","c_3_2","c_3_3"),free=FALSE, values = 0)
Then the mode does not calculate rC between the other 2 parameters is it correct?
4- I am going to keep on working in the saturated model.
Thank you so much!
Log in or register to post comments
In reply to Hi! by JuanJMV
I would like to know the
I'd say, get it from the saturated model, with covariates.
Right. If you really want, you could do a likelihood-ratio test by comparing a model in which the quantity of interest is free versus a model in which it is fixed to the null value (which is usually zero). You can usually fix the quantity to the null value via judicious use of
omxSetParameters(); occasionally, you might need to use an MxConstraint.I think so? Could you post the output that doesn't look right to you?
Log in or register to post comments
In reply to I would like to know the by AdminRobK
Thank you so much! yes of
> fitACE$algebras
$iSD
mxAlgebra 'iSD'
$formula: solve(sqrt(I * V))
$result:
[,1] [,2] [,3]
[1,] 0.79311073 0.0000000 0.00000000
[2,] 0.00000000 1.0937125 0.00000000
[3,] 0.00000000 0.0000000 0.99999978
dimnames: NULL
$A
mxAlgebra 'A'
$formula: a %*% t(a)
$result:
[,1] [,2] [,3]
[1,] 0.81389313 0.46688500 0.54976880
[2,] 0.46688500 0.45165068 0.31551369
[3,] 0.54976880 0.31551369 0.77769329
dimnames: NULL
$C
mxAlgebra 'C'
$formula: c %*% t(c)
$result:
[,1] [,2] [,3]
[1,] 0.16839824 0.15043721 0
[2,] 0.15043721 0.13439190 0
[3,] 0.00000000 0.00000000 0
dimnames: NULL
$E
mxAlgebra 'E'
$formula: e %*% t(e)
$result:
[,1] [,2] [,3]
[1,] 0.607471487 0.141730178 -0.089903038
[2,] 0.141730178 0.249933037 -0.033470557
[3,] -0.089903038 -0.033470557 0.222307147
dimnames: NULL
$V
mxAlgebra 'V'
$formula: A + C + E
$result:
[,1] [,2] [,3]
[1,] 1.58976285 0.75905239 0.45986576
[2,] 0.75905239 0.83597562 0.28204313
[3,] 0.45986576 0.28204313 1.00000043
dimnames: NULL
$rA
mxAlgebra 'rA'
$formula: solve(sqrt(I * A)) %&% A
$result:
[,1] [,2] [,3]
[1,] 1.00000000 0.77006058 0.69102255
[2,] 0.77006058 1.00000000 0.53236893
[3,] 0.69102255 0.53236893 1.00000000
dimnames: NULL
$rC
mxAlgebra 'rC'
$formula: solve(sqrt(I * C)) %&% C
$result: (not yet computed) <0 x 0 matrix>
dimnames: NULL
$rE
mxAlgebra 'rE'
$formula: solve(sqrt(I * E)) %&% E
$result:
[,1] [,2] [,3]
[1,] 1.00000000 0.36373706 -0.24464408
[2,] 0.36373706 1.00000000 -0.14199544
[3,] -0.24464408 -0.14199544 1.00000000
dimnames: NULL
Thanks in advance
Log in or register to post comments
In reply to Thank you so much! yes of by JuanJMV
looks OK
sqrt(I * C)isn't invertible, so no result is returned for MxAlgebra 'rC'. Try something likemxEval(cov2cor(C[1:2,1:2]),fitACE,T)if you want to see the shared-environmental correlation for the first two phenotypes.Log in or register to post comments
In reply to looks OK by AdminRobK
It works perfectly as always!
Thank you so much!! I really appreciate it!
Log in or register to post comments