/***********************************************************************/
/*   SASQuant1.3 - By Gusmini, Wehner, Donaghy                         */
/*      Genetic Estimates for Quantitative Traits                      */
/*                                                                     */
/*    Written in 2005                                                  */
/*    Modified:  07/29/2008  SBD                                       */
/*               Plot, Set, Family, and Gen are now defined as         */ 
/*               numeric variables instead of character variables      */
/*                                                                     */
/*   The user needs to create a dataset called "Data.dat".             */
/*      The dataset should have the following columns:                 */
/*          Plot number, set number, family number, generation number  */
/*          followed by columns of ratings, one per trait measured.    */
/*                                                                     */
/*          Generation numbers should be as follows:                   */
/*               Parent A      = 1                                     */
/*               Parent B      = 2                                     */
/*               (A X B) F1    = 3                                     */
/*               (A X B) F2    = 4                                     */
/*               (F1 X A) BC1A = 5                                     */
/*               (F1 X A) BC1B = 6                                     */
/*                                                                     */
/*   The user needs to modify further statements at the end of the     */
/*   to customize the analysis.                                        */
/*                                                                     */
/*   Abbreviations needed to interpret the outputs are:                */
/*            Tables 1 and 2                                           */
/*               NPa, MPa, VarPa    = N, Mean, and Variance of Parent A */
/*               NPb, MPb, VarPb    = N, Mean, and Variance of Parent B */
/*               NF1, MF1, VarF1    = N, Mean, and Variance of F1      */
/*               NF2, MF2, VarF2    = N, Mean, and Variance of F2      */
/*               NBCa, MBCa, VarBCa = N, Mean, and Variance of BC1A    */
/*               NBCb, MBCb, VarBCb = N, Mean, and Variance of BC1B    */
/*            Table 3                                                  */
/*               VARP = Phenotypic Variance                            */
/*               VARE = Environmental Variance                         */
/*               VARG = Genetic Variance                               */
/*               VARA = Additive Variance                              */
/*               VARD = Dominant Variance                              */
/*               HERN = Narrow-Sense Heritability                      */
/*               HERB = Broad-Sense Heritability                       */
/*            Table 4                                                  */
/*               EF1-EF5 = Number of Effective Factors (Estimates 1-5) */
/*               EFM     = Mean Number of Effective Factors            */
/*               GS05    = Gain from Selection at i= 5%                */
/*               GS10    = Gain from Selection at i=10%                */
/*               GS20    = Gain from Selection at i=20%                */
/*            Table 5                                                  */
/*               m  = Hayman's Overall Mean                            */
/*               a  = Hayman's Additive Component                      */
/*               d  = Hayman's Dominant Component                      */
/*               aa = Hayman's Add.*Add. Epistatic Component           */
/*               ad = Hayman's Add.*Dom. Epistatic Component           */
/*               dd = Hayman's Dom.*Dom. Epistatic Component           */
/*               SEm  = Hayman's Standard Error for 'm'                */
/*               SEa  = Hayman's Standard Error for 'a'                */
/*               SEd  = Hayman's Standard Error for 'd'                */
/*               SEaa = Hayman's Standard Error for 'aa'               */
/*               SEad = Hayman's Standard Error for 'ad'               */
/*               SEdd = Hayman's Standard Error for 'dd'               */
/*               Pm  = Probability of a Larger Student's T for 'm'     */
/*               Pa  = Probability of a Larger Student's T for 'a'     */
/*               Pd  = Probability of a Larger Student's T for 'd'     */
/*               Paa = Probability of a Larger Student's T for 'aa'    */
/*               Pad = Probability of a Larger Student's T for 'ad'    */
/*               Pdd = Probability of a Larger Student's T for 'dd'    */
/*                                                                     */
/*   The number of effective factors is calculated as suggested by:    */
/*               EF1 =     Wright, 1968                                */
/*               EF2 =     Mather & Jinks, 1982                        */
/*               EF3-EF5 = Lande, 1981                                 */
/***********************************************************************/
%macro est_sel;                                                                                                                         
proc sort data=sub;                                                                                                               
      by gen family set;                                                                                                             
data _dst2; set sub;                                                                                                              
     array t &dep;                                                                                                                  
     do over t;                                                                                                                     
        trait=_i_;                                                                                                                  
        y=t;                                                                                                                      
        output;                                                                                                                   
     end;                                                                                                                      
     drop &dep;                                                                                                                 
proc sort; by trait gen family set;                                                                                                     
proc means noprint;                                                                                                               
     by trait gen family set;                                                                                                       
     var y;                                                                                                                      
     output out=_dst2 mean=mean  var=var n=n stderr=stderr;                                                                                     
proc sort data=_dst2;                                                                                                              
     by trait &byvar;                                                                                                            

