!
!  QTL via pihats and covariance matrices
!
#define nvar 1
#define startqpar 3  ! nvar*(nvar+1) +1
#define endqpar   3  ! nvar*(nvar+2)
#define location 1850
Title QTL model with heterogeneity and weights
Data NGroups=2 NInput=15 
 Rectangular File=../q5v1850.dat
 Labels  Pair Locn w(1) w(2) w(3)
 sib1v1 sib1v2
 sib1v3 sib1v4
 sib1v5
 sib2v1 sib2v2
 sib2v3 sib2v4
 sib2v5
 Select sib1v1 sib2v1
  w(1) w(2) w(3) /

 Definition w(1) w(2) w(3) /

Matrices
 A Lower nvar nvar Free  ! background A (or C)
 E Lower nvar nvar Free  ! random envt
 F Full 1 3              ! 0, .5, 1 to compute pihat
 K Full 3 1              ! p(i)'s go here
 Q Full nvar 1 Free      ! qtl effect
 H Full 1 1              ! .5
 M Full 1 nvar Free      ! for mean
 I Unit 1 2              ! To spread mean out for all models & sibs
End Matrices

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

Begin Algebra;
 P = F*K ;                 ! pihats
 V = A*A' + Q*Q' + E*E' ;  ! variance
 U = h@A*A' + p@Q*Q';      ! IBD pihat cov
 Y = V|U_
     U|V ;                 ! IBD pihat covariance matrix
End Algebra;

 Means I@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=qtlpi.fit it=250
 Option multiple
End Group;

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