!MX script for the Longitudinal QTL model from Chapter 12

#Ngroups 4	
#Define nvar 4		! Number of measurements

G1: Calculation
Data calculation
Matrices ;
 X lower nvar nvar free	! Factor loadings Va
 Y full nvar 1 free		! Factor loadings Vq 
 Z lower nvar nvar free	! Factor loadings Ve
 M full 1 nvar free		! Means for sib 1 & 2 for 4 measurements
 R full 1 nvar free		! Age effect for sib 1 & 2.
 S full 1 nvar free		! Sex effect for sib 1 & 2.
 H full 1 1 fix		! .5 matrix
End matrices ;

Matrix H 0.5000;

Begin algebra;
 A = X*X';			! Va
 Q = Y*Y';			! Vq
 E = Z*Z';			! Ve
 K = A%(A+E+Q);		! The relative contribution of component A
 L = Q%(A+E+Q);		! The relative contribution of component Q	
J = E%(A+E+Q);		! The relative contribution of component E
End algebra ;


! Enter starting values.
! Note with MV analysis these might require some accuracy which will also result in a! gain of optimazation speed.
Start 70 M 1 1 1 M 1 1 2 M 1 1 3 M 1 1 4 
Start -2.00 R 1 1 1 R 1 1 2 R 1 1 3 R 1 1 4 
Start -0.20 S 1 1 1 S 1 1 2 S 1 1 3 S 1 1 4 
Start 5.0 X 1 1 1 X 1 2 1 X 1 2 2 X 1 3 1 X 1 3 2 X 1 3 3 X 1 4 1 X 1 4 2 X 1 4 3 X 1 4 4		
Start 5.0 Y 1 1 1 Y 1 2 1 Y 1 3 1 Y 1 4 1
Start 5.0 Z 1 1 1 Z 1 2 1 Z 1 2 2 Z 1 3 1 Z 1 3 2 Z 1 3 3 Z 1 4 1 Z 1 4 2 Z 1 4 3 Z 1 4 4
End

G2: MZs				
Data NI=43
 Missing=-99.000000
 Rectangular File = MZdata.dat	! Data file MZ twins.

Labels						
 family id1 id2 sex1 twzyg totm1
 dosex1 doage1 domsbp1 domsbp1
 nisex1 niage1 nimsbp1 nimsbp1
 dasex1 daage1 damsbp1 damsbp1
 hasex1 haage1 hamsbp1 hamsbp1
 sex2 totm2 dosex2
 doage2 domsbp2 domsbp2
 nisex2 niage2 nimsbp2 nimsbp2
 dasex2 daage2 damsbp2 damsbp2
 hasex2 haage2 hamsbp2 hamsbp2
 ibd0 ibd1 ibd2			! IBD0, IBD1 and IBD2 are fixed at 0, 0, 1.

Select if twzyg < 4;		! Groups 1 & 3 are MZ twins in our data.
Select if twzyg ^= 2;

Select	
dosex1 doage1 hasex1 haage1 dasex1 daage1 nisex1 niage1 domsbp1 hamsbp1 damsbp1 nimsbp1dosex2 doage2 hasex2 haage2 dasex2 daage2 nisex2 niage2 domsbp2 hamsbp2 damsbp2 nimsbp2
ibd1 ibd2;

Definition_variables
dosex1 doage1 dosex2 doage2
hasex1 haage1 hasex2 haage2
dasex1 daage1 dasex2 daage2
nisex1 niage1 nisex2 niage2
ibd1 ibd2;

Matrices = group 1			
 T full 1 8
 V full 1 8
 N full 1 1
 O full 1 1
End matrices ;

Specify T dosex1 hasex1 dasex1 nisex1 dosex2 hasex2 dasex2 nisex2	
Specify V doage1 haage1 daage1 niage1 doage2 haage2 daage2 niage2
Specify N ibd1
Specify O ibd2

Begin algebra			
! Calculate Pihat 
! Note: There are many ways to incorparate this for MZ twins, as it is a constant.
 P = H@N+O;										
End algebra;											 

Means (M|M)+(R|R).V+(S|S).T;
Covariances 
  A+E+Q | P@Q+A _ 
  P@Q+A | A+E+Q;
End

G3: DZs			
 Data NI=43
 Missing=-99.000000
Rectangular File = DZdata.dat	! Data file DZ twins.

Labels
family id1 id2 sex1 twzyg totm1
dosex1 doage1 domsbp1 domsbp1
nisex1 niage1 nimsbp1 nimsbp1
dasex1 daage1 damsbp1 damsbp1
hasex1 haage1 hamsbp1 hamsbp1
sex2 totm2 dosex2
doage2 domsbp2 domsbp2
nisex2 niage2 nimsbp2 nimsbp2
dasex2 daage2 damsbp2 damsbp2
hasex2 haage2 hamsbp2 hamsbp2
ibd0x ibd1x ibd2x		
! These are the IBD estimations for a  given position X on the genome.

Select if twzyg > 1;		! DZ & sibs twins selection
Select if twzyg ^= 3;

Select				
 dosex1 doage1 hasex1 haage1 dasex1 daage1 nisex1 niage1 domsbp1 hamsbp1  
 damsbp1 nimsbp1 dosex2 doage2 hasex2 haage2 dasex2 daage2 nisex2 niage2   
 domsbp2 hamsbp2 damsbp2 nimsbp2 ibd1x ibd2x;
Definition_variables 
 dosex1 doage1 dosex2 doage2
 hasex1 haage1 hasex2 haage2
 dasex1 daage1 dasex2 daage2
 nisex1 niage1 nisex2 niage2
 ibd1x ibd2x;

Matrices = group 1		
! Matrices for definition variables, are recycled for this group.
 T full 1 8											
 V full 1 8
 N full 1 1
 O full 1 1
End matrices ;

Specify T dosex1 hasex1 dasex1 nisex1 dosex2 hasex2 dasex2 nisex2	
Specify V doage1 haage1 daage1 niage1 doage2 haage2 daage2 niage2
Specify N ibd1x
Specify O ibd2x

Begin algebra ;
 P = H@N+O ;
End algebra ;

Means (M|M)+(R|R).V+(S|S).T;
Covariances 
  A+E+Q | P@Q+H@A _
  P@Q+H@A | A+E+Q		! Variance structure Dz + siblings.
 Option Ndecimals = 4
End

G4: Constraint 
Constraint 

Matrices ;
 L computed 4 4 =L1
 P full 1 4 free
End matrices ;

Begin algebra ;
 V=\d2v(L) ;
End Algebra ;

Constraint V=P;
 Option multiple
End

!Constrain % of variance explained by QTL to be equal across time points
EQ P 4 1 1 P 4 1 2 P 4 1 3 P 4 1 4
End

!Test QTL effect
Drop P 4 1 1 P 4 1 2 P 4 1 3 P 4 1 4
End
! Null model can be tested by dropping Q. However, with a genome scan this  
! model is easier to run separately.
! It saves a lot of time, since the null model does not change for each 
! location.