Hi,
I am new to bivariate models and I have tried to adapt the twoACEj.R script by Hermine Maes from an ordinal to a binary variable in the below script. The script runs without errors, but my first question is whether my adaptations are correct?
# Select Binary Variables nth <- 1 # number of thresholds varso <- c('ht') # list of binary variables names nvo <- 1 # number of binary variables ntvo <- nvo*2 # number of total binary variables ordVars <- paste(varso,c(rep(1,nvo),rep(2,nvo)),sep="") # Select Continuous Variables varsc <- c('bmi') # list of continuous variables names nvc <- 1 # number of continuous variables ntvc <- nvc*2 # number of total continuous variables conVars <- paste(varsc,c(rep(1,nvc),rep(2,nvc)),sep="") # Select Variables for Analysis vars <- c('ht','bmi') # list of variables names nv <- nvo+nvc # number of variables ntv <- nv*2 # number of total variables oc <- c(1,0) # 0 for continuous, 1 for ordinal variables selVars <- paste(vars,c(rep(1,nv),rep(2,nv)),sep="") # Select Data for Analysis mzData <- subset(twinData, zyg==1, selVars) dzData <- subset(twinData, zyg==2, selVars) mzDataF <- cbind(mxFactor( x=mzData[,ordVars], levels=c(0:nth)), mzData[,conVars]) dzDataF <- cbind(mxFactor( x=dzData[,ordVars], levels=c(0:nth)), dzData[,conVars]) # Generate Descriptive Statistics apply(mzData,2,myMean) apply(dzData,2,myMean) #?apply(mzData,2,myCor) #?apply(dzData,2,myCor) # Set Starting Values frMV <- c(FALSE, TRUE) # 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(0,22) # start value for means svPa <- .4 # start value for path coefficient svPaD <- vech(diag(svPa,nv,nv)) # start values for diagonal of covariance matrix svPe <- .8 # 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.5 # 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", c('hs','bmi'),sep="_") labTh <- paste(paste("t",1:nth,sep=""),rep(varso,each=nth),sep="_") # ------------------------------------------------------------------------------ # PREPARE MODEL # ACE Model # Create Algebra for expected Mean Matrices meanG <- mxMatrix( type="Full", nrow=1, ncol=ntv, free=frMV, values=svMe, labels=labMe, name="meanG" ) 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="meanG", dimnames=selVars, thresholds="threG", threshnames=ordVars ) expDZ <- mxExpectationNormal( covariance="expCovDZ", means="meanG", dimnames=selVars, thresholds="threG", threshnames=ordVars ) funML <- mxFitFunctionML() # Create Model Objects for Multiple Groups pars <- list(meanG, thinG, inc, threG, matI, invSD, matUnv, matOc, pathA, pathC, pathE, covA, covC, covE, covP, corA, corC, corE) modelMZ <- mxModel( name="MZ", pars, covMZ, expCovMZ, dataMZ, expMZ, funML ) modelDZ <- mxModel( name="DZ", pars, 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,7]","VC[1,9]", "VC[1,11]", "VC[2,8]", "VC[2,10]", "VC[2,12]" )) # Build Model with Confidence Intervals modelACE <- mxModel( "twoACEj", pars, var1, modelMZ, modelDZ, multi, estVC,ciACE)
Subsequently, based on the below output I choose the ACE model over any submodel:
base comparison ep minus2LL df AIC diffLL diffdf p
1 twoACEj 11 32712.800 11178 10356.800 NA NA NA
2 twoACEj twoAEj 8 32722.252 11181 10360.252 9.4523951 3 2.3843227e-02
3 twoACEj twoCEj 8 32940.398 11181 10578.398 227.5982305 3 4.5718550e-49
4 twoACEj twoEj 5 34257.189 11184 11889.189 1544.3890445 6 0.0000000e+00
However, the lower bound of the CI of SA is 0.000, is this correct and does this make sense?
> fitEsts(fitACE)
mean_bmi t1_hs a_1_1 a_2_1 a_2_2 c_1_1 c_2_1 c_2_2 e_1_1 e_2_1 e_2_2
23.1373 2.2241 0.8065 -0.0811 2.8604 0.3302 -1.1363 0.0023 0.4904 -0.0410 1.9097
A A C C E E SA SA SC SC SE SE
VC 0.6505 -0.0654 0.1091 -0.3753 0.2405 -0.0201 0.6505 0.1419 0.1091 0.8144 0.2405 0.0436
VC -0.0654 8.1886 -0.3753 1.2913 -0.0201 3.6485 0.1419 0.6237 0.8144 0.0984 0.0436 0.2779
lbound estimate ubound
twoACEj.VC[1,7] 0.0000 0.6505 0.8507
twoACEj.VC[1,9] 0.0059 0.1091 0.6982
twoACEj.VC[1,11] 0.1067 0.2405 0.4640
twoACEj.VC[2,8] 0.5352 0.6237 0.7075
twoACEj.VC[2,10] 0.0196 0.0984 0.1809
twoACEj.VC[2,12] 0.2580 0.2779 0.2992
Moreover, the below CI details seem to "report a status code "nonzero gradient/red for all CIs. I am not sure how to interpret this?
CI details:
parameter value side fit mean_bmi t1_hs a_1_1 a_2_1 a_2_2 c_1_1
1 twoACEj.VC[1,7] 2.4850239e-05 lower 32716.641 23.140124 2.2065754 0.004985004 0.108697628 2.8455448 0.811215832
2 twoACEj.VC[1,7] 8.5065598e-01 upper 32716.650 23.128152 2.2105490 0.922310133 -0.299920357 2.8246636 0.157758943
3 twoACEj.VC[1,9] 5.8810741e-03 lower 32716.644 23.120816 2.2181279 0.882108184 -0.458473930 2.8182535 0.076688161
4 twoACEj.VC[1,9] 6.9816923e-01 upper 32716.673 23.138658 2.2183166 0.018547226 -2.493000405 1.4130896 0.835565224
5 twoACEj.VC[1,11] 1.0665508e-01 lower 32716.648 23.135947 2.2007855 0.867513790 -0.120282731 2.8472254 0.375186283
6 twoACEj.VC[1,11] 4.6401224e-01 upper 32716.658 23.136136 2.2373873 0.573177370 -0.078829212 2.8644263 0.455475777
7 twoACEj.VC[2,8] 5.3520770e-01 lower 32716.646 23.138107 2.2231490 0.826606517 -0.066884929 2.6534599 0.281284034
8 twoACEj.VC[2,8] 7.0746664e-01 upper 32716.645 23.131552 2.2216208 0.739852906 -0.257968239 3.0366597 0.479652614
9 twoACEj.VC[2,10] 1.9566469e-02 lower 32716.645 23.130503 2.2220530 0.731951748 -0.280816975 3.0194490 0.488947391
10 twoACEj.VC[2,10] 1.8086435e-01 upper 32716.647 23.138985 2.2244832 0.799906860 -0.053095311 2.6695519 0.313275027
11 twoACEj.VC[2,12] 2.5798502e-01 lower 32716.643 23.138385 2.2240622 0.794925301 -0.082005655 2.9456948 0.347706223
12 twoACEj.VC[2,12] 2.9923245e-01 upper 32716.643 23.136116 2.2281715 0.782260737 -0.081690070 2.7667635 0.358493602
c_2_1 c_2_2 e_1_1 e_2_1 e_2_2 method diagnostic statusCode
1 -0.54132925 1.0148719e+00 0.58472646 -0.0458843144 1.9226946 neale-miller-1997 success nonzero gradient/red
2 -1.18976326 4.9083101e-04 0.35278341 0.0193929090 1.9121737 neale-miller-1997 success nonzero gradient/red
3 -0.85016056 7.6804313e-01 0.46476239 0.0483915125 1.9088932 neale-miller-1997 success nonzero gradient/red
4 -0.47162022 1.0220297e+00 0.54907811 -0.0425895786 1.9092756 neale-miller-1997 success nonzero gradient/red
5 -0.93157210 6.9198941e-01 0.32658089 -0.0382025425 1.9088138 neale-miller-1997 success nonzero gradient/red
6 -0.84669292 7.4380220e-01 0.68118617 -0.0357496784 1.9094127 neale-miller-1997 success nonzero gradient/red
7 -1.37263512 6.8850457e-01 0.48744340 -0.0423087655 1.9386693 neale-miller-1997 success nonzero gradient/red
8 -0.54445165 -7.7042247e-05 0.47175324 -0.0136392463 1.8825139 neale-miller-1997 success nonzero gradient/red
9 -0.50536961 -2.7393269e-05 0.47452828 -0.0095315172 1.8977468 neale-miller-1997 success nonzero gradient/red
10 -1.27158725 8.8090563e-01 0.51186716 -0.0442595747 1.9251761 neale-miller-1997 success nonzero gradient/red
11 -1.08633458 2.5297820e-01 0.49718742 -0.0415448035 1.8574297 neale-miller-1997 success nonzero gradient/red
12 -1.04397757 5.3634622e-01 0.50947280 -0.0396417596 1.9642382 neale-miller-1997 success nonzero gradient/red
I appreciate your help!