! Boulder 2003
! multivariate Mx script (cholesky decomposition) for the analysis of three continuous phenotypes
! using raw data input file
! MZM, DZM, MZF, DZF, DOSMF, DOSFM twins
! Evaluated models: 
! ACE males, ACE females, DOS A-correlation = DZ A-correlation
! Mx homepage: http://views.vcu.edu/mx

#define nvar 3				! For every instance of the word nvar in the script, 
							! Mx reads 1. This corresponds to one phenotype to be analysed


G1: calculation group
 Data Calc NGroups=8

  Begin matrices;
     X Lower nvar nvar Free		! MALES additive genetic path coefficient (standard deviation)
     Y Lower nvar nvar FIX	! Free	! assume no common environment ! MALES common environmental path coefficient (standard deviation)
     Z Lower nvar nvar Free		! MALES unique environmental path coefficient (standard deviation) 

     S Lower nvar nvar free			 ! FEMALES additive genetic path coefficient (standard deviation)
     T Lower nvar nvar FIX ! free	! assume no common environment ! FEMALES common environmental path coefficient (standard deviation)
     U Lower nvar nvar  free		! FEMALES unique environmental path coefficient (standard deviation) 

     H Full 1 1	fix					! .5 for dz cov A
     F Full 1 1	fix		!free		! .5 dos cov A
       End matrices;

 Matrix H .5				! Matrix H is fixed and contains .5
 Matrix F .5				! Matrix H is fixed and contains .5

  Begin algebra;
     A= X*X' ;			  	! MALES additive genetic variance
     C= Y*Y' ;   			! MALES common environmental variance
     E= Z*Z' ;				! MALES unique environmental variance

     K= S*S' ;			  	! FEMALES additive genetic variance
     L= T*T' ;   			! FEMALES common environmental variance
     M= U*U' ;				! FEMALES unique environmental variance

    

    End algebra;

st 0.1 all
 st 1.5 X 1 1 1 X 1 2 1 X 1 3 1			! provide starting values
! st 1.5 Y 1 1 1 Y 1 2 1 Y 1 3 1
 st 1.5 Z 1 1 1 Z 1 2 2 Z 1 3 3
 st 1.5 S 1 1 1 S 1 2 1 S 1 3 1			
! st 1.5 T 1 1 1 T 1 2 1 T 1 3 1
 st 1.5 U 1 1 1 U 1 2 2 U 1 3 3


! st 0.5 F 1 1 1
! Bound -0.5 0.5 F 1 1 1

end

G2: MZM twins, datagroup
Data NInput_vars= 14				! nr of inputvars per family
 Missing=-1.00						! missing values = -1.00
 REctangular file=all.dat				! read in raw datafile, must be in same directory

 LABELS     							! tell Mx what is found in the datafile
   nih_id country zygot6 rear ageA sexA ldllnA apoblnA apoelnA ageB sexB ldllnB apoblnB apoelnB

 Select if country = 1 ; 				! select Dutch twins, younger sample (Dorret)
 Select if  zygot6 = 1 ; 				! select MZM's
 Select   ldllnA apoblnA apoelnA  ldllnB apoblnB apoelnB ;		! select nvar phenotypes per twin member
  Matrices = group 1 
G Full 1 6  free				!  mean phenotypes
end matrices
  Means G ;						! model for means 

     Covariances						! model for MZM variance/covariances
	A+C+E	| A+C	 _
	A+C	| A+C+E	;
 st 15 G 1 1 - G 1 nvar
Options RSiduals
End