data _dst2a;                                                                                                                       
     set _dst2;                                                                                                                     
            if gen=1 then do;MPa =mean;VarPa =var;SEPa =stderr;NPa =n;end;                                                                            
            if gen=2 then do;MPb =mean;VarPb =var;SEPb =stderr;NPb =n;end;                                                                            
            if gen=3 then do;MF1 =mean;VarF1 =var;SEF1 =stderr;NF1 =n;end;                                                                            
            if gen=4 then do;MF2 =mean;VarF2 =var;SEF2 =stderr;NF2 =n;end;                                                                            
            if gen=5 then do;MBCa=mean;VarBCa=var;SEBCa=stderr;NBCa=n;end;                                                                            
            if gen=6 then do;MBCb=mean;VarBCb=var;SEBCb=stderr;NBCb=n;end;                                                                            
proc means noprint data=_dst2a;                                                                                                    
      by trait &byvar;                                                                                                            
      var mpa mpb mf1 mf2 mbca mbcb                                                                                               
          varpa varpb varf1 varf2 varbca varbcb                                                                                               
          sepa sepb sef1 sef2 sebca sebcb                                                                                             
          npa npb nf1 nf2 nbca nbcb;                                                                                            
      output out=_dst2b                                                                                                   
              mean= ;                                                                            
data _dst2c;                                                                                                                       
     set _dst2b;                                                                                                                     
            VarA=2*varf2-varbca-varbcb;                                                                                                       
            VarE=(varpa+varpb+2*varf1)/4;                                                                                                     
            VarP=varf2;                                                                                                                   
            VarG=varp-vare;                                                                                                             
            VarD=varg-vara;                                                                                                             
            HerN=vara/varp;                                                                                                             
            HerB=varg/varp;                                                                                                             
            h=(mf1-mpa)/(mpb-mpa);
            EF1=(mpb-mpa)**2*(1.5-2*h*(1-h))/(8*(varf2-.25*(varpa+varpb+2*varf1)));                                
            EF2=.5*(mpb-mpa)**2/(2*varf2-(varbca+varbcb));                                               
            EF3=(mpb-mpa)**2/(8*(varf2-.25*(varpa+varpb+2*varf1)));                                       
            EF4=(mpb-mpa)**2/(8*(2*varf2-(varbca+varbcb)));                                             
            EF5=(mpb-mpa)**2/(8*(varbca+varbcb-(varf1+.5*varpa+.5*varpb)));                                       
            EFm=mean(ef1,ef2,ef3,ef4,ef5);                                                                                              
            GS05=2.06*hern*sqrt(varf2);                                                                    
            GS10=1.76*hern*sqrt(varf2);                                                                    
            GS20=1.40*hern*sqrt(varf2);                                                                    
            m =mf2;                                                                                      
            a =mbca-mbcb;                                                                                                               
            d =(-0.5*mpa)-(0.5*mpb)+mf1-(4*mf2)+(2*mbca)+(2*mbcb);                                                                      
            aa=-(4*mf2)+(2*mbca)+(2*mbcb);                                                                                              
            ad=(-0.5*mpa)+(0.5*mpb)+mbca-mbcb;                                                                                          
            dd=mpa+mpb+(2*mf1)+(4*mf2)-(4*mbca)-(4*mbcb);                                                                               
			DFm=nf2-1;                                                                                                                  
            DFa=mean(nbca,nbcb)-1;                                                                                                      
            DFd=nf2-1;                                                                                                                    
            DFaa=mean(nf2,nbca,nbcb)-1;                                                                                                 
            DFad=mean(nbca,nbcb)-1;                                                                                                                   
            DFdd=mean(nf2,nbca,nbcb)-1;                                                                                                                  
            SEm =sef2;                                                     
			SEa =sebca+sebcb;                                                                                                                                                   
            SEd =(0.5*sepa)+(0.5*sepb)+sef1+(4*sef2)+(2*sebca)+(2*sebcb);                                                                                                          
            SEaa=(4*sef2)+(2*sebca)+(2*sebcb);                                                                                                                                  
            SEad=(0.5*sepa)+(0.5*sepb)+sebca+sebcb;                                                                                                                              
            SEdd=sepa+sepb+(2*sef1)+(4*sef2)+(4*sebca)+(4*sebcb);                                                                                                                                         
            Tm=m/SEm;                                                                                                                   
            Ta=a/SEa;                                                                                                                   
            Td=d/SEd;                                                                                                                   
            Taa=aa/SEaa;                                                                                                                
            Tad=ad/SEad;                                                                                                                
            Tdd=dd/SEdd;                                                                                                                
            Pm=(1-probt(abs(Tm),DFm))*2;                                                                                             
            Pa=(1-probt(abs(Ta),DFa))*2;                                                                                             
            Pd=(1-probt(abs(Td),DFd))*2;                                                                                             
            Paa=(1-probt(abs(Taa),DFaa))*2;                                                                                          
            Pad=(1-probt(abs(Tad),DFad))*2;                                                                                          
            Pdd=(1-probt(abs(Tdd),DFdd))*2;                                                                                          
                                                                                                                                        
      format NPa NPb NF1 NF2 NBCa NBCb 4.0;                                                                                             
      format MPa MPb MF1 MF2 MBCa MBCb 4.1;                                                                                             
      format VarPa VarPb VarF1 VarF2 VarBCa VarBCb
			 VARP VARE VARG VARA VARD HERB HERN                                                               
             m  SEm a  SEa d  SEd
			 aa SEaa ad  SEad dd  SEdd 6.2;                                                
      format EF1 EF2 EF3 EF4 EF5 EFM GS05 GS10 GS20 5.1;                                                                                
      format Pm Pa Pd  Paa Pad Pdd 6.4;                                                                               
      format Tm Ta Td Taa Tad Tdd 5.2;                             
