!
!  QTL via pi-hats using within-pair differences 
!
#define nvar 1
#define startqpar 2  ! nvar*(nvar+1)/2 +1
#define endqpar   2  ! nvar*(nvar+3)/2
#define location 1850
Title QTL model with heterogeneity and weights
Data NGroups=2 NInput=6
 Rectangular File=diff.dat
 Labels  Pair locus p(0) p(1) p(2) diff
 Select   Diff p(0) p(1) p(2) /
 Definition p(0) p(1) p(2) /

Matrices
 E Lower nvar nvar Free  ! random envt
 F Full 1 3              ! 0 .5 1
 K Full 3 1              ! p(i)'s go here
 Q Full nvar 1 Free      ! qtl effect
 H Full 1 1              ! .5
 J Unit 1 1              ! 1.0
 M Full 1 nvar Free      ! for mean
 T Full 1 1              ! 2.0
End Matrices

 Matrix H .5
 Specify K -1 -2 -3;
 Matrix E 1 
 Matrix F 0 .5 1
 Matrix T 2
 Matrix Q 1 
 Start 0 M 1 1 to M 1 nvar

Begin Algebra;
 P = F*K ;                  ! individual pihats
 Y = T.(E*E' + (J-P)@(Q*Q')) ;  ! within pair IBD 0 variance
End Algebra;

Means M ;
Covariance Y ;

 Bound .01 100 E 1 1
 Option rs it=250  
End

Group 2 - calculate diff in fit
Calculation
Matrices
 A full 1 1 =%f1
 B full 1 1
 L full 1 1
End Matrices

 Matrix B 0
 Matrix L location

Begin Algebra;
 C = L | A - B;
End Algebra;

 Label row C fitdiff

 Option mxa=qtldiffpi.fit it=250
 Option multiple
End Group;

! Eliminate effects of QTL
 Matrix 2 B File=qtldiffpi.fit
 Drop @.0 startqpar
Exit