!  Model: recursive model of VIQ, GPA, and CPL
!  with trivariate Cholesky decomposition of A, C, D and E
!  Groups: all (MZ, DZ, FS, HS, CO, NR)
!  Data: single entry covariances
!  Transfo: standardized within categories of RASEX_1$
!  Submodels: BACE and all submodels
!  Restrict ACE to make B identified
!  Formula for K (standardized coefficients)  corrected
!  Formula for O (standardized B) corrected
!  Mating: randomg (k=.5)
!  Bound all in [0 10]

#ngroups 8
#define nvar 3                   ! number of variables
#define twonvar 6                ! two times nvar

Group 1 - Model parameters
  Calculation
  Begin matrices;
    B  sdiag nvar nvar free      ! phenotypic regression coefficients
    X  lower nvar nvar free      ! genetic structure
    Y  lower nvar nvar free     ! common environmental structure
    Z  lower nvar nvar free      ! specific environmental structure
    V  lower nvar nvar      ! dominance structure
    I  iden nvar nvar
    H full 1 1
    Q full 1 1
    T full 1 1
  End matrices;
    Matrix H   .5
    Matrix Q   .25
    Matrix T   .125
    Start .5 all
    Bound 0 10 B 2 1 B 3 1 B 3 2 
  Begin algebra;
    A = X*X'   ;
    C = Y*Y'  ;
    D = V*V'   ;
    E = Z*Z'  ;
    W = (I - B)~  ;
  End algebra;
  Options NO_Output
End

Group 2 - MZ twins
  Data NInputvars=twonvar NObs=170
  Labels VIQ_1 GPA_1 CPL_1 VIQ_2 GPA_2 CPL_2
  CMatrix
    1.0598
    0.2741  0.9223
    0.2975  0.3609  0.9899
    0.8047  0.3193  0.2565  1.1652
    0.1827  0.6173  0.3118  0.3239  0.9489
    0.3339  0.3472  0.6587  0.4035  0.3823  0.9976
  Matrices = Group 1
  Covariance
    W*(A+C+D+E)*W'  | W*(A+C+D)*W'_
    W*(A+C+D)*W'      | W*(A+C+D+E)*W' ;
  Options NO_Output
End

Group 3 - DZ twins
  Data NInputvars=twonvar NObs=290
  Labels VIQ_1 GPA_1 CPL_1 VIQ_2 GPA_2 CPL_2
  CMatrix
  0.8260
  0.2022  0.8680
  0.1474  0.2404  0.8911
  0.3204  0.1219  0.1269  0.9832
  0.0408  0.2943  0.0942  0.2352  0.9081
  0.0270  0.0489  0.2120  0.1146  0.2369  0.7247
  Matrices = Group 1
  Covariance
    W*(A+C+D+E)*W'   | W*(H@A+Q@D+C)*W'_
    W*(H@A+Q@D+C)*W'   | W*(A+C+D+E)*W' ;
! By using Kronecker product H@A each element of A is multiplied by .5
  Options NO_Output
End

Group 4 - FS Full siblings
  Data NInputvars=twonvar NObs=702
  Labels VIQ_1 GPA_1 CPL_1 VIQ_2 GPA_2 CPL_2
  CMatrix
  1.0769
  0.2980  0.9499
  0.2181  0.3671  0.9807
  0.4152  0.1797  0.1422  0.9484
  0.1605  0.3326  0.2299  0.2441  0.8966
  0.1346  0.2156  0.3218  0.2402  0.3766  0.9580
  Matrices = Group 1
  Covariance
W*(A+C+D+E)*W'   | W*(H@A+Q@D+C)*W'_
    W*(H@A+Q@D+C)*W'   | W*(A+C+D+E)*W' ;
  Options NO_Output
End

Group 5 - HS Half siblings
  Data NInputvars=twonvar NObs=242
  Labels VIQ_1 GPA_1 CPL_1 VIQ_2 GPA_2 CPL_2
  CMatrix
  0.9014
  0.1846  1.1316
  0.2399  0.4942  1.1443
  0.2967  0.0984  0.1197  1.0185
 -0.1058  0.3188  0.1181  0.1388  1.1652
  0.1773  0.0880  0.2229  0.2092  0.3244  1.0395
  Matrices = Group 1
  Covariance
    W*(A+C+D+E)*W'   | W*(Q@A+C)*W'_
    W*(Q@A+C)*W'   | W*(A+C+D+E)*W' ;
  Options NO_Output
End

Group 6 - CO Cousins
  Data NInputvars=twonvar NObs=105
  Labels VIQ_1 GPA_1 CPL_1 VIQ_2 GPA_2 CPL_2
  CMatrix
  0.9116
  0.1656  1.0338
  0.0519  0.1879  0.8081
  0.2807  0.1073  0.0055  0.6902
  0.0779  0.0949 -0.0107  0.1433  0.8142
  0.2368  0.2193  0.1130  0.2344  0.2101  1.0828
  Matrices = Group 1
  Covariance
    W*(A+C+D+E)*W'   | W*(T@A+C)*W'_
    W*(T@A+C)*W'   | W*(A+C+D+E)*W' ;
  Options NO_Output
