!MX script for a univariate mixture distribution QTL model from Chapter 10

! Estimate QTL, additive genetic and non-shared environmental components 
! AEQ model 
! mixture distribution approach

#loop $nr 0 260 2
#define nvar 1
#define ndef 2
#ngroups 1

DZ / sib
Data NInput=10 NModel=3 maxrec=4000
 Missing =-1.00
 Rectangular File=ex.dat

 Labels
  fam age_1 sex_1 trait_1 age_2 sex_2 trait_2
  ibd0_0  ibd1_0  ibd2_0
  ibd0_2  ibd1_2  ibd2_2
  ibd0_4  ibd1_4  ibd2_4
  ibd0_6  ibd1_6  ibd2_6
  ibd0_8  ibd1_8  ibd2_8
  ibd0_10 ibd1_10 ibd2_10
  ! there are 400 variables in the entire list we omit some here for didactic    
  ! purposes
  ibd0_250  ibd1_250  ibd2_250
  ibd0_252  ibd1_252  ibd2_252
  ibd0_254  ibd1_254  ibd2_254
  ibd0_256  ibd1_256  ibd2_256
  ibd0_258  ibd1_258  ibd2_258
  ibd0_260  ibd1_260  ibd2_260

Select  age_1 sex_1 trait_1 age_2 sex_2 trait_2 ibd0_$nr ibd1_$nr ibd2_$nr ;
Definition age_1 sex_1 age_2 sex_2 ibd0_$nr ibd1_$nr ibd2_$nr ;

Matrices ;
 X Lower nvar nvar free        ! Additive A
 Z Lower nvar nvar free        ! Unshared E
 V Full  nvar nvar free        ! Qtl Q
 M Full  1    nvar free        ! means
 B Full  ndef nvar Free        ! estimated effects of the covariates
 N Full  1    ndef  fix        ! observed covariates for twin 1
 O Full  1    ndef  fix        ! observed covariates for twin 2
 H Full  1    1     fix        ! will contain .5
 K Full  3    1     fix        ! IBD probabilities (from Merlin)
 U Unit  3    1     fix
End Matrices ;

!now contains the first definition variable
Specify K ibd0_$nr ibd1_$nr ibd2_$nr ;   
Specify N age_1 sex_1
Specify O age_2 sex_2
Matrix H 0.5                    !now contains 0.5

st 9 X 1 1 1
st 1 V 1 1 1
st 3 Z 1 1 1
st 100 M 1 1 1
st 2 B 1 1 1 B 1 2 1

Begin Algebra ;
 A = X*X';
 E = Z*Z';
 Q = V*V';
 S = A%(A+E+Q) |E%(A+E+Q)| Q%(A+E+Q) ;
 T = A|Q|E;
End Algebra ;

Means   U@(M+N*B) | U@(M+O*B);
Covariances  
 A+E+Q   | H@A     _
 H@A     | A+E+Q   _       ! IBD 0 Covariance matrix
 A+E+Q   | H@A+H@Q _
 H@A+H@Q | A+E+Q   _       ! IBD 1 Covariance matrix
 A+E+Q   | H@A+Q   _
 H@A+Q   | A+E+Q ;         ! IBD 2 Covariance matrix

Weights K;                 ! IBD probabilities
 Option  nd=4               ! request 4 decimal places in output
 OPtion RS                  ! request residuals
 Option Sub = 8137.729, 969
Exit
#end loop