#define nvar 2
G1: DZ twins  
Data NGroups=1 NInput=13 NModel=1
 Rectangular File=basic2.dat

 Labels
famid sex1 sex2 zyg1 zyg2 age1 age2 
lffqf1 lffqf2 adh2a_1 adh2b_1 adh2a_2 adh2b_2 

!select if zyg1 > 2 !DZ's
select lffqf1 lffqf2 adh2a_1 adh2b_1 adh2a_2 adh2b_2 age1 age2 sex1 sex2  /

 Definition  adh2a_1 adh2b_1 adh2a_2 adh2b_2 age1 age2 sex1 sex2   / 

Matrices
 A Lower nvar nvar Free  	! Additive genetic influences 
 F Full 1 1 			! adh_t1_a
 G Full 1 1			! adh_t1_b
 H Full 1 1  			! adh_t1_a
 I Full 1 1			! adh_t1_b

 J Unit 1 1

 B Full 1 1  Free		! Beta for adh use for overall association
 L Full 2 2			! Age & Sex of twins
 M Full 1 2  Free		! b2, Age regression beta and male (=1) difference
 N Full 1 1  Free		! Overall mean for twin 

 C Full 1 2  			! how many classes for the mark
 V Full 9 1			! Disequilibrim/Candidate Gene Approch
 X Full 9 1			! V is the Mean, X is the 1/2 sib diff

 W Full 2 1  !Free               ! ab and aw parameters
 End Matrices
!drop @1 7
 Matrix C 1 0  
!Specify V 101 102 103 102 104 105 103 105 106
 Matrix V 
1
.5 
0 
.5
0
-.5
0 
-.5
-1
Matrix X
0
.5
1
-.5
0
.5
-1
-.5
0

 Specify F -1 ;
 Specify G -2 ;
 Specify H -3 ;
 Specify I -4 ;
 Specify L -5 -6 -7 -8 ;
Specify M 0 300  ! ensures age not included abd sex deviation 
                 !estimated as param 300
Begin Algebra;
 K = C* (((J-\sum(F))_F)| ((J-\sum(G))_G)); 
 D = C* (((J-\sum(H))_H)| ((J-\sum(I))_I));
 ! ugly (J-\sum(F))_F) gets the allele 2 status then others
 ! K is now 1 or 0 indicator for allele 1 2 sib 1
 ! D is same for sib 2
 E = \sum(K)+ \sum(K)+  \sum(K)+ \sum(D) + J;
 Z = \part((V|X),(E|J|E|J+J))*W | \part((V|-X),(E|J|E|J+J))*W ;
End Algebra;

Means (N + B*F + B*G | N + B*H + B*I)  + M*L + Z;
Covariance A*A';
 


Start .9 A 1 1 1 - A 1 2 2

!Start 0.0 B 1 1
!Start -.1 M 1 1
Start .9 N 1 1   
!Bound  -5 5 B 1 1 
!Bound -.5 .5 M 1 1
!Bound 10 150 N 1 1  
!INtervals B 1 1 
Option  nd=4 it=5000
OPtion RS Multiple
!Option issat
End

save dseq5.mxs
End

!drop allelic deviation - beta, matrix B
get dseq5.mxs
drop @.0 4
End

!drop sex deviation
get dseq5.mxs
drop @.0 300
End

!get dseq5.mxs
!drop @.0 6
!End

!get dseq5.mxs
!drop @.0 7
!End

!get dseq5.mxs
!eq W 1 1 1 W 1 2 1
!End

!get dseq5.mxs
!drop @.0 6 7
!End