End

Group 7 - NR Unrelated siblings
  Data NInputvars=twonvar NObs=174
  Labels VIQ_1 GPA_1 CPL_1 VIQ_2 GPA_2 CPL_2
  CMatrix
   0.8979
   0.2526  0.8364
   0.1492  0.2947  0.7491
   0.0553 -0.0561  0.0792  0.8613
  -0.0878  0.0676  0.1349  0.2162  0.8480
   0.0129 -0.0070  0.1709  0.1496  0.1724  1.0792
  Matrices = Group 1
  Covariance
    W*(A+C+D+E)*W'  | W*(C)*W'_
    W*(C)*W'      | W*(A+C+D+E)*W' ;
  Options NO_Output
End

Group 8 - Calculate standardized solution, etc.
  Calculation
  Matrices = Group 1  
  Begin algebra;
! Next calculate A, C, E as proportions of total predicted
! covariance matrix S = A+C+E
    S = W*(A+C+D+E)*W'  ;
    K = A%S | C%S | D%S | E%S  ;
! Next calculate genetic, shared environmental,
! and unshared environmental correlations
  !  L = \stnd(A) | \stnd(C) | \stdn[D] |  \stnd(E) ;
! Next calculate standardized path coefficients
! For paths from latent variables just divide by predicted SD of head
     M = (\sqrt(S.I))~*X | (\sqrt(S.I))~*Y | (\sqrt(S.I))~*V | (\sqrt(S.I))~*Z ;
! For standardized coefficients for B also multiply by predicted SD of tail 
     O = (\sqrt(S.I))~*B*(\sqrt(S.I))  ; 
! Next calculate squared standardized paths = components
! of heritabilities, environmentalities, and specificities 
     N = M.M ;
  End algebra;
!  Labels Columns K A1 A2 A3 C1 C2 C3 E1 E2 E3
!  Labels Columns L A1 A2 A3 C1 C2 C3 E1 E2 E3
!  Labels Columns M A1 A2 A3 C1 C2 C3 E1 E2 E3
!  Labels Columns N A1 A2 A3 C1 C2 C3 E1 E2 E3
!  Labels Rows K VIQ GPA CPL
!  Labels Rows L VIQ GPA CPL
!  Labels Rows M VIQ GPA CPL
!  Labels Rows N VIQ GPA CPL
!  Intervals K 8 1 1 K 8 2 2 K 8 3 3 
!  Intervals K 8 1 4 K 8 2 5 K 8 3 6 
!  Intervals K 8 1 7 K 8 2 8 K 8 3 9
!  Intervals O 8 2 1 O 8 3 1 O 8 3 2 
!  Options THard=10
  Options NDecimals=4
  Options RSiduals Multiple Issat
End

Save mx0204c.mxs

! Drop A, C2, C3, off diag E -> BC1Ed
! Maximal identifiable model with B paths
get mx0204c.mxs
Drop 4 to 9
Drop 12 14 15
Drop 17 19 20
Options Multiple Issat
End

! Drop A, B, C2, C3, off diag E -> C1Ed
! To check if B identified
get mx0204c.mxs
Drop 1 2 3
Drop 4 to 9
Drop 12 14 15
Drop 17 19 20
Options Multiple
End

! Drop B -> ACE
get mx0204c.mxs
Drop 1 2 3
Options Multiple Issat
End

! Drop B & C -> AE
get mx0204c.mxs
Drop 1 2 3
Drop 10 to 15
Options Multiple
End

! Drop B & A -> CE
get mx0204c.mxs
Drop 1 2 3
Drop 4 to 9
Options Multiple
End

! Drop B & A2 A3 -> A1CE
get mx0204c.mxs
Drop 1 2 3
Drop 6 8 9 
Options Multiple
End

! Drop B & off diag A -> AdCE
get mx0204c.mxs
Drop 1 2 3
Drop 5 7 8
Options Multiple
End

! Drop B & off diag C -> ACdE
get mx0204c.mxs
Drop 1 2 3
Drop 11 13 14
Options Multiple
End

! Drop B & C3 -> AC12E
get mx0204c.mxs
Drop 1 2 3
Drop 15
Options Multiple
End

! Drop B & C2 C3  -> AC1E
get mx0204c.mxs
Drop 1 2 3
Drop 12 14 15
Options Multiple
End

! Drop B & C2 C3 & off diag E  -> AC1Ed
!  Test against ACE
get mx0204c.mxs
Drop 1 2 3
Drop 12 14 15
Drop 17 19 20
Options Multiple
End

! Drop B & C2 C3  -> AC1E
get mx0204c.mxs
Drop 1 2 3
Drop 12 14 15
!  Define this model as saturated
Options Multiple Issat
End

! Drop B & C2 C3 & off diag E  -> AC1Ed
! Test against AC1E
get mx0204c.mxs
Drop 1 2 3
Drop 12 14 15
Drop 17 19 20
Options Multiple
End