! Genetic Cholesky Model -- RAW DATA

#define nvar=2	          !number of examined phenotypes/ages
#define nsib 2		! number of sibs (twin = 2)

G1: Model parameters 
Data Calc NGroups=4
Begin Matrices;
X Lower nvar nvar Free    ! genetic structure
Y Lower nvar nvar  Free   ! common environmental structure
Z Lower nvar nvar  Free   ! specific environmental structure
!B Full 1 nvar  Free            ! means refers to the number of variables selected
P Symm 4 4		! interaction parameters
I Iden	4 4
End Matrices;
Specify P	
0
0 0
91 0 0
0 92 0 0

 st .5  X 1 1 1 Y 1 1 1 Z 1 1 1 		! provide starting values sqrt(var/3) 
 st .5  X 1 2 2 Y 1 2 2 Z 1 2 2

! st .5  X 1 2 2 Y 1 2 2 Z 1 2 2
! st .1 B 1 1 1  					! based on means from SAS
! st .5 B 1 1 2

Begin Algebra;
A= X*X' ;   
C= Y*Y' ;  
E= Z*Z' ;   
W= (I - P)~;
End Algebra;


End


G2: MZ TWINS
Data NInput_vars=62						! number of variables in data file
 Missing=-999999							! missing value = -1.00
 Rectangular File=P_factor.dat					! read in raw data file
 LABELS
	IFAM
zygosity zygDich zyg sex_1 sex_2
 P3_1 INT3_1 EXT3_1 P5_1
INT5_1 EXT5_1 P6_1 INT6_1 EXT6_1
P9_1 INT9_1 EXT9_1 P3_2 INT3_2 EXT3_2 P5_2 INT5_2
EXT5_2 P6_2 INT6_2 EXT6_2 P9_2 INT9_2
EXT9_2 P3R_1 P5R_1 P6R_1 P9R_1 P3R_2 P5R_2
P6R_2 P9R_2 P3res_1 P5res_1 P6res_1 P9res_1 EXT3res_1
EXT5res_1 EXT6res_1 EXT9res_1 INT3res_1 INT5res_1 INT6res_1 INT9res_1
P3res_2 P5res_2 P6res_2 P9res_2 EXT3res_2 EXT5res_2
EXT6res_2 EXT9res_2 INT3res_2 INT5res_2 INT6res_2 INT9res_2 		! Tell mx what is found in the datafile

Select if zygDich = 1;   						! select MZ twins
Select  EXT3_1 EXT5_1 EXT3_2 EXT5_2; 				! select a phenotype 
Matrices= Group 1
B Full 4 1 Free
Means B/
Covariances W &( A+C+E | A+C _
            A+C   | A+C+E) ;
Option RSidual
End 


G3:  DZ TWIN PAIRS
Data NInput_vars=62						! number of variables in data file
 Missing=-999999							! missing value = -1.00
 Rectangular File=P_factor.dat					! read in raw data file
 LABELS
	IFAM
zygosity zygDich zyg sex_1 sex_2
 P3_1 INT3_1 EXT3_1 P5_1
INT5_1 EXT5_1 P6_1 INT6_1 EXT6_1
P9_1 INT9_1 EXT9_1 P3_2 INT3_2 EXT3_2 P5_2 INT5_2
EXT5_2 P6_2 INT6_2 EXT6_2 P9_2 INT9_2
EXT9_2 P3R_1 P5R_1 P6R_1 P9R_1 P3R_2 P5R_2
P6R_2 P9R_2 P3res_1 P5res_1 P6res_1 P9res_1 EXT3res_1
EXT5res_1 EXT6res_1 EXT9res_1 INT3res_1 INT5res_1 INT6res_1 INT9res_1
P3res_2 P5res_2 P6res_2 P9res_2 EXT3res_2 EXT5res_2
EXT6res_2 EXT9res_2 INT3res_2 INT5res_2 INT6res_2 INT9res_2 		! Tell mx what is found in the datafile


Select if zygDich = 0;   						! select DZ twins
Select  EXT3_1 EXT5_1 EXT3_2 EXT5_2; 				! select a phenotype 

Matrices= Group 1
H Full 1 1
B Full 4 1 Free
Means B/
Covariances  W &( A+C+E   | H@A+C _
            H@A+C | A+C+E) ;
Matrix H .5
Option RSidual
Options NDecimals=4
End


G4:  Calculation of standardized solution 
Calculation
Matrices= Group 1
I Iden nvar nvar
Begin Algebra;
V= A+C+E;
G=  \sqrt(I.V)~*X; !!standardized genetic parameters
K=  \sqrt(I.V)~*Y; !standardized c parameters
L=  \sqrt(I.V)~*Z; !standardized e parameters
D= \SQRT(\V2D(\D2V(A+C+E))); !PHENOTYPIC STAND DEV MATRIX !
M= D~*A*D~; !h2 on diag 
N= D~*C*D~; !c2 on diag 
O= D~*E*D~; !e2 on diag 
Q = \sqrt( I . A )~ * A * \sqrt( I . A )~;  ! Genetic correlations same as ! Q= \STND(A); !rg
R = \sqrt( I . C )~ * C * \sqrt( I . C )~; ! Shared environment correlations same as !R= \STND(C); !rc
S = \sqrt( I . E )~ * E * \sqrt( I . E )~; ! Nonshared environment correlations same as !S= \STND(E); !re
U = A%V _ C%V _ E%V;		! bivariate heritability and enviromentalities
End Algebra;
Label Col G a_stand a_stand
Label Col K c_stand c_stand
Label Col L e_stand e_stand
Label Col M H2_1 H2_2
Label Col N C2_1 C2_2
Label Col O E2_1 E2_2
INTERVALS M 4 1 1 M 4 2 2
INTERVALS N 4 1 1 N 4 2 2 
INTERVALS O 4 1 1 O 4 2 2 
INTERVALS G 4 1 1 G 4 2 1 G 4 2 2 
INTERVALS K 4 1 1 K 4 2 1  K 4 2 2 
INTERVALS L 4 1 1 L4 2 1  L 4 2 2 
INTERVALS Q 4 2 1
INTERVALS R 4 2 1
INTERVALS S 4 2 1
Intervals P 4 3 1 P 4 4 2
option ndecimals=2
Option sat=9654.101,4615 ! sat model
Options Multiple ISSub
End

save ace.mxs

! drop SI
get ace.mxs
drop 91 92 
end

! drop c on time 1
!get ace.mxs
!drop 4 5
!end

!drop c on time 2
!get ace.mxs
!drop 5 6
!end

!drop a time 1
!get ace.mxs
!drop  1 2
!end

!drop a time 2
!get ace.mxs
!drop 2 3
!end

!drop all a
!get ace.mxs
!drop 1 2 3
!end

! drop all c
get ace.mxs
drop 4 5 6
end

! drop both a and c
!get ace.mxs
!drop 1 2 3 4 5 6
!end