! Combined Linkage and Association analyses, parental genotypes not available
! Fulker et al., 1999 ; Posthuma et al., 2004

#define n 3				! number of alleles is 3, coded 1, 2, 3
#define nvar 1				! one phenotype per individual 
#ngroups 5				! 5 Mx groups
#define ndef 2				! 2 covariates (definition variables - age and sex)

G1: calculation group between and within effects

Data Calc 

 Begin matrices;
  A Full 1 n  free			! additive allelic effects within
  C Full 1 n  free			! additive allelic effects between
  D Sdiag n n free 			! dominance deviations within
  F Sdiag n n free 			! dominance deviations between
  I Unit 1 n 				! one's
 End matrices;

Specify A 100 101 102
Specify C 200 201 202
Specify D 800 801 802
Specify F 900 901 902

 Begin algebra;
  K = (A'@I) + (A@I') ;			! Within effects, additive
  L = D + D' ;				! Within effects, dominance
  W = K+L ;				! Within effects - additive and dominance in one matrix

  M = (C'@I) + (C@I') ;			! Between effects, additive
  N = F + F' ;				! Between effects, dominance
  B = M+N ;				! Between effects - additive and dominance in one matrix
 End algebra ; 

st .2 all
end

G2: datagroup: sibship size two, MZ
 Data NInput=124
 Missing =-99.00              
 Rectangular File=example.dat
 Labels 
   famnr  zygos
   fata1 fata2 mota1 mota2		! fatherallele1 fatherallele2 motherallele1 mothereallele2
   pheno1 age1 sex1 tw1a1 tw1a2		! phenotypetwin1 agetwin1 sextwin1 twin1allele1 twin1allele2
   pheno2 age2 sex2 tw2a1 tw2a2		! phenotypetwin2 agetwin2 sextwin2 twin2allele1 twin2allele2

 !below probability ibd0(z0), ibd1(z1) and ibd2(z2) for positions 0 to 35

   z0_0 z1_0 z2_0 z0_1 z1_1 z2_1 z0_2 z1_2 z2_2 z0_3 z1_3 z2_3 z0_4 z1_4 z2_4 
   z0_5 z1_5 z2_5 z0_6 z1_6 z2_6 z0_7 z1_7 z2_7 z0_8 z1_8 z2_8 z0_9 z1_9 z2_9 
   z0_10 z1_10 z2_10 z0_11 z1_11 z2_11 z0_12 z1_12 z2_12 z0_13 z1_13 z2_13 z0_14 z1_14 z2_14
   z0_15 z1_15 z2_15 z0_16 z1_16 z2_16 z0_17 z1_17 z2_17 z0_18 z1_18 z2_18 z0_19 z1_19 z2_19
   z0_20 z1_20 z2_20 z0_21 z1_21 z2_21 z0_22 z1_22 z2_22 z0_23 z1_23 z2_23 z0_24 z1_24 z2_24
   z0_25 z1_25 z2_25 z0_26 z1_26 z2_26 z0_27 z1_27 z2_27 z0_28 z1_28 z2_28 z0_29 z1_29 z2_29
   z0_30 z1_30 z2_30 z0_31 z1_31 z2_31 z0_32 z1_32 z2_32 z0_33 z1_33 z2_33 z0_34 z1_34 z2_34
   z0_35 z1_35 z2_35

