!  Power calculation for IBD statuses sibships size 2

#define nsib 2
#define nsibm1 1
#define nvar 1
#define nvarnsib 2

G1: Model parameters
Calc NGroups=4
Begin Matrices;
 U Unit nsib nsib       ! for model for C 
 I Iden nsib nsib       ! for model for E 
 W Lower nvar nvar      ! QTL path
 X Lower nvar nvar      ! A path
 Y Lower nvar nvar      ! C path
 Z Lower nvar nvar      ! E path
 H Stan nsib nsib       ! for A of sib pairs,
End Matrices;

 Value .5 H 2 1 to H nsib nsibm1
 Matrix X 0.547722557505166 ! sqrt of .3 residual A 
 Matrix Z 0.774596669241483 ! sqrt of .6 residual E
 Matrix W 0.316227766016838 ! sqrt of .1 QTL variance

Begin Algebra;
 A= X*X';
 C= Y*Y';
 E= Z*Z';
 Q= W*W';
End Algebra:
End


IBD 2 pairs
Data NInput_vars=nvarnsib NObservations=80
Labels sib1 sib2
CMatrix Symmetric File=dummy.cov
Matrices= Group 1
 S Stan nsib nsib       ! for pihat        
End Matrices;
 Matrix S 1 
Cov H@A + U@C + I@E + S@Q /
Options RSidual
Option mx%e=2.cov
End 

IBD 1 pairs
Data NInput_vars=nvarnsib NObservations=80
Labels sib1 sib2
CMatrix Symmetric File=dummy.cov
Matrices= Group 1
 S Stan nsib nsib       ! for pihat        
End Matrices;
 Matrix S .5 
Cov H@A + U@C + I@E + S@Q /
Options RSidual
Option mx%e=1.cov
End 

IBD 0 pairs
Data NInput_vars=nvarnsib NObservations=80
Labels sib1 sib2
CMatrix Symmetric File=dummy.cov
Matrices= Group 1
 S Stan nsib nsib       ! for pihat        
End Matrices;
 Matrix S 0 
Cov H@A + U@C + I@E + S@Q /
Options RSidual
Option mx%e=0.cov
End 

!  Power calculation for IBD statuses sibships size 4

#define nsib 2
#define nsibm1 1
#define nvar 1
#define nvarnsib 2

G1: Model parameters
Calc NGroups=4
Begin Matrices;
 U Unit nsib nsib       ! for model for C 
 I Iden nsib nsib       ! for model for E 
 W Lower nvar nvar Free     ! QTL path
 X Lower nvar nvar Free     ! A path
 Y Lower nvar nvar          ! C path
 Z Lower nvar nvar Free     ! E path
 H Stan nsib nsib       ! for A of sib pairs,
End Matrices;

 Value .5 H 2 1 to H nsib nsibm1
 Matrix X .5
 Matrix Z .5
 Matrix W .5

Begin Algebra;
 A= X*X';
 C= Y*Y';
 E= Z*Z';
 Q= W*W';
End Algebra:
End


IBD 2 pairs
Data NInput_vars=nvarnsib NObservations=1000
Labels sib1 sib2
CMatrix Symmetric File=2.cov
Matrices= Group 1
 S Stan nsib nsib       ! for pihat        
End Matrices;
 Matrix S 1 
Cov H@A + U@C + I@E + S@Q /
Options RSidual
End 

IBD 1 pairs
Data NInput_vars=nvarnsib NObservations=2000
Labels sib1 sib2
CMatrix Symmetric File=1.cov
Matrices= Group 1
 S Stan nsib nsib       ! for pihat        
End Matrices;
 Matrix S .5 
Cov H@A + U@C + I@E + S@Q /
Options RSidual
End 

IBD 0 pairs
Data NInput_vars=nvarnsib NObservations=1000
Labels sib1 sib2
CMatrix Symmetric File=0.cov
Matrices= Group 1
 S Stan nsib nsib       ! for pihat        
End Matrices;
 Matrix S 0 
Cov H@A + U@C + I@E + S@Q /
Options RSidual Multiple
End 

Drop W 1 1 1
Option power=2.0166445E-4 1
End