%macro FixedBiplot(data,env_n,gen_n,rep_n,var,env,gen,rep,gencod);

/**************TO GET RESIDUALS FROM THE ADDITIVE MODEL  *************/

proc mixed data=&DATA info noitprint noclprint noinfo;
  class &GEN &ENV &REP;
  model &VAR = &ENV &GEN /outpred=resid1;
  random &REP(&ENV);
  ods listing exclude outpred;
  id &GEN &ENV &REP;


/************* COMPUTE AVERAGE OF RESIDUALS TO OBTAIN THE E MATRIX ******/

proc sort data=resid1;
       by &ENV &GEN ;
proc means mean noprint;
       by &ENV &GEN ;
       var resid;
output out=resid2 mean=res;

/************************************************************************/
/******************  START IML FOR RESIDUAL ANALYSIS  *******************/
/************************************************************************/

proc iml;
    use resid2 var{&ENV &GEN res};
    read all var{res} into r;

    e=shape(r,&ENV_N);

    call svd(envvec,val1,genvec,e);
    val=val1#val1*&REP_N;

/******************* PRINT PCA RESULTS  *********************************/

    propeig= val/val[+];
    propacum=propeig;
    do i=2 to nrow(val);
       propacum[i]=propacum[i-1]+propeig[i];
    end;

    print 'Proportions of each component', val propeig propacum;


    labe='env_1':'env_25';
    labg='gen_1':'gen_25';


    create vece from envvec (|colname=labe|);
    append from envvec;


    create vecg from genvec (|colname=labg|);
    append from genvec;


/*** CREATE DATASET FOR COMPUTING PREDICTED VALUES ***********************/

   matcomb = envvec` || genvec` || val1;
   label1="env1":"env&ENV_N" ;
   label2="gen1":"gen&GEN_N";
   label3={"lambda"};
   label=label1||label2||label3;
   create matrix from matcomb (|colname=label|);
   append from matcomb;

quit;



/*** LR TESTS TO SELECT NUMBER OF SIGNIFICANT MULTIPLICATIVE TERMS ******/

proc mixed noitprint noclprint data=&DATA method=ML noinfo;
class &ENV &GEN &REP;
model &VAR = &ENV &GEN /outpredm=pre0;
random &REP(&ENV) ;
ods listing exclude outpredm;
ods output FitStatistics=pca0;
id &ENV &GEN &REP &VAR;

data pre0; set pre0;check=1;
data matrix1; set matrix;if _n_=1;check=1;


data pred1;
merge pre0 matrix1;
by check;
array environ {&ENV_N} env1-env&ENV_N;
array genotyp {&GEN_N} gen1-gen&GEN_N;
keep &VAR &ENV &GEN &REP PRED RESID check;
PRED = PRED+lambda*environ(&ENV)*genotyp(&GEN);
RESID= &VAR - PRED;

data matrix1; set matrix;if _n_=2;check=1;

data pred2;
merge pred1 matrix1;
by check;
array environ{&ENV_N} env1-env&ENV_N;
array genotyp{&GEN_N} gen1-gen&GEN_N;
keep &VAR &ENV &GEN &REP PRED RESID check;
PRED  = PRED+lambda*environ(&ENV)*genotyp(&GEN);
RESID = &VAR - PRED;

data matrix1; set matrix;if _n_=3;check=1;

data pred3;
merge pred2 matrix1;
by check;
array environ{&ENV_N} env1-env&ENV_N;
array genotyp{&GEN_N} gen1-gen&GEN_N;
keep &VAR &ENV &GEN &REP PRED RESID check;
PRED = PRED+lambda*environ(&ENV)*genotyp(&GEN);
RESID = &VAR - PRED;

proc mixed noitprint noclprint data=pred1 method=ml noinfo;
class &ENV &GEN &REP;
model RESID= /noint;
random &REP(&ENV);
ods output FitStatistics=pca1;
run;

proc mixed noitprint noclprint data=pred2 method=ml noinfo;
class &ENV &GEN &REP;
model RESID= /noint;
random &REP(&ENV);
ods output FitStatistics=pca2;
run;

proc mixed noitprint noclprint data=pred3 method=ml noinfo;
class &ENV &GEN &REP;
model RESID = /noint;
random &REP(&ENV);
ods output FitStatistics=pca3;
run;

proc mixed noitprint noclprint data=&DATA method=ml noinfo;
class &ENV &GEN &REP;
model &VAR = &ENV|&GEN;
random &REP(&ENV);
ods output FitStatistics=pcall;
run;


data pca0; set pca0; if descr='-2 Log Likelihood'; pca0=value; keep pca0;
run;


data pca1; set pca1; if descr='-2 Log Likelihood'; pca1=value; keep pca1;
run;


data pca2; set pca2; if descr='-2 Log Likelihood'; pca2=value; keep pca2;
run;


data pca3; set pca3; if descr='-2 Log Likelihood'; pca3=value; keep pca3;
run;

data pcall; set pcall; if descr='-2 Log Likelihood'; pcall=value; keep pcall;
run;

data likratio; merge pca0 pca1 pca2 pca3  pcall;
pca=0; chisq=pca0-pcall; df=(&ENV_N-1)*(&GEN_N-1);
       p_value=1-probchi(chisq,df); output;