Select if zygos < 3 ; 			! MZ, zygos is coded as 1=MZM, 2=MZF, 3=DZM, 4=DZF, 5=DOS

 Select 	pheno1 pheno2
		tw1a1 tw1a2
		tw2a1 tw2a2 
		age1 sex1
		age2 sex2;

  Definition_variables 
		tw1a1 tw1a2
		tw2a1 tw2a2 
		age1 sex1
		age2 sex2;

  Begin Matrices;
	X Lower nvar nvar Free  	! Additive A
	Z Lower nvar nvar Free  	! Unshared E
	D Lower nvar nvar Free  	! Qtl Q
  
	G Full 1 nvar	 	Free	! grand mean 
	F Full 1 ndef	Fix		! will contain age1  sex1 
	I Full 1 ndef Fix		! will contain age2  sex2 
	R Full nvar ndef Free		! parameter estimates age and sex 	
  
	B Computed n n  = B1 	! spurious and genuine genotypic effects (between)
	W Computed n n  = W1	! genuine genotypic effects (within)
   
	K Full 1 4 Fix			! Will contain first  and second allele of twin1 
	L Full 1 4 Fix        		! Will contain first  and second allele of twin2 
	S Full 1 1 Fix			! Will contain 2  (for two individuals per family)
  End Matrices;

 Matrix S 2				! Now contains 2
 
 Matrix K 1 1 1 1
 Matrix L 1 1 1 1

 Specify K tw1a1 tw1a2 tw1a1 tw1a2	! allele1twin1   allele2twin1   allele1twin1   allele2twin1 , used for \part
 Specify L tw2a1 tw2a2 tw2a1 tw2a2	! allele1twin2   allele2twin2   allele1twin2   allele2twin2 , used for \part
 
 Specify F age1 sex1
 Specify I age2 sex2
 
 Begin Algebra;
  V = (\part(B,K) + \part(B,L) ) %S ; 	! Calculates sib genotypic mean (= Between effects)
  C = (\part(W,K) + \part(W,L) ) %S ;	! Calculates sib genotypic mean, used to derive deviation from this mean below (Within effects)
  A = X*X';
  E = Z*Z';
  Q = D*D';
 End Algebra;

 Means G + F*R '+ V + (\part(W,K)-C) | G + I*R' + V +(\part(W,L)-C);  
! grand mean + age and sex effects + between allele effects + within allele effects

 Covariance 
	 A+E+Q| A+Q_
         A+Q | A+E+Q ;

st 2 X 2 1 1  
st 1 Z 2 1 1 
st .5 D 2 1 1
st 10 G  2 1 1  
st -.2 R 2 1 1 
st 1.5 R 2 1 2
 Option  nd=4       			 ! request 4 decimal places in output
 OPtion RS 	 			 ! request residuals 
End

G3: datagroup: sibship size two, DZ
 Data NInput=124
 Missing =-99.00              
 Rectangular File=example.dat
 Labels 
   famnr  zygos
   fata1 fata2 mota1 mota2		! fatherallele1 fatherallele2 motherallele1 mothereallele2
   pheno1 age1 sex1 tw1a1 tw1a2		! phenotypetwin1 agetwin1 sextwin1 twin1allele1 twin1allele2
   pheno2 age2 sex2 tw2a1 tw2a2		! phenotypetwin2 agetwin2 sextwin2 twin2allele1 twin2allele2

 !below probability ibd0(z0), ibd1(z1) and ibd2(z2) for positions 0 to 35

   z0_0 z1_0 z2_0 z0_1 z1_1 z2_1 z0_2 z1_2 z2_2 z0_3 z1_3 z2_3 z0_4 z1_4 z2_4 
   z0_5 z1_5 z2_5 z0_6 z1_6 z2_6 z0_7 z1_7 z2_7 z0_8 z1_8 z2_8 z0_9 z1_9 z2_9 
   z0_10 z1_10 z2_10 z0_11 z1_11 z2_11 z0_12 z1_12 z2_12 z0_13 z1_13 z2_13 z0_14 z1_14 z2_14
   z0_15 z1_15 z2_15 z0_16 z1_16 z2_16 z0_17 z1_17 z2_17 z0_18 z1_18 z2_18 z0_19 z1_19 z2_19
   z0_20 z1_20 z2_20 z0_21 z1_21 z2_21 z0_22 z1_22 z2_22 z0_23 z1_23 z2_23 z0_24 z1_24 z2_24
   z0_25 z1_25 z2_25 z0_26 z1_26 z2_26 z0_27 z1_27 z2_27 z0_28 z1_28 z2_28 z0_29 z1_29 z2_29
   z0_30 z1_30 z2_30 z0_31 z1_31 z2_31 z0_32 z1_32 z2_32 z0_33 z1_33 z2_33 z0_34 z1_34 z2_34
   z0_35 z1_35 z2_35

