! ACE Script for dichotomised variables (twin1/twin2)
! 2 categories: 1 Threshold per variable, one for twin1 and one for twin2
! 0 = bottom 88% 	(non-agressive)
! 1 = top 12%	(agressive)	 
! Thresholds are EQ across twins within each groups and across zygosity group. 

! THRESHOLD MATRIX 1x2 :
! T(1,1) = threshold twin1
! T(1,2) = threshold twin2


#define nvar 2      		! number of variables
#define nthresh 1 		! number of thresholds

G1: Matrices for the Genetic Analysis of the Binary data
Calculation NGroups=4
Begin Matrices;
X FULL 1 1 FREE			! Additive genetic path coefficient
Y FULL 1 1 FREE          		! Common environmental path coefficient
Z FULL 1 1 FREE			! Unique environmental path coefficient
End matrices;
Begin Algebra;
A=X*X' ;				!Additive genetic variance (path X squared)
C=Y*Y' ;				!Common Environmental variance(path Y squared)
E=Z*Z' ;				!Unique Environmental variance(path Z squared)
End Algebra;
start .6 all				!starting value for X, Y, Z
Interval @95 A 1 1 C 1 1 E 1 1
End 

G2: Data and model for MZ pairs
DAta NInput_vars=3  
Missing=.
Ordinal File=bin.dat
Labels 
zyg bin1 bin2
SELECT IF zyg = 1
SELECT bin1 bin2 /

Matrices = group 1
T FULL nthresh nvar FREE		!thresh tw1, thresh tw2
L LOW nthresh nthresh
end matrices;

Value 1 L 1 1 to L maxthresh maxthresh 	! initialize L

COVARIANCE
 ( A + C + E | A + C _
      A + C     | A + C + E ) /		!Predicted Correlation matrix for MZ pairs

Thresholds  L*T ;			!to ensure t1>t2>t3 etc.......


Matrix T   
1.2 1.2 			!starting values for the thresholds
EQ T 1 1 1 T 1 1 2		!constrain TH to be equal across twins (twin1/twin2) 

Option rs
End
 
G3: Data and model for DZ pairs
DAta NInput_vars=3  
Missing=.
Ordinal File=bin.dat
Labels 
zyg bin1 bin2
SELECT IF zyg = 2
SELECT bin1 bin2 /

Matrices = group 1
T FULL maxthresh nvar =T2		!thresh tw1, thresh tw2
L LOW maxthresh maxthresh =L2
H FULL 1 1			!.5
end matrices;

COVARIANCE
 ( A + C + E | H@A + C _
   H@A + C   | A + C + E ) /
				!Predicted Correlation matrix for DZ pairs

Thresholds  L*T ;			!to ensure t1>t2>t3 etc.......

MA H .5
Option func=1.E-10  		!function precision is less than usual
Option rs
End 

G4: CONSTRAIN VARIANCES OF OBSERVED VARIABLES TO 1
CONSTRAINT
Matrices = Group 1
I UNIT 1 1
CO A+C+E= I /
Options RS ISSAT Multiple
End Group;

Save ace.mxs        !this saves the full model so you can test reduced models

Get ace.mxs
drop 1          !this drops the A term - tests CE model
end group;

get ace.mxs
drop 2          !this drops the C term - tests AE model
end group;

get ace.mxs
drop 1 2           !this drops A and C from ACE - tests E model
end group;