!
!  QTL via weighted likelihood of three 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 NModel=3
 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
 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 3 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 
 Start 0 M 1 1 to M 1 nvar

Begin Algebra;
 V = A*A' + Q*Q' + E*E' ;  ! variance
 U = h@A*A' ;              ! IBD 0 cov (=non-qtl genetic cov)
 X = U + h@Q*Q' ;          ! IBD 1 cov
 W = U + Q*Q' ;            ! IBD 2 cov
 Y = V|U_
     U|V_      ! IBD 0 matrix
     V|X_
     X|V_      ! IBD 1 matrix
     V|W_
     W|V ;     ! IBD 2 matrix
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=qtl.fit samecov it=250
 Option multiple
End Group;

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