Select if zygos > 2 ; 			! DZ, zygos is coded as 1=MZM, 2=MZF, 3=DZM, 4=DZF, 5=DOS


 Select 	pheno1 pheno2
		tw1a1 tw1a2
		tw2a1 tw2a2 
		age1 sex1
		age2 sex2
		z1_20 z2_20;		! this is for position 20 only, but can also be used in batch mode

  Definition_variables 
		tw1a1 tw1a2
		tw2a1 tw2a2 
		age1 sex1
		age2 sex2
		z1_20 z2_20;

 Begin Matrices;
	X Lower nvar nvar = X2  	! Additive A
	Z Lower nvar nvar = Z2  	! Unshared E
	D Lower nvar nvar = D2  	! Qtl Q
  
	G Full 1 nvar	  = G2		! grand mean 
	F Full 1 ndef			! Will contain age1 sex1   
	I Full 1 ndef 			! Will contain age2 sex2
	R Full nvar ndef  = R2		! parameter estimates age and sex
  
	B Computed n n    = B1 		! spurious and genuine genotypic effects (between)
	W Computed n n    = W1		! genuine genotypic effects (within)
   
	K Full 1 4 Fix			! Will contain first  and second allele of twin1 
	L Full 1 4 Fix        		! Will contain first  and second allele of twin2 
	S Full 1 1 Fix			! Will contain 2  (for two individuals per family)

	T Full 1 1 fix			! Will contain probability ibd=1
	U Full 1 1 fix 			! Will contain probability ibd=2
	H Full 1 1 fix			! .5	

 End Matrices;

 Matrix S 2					
 Matrix H 0.5


 Matrix K 1 1 1 1
 Matrix L 1 1 1 1

 Specify K tw1a1 tw1a2 tw1a1 tw1a2	! allele1twin1   allele2twin1   allele1twin1   allele2twin1 , used for \part
 Specify L tw2a1 tw2a2 tw2a1 tw2a2	! allele1twin2   allele2twin2   allele1twin2   allele2twin2 , used for \part
 
 Specify F age1 sex1
 Specify I age2 sex2
 
 Specify T z1_20
 Specify U z2_20

 
 Begin Algebra;
  V = (\part(B,K) + \part(B,L) ) %S ; 	! Calculates sib genotypic mean (= Between effects)
  C = (\part(W,K) + \part(W,L) ) %S ;	! Calculates sib genotypic mean, used to derive deviation from this mean below (Within effects)
  A = X*X';
  E = Z*Z';
  Q = D*D';
  J = (H*T) + U ;			! calculates pi-hat in matrix J
 End Algebra;

 Means G + F*R' + V + (\part(W,K)-C) | G + I*R' + V + (\part(W,L)-C) ;

 Covariance 
	 A+E+Q| H@A+J@Q_
         H@A+J@Q | A+E+Q ;

Options rs
 End

Constrain sum additive allelic within effects = 0 
Constraint ni=1
Begin Matrices; 
  A full 1 n = A1	
  O zero 1 1  
End Matrices;
Begin algebra;
  B = \sum(A) ;
End Algebra;
Constraint O = B ;
end

Constrain sum additive allelic between effects = 0 
Constraint ni=1
Begin Matrices;
 C full 1 n = C1	
 O zero 1 1   
End Matrices;
Begin algebra;
 B = \sum(C) ;
End Algebra;
Constraint O = B ;
option multiple issat ! this is saturated model for submodel comparison
end


!test for linkage in presence of full association
Drop D 2 1 1 
end


!Test for population stratification: between effects = within effects
Specify 1 A 100 101 102
Specify 1 C 100 101 202			!note: last parameter need not be equal because of constraint group
Specify 1 D 800 801 802
Specify 1 F 800 801 802
end

!Test for presence of dominance 
Specify 1 A 100 101 102
Specify 1 C 100 101 202
Specify 1 D 800 801 802
Specify 1 F 800 801 802
Drop @0 800 801 802
end

!Test for presence of full association
Specify 1 A 100 101 102
Specify 1 C 100 101 202
Specify 1 D 800 801 802
Specify 1 F 800 801 802
Drop @0 800 801 802 100 101	!note: parameters 102 and 202 need not be dropped because of constraint group
end

! test for linkage in absence of association, compare with -2ll from one model above
Free D 2 1 1 
end