! Bivariate Genetic Cholesky Model CCC
! Simulated ordinal data 

#ngroups 4

#define nthresh1 1
#define nthresh2 1
#define nvar 4

Group 1: set up model
 Calculation
 Begin Matrices;
  B Full nvar nvar Free    ! causal pathways
  J Iden 2 2
  K Iden 4 4
  X diag nvar nvar Free    ! genetic structure
  Y diag nvar nvar Free    ! common environmental structure
  Z diag nvar nvar Free    ! specific environmental structure
  W Lower nvar nvar        ! dominance structure (set to zero)
  T Full nthresh2 8 Free
 End Matrices;

 Bound .0 2 X 1 1 X 2 2 X 3 3 X 4 4
 Bound .0 2 Y 1 1 Y 2 2 Y 3 3 Y 4 4
 Bound .1 2 Z 1 1 Z 2 2 Z 3 3 Z 4 4

 Matrix T 
  0 0 0 0 0 0 0 0
 Bound -3 3 T 1 1 - T 1 8

 Specify B
  0   0 104 0 
  101 0 105 0
  102 0 0   0
  103 0 106 0
 Bound -.99 .99 B 1 1 to B 4 4

 Labels Col B IS DS IC DC
 Labels Row B IS DS IC DC

 Matrix X .8 .7071  .6 .5
 Matrix Z .6 .5 .8 .7071
 Matrix B
  0 0 .3 0 
  .5 0 .3 0 
  .3 0 0 0
  .3 0 .5 0

 Begin Algebra;
  A= X*X' ;
  C= Y*Y' ;
  E= Z*Z' ;
  D= W*W' ;
  F= (J@(K-B))~  ;
 End Algebra;
End

Group 2: Fit model to MZ twin pairs
 Data Ninput=10
 Labels  zyg agea  nicusea nicpc2a canusea canabua
               nicuseb nicpc2b canuseb canabub
 Ordinal file=ffinits.rec
 select if zyg = 1
 select nicusea nicpc2a canusea canabua nicuseb nicpc2b canuseb canabub ;

 Begin Matrices= Group 1;

 Covariances F&(A+C+E+D | A+C+D _
                A+C+D   | A+C+E+D) /
 Thresholds T ;
 Option RSidual
End 

Group 3: Fit model to DZ twin pairs
 Data Ninput=10
 Labels  zyg agea  nicusea nicpc2a canusea canabua
               nicuseb nicpc2b canuseb canabub
 Ordinal file=ffinits.rec
 Select if zyg = 2
 Select nicusea nicpc2a canusea canabua nicuseb nicpc2b canuseb canabub ;

 Begin Matrices= Group 1;
  H Full 1 1
  Q Full 1 1
 End Matrices;
 Matrix H .5
 Matrix Q .25
 Covariances F&(A+C+E+D   | H@A+C+Q@D _
                H@A+C+Q@D | A+C+E+D) /
 Thresholds T ;
 Option RSidual
 Options NDecimals=4
 Option func=1.e-8
End

G4: Constrain variances
 Constraint NI=1
 Begin Matrices;
  U unit 1 4
  E symm 8 8 = %e2
  Z iz 8 4
 End Matrices;

 Constraint U = \d2v(E) * Z ;
 Option
End