proc print uniform data=_dst2c noobs;                                                                                              
     by trait;                                                                                                                      
     var &byvar NPa MPa NPb MPb NF1 MF1 NF2 MF2 NBCa MBCa NBCb MBCb;                                                             
     title5 "TABLE 1. Number of Observations and Generation Means by &byvar";                                                    
proc print uniform data=_dst2c noobs;                                                                                              
     by trait;                                                                                                                      
     var &byvar VarPa VarPb VarF1 VarF2 VarBCa VarBCb;                                                                                       
     title5 "TABLE 2. Generation Variances by &byvar";  
proc print uniform data=_dst2c noobs;                                                                                              
     by trait;                                                                                                                      
     var &byvar VarP VarE VarG VarA VarD HerB HerN;                                                                              
     title5 "TABLE 3. Genetic Variances and Heritability by &byvar";                                                               
proc print uniform data=_dst2c noobs;                                                                                              
     by trait;                                                                                                                      
     var &byvar EF1 EF2 EF3 EF4 EF5 EFM GS05 GS10 GS20;                                                                          
     title5 "TABLE 4. Effective Factors and Gain from Selection by &byvar";
proc print uniform data=_dst2c noobs;                                                                                              
     by trait;                                                                                                                      
     var &byvar m  SEm  Pm
            	a  SEa  Pa
				d  SEd  Pd;                                            
     title5 "TABLE 5. Hayman's Main Gene Effects by &byvar";
proc print uniform data=_dst2c noobs;                                                                                              
     by trait;                                                                                                                      
     var &byvar aa SEaa Paa
				ad SEad Pad
				dd SEdd Pdd;                                            
     title5 "TABLE 6. Hayman's Epistatic Gene Effects by &byvar";
proc datasets library=work nolist;
   delete _dst2 _dst2a _dst2b _dst2c / memtype=data;
run;
quit;
%mend est_sel;                                                                                                                          
title  'SASQuant1.1 - By Gusmini, Wehner, Donaghy';                                                                                     
title3 'Genetic Estimates for Quantitative Traits';                                                                                     
title4 ' ';                                                                                                                             
options pageno=1 nodate;

/***********************************************************************/
/*   The user needs to modify the following lines as indicated below:  */
/*                                                                     */
/*        %let dep=trait1 trait2;                                      */
/*        input plot  set  family  gen  trait1 trait2;                 */
/*                 Modify the list traits to be analyzed.              */
/*                                                                     */
/*        infile 'Data.dat' missover;                                  */
/*                 Modify the name of the data file, if different.     */
/*                                                                     */
/*        %let byvar=set family;                                       */
/*                 Modify the list of blocking factors to be used      */
/*                 for the analysis as follows:                        */
/*                      Set        = blocking families by set          */
/*                      Family     = blocking sets by family           */
/*                      Set Family = no blocking                       */
/***********************************************************************/

%let dep=trait1 trait2;                                                                                                                 
%let byvar=Set Family;                                                                                                                                                                        
data dst1;                                                                                                                              
     infile 'data.dat' missover;                                                                                                 
     input Plot  Set  Family  Gen  trait1 trait2;
 
proc sort data=dst1;                                                                                                                    
     by &byvar;
data sub; set dst1;

/***********************************************************************/
/*   The user needs to modify the following lines as indicated below:  */
/*                                                                     */
/*   To analyze data for a sub-set of the original dataset:            */
/*        1. Uncomment the IF FAMILY IN... statement by deleting the *.   */
/*               The * will be needed to analyze the entire dataset.   */
/*        2. List between () the families to be analyzed.              */
/*               Use the format: ('family number','family number')     */
/***********************************************************************/

   * if family in(17,18,19); 

                                                                                                                                        
%est_sel;                                                                                                                               
run;