!MX script for a univariate pi-hat QTL model suitable for use with larger sibships from Chapter 10

! Estimate QTL, additive genetic and non-shared environmental components 
! AEQ model 
! pi-hat approach, multiple siblings

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

DZ / sib
Data NInput=32
 Missing =-1.00
 Rectangular File=chr2_$nr.rec
 
 Labels
 fam id1 id2 id3 id4 id5 p12 p13 p14 p15 p23 p24 p25 p34 p35 p45 twin_add
 age_1 age_2 age_3 age_4 age_5
 sex_1 sex_2 sex_3 sex_4 sex_5
 trait_1 trait_2 trait_3 trait_4 trait_5

Select
 p12 p13 p14 p15 p23 p24 p25 p34 p35 p45 twin_add
 age_1 age_2 age_3 age_4 age_5
 sex_1 sex_2 sex_3 sex_4 sex_5
 trait_1 trait_2 trait_3 trait_4 trait_5 ;
Definiton_variables
 p12 p13 p14 p15 p23 p24 p25 p34 p35 p45 twin_add
 age_1 age_2 age_3 age_4 age_5
 sex_1 sex_2 sex_3 sex_4 sex_5 ;

Matrices ;
 X Lower nvar nvar Free    ! genetic structure
 Z Lower nvar nvar Free    ! specific environmental structure
 V Full  nvar 1 Free       ! QTL effect
 M Full  1 nvar Free       ! means
 I Unit  nsib 1
 B Full  ndef nvar Free    ! estimated effects of the covariates
 N Full  nsib ndef         ! observed covariates for all sibs
 F stand nsib nsib         ! A coefficents
 P stand nsib nsib         ! Q coefficents - pihats
 R ident nsib nsib         ! E coefficents
End Matrices ;

Specify N age_1 sex_1 age_2 sex_2 age_3 sex_3 age_4 sex_4 age_5 sex_5
Specify P p12 p13 p23 p14 p24 p34 p15 p25 p35 p45 ;
Specify F
twin_add
0 0
0 0 0
0 0 0 0
Value .5 F 3 1 - F 5 4 ;

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 Albegra ;

Means (I@M)+N*B ;
Covariances
 F@A + R@E + P@Q ;

 Option nd=4                   ! request 4 decimal places in output
 OPtion RS                      ! request residuals
 Exit
#end loop