!
! Mx script for causal conditional model
!

#ngroup 6

Group 1 Compute MZ Correlations
  Calculation
 Begin Matrices;
  A Di 2 2 Free
  C Di 2 2 Free
  E Di 2 2 Free
  B Fu 2 2
  I Id 2 2 
 End Matrices; 

  Specify B  ! causal parameter from initiation to progression
   0 0
   7 0

! starting values
  Matrix A .7 .7
  Matrix C .5 .5
  Matrix E .5 .5
  Matrix B 0 0 .4 0

! parameter bounds
 Bound 0 1 A 1 1 A 2 2 C 1 1 C 2 2
 
 Begin algebra;
  X = A*A' ;
  Y = C*C' ;
  Z = E*E' ;
  R =       (I@((I-B)~)) & ( X+Y+Z| X+Y _
                               X+Y| X+Y+Z);
 End algebra;

End group

Group 2 DZ Correlation matrix
 Calculation
 Begin Matrices=Group 1;
  H Fu 1 1 ! .5
 End Matrices;

 Matrix H .5

 Begin algebra;
  R =            (I@((I-B)~)) &( X+Y+Z| h@X+Y _
                                h@X+Y| X+Y+Z);
 End algebra;
End

Fit model to MZ data with user-defined fit function (ML)
 Data Ni=1 No=1
 Begin Matrices;
  d full 1 1 ! two
  i zero 1 1
  n full 1 1 ! scalar 2.0
  o full 9 1
  r computed =R1 ! correlation matrix A1B1A2B2
  t full 1 4 ! thresholds abab
  w zero 1 4 ! means
  z unit 1 1
 End matrices;

 matrix d 2
 matrix n 2
 matrix o ! non, ex, current cell frequencies
  214 53 5 
  55 117 17 
  1 20 18

! mnor function takes matrices with 4 more rows than columns.
! first n (=4) rows are correlation matrix
! row n+1 is mean vector
! row n+2 is upper thresholds
! row n+3 is lower threshold
! row n+4 is indicator, 0 = integrate from -infinity to upper threshold
!                       1 = integrate from lower threshold to +infinity
!
 Begin algebra ;
  e = -o.\ln
       (\mnor(r_w_t_t_(i|i|i|i)) +
        \mnor(r_w_t_t_(i|i|i|z)) +
        \mnor(r_w_t_t_(i|z|i|i)) +
        \mnor(r_w_t_t_(i|z|i|z)) _

        \mnor(r_w_t_t_(i|i|z|i)) +
        \mnor(r_w_t_t_(i|z|z|i)) _

        \mnor(r_w_t_t_(i|i|z|z)) +
        \mnor(r_w_t_t_(i|z|z|z)) _

        \mnor(r_w_t_t_(z|i|i|i)) +
        \mnor(r_w_t_t_(z|i|i|z)) _

        \mnor(r_w_t_t_(z|i|z|i)) _

        \mnor(r_w_t_t_(z|i|z|z)) _

        \mnor(r_w_t_t_(z|z|i|i)) +
        \mnor(r_w_t_t_(z|z|i|z)) _

        \mnor(r_w_t_t_(z|z|z|i)) _

        \mnor(r_w_t_t_(z|z|z|z)) );
 End algebra ;

 Compute d. \sum(e) ;
 Option user rs
End

Fit model to DZ data with user-defined fit function (ML)
 Data Ni=1 No=1
 Begin Matrices = Group 3;
  ! re-declare o and r as they are different for DZ's
  o full 9 1
  r comp =R2 ! correlation matrix A1B1A2B2
 End matrices;

 Specify t 10 11 10 11 ! equate thresholds for twin 1 & twin 2 (and MZ/DZ)
 Matrix t .2 .0 .2 .0
 Matrix o 
  100 52 3 
  45 82 16 
  6 18 4
 Bound -2 3 10 11

 Begin algebra ;
  e = -(o). \ln 
       (\mnor(r_w_t_t_(i|i|i|i)) +
        \mnor(r_w_t_t_(i|i|i|z)) +
        \mnor(r_w_t_t_(i|z|i|i)) +
        \mnor(r_w_t_t_(i|z|i|z)) _

        \mnor(r_w_t_t_(i|i|z|i)) +
        \mnor(r_w_t_t_(i|z|z|i)) _

        \mnor(r_w_t_t_(i|i|z|z)) +
        \mnor(r_w_t_t_(i|z|z|z)) _

        \mnor(r_w_t_t_(z|i|i|i)) +
        \mnor(r_w_t_t_(z|i|i|z)) _

        \mnor(r_w_t_t_(z|i|z|i)) _

        \mnor(r_w_t_t_(z|i|z|z)) _

        \mnor(r_w_t_t_(z|z|i|i)) +
        \mnor(r_w_t_t_(z|z|i|z)) _

        \mnor(r_w_t_t_(z|z|z|i)) _

        \mnor(r_w_t_t_(z|z|z|z)) );
 End algebra ;

 Compute d. \sum(e);
 Option user 
End

Group 5 constrain variances to 1
 Constraint NI=1
 Begin Matrices = Group 1;
  U Unit 1 2
  z izero 4 2
 End Matrices;
 Constraint U = \d2v(r)*z;
End

Group 6 - standardize estimates
 Calculation
 Begin Matrices = Group 1;
  I Id 2 2
  J iden 4 4
 End Matrices;

 Begin algebra;
  K = I@((I-B)~);
  L =  \v2d(\sqrt(\d2v((I@((I-B)~))*( X+Y+Z| X+Y _
                                X+Y| X+Y+Z)* (I@((I-B)~)'))));
  M = L~*K*L;
  R = \d2v(X+Y+Z);
  S = ((\d2v(X))%R)_((\d2v(Y))%R)_((\d2v(Z))%R);
 End algebra;

Labels row S
 A_A A_B C_A C_B E_A E_B
Labels row L
 MZT1Abeta MZT1Bbeta
 MZT1Abeta MZT1Bbeta

 option func=1.e-10    ! function precision for optimization
 option df=18          ! adjust df
 option nd=4           ! 4 decimal places
 option eps=.00000001  ! integration precision for mnor
 option th=-2          ! retry optimization from final point twice
end