! Sib pair ACE data

#define nvar 1
#define nvarx2 2
#define nallele 2
#ngroups 2

ACE data
 Data NInput=13                ! 13 input variables
 missing=-99.989998
 Rectangular File=marker1.dat ! read data from file
 Labels locus pednum id1 id2 pibd0 pibd1 pibd2 ace1 ace2 a1s1 a2s1 a1s2 a2s2
 ! aisj is allele i of sibling j
! select only those sib pairs with no missing marker genotypes
 select if a1s1 ^= 0
 select if a2s1 ^= 0
 select if a1s2 ^= 0
 select if a2s2 ^= 0
! Phenotype sib 1, phenotype sib 2, probabilities of sharing 0, 1, and 2
! alleles ibd, genotype of sib 1, genotype of sib 2
 Select ace1 ace2 pibd0 pibd1 pibd2 a1s1 a2s1 a1s2 a2s2 ; 
! use the genotypic information to define the mean and covariance models
 Definition_variables pibd0 pibd1 pibd2 a1s1 a2s1 a1s2 a2s2 ;


 Begin Matrices;
  A Lower nvar nvar Free  	! Polygenic component
  z lower nvar nvar Free        ! Unique environmental component
  q lower nvar nvar free        ! Linkage component
  b full 3 1                    ! IBD probabilities (from Genehunter or 
                                ! Mapmaker/SIBS)
  d full 1 3                    ! will contain 0, .5, 1
  c fu 1 1                      ! constant to contain .5
! association specification
  E Full 1 4 			! First  allele of sib 1
  F Full 1 4                    ! Second allele of sib 1
  G Full 1 4 			! First  allele of sib 2
  H Full 1 4                    ! Second allele of sib 2

  L Full 1 4                    ! 0 1 0 1
  m fu 1 1 free                 ! grand mean
  U Unit 4 1                    ! Unit matrix
 
  W Full nallele 2  Free        ! a_b and a_w parameters, rows are alleles,
                                ! a_b is first column, a_w second column
 End Matrices;

! fix the constants
 matrix d 0 .5 1
 matrix c .5
! starting values for matrices used in linkage component of model
 Matrix A  .5
 Matrix z  .7
 Matrix q  .1
! matrices used for association component of model
! which take genotypes in columns 1 and 3, with columns 2 and 4
! fixed to 1.
 Matrix E 1 1 1 1
 Matrix F 1 1 1 1
 Matrix G 1 1 1 1
 Matrix H 1 1 1 1
! L is a constant
 Matrix L 0 1 0 1

 matrix m 0   ! starting value for estimate of the overall mean

! put ibd probabilities in b
 Specify b pibd0 pibd1 pibd2

! put sib genotypes into matrices
 Specify E a1s1 0 a1s1 0
 Specify F a2s1 0 a2s1 0
 Specify G a1s2 0 a1s2 0
 Specify H a2s2 0 a2s2 0

 Begin Algebra;
! compute the between effect, which is the sum of all the alleles
! in the sibship, obtained by pulling out the appropriate row of
! the first column of the W matrix, which contains the between
! effects for alleles 1 and 2, and summing those effects
  I = \part(W,E) + \part(W,F) + \part(W,G) + \part(W,H) ;
! compute the within effect, 
  J = \part(W,(E+L)) + \part(W,(F+L)) - \part(W,(G+L)) - \part(W,(H+L)) ;
! compute pihat from ibd probabilities
  p = d*b;
 End Algebra;

 Means (M + I + J | M + I - J);
 Covariance A*A' + z*z' + q*q' | c@a*a' + p@q*q' _
            c@a*a' + p@q*q' | a*a' + z*z' + q*q' ;

 Option  nd=4        ! request 4 decimal places in output
 OPtion RS Multiple  ! request residuals, multiple fit
 Option issat        ! this is saturated model for submodel comparison
End

Constrain sum of weights on diffs to equal 0
 Constraint NI=2
 Begin Matrices = Group 1;
  O Zero 1 2
  P Unit 1 nallele
 End Matrices;
 
 Constrain P*W = O ;
 Option multiple
 option issat
End                                                                            

! equate a_b to a_w
Specify 1 w
  5 5
  6 6
option df=-1
End

save combined.mxs
option issat
end

! drop both a_b and a_w against a_b = a_w
Drop w 1 1 1 to w 2 2 2
option df=-2
End

! get the a_b = a_w result
! and then drop linkage parameter
get combined.mxs

drop q 1 1 1
end