G3: DZM twins, datagroup
Data NInput_vars= 14				! nr of inputvars per family
 Missing=-1.00						! missing values = -1.00
 REctangular file=all.dat				! read in raw datafile, must be in same directory

 LABELS     							! tell Mx what is found in the datafile
   nih_id country zygot6 rear ageA sexA ldllnA apoblnA apoelnA ageB sexB ldllnB apoblnB apoelnB

 Select if country = 1 ; 				! select Dutch twins, younger sample (Dorret)
 Select if  zygot6 = 2 ; 				! select DZM's
 Select   ldllnA apoblnA apoelnA  ldllnB apoblnB apoelnB ;				! select nvar phenotypes per twin member
 Matrices = group 1 
    G Full 1 6  free				!  mean phenotypes
 end matrices;

    Means G ;			! model for means
    
    Covariances						! model for DZM variance/covariances
	A+C+E	| H@A+C	 _
	H@A+C	| A+C+E	;
 st 15 G 1 1 - G 1 nvar
 Options RSiduals
 End

G4: MZF twins, datagroup
 Data NInput_vars= 14				! nr of inputvars per family
 Missing=-1.00						! missing values = -1.00
 REctangular file=all.dat				! read in raw datafile, must be in same directory

 LABELS     							! tell Mx what is found in the datafile
   nih_id country zygot6 rear ageA sexA ldllnA apoblnA apoelnA ageB sexB ldllnB apoblnB apoelnB

 Select if country = 1 ; 				! select Dutch twins, younger sample (Dorret)
 Select if  zygot6 = 3 ; 				! select MZF's
 Select   ldllnA apoblnA apoelnA  ldllnB apoblnB apoelnB ;			! select nvar phenotypes per twin member
  Matrices = group 1 
    G Full 1 6  free				!  mean phenotypes
 end matrices;

    Means G ;			! model for means, assuming mean twin 1 = mean twin 2, and DZ mean=MZ mean
    
    Covariances						! model for MZF variance/covariances
	K+L+M	| K+L	 _
	K+L	| K+L+M	;
 st 15 G 1 1 - G 1 nvar
 Options RSiduals
 
End

G5: DZF twins, datagroup
 Data NInput_vars= 14				! nr of inputvars per family
 Missing=-1.00						! missing values = -1.00
 REctangular file=all.dat				! read in raw datafile, must be in same directory

 LABELS     							! tell Mx what is found in the datafile
   nih_id country zygot6 rear ageA sexA ldllnA apoblnA apoelnA ageB sexB ldllnB apoblnB apoelnB

 Select if country = 1 ; 				! select Dutch twins, younger sample (Dorret)
 Select if  zygot6 = 4 ; 				! select DZF's
  Select   ldllnA apoblnA apoelnA  ldllnB apoblnB apoelnB ;				! select nvar phenotypes per twin member,
 Matrices = group 1 
   G Full 1 6  free				!  mean phenotypes
 end matrices;
    Means G ;					! model for means,
    
    Covariances						! model for DZF variance/covariances
	K+L+M	| H@K+L	 _
	H@K+L	| K+L+M	;
st 15 G 1 1 - G 1 nvar
 Options RSiduals
End

G6: DOS-MF twins, datagroup
Data NInput_vars= 14				! nr of inputvars per family
 Missing=-1.00						! missing values = -1.00
 REctangular file=all.dat				! read in raw datafile, must be in same directory

 LABELS     							! tell Mx what is found in the datafile
   nih_id country zygot6 rear ageA sexA ldllnA apoblnA apoelnA ageB sexB ldllnB apoblnB apoelnB

 Select if country = 1 ; 				! select Dutch twins, younger sample (Dorret)
 Select if  zygot6 = 5 ; 				! select DOSmf's
 Select  ldllnA apoblnA apoelnA  ldllnB apoblnB apoelnB ;				! select nvar phenotypes per twin member 
 Matrices = group 1 
     G Full 1 6  free				!  mean phenotypes
 end matrices;

    Means G ;			! model for means, assuming mean twin 1 = mean twin 2, and DZ mean=MZ mean
    
    Covariances						! model for DOSMF variance/covariances
	A+C+E	| F@(X*S')+(Y*T')	 _
	F@(S*X')+(T*Y')	| K+L+M	;
st 15 G 1 1 - G 1 nvar
 Options RSiduals
 
End


