You are here

Trivariate Model

8 posts / 0 new
Last post
JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Trivariate Model

Hi Everyone,

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

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
I have managed to make it

I have managed to make it work but I'm not sure because I have to put 2 thresholds otherwise it does not work, but my ordinal variable is dichotomous. is it correct?

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.

AdminNeale's picture
Offline
Joined: 03/01/2013 - 14:09
Only one threshold needed, check other issues

Hi Juan

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

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
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.

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
ordinal variable

The covariates look OK to me.

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)
JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Thank you! I have made the

Thank you! I have made the changes but I have a warning message in this line:

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!!

JuanJMV's picture
Offline
Joined: 07/20/2016 - 13:13
Sorry I forgot to say that I

Sorry I forgot to say that I use my own data

Thanks

AdminRobK's picture
Offline
Joined: 01/24/2014 - 12:15
warning message
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")

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