!
!  QTL via weighted likelihood of 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 NModel=3
 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
 K Full 3 1              ! p(i)'s go here
 Q Full nvar 1 Free      ! qtl effect
 H Full 1 1              ! .5
 I Unit 3 1              ! To spread mean out for all models & sibs
 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 T 2
 Matrix Q 1 
 Start 0 M 1 1 to M 1 nvar

Begin Algebra;
 U = T. E*E' ;
 V = U + T.(Q*Q') ;  ! within pair IBD 0 variance
 W = U + Q*Q' ;      ! within pair IBD 1 variance
 X = U ;             ! within pair IBD 2 variance
 Y = V_
     W_
     X ;
End Algebra;

Means I@M ;
Covariance Y ;
Weights K ;

 Bound .01 100 E 1 1
 Option rs samecov it=250  
! option samecov tells it not to re-invert cov matrix as it doesn't change
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=qtldiff.fit samecov it=250
 Option multiple
End Group;

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