G7: DOSFM twins, datagroup
 Data NInput_vars= 14				! nr of inputvars per family
 Missing=-1.00						! missing values = -1.00
 REctangular file=all.dat				! read in raw datafile, must be in same directory

 LABELS     							! tell Mx what is found in the datafile
   nih_id country zygot6 rear ageA sexA ldllnA apoblnA apoelnA ageB sexB ldllnB apoblnB apoelnB

 Select if country = 1 ; 				! select Dutch twins, younger sample (Dorret)
 Select if  zygot6 = 6 ; 				! select DOSfm's
  Select   ldllnA apoblnA apoelnA ldllnB apoblnB apoelnB ;		
		! select nvar phenotypes per twin member
  Matrices = group 1 
G Full 1 6  free				!  mean phenotypes
    end matrices;


    Means G ;			! model for means, assuming mean twin 1 = mean twin 2, and DZ mean=MZ mean
    
    Covariances						! model for DZ variance/covariances
	K+L+M	| F@(S*X')+(T*Y')	 _
	F@(X*S')+(Y*T')	| A+C+E	;
st 15 G 1 1 - G 1 nvar
 Options RSiduals
 !Option Multiple issat  				! allows comparison of submodels with present model (ACE KLM)
 !Interval @95 V 1 1 1 V 1 1 2 V 1 1 3 V 1 2 1 V 1 2 2 V 1 2 3 P 1 1 1 P 1 2 1 G 1 1 1 F 1 1 1 

End

G8: Calculate genetic and environmental correlations and standardized variance-covariance components
 Calculation
 Begin Matrices
A comp nvar nvar = A1
C comp nvar nvar = C1
E comp nvar nvar = E1
K comp nvar nvar = K1
L comp nvar nvar = L1
M comp nvar nvar = M1   
I iden nvar nvar
 End Matrices;
 Begin Algebra;


!men: calculation of correlations using more elaborated approach (no error messages in output) 
  R=\sqrt(I.A)~*A*\sqrt(I.A)~ ;	! genetic correlations
  S=\sqrt(I.C)~*C*\sqrt(I.C)~ ;	! shared environmental correlations
  T=\sqrt(I.E)~*E*\sqrt(I.E)~ ;	! unique environmental correlations

! women: calculation of correlations using \stnd function (mx will give error messages in output, but this is okay)
  J = \stnd(K) ;			! standardize additive genetic variance/covariance matrix females to obtain genetic correlations between variables
  N = \stnd(L) ;			! standardize common environmental variance/covariance matrix females to obtain com. env. correlations between variables
  O = \stnd(M) ;			! standardize unique environmental variance/covariance matrix females to obtain unique env. correlations between variables
  
!standardized variance and covariance components 
V = 
     A%(A+C+E) | C%(A+C+E) | E%(A+C+E) _	    !MALES FIRST ROWs standardized variance/covariance components
     K%(K+L+M) | L%(K+L+M) | M%(K+L+M) ;	    !FEMALES SECOND ROWs standardized variance/covariance components

!all correlations in 1 matrix
  W = R | S | T _
	   J | N | O ;
 End Algebra;
  Label Col W r_a r_a r_a r_c r_c r_c r_e r_e r_e
 Label Row W    ldlmen apobmen apoemen ldlwomen apobwomen apoewomen
Label Col V h2_ldl h2_ApoB  h2_ApoE  c2_ldl c2_ApoB  c2_ApoE e2_ldl e2_ApoB  e2_ApoE    
 Label Row V    men men men women women women
 Option Iterations=5000 
 Option NDecimals=4 ! THard=-2
 Option Multiple
 Option Issat		! compare fit of submodels to full model
End

! FIX DOS COVARIANCE FOR A AT 0.5, WHICH EQUALS THE DZ COVARIANCE
! THIS PROVIDES A TEST OF WHETHER THE SAME GENES ARE EXPRESSED IN MALES AND FEMALES
!DROP @.5 F 1 1 1 
!end