> require(OpenMx) > mydatafile=read.csv("mydataanthro.csv", header=T, na.strings=".") > summary(mydatafile) Zygosity Gender1 Gender2 Age1 Age2 Min. :0.0000 Min. :0.0000 Min. :0.0000 Min. :19.00 Min. :19.00 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:24.00 1st Qu.:24.00 Median :0.0000 Median :1.0000 Median :1.0000 Median :33.00 Median :33.00 Mean :0.3478 Mean :0.5652 Mean :0.5652 Mean :35.52 Mean :35.52 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:1.0000 3rd Qu.:45.00 3rd Qu.:45.00 Max. :1.0000 Max. :1.0000 Max. :1.0000 Max. :65.00 Max. :65.00 height1 height2 waist1 waist2 hip1 Min. :1.531 Min. :1.435 Min. : 62.00 Min. : 64.00 Min. : 77.50 1st Qu.:1.631 1st Qu.:1.640 1st Qu.: 74.50 1st Qu.: 74.50 1st Qu.: 92.50 Median :1.679 Median :1.686 Median : 79.00 Median : 80.00 Median : 98.50 Mean :1.698 Mean :1.700 Mean : 80.16 Mean : 79.72 Mean : 98.71 3rd Qu.:1.758 3rd Qu.:1.758 3rd Qu.: 86.00 3rd Qu.: 85.00 3rd Qu.:103.50 Max. :1.925 Max. :1.899 Max. :110.00 Max. :101.00 Max. :121.00 hip2 whr1 whr2 bodyfat1 bodyfat2 Min. : 77.50 Min. :0.6735 Min. :0.6733 Min. : 6.00 Min. : 7.90 1st Qu.: 93.00 1st Qu.:0.7556 1st Qu.:0.7632 1st Qu.:17.70 1st Qu.:18.60 Median : 97.00 Median :0.8135 Median :0.7960 Median :25.00 Median :26.10 Mean : 97.38 Mean :0.8127 Mean :0.8208 Mean :26.49 Mean :25.78 3rd Qu.:102.00 3rd Qu.:0.8650 3rd Qu.:0.8854 3rd Qu.:35.30 3rd Qu.:31.30 Max. :115.00 Max. :1.0138 Max. :1.0000 Max. :52.00 Max. :45.80 bodymass1 bodymass2 bmi1 bmi2 sys1 Min. : 50.52 Min. :48.84 Min. :17.90 Min. :18.13 Min. : 98.0 1st Qu.: 62.59 1st Qu.:61.30 1st Qu.:22.10 1st Qu.:22.00 1st Qu.:113.2 Median : 69.27 Median :70.01 Median :23.79 Median :23.95 Median :122.0 Mean : 70.55 Mean :69.88 Mean :24.46 Mean :24.09 Mean :121.7 3rd Qu.: 78.16 3rd Qu.:75.53 3rd Qu.:26.72 3rd Qu.:26.09 3rd Qu.:130.0 Max. :101.27 Max. :98.32 Max. :35.61 Max. :30.69 Max. :146.0 sys2 Min. : 93.67 1st Qu.:111.60 Median :120.60 Mean :121.24 3rd Qu.:127.00 Max. :164.00 > mycols=colnames(mydatafile) > VarsA <- 'bmi1' > VarsB <- 'bmi2' > Strtpaths <- 0.6 > StrtMean <- 20 > selVars=c(VarsA,VarsB) > print(selVars) [1] "bmi1" "bmi2" > cat("\n") > mzData<-data.frame(subset(mydatafile, Zygosity==0, c(VarsA, VarsB, 'Age1','Age2', 'Gender1', 'Gender2'))) > dzData<-data.frame(subset(mydatafile, Zygosity==1, c(VarsA, VarsB, 'Age1','Age2', 'Gender1', 'Gender2'))) > datasetname=c(mydatafile,mycols) > datasetname $Zygosity [1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 [45] 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 $Gender1 [1] 1 0 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 1 0 0 0 [45] 1 0 0 1 1 1 1 1 0 1 1 0 0 1 0 0 0 0 0 0 1 1 1 1 1 $Gender2 [1] 1 0 1 1 1 1 1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 1 1 1 0 0 0 0 1 0 0 0 [45] 1 0 0 1 1 1 1 1 0 1 1 0 0 1 0 0 0 0 0 0 1 1 1 1 1 $Age1 [1] 37 23 20 49 34 19 60 45 50 60 33 44 30 42 31 22 46 46 46 24 60 47 21 28 28 21 24 45 44 [30] 35 30 25 25 50 41 21 36 33 42 28 45 52 38 19 43 33 34 19 59 19 44 33 20 46 53 65 22 40 [59] 22 31 21 19 26 23 25 37 20 33 65 $Age2 [1] 37 23 20 49 34 19 60 45 50 60 33 44 30 42 31 22 46 46 46 24 60 47 21 28 28 21 24 45 44 [30] 35 30 25 25 50 41 21 36 33 42 28 45 52 38 19 43 33 34 19 59 19 44 33 20 46 53 65 22 40 [59] 22 31 21 19 26 23 25 37 20 33 65 $height1 [1] 1.685 1.784 1.678 1.680 1.700 1.723 1.665 1.620 1.810 1.531 1.664 1.608 1.665 1.672 [15] 1.677 1.641 1.624 1.651 1.647 1.620 1.575 1.716 1.705 1.820 1.631 1.850 1.859 1.734 [29] 1.720 1.670 1.755 1.695 1.665 1.632 1.603 1.566 1.820 1.758 1.817 1.637 1.614 1.820 [43] 1.781 1.764 1.638 1.723 1.834 1.738 1.571 1.629 1.672 1.679 1.780 1.600 1.575 1.725 [57] 1.802 1.631 1.827 1.823 1.624 1.925 1.728 1.799 1.692 1.746 1.623 1.600 1.598 $height2 [1] 1.688 1.735 1.687 1.665 1.701 1.758 1.652 1.622 1.835 1.521 1.668 1.614 1.664 1.667 [15] 1.692 1.644 1.619 1.671 1.638 1.615 1.615 1.735 1.695 1.812 1.627 1.840 1.843 1.773 [29] 1.708 1.654 1.743 1.704 1.673 1.600 1.620 1.592 1.839 1.768 1.833 1.686 1.632 1.834 [43] 1.793 1.766 1.612 1.740 1.842 1.655 1.599 1.686 1.639 1.693 1.838 1.645 1.625 1.669 [57] 1.899 1.725 1.756 1.763 1.732 1.892 1.789 1.752 1.640 1.657 1.645 1.645 1.435 $waist1 [1] 70.5 77.0 75.5 63.0 74.5 65.6 78.0 90.0 85.5 87.5 66.5 85.0 63.5 71.0 [15] 70.5 62.5 71.0 77.0 78.0 99.0 84.5 77.5 66.0 87.0 75.0 89.5 81.5 110.0 [29] 91.5 79.0 75.5 82.0 89.9 104.0 77.0 66.0 80.0 79.0 96.0 74.0 74.0 82.0 [43] 81.0 76.5 93.0 81.0 73.0 87.0 77.5 85.5 89.5 90.0 88.0 75.0 70.5 86.0 [57] 80.0 94.0 78.5 86.5 77.0 84.0 84.5 83.0 70.0 85.0 75.0 62.0 86.0 $waist2 [1] 75.0 73.0 75.0 64.0 75.5 74.5 81.0 92.0 86.5 87.0 68.0 78.0 83.0 73.5 [15] 70.0 65.0 73.0 75.5 80.0 79.5 87.5 76.5 67.5 92.0 77.5 91.0 80.5 101.0 [29] 89.0 78.0 80.0 85.0 89.0 89.5 85.0 66.5 82.0 82.0 100.0 74.0 78.5 88.0 [43] 86.0 79.0 81.0 83.0 80.0 68.5 80.0 78.0 82.0 79.0 87.5 71.0 74.5 91.5 [57] 80.5 72.5 85.0 80.0 78.5 85.5 86.0 82.5 74.0 77.5 73.5 64.0 71.5 $hip1 [1] 103.5 94.0 99.5 87.5 96.5 96.5 100.5 100.0 98.5 100.5 98.5 112.5 92.0 88.0 [15] 95.0 92.5 94.0 98.5 105.5 121.0 110.0 98.5 98.0 106.0 77.5 108.5 103.0 108.5 [29] 96.5 90.0 88.0 91.0 101.5 109.5 94.0 89.0 101.5 94.0 116.0 86.0 97.0 97.0 [43] 92.0 91.0 109.0 86.0 93.0 108.0 101.0 116.5 105.5 108.0 110.0 101.0 94.5 91.0 [57] 91.0 103.5 96.5 100.0 87.0 104.0 91.5 101.0 93.5 109.0 102.0 88.0 102.0 $hip2 [1] 103.5 91.5 101.5 85.5 96.5 106.0 101.0 97.0 97.5 99.0 101.0 103.5 105.0 95.5 [15] 94.0 89.0 95.5 97.0 104.0 113.0 115.0 96.5 100.0 95.0 77.5 110.0 104.5 104.5 [29] 98.0 88.0 84.0 86.0 100.0 102.0 96.0 85.5 104.0 93.0 110.0 88.0 94.5 99.0 [43] 95.0 91.0 102.0 89.0 93.0 96.0 100.5 108.0 105.5 102.0 100.0 95.0 96.5 97.5 [57] 92.0 95.0 103.5 96.5 95.5 92.0 97.5 89.0 97.0 100.0 98.0 92.5 92.0 $whr1 [1] 0.6812 0.8191 0.7588 0.7200 0.7720 0.6798 0.7761 0.9000 0.8680 0.8706 0.6751 0.7556 [13] 0.6902 0.8068 0.7421 0.6757 0.7553 0.7817 0.7393 0.8182 0.7682 0.7868 0.6735 0.8208 [25] 0.9677 0.8249 0.7913 1.0138 0.9482 0.8778 0.8580 0.9011 0.8857 0.9498 0.8191 0.7416 [37] 0.7882 0.8404 0.8276 0.8605 0.7629 0.8454 0.8804 0.8407 0.8532 0.9419 0.7849 0.8056 [49] 0.7673 0.7339 0.8483 0.8333 0.8000 0.7426 0.7460 0.9451 0.8791 0.9082 0.8135 0.8650 [61] 0.8851 0.8077 0.9235 0.8218 0.7487 0.7798 0.7353 0.7045 0.8431 $whr2 [1] 0.7246 0.7978 0.7389 0.7485 0.7824 0.7028 0.8020 0.9485 0.8872 0.8788 0.6733 0.7536 [13] 0.7905 0.7696 0.7447 0.7303 0.7644 0.7784 0.7692 0.7035 0.7609 0.7927 0.6750 0.9684 [25] 1.0000 0.8273 0.7703 0.9665 0.9082 0.8864 0.9524 0.9884 0.8900 0.8775 0.8854 0.7778 [37] 0.7885 0.8817 0.9091 0.8409 0.8307 0.8889 0.9053 0.8681 0.7941 0.9326 0.8602 0.7135 [49] 0.7960 0.7222 0.7773 0.7745 0.8750 0.7474 0.7720 0.9385 0.8750 0.7632 0.8213 0.8290 [61] 0.8220 0.9293 0.8821 0.9270 0.7629 0.7750 0.7500 0.6919 0.7772 $bodyfat1 [1] 32.2 12.9 33.4 14.8 25.0 28.2 35.6 42.6 24.5 35.0 25.8 41.8 25.9 19.7 27.4 24.5 28.1 [18] 35.1 36.2 52.0 45.7 25.6 21.5 20.8 16.0 17.7 16.7 35.3 28.6 15.1 9.5 20.3 30.5 44.4 [35] 20.6 23.0 6.0 14.0 23.7 21.6 33.4 12.9 11.4 21.2 43.2 23.7 12.5 37.1 35.3 39.3 43.8 [52] 40.8 23.7 38.0 36.4 20.1 11.0 39.8 10.4 22.4 10.4 14.7 19.1 17.7 28.4 33.2 38.5 16.4 [69] 35.6 $bodyfat2 [1] 37.6 10.4 31.3 13.1 26.2 34.5 38.6 42.6 23.8 33.7 22.7 31.1 42.9 24.5 25.6 28.6 27.7 [18] 28.9 38.2 44.5 45.8 26.7 24.9 25.1 19.7 18.4 12.2 31.7 23.6 13.5 21.0 26.0 29.7 37.2 [35] 29.9 22.4 7.9 15.4 26.1 16.9 34.0 10.1 18.6 16.0 31.9 18.6 13.2 20.6 36.1 30.6 44.6 [52] 30.5 15.3 29.0 27.6 25.5 10.4 26.3 18.6 13.9 11.5 16.9 19.9 18.0 28.4 34.1 36.8 30.5 [69] 30.4 $bodymass1 [1] 62.755 70.422 63.843 50.518 59.013 58.709 64.821 70.136 78.122 65.569 57.791 [12] 81.528 52.162 52.492 62.275 53.819 57.213 65.634 70.004 93.446 74.539 62.870 [23] 59.359 94.562 59.874 92.980 89.830 101.267 83.198 67.202 71.280 64.382 80.276 [34] 82.719 54.698 53.261 78.160 74.320 96.855 59.809 62.591 73.327 74.941 65.897 [45] 76.679 70.425 67.138 80.666 64.669 80.551 69.274 85.168 91.031 63.282 53.389 [56] 68.269 75.871 73.455 75.832 80.246 64.592 82.226 71.036 82.347 61.669 74.431 [67] 67.959 50.839 68.533 $bodymass2 [1] 67.334 65.300 64.471 50.250 60.182 70.012 66.211 70.437 76.272 61.302 57.448 69.118 [13] 73.100 57.188 61.892 54.152 56.987 64.088 70.741 76.648 80.059 62.089 61.022 92.303 [25] 58.246 94.160 87.157 93.330 77.144 64.314 75.526 70.508 79.744 70.632 58.740 51.876 [37] 82.614 73.262 96.703 64.448 60.821 80.483 78.502 68.300 65.212 73.456 74.991 57.436 [49] 65.073 71.756 70.079 72.342 80.916 58.320 57.572 73.714 75.214 60.249 84.795 74.133 [61] 75.504 98.320 79.119 71.388 65.354 67.660 64.833 60.393 48.845 $bmi1 [1] 22.1029 22.1268 22.6741 17.8990 20.4197 19.7758 23.3823 26.7246 23.8460 27.9736 20.8715 [12] 31.5308 18.8159 18.7768 22.1436 19.9857 21.6932 24.0788 25.8069 35.6066 30.0485 21.3505 [23] 20.4192 28.5479 22.5077 27.1673 25.9934 33.6798 28.1226 24.0962 23.1427 22.4091 28.9572 [34] 31.0574 21.2865 21.7183 23.5962 24.0474 29.3368 22.3187 24.0273 22.1371 23.6261 21.1772 [45] 28.5791 23.7223 19.9604 26.7049 26.2026 30.3549 24.7798 30.2117 28.7309 24.7195 21.5224 [56] 22.9427 23.3650 27.6130 22.7183 24.1463 24.4910 22.1895 23.7898 25.4440 21.5410 24.4155 [67] 25.7994 19.8590 26.8378 $bmi2 [1] 23.6314 21.6927 22.6534 18.1262 20.7997 22.6535 24.2611 26.7731 22.6513 26.4982 20.6482 [12] 26.5329 26.4004 20.5794 21.6189 20.0360 21.7411 22.9522 26.3659 29.3870 30.6948 20.6260 [23] 21.2396 28.1125 22.0035 27.8119 25.6597 29.6896 26.4440 23.5090 24.8600 24.2828 28.4909 [34] 27.5906 22.3823 20.4682 24.4281 23.4377 28.7816 22.6722 22.8356 23.9280 24.4186 21.8998 [45] 25.0956 24.2621 22.1019 20.9695 25.4509 25.2431 26.0874 25.2393 23.9521 21.5519 21.8024 [56] 26.4629 20.8569 20.2475 27.4993 23.8510 25.1695 27.4663 24.7207 23.2572 24.2988 24.6426 [67] 23.9588 22.3180 23.7201 $sys1 [1] 116.6000 131.4000 123.0000 111.0000 111.2000 118.4000 126.4000 110.0000 146.0000 [10] 135.6000 104.8333 135.2000 124.1667 115.0000 112.0000 115.2000 119.6000 101.0000 [19] 125.6667 122.0000 135.4000 122.5000 113.8000 123.6000 122.2500 140.0000 130.0000 [28] 138.2000 138.0000 134.6000 119.8000 123.5000 124.8333 112.4000 142.8000 98.0000 [37] 127.0000 113.2000 129.1667 125.8000 118.6000 110.6000 136.4000 122.2000 102.6000 [46] 133.6000 107.8000 119.2000 137.6000 108.4000 124.0000 121.8333 132.6000 103.8000 [55] 136.6000 142.0000 121.8000 130.2000 121.8000 127.6000 115.5000 110.6000 128.0000 [64] 114.6000 110.4000 116.0000 105.2000 103.0000 114.5000 $sys2 [1] 108.0000 119.8000 125.6000 118.0000 115.2000 111.6000 127.4000 111.1429 142.6000 [10] 114.6000 107.5000 118.8000 123.5000 119.4000 108.6000 111.2000 117.4000 102.0000 [19] 115.8000 116.0000 120.6000 118.2000 108.2000 127.0000 120.8000 125.6000 125.2000 [28] 139.0000 131.0000 133.8000 125.0000 124.6000 133.0000 93.6667 141.2000 103.2000 [37] 124.2000 124.4000 144.2000 129.8000 112.4000 109.6000 132.8000 119.6000 104.6000 [46] 139.1667 124.4000 120.2000 126.0000 106.4000 130.0000 108.2000 144.0000 110.4000 [55] 101.8000 164.0000 126.6667 123.6000 125.5000 126.2000 132.0000 127.5000 125.6000 [64] 116.8000 105.6000 113.4000 113.4000 98.5000 150.5000 [[22]] [1] "Zygosity" [[23]] [1] "Gender1" [[24]] [1] "Gender2" [[25]] [1] "Age1" [[26]] [1] "Age2" [[27]] [1] "height1" [[28]] [1] "height2" [[29]] [1] "waist1" [[30]] [1] "waist2" [[31]] [1] "hip1" [[32]] [1] "hip2" [[33]] [1] "whr1" [[34]] [1] "whr2" [[35]] [1] "bodyfat1" [[36]] [1] "bodyfat2" [[37]] [1] "bodymass1" [[38]] [1] "bodymass2" [[39]] [1] "bmi1" [[40]] [1] "bmi2" [[41]] [1] "sys1" [[42]] [1] "sys2" > colnames(mzData) = c(VarsA, VarsB, 'Age1', 'Age2', 'Gender1', 'Gender2') > colnames(dzData) = c(VarsA, VarsB, 'Age1', 'Age2', 'Gender1', 'Gender2') > summary(mzData) bmi1 bmi2 Age1 Age2 Gender1 Min. :17.90 Min. :18.13 Min. :19.00 Min. :19.00 Min. :0.0000 1st Qu.:21.69 1st Qu.:21.90 1st Qu.:25.00 1st Qu.:25.00 1st Qu.:0.0000 Median :23.38 Median :23.63 Median :36.00 Median :36.00 Median :1.0000 Mean :24.35 Mean :24.16 Mean :36.49 Mean :36.49 Mean :0.5778 3rd Qu.:27.17 3rd Qu.:26.44 3rd Qu.:45.00 3rd Qu.:45.00 3rd Qu.:1.0000 Max. :35.61 Max. :30.69 Max. :60.00 Max. :60.00 Max. :1.0000 Gender2 Min. :0.0000 1st Qu.:0.0000 Median :1.0000 Mean :0.5778 3rd Qu.:1.0000 Max. :1.0000 > summary(dzData) bmi1 bmi2 Age1 Age2 Gender1 Min. :19.86 Min. :20.25 Min. :19.00 Min. :19.00 Min. :0.0000 1st Qu.:22.89 1st Qu.:22.26 1st Qu.:21.75 1st Qu.:21.75 1st Qu.:0.0000 Median :24.45 Median :24.11 Median :32.00 Median :32.00 Median :1.0000 Mean :24.67 Mean :23.96 Mean :33.71 Mean :33.71 Mean :0.5417 3rd Qu.:26.33 3rd Qu.:25.24 3rd Qu.:41.00 3rd Qu.:41.00 3rd Qu.:1.0000 Max. :30.35 Max. :27.50 Max. :65.00 Max. :65.00 Max. :1.0000 Gender2 Min. :0.0000 1st Qu.:0.0000 Median :1.0000 Mean :0.5417 3rd Qu.:1.0000 Max. :1.0000 > colMeans(mzData,na.rm=TRUE) bmi1 bmi2 Age1 Age2 Gender1 Gender2 24.3463289 24.1636711 36.4888889 36.4888889 0.5777778 0.5777778 > colMeans(dzData,na.rm=TRUE) bmi1 bmi2 Age1 Age2 Gender1 Gender2 24.6692375 23.9637583 33.7083333 33.7083333 0.5416667 0.5416667 > cov(mzData,use="complete") bmi1 bmi2 Age1 Age2 Gender1 Gender2 bmi1 16.7449374 10.4727642 11.511411 11.511411 -0.3239989 -0.3239989 bmi2 10.4727642 8.8089407 7.154574 7.154574 -0.3560034 -0.3560034 Age1 11.5114106 7.1545735 140.755556 140.755556 1.5292929 1.5292929 Age2 11.5114106 7.1545735 140.755556 140.755556 1.5292929 1.5292929 Gender1 -0.3239989 -0.3560034 1.529293 1.529293 0.2494949 0.2494949 Gender2 -0.3239989 -0.3560034 1.529293 1.529293 0.2494949 0.2494949 > cov(dzData,use="complete") bmi1 bmi2 Age1 Age2 Gender1 Gender2 bmi1 7.9554556 0.1690636 -3.0705364 -3.0705364 0.4287571 0.4287571 bmi2 0.1690636 4.0587545 -0.4444953 -0.4444953 -0.2608069 -0.2608069 Age1 -3.0705364 -0.4444953 217.2590580 217.2590580 2.3822464 2.3822464 Age2 -3.0705364 -0.4444953 217.2590580 217.2590580 2.3822464 2.3822464 Gender1 0.4287571 -0.2608069 2.3822464 2.3822464 0.2590580 0.2590580 Gender2 0.4287571 -0.2608069 2.3822464 2.3822464 0.2590580 0.2590580 > MZlist=as.vector(as.matrix(mzData)) > MZfinite=MZlist[is.finite(MZlist)] > DZlist=as.vector(as.matrix(dzData)) > DZfinite=DZlist[is.finite(DZlist)] > Nobs=length(MZfinite)=length(DZfinite) > print(head(mzData)) bmi1 bmi2 Age1 Age2 Gender1 Gender2 1 22.1029 23.6314 37 37 1 1 2 22.1268 21.6927 23 23 0 0 3 22.6741 22.6534 20 20 1 1 4 17.8990 18.1262 49 49 1 1 5 20.4197 20.7997 34 34 1 1 6 19.7758 22.6535 19 19 1 1 > print(head(dzData)) bmi1 bmi2 Age1 Age2 Gender1 Gender2 46 23.7223 24.2621 33 33 0 0 47 19.9604 22.1019 34 34 0 0 48 26.7049 20.9695 19 19 1 1 49 26.2026 25.4509 59 59 1 1 50 30.3549 25.2431 19 19 1 1 51 24.7798 26.0874 44 44 1 1 > ################# > #Saturated model > ################# > mytwinSatModel <- mxModel("twinSat", + mxModel("MZ", + mxMatrix(type="Full", nrow=1,ncol= 2,free=TRUE,values= c(0,0),label= "mean", name="expMean"), + mxMatrix(type="Full", nrow=1, ncol=2, free=TRUE, values= 0.1, label=c("betaGender "), name="bGender"), + mxMatrix(type="Full", nrow=1, ncol=2, free=TRUE, values= 0.1, label=c("betaAge"), name="bAge"), + mxMatrix("Lower", nrow= 2, ncol=2,free=TRUE,values=.5,name="CholMZ"), + mxAlgebra(expression= CholMZ %*% t(CholMZ), name="expCovMZ"), + mxMatrix(type="Full",nrow=2,ncol=2,free=F,label=c("data.Age1", "data.Age2"), name="Age"), + mxMatrix(type="Full",nrow=2,ncol=2,free=F,label=c("data.Gender1","data.Gender2"), name="Gender"), + meanAge= mxAlgebra(bAge%*%Age, name= "AgeR"), + meanGender = mxAlgebra (bGender%*%Gender, name= "GenderR"), + mxAlgebra( expression=expMean +AgeR + GenderR , name="expMeanMZ"), + mxData(mzData, type="raw"), + mxExpectationNormal("expCovMZ", "expMeanMZ", dimnames=selVars) , + mxFitFunctionML()), + + mxModel("DZ", + mxMatrix(type="Full", nrow=1,ncol= 2,free=TRUE,values= c(0,0),label= "mean", name="expMean"), + mxMatrix(type="Full", nrow=1, ncol=2, free=TRUE, values= 0.1, label=c("betaAge"), name="bAge"), + mxMatrix(type="Full", nrow=1, ncol=2, free=TRUE, values= 0.1, label=c("betaGender"), name="bGender"), + mxMatrix("Lower", nrow= 2, ncol=2,free=TRUE,values=.5,name="CholDZ"), + mxAlgebra(expression= CholDZ %*% t(CholDZ), name="expCovDZ"), + mxMatrix(type="Full",nrow=2,ncol=2,free=F,label=c("data.Age1", "data.Age2"), name="Age"), + mxMatrix(type="Full",nrow=2,ncol=2,free=F,label=c("data.Gender1", "data.Gender2"), name="Gender"), + meanAge= mxAlgebra(bAge%*%Age, name="AgeR"), + meanGender = mxAlgebra (bGender%*%Gender, name="GenderR"), + mxAlgebra( expression=expMean + AgeR + GenderR, name="expMeanDZ"), + mxData(dzData, type="raw"), + mxExpectationNormal("expCovDZ", "expMeanDZ", dimnames=selVars) , + mxFitFunctionML()), + mxFitFunctionMultigroup( c("MZ.fitfunction", "DZ.fitfunction"))) > mytwinSatFit <- mxRun(mytwinSatModel) #The mxRun command evaluates the model. Running twinSat with 10 parameters > LL_Sat <- mxEval(objective, mytwinSatFit) > print(summary(mxRun(mytwinSatModel))) Running twinSat with 10 parameters Summary of twinSat free parameters: name matrix row col Estimate Std.Error A 1 mean MZ.expMean 1 1 23.983276148 0.74038249 2 betaGender MZ.bGender 1 1 -0.617519874 0.34663600 3 betaAge MZ.bAge 1 1 0.009317255 0.01066006 4 MZ.CholMZ[1,1] MZ.CholMZ 1 1 -3.959701616 0.42330905 5 MZ.CholMZ[2,1] MZ.CholMZ 2 1 -2.402890090 0.34083105 6 MZ.CholMZ[2,2] MZ.CholMZ 2 2 -1.454524462 0.15365570 7 betaGender DZ.bGender 1 1 -0.305461946 0.34768174 8 DZ.CholDZ[1,1] DZ.CholDZ 1 1 2.915439768 0.45029661 9 DZ.CholDZ[2,1] DZ.CholDZ 2 1 0.107348756 0.41621054 10 DZ.CholDZ[2,2] DZ.CholDZ 2 2 1.949847149 0.29412203 observed statistics: 138 estimated parameters: 10 degrees of freedom: 128 fit value ( -2lnL units ): 632.6164 number of observations: 69 Information Criteria: | df Penalty | Parameters Penalty | Sample-Size Adjusted AIC: 376.61640 652.6164 NA BIC: 90.65077 674.9575 643.4627 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: 2016-07-05 10:10:12 Wall clock time (HH:MM:SS.hh): 00:00:00.32 optimizer: SLSQP OpenMx version number: 2.5.2 Need help? See help(mxSummary) > DF_Sat=Nobs-nrow(mytwinSatFit@output$standardErrors) > > ################ > #Run ACE model > ################ > #start value for paths (Variance) is (sqrt(variance/3)) > twinACEModel <- mxModel("twinACE", + mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=Strtpaths, label="a", name="X" ), + mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=Strtpaths, label="c", name="Y" ), + mxMatrix( type="Full", nrow=1, ncol=1, free=TRUE, values=Strtpaths, label="e", name="Z" ), + mxAlgebra( expression=X %*% t(X), name="A" ), + mxAlgebra( expression=Y %*% t(Y), name="C" ), + mxAlgebra( expression=Z %*% t(Z), name="E" ), + mxMatrix(type="Full", nrow=1,ncol= 2,free=TRUE,values= StrtMean,label= "mean", name="expMean"), + mxMatrix(type="Full", nrow=1, ncol=2, free=TRUE, values= 0.1, label=c("betaAge"), name="bAge"), + mxMatrix (type="Full", nrow=1, ncol=2, free=TRUE, values= 0.1, label=c("betaGender"), name="bGender"), + mxAlgebra(expression=rbind(cbind(A+C+E , A+C),cbind(A+C, A+C+E)),name="expCovMZ"), + mxAlgebra(expression=rbind(cbind(A+C+E,0.5%x%A+C),cbind(0.5%x%A+C,A+C+E)), name="expCovDZ"), + + mxModel("MZ", + mxData( observed=mzData, type="raw" ), + mxMatrix(type="Full",nrow=2,ncol=2,free=F, label=c("data.Age1", "data.Age2"), name="Age"), + mxMatrix(type="Full",nrow=2,ncol=2,free=F, label=c("data.Gender1", "data.Gender2"), name="Gender"), + meanAge= mxAlgebra(twinACE.bAge%*%Age, name="AgeR"), + meanGender = mxAlgebra (twinACE.bGender%*%Gender, name="GenderR"), + mxAlgebra( expression=twinACE.expMean + GenderR + AgeR, name="expMeanMZ"), + mxExpectationNormal( covariance="twinACE.expCovMZ", means="expMeanMZ", dimnames=selVars), + mxFitFunctionML()), + + mxModel("DZ", + mxData( observed=dzData, type="raw" ), + mxMatrix(type="Full",nrow=2,ncol=2,free=F, label=c("data.Age1", "data.Age2"), name="Age"), + mxMatrix(type="Full",nrow=2,ncol=2,free=F, label=c("data.Gender1", "data.Gender2"), name="Gender"), + meanAge= mxAlgebra(twinACE.bAge%*%Age, name="AgeR"), + meanGender = mxAlgebra (twinACE.bGender%*%Gender, name="GenderR"), + mxAlgebra( expression=twinACE.expMean + GenderR + AgeR, name="expMeanDZ"), + mxExpectationNormal( covariance="twinACE.expCovDZ", means="expMeanDZ", dimnames=selVars), + mxFitFunctionML()), + mxFitFunctionMultigroup( c("MZ.fitfunction", "DZ.fitfunction"))) > twinACEFit <- mxRun(twinACEModel) Running twinACE with 6 parameters > print(summary(twinACEFit)) Summary of twinACE free parameters: name matrix row col Estimate Std.Error A 1 a X 1 1 2.697179e+00 0.25363808 2 c Y 1 1 2.603630e-06 1.58913890 3 e Z 1 1 1.498154e+00 0.15347292 4 mean expMean 1 1 2.345935e+01 0.97588689 5 betaAge bAge 1 1 1.925819e-02 0.01341096 6 betaGender bGender 1 1 -4.807598e-01 0.35056534 observed statistics: 138 estimated parameters: 6 degrees of freedom: 132 fit value ( -2lnL units ): 659.3252 number of observations: 69 Information Criteria: | df Penalty | Parameters Penalty | Sample-Size Adjusted AIC: 395.3252 671.3252 NA BIC: 100.4232 684.7299 665.833 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: 2016-07-05 10:10:12 Wall clock time (HH:MM:SS.hh): 00:00:00.16 optimizer: SLSQP OpenMx version number: 2.5.2 Need help? See help(mxSummary) > > ########################################################## > #Summary statistics > ########################################################## > DF_ACE=Nobs-nrow(twinACEFit@output$standardErrors) > LL_ACE <- mxEval(objective, twinACEFit) > mychi_ACE= LL_ACE - LL_Sat > mychi_DF_ACE=DF_ACE-DF_Sat > mychi_p_ACE=1-pchisq(mychi_ACE,mychi_DF_ACE) > expMZcov_ACE <- mxEval(expCovMZ, twinACEFit) > expDZcov_ACE <- mxEval(expCovDZ, twinACEFit) > expMeans_ACE <- mxEval(expMean, twinACEFit) > A_ACE <- mxEval(a*a, twinACEFit) > C_ACE <- mxEval(c*c, twinACEFit) > E_ACE <- mxEval(e*e, twinACEFit) > V <- (A_ACE+C_ACE+E_ACE) > a2_ACE <- A_ACE/V > c2_ACE <- C_ACE/V > e2_ACE <- E_ACE/V > ACE_mySE=round(twinACEFit@output$standardErrors,3) > ACE_myest=round(twinACEFit@output$estimate,3) > ACE_mylower=round(ACE_myest-1.96*ACE_mySE,3) > ACE_myupper=round(ACE_myest+1.96*ACE_mySE,3) > > ################ > #Run AE model > ################ > twinAEModel <- mxModel(twinACEModel, name = "twinACE", + mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=0, label="c", name="Y")) > twinAEFit <- mxRun(twinAEModel) Running twinACE with 5 parameters > > ########################################################## > #Summary statistics > ########################################################## > LL_AE <- mxEval(objective, twinAEFit) > DF_AE=Nobs-nrow(twinAEFit@output$standardErrors) > mychi_AE= LL_AE - LL_Sat #subtract LL for Saturated model from LL for AE > mychi_DF_AE=DF_AE-DF_Sat #subtract DF for Saturated model from DF for AE > mychi_p_AE=1-pchisq(mychi_AE,mychi_DF_AE)# compute chi square probability > #expected values for covs and means can be found in mxEval(expCovMZ, twinAEFit) > A_AE <- mxEval(a*a, twinAEFit) > C_AE <- mxEval(c*c, twinAEFit) > E_AE <- mxEval(e*e, twinAEFit) > V <- (A_AE+C_AE+E_AE) > a2_AE <- A_AE/V > c2_AE <- C_AE/V > e2_AE <- E_AE/V > mychiAEdiff=mychi_AE-mychi_ACE > myDFAEdiff=mychi_DF_AE-mychi_DF_ACE > mychiAEdiff_p=1-pchisq(mychiAEdiff,myDFAEdiff) > AE_mySE=round(twinAEFit@output$standardErrors,3) > AE_myest=round(twinAEFit@output$estimate,3) > AE_mylower=round(AE_myest-1.96*AE_mySE,3) > AE_myupper=round(AE_myest+1.96*AE_mySE,3) > print(summary(mxRun(twinAEModel))) Running twinACE with 5 parameters Summary of twinACE free parameters: name matrix row col Estimate Std.Error A 1 a X 1 1 -2.6971796 0.2536364 2 e Z 1 1 1.4981540 0.1534725 3 mean expMean 1 1 23.4593494 0.9758800 4 betaAge bAge 1 1 0.0192582 0.0134109 5 betaGender bGender 1 1 -0.4807605 0.3505241 observed statistics: 138 estimated parameters: 5 degrees of freedom: 133 fit value ( -2lnL units ): 659.3252 number of observations: 69 Information Criteria: | df Penalty | Parameters Penalty | Sample-Size Adjusted AIC: 393.32525 669.3252 NA BIC: 96.18908 680.4958 664.7484 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: 2016-07-05 10:10:13 Wall clock time (HH:MM:SS.hh): 00:00:00.17 optimizer: SLSQP OpenMx version number: 2.5.2 Need help? See help(mxSummary) > > ################ > #Run CE model > ################ > twinCEModel <- mxModel(twinACEModel, name = "twinACE", + mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=0, label="a", name="X")) > twinCEFit <- mxRun(twinCEModel) Running twinACE with 5 parameters > DF_CE=Nobs-nrow(twinCEFit@output$standardErrors) > LL_CE <- mxEval(objective, twinCEFit) > mychi_CE= LL_CE - LL_Sat > mychi_DF_CE=DF_CE-DF_Sat > mychi_p_CE=1-pchisq(mychi_CE,mychi_DF_CE) > A_CE <- mxEval(a*a, twinCEFit) > C_CE <- mxEval(c*c, twinCEFit) > E_CE <- mxEval(e*e, twinCEFit) > V <- (A_CE+C_CE+E_CE) > a2_CE <- A_CE/V > c2_CE <- C_CE/V > e2_CE <- E_CE/V > mychiCEdiff=mychi_CE-mychi_ACE > myDFCEdiff=mychi_DF_CE-mychi_DF_ACE > mychiCEdiff_p=1-pchisq(mychiCEdiff,myDFCEdiff) > CE_mySE=round(twinCEFit@output$standardErrors,3) > CE_myest=round(twinCEFit@output$estimate,3) > CE_mylower=round(CE_myest-1.96*CE_mySE,3) > CE_myupper=round(CE_myest+1.96*CE_mySE,3) > print(summary(mxRun(twinCEModel))) Running twinACE with 5 parameters Summary of twinACE free parameters: name matrix row col Estimate Std.Error A 1 c Y 1 1 -2.49854038 0.27902260 2 e Z 1 1 1.87424076 0.15954573 3 mean expMean 1 1 23.31576975 1.01207046 4 betaAge bAge 1 1 0.02223037 0.01386036 5 betaGender bGender 1 1 -0.54731311 0.35797566 observed statistics: 138 estimated parameters: 5 degrees of freedom: 133 fit value ( -2lnL units ): 669.6199 number of observations: 69 Information Criteria: | df Penalty | Parameters Penalty | Sample-Size Adjusted AIC: 403.6199 679.6199 NA BIC: 106.4838 690.7905 675.0431 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: 2016-07-05 10:10:14 Wall clock time (HH:MM:SS.hh): 00:00:00.16 optimizer: SLSQP OpenMx version number: 2.5.2 Need help? See help(mxSummary) > > ################ > #Run E model > ################ > twinEModel <- mxModel(twinAEModel, name = "twinACE", mxMatrix(type="Full", nrow=1, ncol=1, free=FALSE, values=0, label="a", name="X")) > twinEFit <- mxRun(twinEModel) Running twinACE with 4 parameters > DF_E=Nobs-nrow(twinEFit@output$standardErrors) > LL_E <- mxEval(objective, twinEFit) > mychi_E= LL_E - LL_Sat > mychi_DF_E=DF_E-DF_Sat > mychi_p_E=1-pchisq(mychi_E,mychi_DF_E) > A_E <- mxEval(a*a, twinEFit) > C_E <- mxEval(c*c, twinEFit) > E_E <- mxEval(e*e, twinEFit) > V <- (A_E+C_E+E_E) > a2_E <- A_E/V > c2_E <- C_E/V > e2_E <- E_E/V > mychiEdiff=mychi_E-mychi_ACE > myDFEdiff=mychi_DF_E-mychi_DF_ACE > mychiEdiff_p=1-pchisq(mychiEdiff,myDFEdiff) > E_mySE=round(twinEFit@output$standardErrors,3) > E_myest=round(twinEFit@output$estimate,3) > E_mylower=round(E_myest-1.96*E_mySE,3) > E_myupper=round(E_myest+1.96*E_mySE,3) > print(summary(mxRun(twinEModel))) Running twinACE with 4 parameters Summary of twinACE free parameters: name matrix row col Estimate Std.Error A 1 e Z 1 1 3.12337683 0.18800534 2 mean expMean 1 1 23.31576575 0.79031518 3 betaAge bAge 1 1 0.02223044 0.01082342 4 betaGender bGender 1 1 -0.54731492 0.27954708 observed statistics: 138 estimated parameters: 4 degrees of freedom: 134 fit value ( -2lnL units ): 705.9675 number of observations: 69 Information Criteria: | df Penalty | Parameters Penalty | Sample-Size Adjusted AIC: 437.9675 713.9675 NA BIC: 138.5972 722.9039 710.306 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: 2016-07-05 10:10:15 Wall clock time (HH:MM:SS.hh): 00:00:00.15 optimizer: SLSQP OpenMx version number: 2.5.2 Need help? See help(mxSummary) > ########################## > #Output to compare models > ########################## > myoutput <- rbind(cbind("__________________________","_______________","_________________","______________"), + + cbind("ACE model","A","C","E"), cbind("Unsquared path estimates", ACE_myest[1],ACE_myest[2],ACE_myest[3]), + + cbind("Standard errors",ACE_mySE[1],ACE_mySE[2],ACE_mySE[3]), + + cbind("Lower 95% CI",ACE_mylower[1], ACE_mylower[2], ACE_mylower[3]), cbind("Upper 95% CI",ACE_myupper[1], ACE_myupper[2], ACE_myupper[3]), cbind(".",".",".","."), + + cbind("Unstandardized variance comps",round(A_ACE,3), round(C_ACE,3), round(E_ACE,3)), cbind("Standardized variance comps",round(a2_ACE,3),round(c2_ACE,3),round(e2_ACE,3)), + + cbind(".","chisq","DF","p"), cbind("saturated vs. ACE",round(mychi_ACE,3), mychi_DF_ACE,round(mychi_p_ACE,3)), cbind("__________________________","_______________","_________________","______________"), + + cbind("AE model","A","C","E"), cbind("Unsquared path estimates", AE_myest[1],".", AE_myest[2]), cbind("Standard errors",AE_mySE[1],".",AE_mySE[2]), cbind("Lower 95% CI",AE_mylower[1],".",AE_mylower[2]), cbind("Upper 95% CI",AE_myupper[1],".",AE_myupper[2]), cbind(".",".",".","."), cbind("Unstandardized variance comps",round(A_AE,3), round(C_AE,3), round(E_AE,3)), cbind("Standardized variance comps",round(a2_AE,3), round(c2_AE,3),round(e2_AE,3)), cbind(".","chisq","DF","p"), + + cbind("saturated vs. AE",round(mychi_AE,3), mychi_DF_AE,round(mychi_p_AE,3)), cbind("ACE vs AE",round (mychiAEdiff,3),1,round(mychiAEdiff_p,3)), cbind(".",".",".","."), cbind ("__________________________","_______________","_________________","______________"), + + cbind("CE model","A","C","E"), cbind("Unsquared path estimates",".", CE_myest[1],CE_myest[2]), cbind("Standard errors",".",CE_mySE[1],CE_mySE[2]),cbind("Lower 95% CI",".", CE_mylower[1], CE_mylower[2]), cbind("Upper 95% CI",".", CE_myupper[1], CE_myupper[2]), cbind(".",".",".","."), cbind("Unstandardized variance comps",round(A_CE,3),round(C_CE,3),round(E_CE,3)), + + cbind("Standardized variance comps",round (a2_CE,3),round(c2_CE,3),round(e2_CE,3)), cbind(".","chisq","DF","p"), + + cbind("saturated vs. CE",round(mychi_CE,3), mychi_DF_CE,round(mychi_p_CE,3)), cbind("ACE vs CE",round(mychiCEdiff,3), 1, round(mychiCEdiff_p,3)), cbind ("__________________________","_______________","_________________","______________"), + + cbind("E model","A","C","E"), cbind("Unsquared path estimates",".", E_myest[1],E_myest[2]), cbind("Standard errors",".",E_mySE[1],E_mySE[2]),cbind("Lower 95% CI",".", E_mylower[1], E_mylower[2]), cbind("Upper 95% CI",".", E_myupper[1], E_myupper[2]), cbind(".",".",".","."), cbind("Unstandardized variance comps",round(A_E,3),round(C_E,3),round(E_E,3)), + + cbind("Standardized variance comps",round (a2_E,3),round(c2_E,3),round(e2_E,3)), cbind(".","chisq","DF","p"), + + cbind("saturated vs. E",round(mychi_E,3), mychi_DF_E,round(mychi_p_E,3)), cbind("ACE vs E",round(mychiEdiff,3), 1, round(mychiEdiff_p,3)), cbind ("__________________________","_______________","_________________","______________"),cbind(date(),".",".",".")) > > myoutput2 <- data.frame(myoutput) Warning message: In data.row.names(row.names, rowsi, i) : some row.names duplicated: 2,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,28,29,30,31,32,33,34,35,36,37,38,40,41,42,43,44,45,46,47,48,49,50 --> row.names NOT used > > names(myoutput2)=selVars > > #write.table(myoutput2, "myfileout.txt",sep = "\t",quote=F) > > write.table(myoutput2, sprintf("%s.txt", VarsA) ,sep = "\t",quote=F) > > print(myoutput2) bmi1 bmi2 NA NA 1 __________________________ _______________ _________________ ______________ 2 ACE model A C E 3 Unsquared path estimates 2.697 0 1.498 4 Standard errors 0.254 1.589 0.153 5 Lower 95% CI 2.199 -3.114 1.198 6 Upper 95% CI 3.195 3.114 1.798 7 . . . . 8 Unstandardized variance comps 7.275 0 2.244 9 Standardized variance comps 0.764 0 0.236 10 . chisq DF p 11 saturated vs. ACE 26.709 4 0 12 __________________________ _______________ _________________ ______________ 13 AE model A C E 14 Unsquared path estimates -2.697 . 1.498 15 Standard errors 0.254 . 0.153 16 Lower 95% CI -3.195 . 1.198 17 Upper 95% CI -2.199 . 1.798 18 . . . . 19 Unstandardized variance comps 7.275 0 2.244 20 Standardized variance comps 0.764 0 0.236 21 . chisq DF p 22 saturated vs. AE 26.709 5 0 23 ACE vs AE 0 1 1 24 . . . . 25 __________________________ _______________ _________________ ______________ 26 CE model A C E 27 Unsquared path estimates . -2.499 1.874 28 Standard errors . 0.279 0.16 29 Lower 95% CI . -3.046 1.56 30 Upper 95% CI . -1.952 2.188 31 . . . . 32 Unstandardized variance comps 0 6.243 3.513 33 Standardized variance comps 0 0.64 0.36 34 . chisq DF p 35 saturated vs. CE 37.004 5 0 36 ACE vs CE 10.295 1 0.001 37 __________________________ _______________ _________________ ______________ 38 E model A C E 39 Unsquared path estimates . 3.123 23.316 40 Standard errors . 0.188 0.79 41 Lower 95% CI . 2.755 21.768 42 Upper 95% CI . 3.491 24.864 43 . . . . 44 Unstandardized variance comps 0 0 9.755 45 Standardized variance comps 0 0 1 46 . chisq DF p 47 saturated vs. E 73.351 6 0 48 ACE vs E 46.642 1 0 49 __________________________ _______________ _________________ ______________ 50 Tue Jul 5 10:10:15 2016 . . . > > mymessage="If table not formatted properly, make console screen full size and type myoutput2" > > > > } Error: unexpected '}' in "}" > > >