> Vars <- c("EXT3", "EXT5") > nvar <- length(Vars) # number of variables > tnvar <- nvar*2 # number of variables*max family size > nlower <-nvar*(nvar+1)/2 # number of free elements in a lower matrix nvar*nvar > > selVars <- paste(Vars, c(rep(1,nvar), rep(2,nvar)), sep="_") > > ## using both DZ-S and DZ-O. > mzData <- as.matrix(subset(twin_data, zygDich==1, selVars)) > dzData <- as.matrix(subset(twin_data, zygDich==0, selVars)) > > describe(mzData, skew=F) vars n mean sd min max range se EXT3_1 1 291 0.00 0.68 -1.47 2.16 3.63 0.04 EXT5_1 2 185 0.04 0.70 -1.26 1.93 3.19 0.05 EXT3_2 3 291 0.05 0.69 -1.75 1.97 3.72 0.04 EXT5_2 4 185 0.05 0.67 -1.39 2.05 3.45 0.05 > describe(dzData, skew=F) vars n mean sd min max range se EXT3_1 1 1084 -0.04 0.74 -1.67 2.38 4.05 0.02 EXT5_1 2 761 0.01 0.74 -1.53 2.32 3.85 0.03 EXT3_2 3 1080 0.07 0.73 -1.72 2.31 4.03 0.02 EXT5_2 4 757 0.05 0.77 -1.47 2.32 3.79 0.03 > dim(mzData) [1] 342 4 > dim(dzData) [1] 1257 4 > > # Generate Descriptive Statistics > colMeans(mzData,na.rm=TRUE) EXT3_1 EXT5_1 EXT3_2 EXT5_2 0.002611684 0.039275676 0.052938144 0.049670270 > colMeans(dzData,na.rm=TRUE) EXT3_1 EXT5_1 EXT3_2 EXT5_2 -0.04113100 0.00911958 0.07193519 0.05295112 > cov(mzData,use="complete") EXT3_1 EXT5_1 EXT3_2 EXT5_2 EXT3_1 0.4446349 0.19759141 0.18652768 0.1345051 EXT5_1 0.1975914 0.48636699 0.09653651 0.1969580 EXT3_2 0.1865277 0.09653651 0.41469496 0.1596745 EXT5_2 0.1345051 0.19695805 0.15967450 0.4287291 > cov(dzData,use="complete") EXT3_1 EXT5_1 EXT3_2 EXT5_2 EXT3_1 0.55458872 0.29376275 -0.09376566 -0.08926268 EXT5_1 0.29376275 0.54333388 -0.09395820 -0.05508327 EXT3_2 -0.09376566 -0.09395820 0.52783154 0.26437616 EXT5_2 -0.08926268 -0.05508327 0.26437616 0.57903044 > cor(mzData,use="complete") EXT3_1 EXT5_1 EXT3_2 EXT5_2 EXT3_1 1.0000000 0.4248975 0.4343873 0.3080671 EXT5_1 0.4248975 1.0000000 0.2149538 0.4313206 EXT3_2 0.4343873 0.2149538 1.0000000 0.3786864 EXT5_2 0.3080671 0.4313206 0.3786864 1.0000000 > cor(dzData,use="complete") EXT3_1 EXT5_1 EXT3_2 EXT5_2 EXT3_1 1.0000000 0.53515282 -0.1733048 -0.15751928 EXT5_1 0.5351528 1.00000000 -0.1754501 -0.09820545 EXT3_2 -0.1733048 -0.17545011 1.0000000 0.47821581 EXT5_2 -0.1575193 -0.09820545 0.4782158 1.00000000 > > nv <- length(Vars) # number of variables > #tnvar <- nvar*2 # number of variables*max family size > #nlower <-nvar*(nvar+1)/2 # number of free elements in a lower matrix nvar*nvar > > > # Set Starting Values > meVals <- c(0.5,0.5) # start value for means > meanLabs <- paste(Vars,"mean",sep="") # create labels for means > paVal <- .6 # start value for path coefficient > paLBo <- .0001 # start value for lower bounds > paValD <- vech(diag(paVal,nv,nv)) # start values for diagonal of covariance matrix > paLBoD <- diag(paLBo,nv,nv) # lower bounds for diagonal of covariance matrix > paLBoD[lower.tri(paLBoD)] <- -10 # lower bounds for below diagonal elements > paLBoD[upper.tri(paLBoD)] <- NA # lower bounds for above diagonal elements > > # Create Labels for Lower Triangular Matrices > aLabs <- paste("a",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") > cLabs <- paste("c",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") > eLabs <- paste("e",rev(nv+1-sequence(1:nv)),rep(1:nv,nv:1),sep="_") > > # Matrices declared to store a, c, and e Path Coefficients > pathA <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, + values=paValD, labels=aLabs, name="a" ) > pathC <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, + values=paValD, labels=cLabs, name="c" ) > pathE <- mxMatrix( type="Lower", nrow=nv, ncol=nv, free=TRUE, + values=paValD, labels=eLabs, name="e" ) > > # Matrices generated to hold A, C, and E computed Variance Components > covA <- mxAlgebra( expression=a %*% t(a), name="A" ) > covC <- mxAlgebra( expression=c %*% t(c), name="C" ) > covE <- mxAlgebra( expression=e %*% t(e), name="E" ) > > # Algebra to compute total variances and standard deviations (diagonal only) > covP <- mxAlgebra( expression=A+C+E, name="V" ) > matI <- mxMatrix( type="Iden", nrow=nv, ncol=nv, name="I") > invSD <- mxAlgebra( expression=solve(sqrt(I*V)), name="iSD") > > # Algebra to obtain solve(I-B) for sib interaction > bMat <- mxMatrix(type="Full", nrow=2*nv, ncol=2*nv, free=c(F,T,T,F), labels=c(NA,"b","b",NA), name="bMat") > iMat <- mxMatrix(type="Iden", nrow=2*nv, ncol=2*nv, name="iMat") > sibInt <- mxAlgebra(solve(iMat-bMat), name="sibInt") > > # Algebra for expected Mean and Variance/Covariance Matrices in MZ & DZ twins > meanG <- mxMatrix( type="Full", nrow=1, ncol=nv, free=TRUE, + values=meVals, labels=meanLabs, name="Mean" ) > meanT <- mxAlgebra( expression= cbind(Mean,Mean), name="expMean" ) > covMZ <- mxAlgebra( expression= sibInt %&% ( + rbind( cbind(A+C+E , A+C), + cbind(A+C , A+C+E))), name="expCovMZ" ) > covDZ <- mxAlgebra( expression= sibInt %&% ( + rbind( cbind(A+C+E , 0.5%x%A+C), + cbind(0.5%x%A+C , A+C+E))), name="expCovDZ" ) > > # Standeraized parameters and CI (getting CIs for squared path estimates) > sta <- mxAlgebra((iSD %*% a)*(iSD %*% a), name="sta") > stc <- mxAlgebra((iSD %*% c)*(iSD %*% c), name="stc") > ste <- mxAlgebra((iSD %*% e)*(iSD %*% e), name="ste") > StandardizedA <- mxAlgebra(A/V, name="StandardizedA") > StandardizedC <- mxAlgebra(C/V, name="StandardizedC") > StandardizedE <- mxAlgebra(E/V, name="StandardizedE") > > CI <- mxCI(c('sta', 'stc', 'ste')) > > # Data objects for Multiple Groups > dataMZ <- mxData( observed=mzData, type="raw" ) > dataDZ <- mxData( observed=dzData, type="raw" ) > > # Objective objects for Multiple Groups > objMZ <- mxExpectationNormal( covariance="expCovMZ", means="expMean", + dimnames=selVars ) > objDZ <- mxExpectationNormal( covariance="expCovDZ", means="expMean", + dimnames=selVars ) > funML <- mxFitFunctionML() > > # Combine Groups > pars <- list( pathA, pathC, pathE, covA, covC, covE, covP, + matI, invSD, meanG, meanT, sta, stc, ste, + StandardizedA, StandardizedC, StandardizedE, bMat, iMat, sibInt) > modelMZ <- mxModel( pars, covMZ, dataMZ, objMZ, funML, name="MZ" ) > modelDZ <- mxModel( pars, covDZ, dataDZ, objDZ, funML, name="DZ" ) > fitML <- mxFitFunctionMultigroup(c("MZ.fitfunction","DZ.fitfunction") ) > CholAceModel <- mxModel( "CholACE", pars, modelMZ, modelDZ, fitML, CI) > > # Run Cholesky Decomposition ACE model > CholAceFit <- mxRun(CholAceModel, intervals = TRUE) Running CholACE with 12 parameters [8/2] MxComputeGradientDescent(SLSQP) evaluations 7774 fit 2.16279e-10 change -1.462e-08 > CholAceSumm <- summary(CholAceFit, verbose=T) > CholAceSumm Summary of CholACE data: $MZ.data EXT3_1 EXT5_1 EXT3_2 EXT5_2 Min. :-1.47100 Min. :-1.26300 Min. :-1.74500 Min. :-1.39400 1st Qu.:-0.51550 1st Qu.:-0.51900 1st Qu.:-0.44650 1st Qu.:-0.46000 Median : 0.02500 Median :-0.07200 Median : 0.05100 Median : 0.05300 Mean : 0.00261 Mean : 0.03928 Mean : 0.05294 Mean : 0.04967 3rd Qu.: 0.46900 3rd Qu.: 0.49500 3rd Qu.: 0.54950 3rd Qu.: 0.50600 Max. : 2.15800 Max. : 1.92900 Max. : 1.97500 Max. : 2.05200 NA's :51 NA's :157 NA's :51 NA's :157 $DZ.data EXT3_1 EXT5_1 EXT3_2 EXT5_2 Min. :-1.66500 Min. :-1.5320 Min. :-1.71600 Min. :-1.471 1st Qu.:-0.65325 1st Qu.:-0.5950 1st Qu.:-0.48200 1st Qu.:-0.580 Median :-0.09400 Median :-0.0680 Median : 0.03850 Median : 0.025 Mean :-0.04113 Mean : 0.0091 Mean : 0.07194 Mean : 0.053 3rd Qu.: 0.47750 3rd Qu.: 0.5010 3rd Qu.: 0.56400 3rd Qu.: 0.575 Max. : 2.38100 Max. : 2.3170 Max. : 2.31000 Max. : 2.323 NA's :173 NA's :496 NA's :177 NA's :500 free parameters: name matrix row col Estimate Std.Error A lbound ubound 1 a_1_1 a 1 1 2.966120e-01 0.049463510 2 a_2_1 a 2 1 9.537415e-02 0.086803877 3 a_2_2 a 2 2 2.927853e-01 0.043619444 4 c_1_1 c 1 1 1.524232e-07 0.056435773 5 c_2_1 c 2 1 3.196537e-07 0.089604287 6 c_2_2 c 2 2 -1.189289e-07 0.088790143 7 e_1_1 e 1 1 6.795622e-01 0.020517713 8 e_2_1 e 2 1 3.809450e-01 0.032770733 9 e_2_2 e 2 2 5.765765e-01 0.020993705 10 EXT3mean Mean 1 1 1.999969e-02 0.014146728 11 EXT5mean Mean 1 2 3.097525e-02 0.016904570 12 b bMat 2 1 -2.584269e-02 0.009168575 confidence intervals: lbound estimate ubound note CholACE.sta[1,1] 5.921499e-02 1.600243e-01 0.26161707 CholACE.sta[2,1] 5.287978e-34 1.589197e-02 0.10463054 CholACE.sta[1,2] 0.000000e+00 0.000000e+00 0.00000000 !!! CholACE.sta[2,2] 5.048442e-02 1.497666e-01 0.23156572 CholACE.stc[1,1] 8.248534e-30 4.225814e-14 0.02157557 CholACE.stc[2,1] 2.111290e-34 1.785155e-13 0.04828541 CholACE.stc[1,2] 0.000000e+00 0.000000e+00 0.00000000 !!! CholACE.stc[2,2] 9.528482e-32 2.471103e-14 0.04813732 CholACE.ste[1,1] 7.383833e-01 8.399757e-01 0.94078648 CholACE.ste[2,1] 1.731881e-01 2.535366e-01 0.32707079 CholACE.ste[1,2] 0.000000e+00 0.000000e+00 0.00000000 !!! CholACE.ste[2,2] 5.003981e-01 5.808048e-01 0.66696672 CI details: parameter side value fit diagnostic statusCode method 1 CholACE.sta[1,1] lower 5.921499e-02 9791.776 success OK neale-miller-1997 2 CholACE.sta[1,1] upper 2.616171e-01 9791.775 success OK neale-miller-1997 3 CholACE.sta[2,1] lower 5.287978e-34 9791.771 success OK neale-miller-1997 4 CholACE.sta[2,1] upper 1.046305e-01 9791.776 success OK neale-miller-1997 5 CholACE.sta[1,2] lower 0.000000e+00 9791.775 success OK neale-miller-1997 6 CholACE.sta[1,2] upper 0.000000e+00 9791.775 success OK neale-miller-1997 7 CholACE.sta[2,2] lower 5.048442e-02 9791.776 success OK neale-miller-1997 8 CholACE.sta[2,2] upper 2.315657e-01 9791.775 success OK neale-miller-1997 9 CholACE.stc[1,1] lower 8.248534e-30 9791.773 success OK neale-miller-1997 10 CholACE.stc[1,1] upper 2.157557e-02 9791.776 success OK neale-miller-1997 11 CholACE.stc[2,1] lower 2.111290e-34 9791.768 success OK neale-miller-1997 12 CholACE.stc[2,1] upper 4.828541e-02 9791.776 success OK neale-miller-1997 13 CholACE.stc[1,2] lower 0.000000e+00 9791.775 success OK neale-miller-1997 14 CholACE.stc[1,2] upper 0.000000e+00 9791.775 success OK neale-miller-1997 15 CholACE.stc[2,2] lower 9.528482e-32 9791.746 success OK neale-miller-1997 16 CholACE.stc[2,2] upper 4.813732e-02 9791.776 success OK neale-miller-1997 17 CholACE.ste[1,1] lower 7.383833e-01 9791.775 success OK neale-miller-1997 18 CholACE.ste[1,1] upper 9.407865e-01 9791.776 success OK neale-miller-1997 19 CholACE.ste[2,1] lower 1.731881e-01 9791.776 success OK neale-miller-1997 20 CholACE.ste[2,1] upper 3.270708e-01 9791.775 success OK neale-miller-1997 21 CholACE.ste[1,2] lower 0.000000e+00 9791.775 success OK neale-miller-1997 22 CholACE.ste[1,2] upper 0.000000e+00 9791.775 success OK neale-miller-1997 23 CholACE.ste[2,2] lower 5.003981e-01 9791.776 success OK neale-miller-1997 24 CholACE.ste[2,2] upper 6.669667e-01 9791.775 success OK neale-miller-1997 a_1_1 a_2_1 a_2_2 c_1_1 c_2_1 c_2_2 e_1_1 1 -0.1791799 4.482190e-02 0.2628680 -3.144940e-05 5.638027e-05 -6.870275e-05 0.7141981 2 0.3831676 1.774438e-01 0.2993529 -1.656143e-07 -6.236603e-06 5.110213e-09 0.6437193 3 0.2877474 1.757501e-17 0.2929738 1.545152e-07 3.209140e-07 -1.171140e-07 0.6807944 4 0.3350091 2.475169e-01 0.2941886 -2.335216e-04 2.211207e-04 -3.463375e-05 0.6679825 5 0.2952152 8.147860e-02 0.3043934 -1.147170e-02 -4.031145e-02 3.576713e-02 0.6786527 6 0.3172594 8.190361e-02 0.2896685 -5.436502e-03 -3.555027e-02 -5.198509e-03 0.6619769 7 0.2244064 -8.786803e-02 -0.1684988 9.883656e-04 -4.543656e-02 -9.976699e-03 0.7046593 8 0.3004123 1.241346e-01 0.3660612 6.899049e-08 -7.754691e-08 -2.661873e-07 0.6782227 9 0.3168946 7.312977e-02 0.3234478 2.143243e-15 2.551292e-02 4.467853e-03 0.6756210 10 0.2627329 1.154554e-01 0.2850676 1.089923e-01 -1.534000e-02 -2.291077e-09 0.6853344 11 0.2838279 1.061311e-01 0.3237822 -1.160309e-02 1.111701e-17 4.682493e-02 0.6924718 12 0.2967020 9.602793e-02 0.2026246 -5.949076e-03 1.664170e-01 6.655811e-06 0.6794090 13 0.2762674 1.454386e-01 0.3058666 4.007614e-02 3.609052e-02 2.705917e-02 0.6954754 14 0.3090294 9.903026e-02 0.3060152 -5.724688e-03 3.864484e-02 -3.486776e-02 0.6800429 15 0.3186726 1.244519e-01 0.2869771 -2.426222e-02 -1.476362e-02 -2.366910e-16 0.6840166 16 0.2966497 9.032479e-02 0.2048518 -1.778062e-05 1.009606e-04 -1.661784e-01 0.6795304 17 0.3831656 1.774086e-01 0.2993341 -5.495728e-08 6.342705e-07 -2.230075e-09 0.6437167 18 -0.1791830 4.487373e-02 0.2628034 -2.672434e-05 2.953762e-05 1.678011e-05 0.7142200 19 0.3375149 2.319317e-01 0.2975769 3.078347e-06 1.239982e-06 6.946254e-10 0.6632122 20 0.2478249 -8.767609e-02 -0.2005438 -1.902694e-04 -2.933143e-03 -2.222062e-03 0.6999265 21 0.3149549 1.047724e-01 0.2699546 1.109786e-03 5.293696e-04 -1.038720e-02 0.6660256 22 0.3343451 6.918420e-02 0.2859929 1.045955e-02 -9.728542e-03 1.683581e-02 0.6536515 23 0.3036929 1.182761e-01 0.3603478 -1.336972e-05 1.861870e-06 2.238908e-07 0.6803705 24 0.2889077 7.935302e-02 -0.2064947 1.172570e-05 -3.214901e-05 6.465120e-10 0.6786534 e_2_1 e_2_2 EXT3mean EXT5mean b 1 0.4041408 0.5754971 0.019219579 0.03034892 -0.01795438 2 0.3545583 0.5779482 0.020876467 0.03170443 -0.03366711 3 0.4053475 0.5779094 0.019908631 0.03120714 -0.02525442 4 0.3254629 0.5760164 0.020631019 0.03108662 -0.03338358 5 0.3742858 0.5716828 0.004894757 0.00706292 -0.02842985 6 0.3702606 0.5636224 0.028834157 0.04367339 -0.01978428 7 0.4249016 0.5861486 0.021920700 0.03345703 -0.01832472 8 0.3670117 0.5427386 0.020213303 0.03110817 -0.02770755 9 0.3988686 0.5533152 0.011161991 0.01244337 -0.02323148 10 0.3772107 0.5788541 0.019983059 0.03092931 -0.02588472 11 0.3933835 0.5588606 0.010580137 0.01084321 -0.02834827 12 0.3822836 0.5911409 0.019919472 0.03103651 -0.02594087 13 0.3716117 0.5844778 0.027490138 0.02700470 -0.02826214 14 0.4114676 0.5687822 0.015638145 0.02361725 -0.02828929 15 0.3858504 0.5839576 0.016446320 0.05080298 -0.03331515 16 0.3834448 0.5906851 0.019909286 0.03107869 -0.02590674 17 0.3545600 0.5779713 0.020894397 0.03170790 -0.03366154 18 0.4041475 0.5755160 0.019248252 0.03032590 -0.01795693 19 0.3165707 0.5797400 0.020449026 0.03110263 -0.03192015 20 0.4314573 0.5788657 0.015152692 0.03397242 -0.01936569 21 0.3614071 0.5915371 0.013303120 0.04377767 -0.03188161 22 0.3748656 0.5779067 0.025937276 0.03746548 -0.02277717 23 0.3845620 0.5405483 0.020406589 0.03119117 -0.02883172 24 0.3725405 0.6131510 0.019570003 0.03075231 -0.02238281 Model Statistics: | Parameters | Degrees of Freedom | Fit (-2lnL units) Model: 12 4622 9787.934 Saturated: NA NA NA Independence: NA NA NA Number of observations/statistics: 1599/4634 condition number of the information matrix: 248.1277 maximum absolute gradient: 0.003136516 ( e_1_1 ) chi-square: χ² ( df=NA ) = NA, p = NA Information Criteria: | df Penalty | Parameters Penalty | Sample-Size Adjusted AIC: 543.934 9811.934 9812.131 BIC: -24309.178 9876.460 9838.338 CFI: NA TLI: 1 (also known as NNFI) RMSEA: 0 [95% CI (NA, NA)] Prob(RMSEA <= 0.05): NA To get additional fit indices, see help(mxRefModels) timestamp: 2021-06-08 15:30:57 frontend time: 0.1499999 secs backend time: 0.7500012 secs independent submodels time: 0 secs cpu time: 0.900001 secs Wall clock time: 0.900001 secs OpenMx version number: 2.18.1 Need help? See help(mxSummary) >