pca=1; chisq=pca1-pcall; df=(&ENV_N-1)*(&GEN_N-1)-(&ENV_N+&GEN_N-3);
       p_value=1-probchi(chisq,df); output;
pca=2; chisq=pca2-pcall; df=(&ENV_N-1)*(&GEN_N-1)-(2*&ENV_N+2*&GEN_N-6);
       p_value=1-probchi(chisq,df); output;
pca=3; chisq=pca3-pcall; df=(&ENV_N-1)*(&GEN_N-1)-(3*&ENV_N+3*&GEN_N-11);
       p_value=1-probchi(chisq,df); output;

proc print data=likratio;
var pca chisq df p_value;
Title 'Goodness of fit LRTs for adding PCA terms';
run;



/***********  PREPARE AND GENERATE BIPLOTS  ************************************/

proc sort data=&DATA;
    by &ENV;

proc means data=&DATA noprint;
    var &VAR;
    by &ENV;id &gencod;
    output out=enva mean=ydoteye;

proc sort data=&DATA;
    by &GEN;

proc means data=&DATA noprint;
    var &VAR;
    by &GEN; id &gencod;
    output out=gena mean=ydoteye;

proc sort data=&DATA;
    by &GEN &ENV;
proc means data=&DATA noprint;
    var &VAR;
    by &GEN &ENV; id &gencod;
    output out=means mean=;

data enva;
    merge enva vece;
    PC1=env_1;
    PC2=env_2;
    PC3=env_3;
proc print;
data gena;
    merge gena vecg;
    PC1=gen_1;
    PC2=gen_2;
    PC3=gen_3;

data envanno(keep=xsys ysys x y color function position size text style);
    length text $ 8;
    set enva;
    text=&ENV;
    style = 'SWISSB';
    xsys='2'; ysys='2'; color='blue'; position='5'; function='label';
    size=1.5;
    x=PC1;
    y=PC2;
data genanno(keep=xsys ysys x y color function position size text style);
    length text $ 8;
    set gena;
    text=&GENCOD;
    style = 'ZAPFB';
    xsys='2'; ysys='2'; color='red'; position='5'; function='label';
    size=1;
    x=PC1;
    y=PC2;
data vecann1;
    set envanno genanno;


data envanno(keep=xsys ysys x y color function position size text style);
    length text $ 8;
    set enva;
    text=&ENV;
    style = 'SWISSB';
    xsys='2'; ysys='2'; color='blue'; position='5'; function='label';
    size=1.5;
    x=PC1;
    y=PC3;
data genanno(keep=xsys ysys x y color function position size text style);
    length text $ 8;
    set gena;
    text=&GENCOD;
    style = 'ZAPFB';
    xsys='2'; ysys='2'; color='red'; position='5'; function='label';
    size=1;
    x=PC1;
    y=PC3;
data vecann13;
    set envanno genanno;



data envanno(keep=xsys ysys x y color function position size text style);
    length text $ 8;
    set enva;
    text=&ENV;
    style = 'SWISSB';
    xsys='2'; ysys='2'; color='blue'; position='5'; function='label';
    size=1.5;
    x=ydoteye;
    y=PC1;
data genanno(keep=xsys ysys x y color function position size text style);
    length text $ 8;
    set gena;
    text=&GENCOD;
    style = 'ZAPFB';
    xsys='2'; ysys='2'; color='red'; position='5'; function='label';
    size=1;
    x=ydoteye;
    y=PC1;
data vecann21;
    set envanno genanno;


data envanno(keep=xsys ysys x y color function position size text style);
    length text $ 8;
    set enva;
    text=&ENV;
    style = 'SWISSB';
    xsys='2'; ysys='2'; color='blue'; position='5'; function='label';
    size=1.5;
    x=ydoteye;
    y=PC2;
data genanno(keep=xsys ysys x y color function position size text style);
    length text $ 8;
    set gena;
    text=&GENCOD;
    style = 'ZAPFB';
    xsys='2'; ysys='2'; color='red'; position='5'; function='label';
    size=1;
    x=ydoteye;
    y=PC2;
data vecann22;
    set envanno genanno;


data vectors;
    set enva gena;

title '                                                    ';
proc gplot data = vectors;
    symbol1 v=none i=none color=white;
    plot PC1*ydoteye=1 PC1*ydoteye=1/anno=vecann21 overlay vref=0 href=0;
    *title 'First Multiplicative Component vs. Mean Response';
run;

proc gplot data=vectors;
    symbol1 v=none i=none color=white;
    plot PC2*PC1=1 PC2*PC1=1/anno=vecann1 overlay vref=0 href=0 ;
    *title1 'Second vs. First Multiplicative Component Scores';
run;

proc gplot data=vectors;
    symbol1 v=none i=none color=white;
    plot PC3*PC1=1 PC3*PC1=1/anno=vecann13 overlay vref=0 href=0 ;
    *title1 'Third vs. First Multiplicative Component Scores';
run;

/******************  Show GEI table ******************************************/
title '                                                    ';
proc iml;
    use means var{&ENV &GEN &VAR};
    read all var{&VAR} into y;
    e=shape(y,&GEN_N);
    print 'Matrix subjected to SVD';
    print e;
quit;

%mend;