%macro res_plot(ds);

 /************************************************************
 /  Plots Y-residuals vs. predicted values for each PLS      /
 /      component.                                           /
 /  Variable:                                                /
 /      DS -        The input data set:  Must at least       /
 /                  contain variables N for observation      /
 /                  numbers, as well as predicted values and /
 /                  residuals and should not contain missing /
 /                  values.                                  /
 /  Required Global Variables:                               /
 /      %RESNAME -  Prefix given to Y-residuals when OUT=    /
 /                  data set was defined in PROC PLS.        /
 /      %PREDNAME - Prefix given to predictions when OUT=    /
 /                  data set was defined in PROC PLS.        /
 /      %NUM_Y    - Number of response variables Y.          /
 ************************************************************/

data _NULL_; set &ds;
  call symput('max_n',n);
run;

%do i=1 %to &num_y;

  axis1 label=(angle=270 rotate=90 "Residual")
        major=(number=5) minor=none;

  axis2 label=("Prediction for Response &i") minor=none;

  data res_anno;     *** Annotation Data Set for Plot ***;
    length text $ %length(&max_n);
    retain function 'label' position '5' hsys '3' xsys '2' ysys '2' ;
    set &ds;
    text=%str(n);  x=&predname&i;  y=y&resname&i;
  run;

  proc gplot data=&ds;
    plot y&resname&i*&predname&i/anno=res_anno vaxis=axis1
                                 haxis=axis2 vref=0 lvref=2 frame;
    symbol1 v=none i=none;
  run;

%end;

%mend;

%macro nor_plot(ds);

 /************************************************************
 /  Plots Y-residuals vs. Normal quantiles for each PLS      /
 /  component.                                               /
 /  Variable:                                                /
 /      DS -       The input data set:  Must at least        /
 /                 contain variables for observation         /
 /                 numbers, predicted values and residuals   /
 /                 and should not contain missing values.    /
 /  Required Global Variables:                               /
 /      %RESNAME - Prefix given to Y-residuals when OUT=     /
 /                 data set was defined in PROC PLS.         / 
 /      %NUM_Y   - Number of response variables Y.           /
 ************************************************************/

data ds; set &ds;
run;

data _NULL_; set &ds;
  call symput('max_n',n);
run;


%do i=1 %to &num_y;
  data ds; set ds;
    if y&resname&i=. then delete;
  run;
%end;

data _NULL_; set ds;
  call symput('numobs',_N_);
run;

%do i=1 %to &num_y;

  proc sort data=ds; by y&resname&i;

  /***********************************************************
  /  Calculate the expected values under normality for each  /
  /  residual.                                               /
  ***********************************************************/

  data resid&i; set ds(keep=n y&resname&i);
    v=(_n_ - 0.375)/(&numobs+0.25);
    z=probit(v);
  run;

  axis1 label=(angle=270 rotate=90 "Y&i Residual")
        major=(number=5) minor=none;

  axis2 label=('Normal Quantile') minor=none;

  data nor_anno;             *** Annotation Data Set for Plot ***;
    length text $ %length(&max_n);
    retain function 'label' position '5' hsys '3' xsys '2' ysys '2' ;
    set resid&i;
    text=%str(n); x=z; y=y&resname&i;
  run;

  proc gplot data=resid&i;
    plot y&resname&i*z/anno=nor_anno vaxis=axis1 haxis=axis2
                      frame;
    symbol1 v=none i=none;
  run;

%end;

%mend;

%macro plot_scr(ds,
                max_lv=&lv);

 /************************************************************
 /  Plots the Y-scores vs. the corresponding X-scores for    /
 /  each PLS component.                                      /
 /  Variables:                                               /
 /      DS -        The data set containing the scores and a /
 /                  variable N containing the observation    /
 /                  numbers.                                 /
 /      MAX_LV -    Number of the last PLS component to have /
 /                  its scores plotted.                      / 
 /  Required Global Variables:                               /
 /      %XSCRNAME - Prefix given to X-scores when OUT= data  /
 /                  set was defined.                         /
 /      %YSCRNAME - Prefix given to Y-scores when OUT= data  /
 /                  set was defined in PROC PLS.             /
 ************************************************************/

data dsout; set &ds;      *** Uses nonmissing observations ***;
  if n ^= .;
run;

data _NULL_; set &ds;
  call symput('max_n',n);
run;


%do i=1 %to &max_lv;

  data pltanno;           *** Annotation Data Set for Plot ***;
    length text $ %length(&max_n);
    retain function 'label' position '5' hsys '3' xsys '2' ysys '2' ;
    set dsout;
    text=%str(n); x=&xscrname&i; y=&yscrname&i;
  run;

  axis1 label=(angle=270 rotate=90 "Y score &i")
        major=(number=5) minor=none;

  axis2 label=("X-score &i") minor=none;

  proc gplot data=dsout;
    plot &yscrname&i*&xscrname&i/anno=pltanno vaxis=axis1
                                 haxis=axis2 frame;
    symbol1 v=none i=none;
  run;

%end;
%mend plot_scr;

%macro plotxscr(ds,
                max_lv=&lv);

 /************************************************************
 /  Plots X-scores for a given number of PLS components      /
 /  vs. those of the preceding PLS component.                /
 /  Variables:                                               /
 /      DS -        The data set containing the scores and a /
 /                  variable N containing the observation    /
 /                  numbers.                                 /
 /      MAX_LV -    Number of the last PLS component to have /
 /                  its scores plotted.                      /
 /  Required Global Variables                                /
 /      %XSCRNAME - Prefix given to X-scores when OUT= data  /
 /                  set was defined in PROC PLS.             /
 ************************************************************/

data dsout; set &ds;
  if n ^= .;               *** Uses nonmissing observations ***;
run;

data _NULL_; set &ds;
  call symput('max_n',n);
run;

%do i=1 %to %eval(&max_lv-1);

  %let j=%eval(&i+1);

  data pltanno;            *** Annotation Data Set for Plot ***;
    length text $ %length(&max_n);
    retain function 'label' position '5' hsys '3' xsys '2' ysys '2' ;
    set dsout;
    text=%str(n); x=&xscrname&i; y=&xscrname&j;
  run;

  axis1 label=(angle=270 rotate=90 "X score &j")
        major=(number=5) minor=none;

  axis2 label=("X-score &i") minor=none;

  proc gplot data=dsout;
    plot &xscrname&j*&xscrname&i/anno=pltanno vaxis=axis1
                                 haxis=axis2 frame;
    symbol1 v=none i=none;
  run;

%end;

%mend plotxscr;

%macro get_wts(dsoutmod,
               dsxwts=xwts);

 /***********************************************************
 /  Gets X-weights w from OUTMODEL data set:                /
 /    1.  Gets appropriate section of OUTMODEL data set.    /
 /    2.  Outputs this data set as DSXWTS1 (will be used    /
 /        in VIP calculation.)                              /
 /    3.  Transposes the data set so the w's are the        /
 /        variables.                                        /
 /    4.  Renames the columns to w1 - wA, where A is the    /
 /        number of PLS components LV in the final model.   /
 /  Variables:                                              /
 /      DSOUTMOD - Name of the OUTMODEL data set generated  /
 /                 by proc PLS.                             /
 /      DSXWTS -   Name of the data set containing the      /
 /                 X-weights as variables that is output    /
 /                 by this macro.                           /
 /  Required Global Variable:                                /
 /      %XVARS -   Macro variable containing the names of    /
 /                 the X-variables in a string with space    /
 /                 delimiters.                               / 
 ************************************************************/

data &dsxwts; set &dsoutmod(keep=_TYPE_ _LV_ &xvars);
  if _TYPE_='WB' then output;

proc transpose data=&dsxwts out=&dsxwts; run;

data &dsxwts; set &dsxwts;
  if _NAME_='_LV_' then delete;
  n=_n_-1;
run;

%do i=1 %to &lv;

  data &dsxwts; set &dsxwts;
    rename col&i=w&i;
  run;

%end;

%mend;

%macro plot_wt(ds,
               max_lv=&lv);

 /************************************************************
 /  Plots X-weights for a given number of PLS components     /
 /  vs. those of the preceding PLS component.                /
 /  Variables:                                               /
 /      DS -       Name of the data set containing the       /
 /                 weights as variables w1-wA, where A=LV,   /
 /                 the number of PLS components, and a       /
 /                 character variable _NAME_ containing the  /
 /                 X-variable names.                         /
 /      MAX_LV -   Number of the last PLS component to have  /
 /                 its weights plotted.                      /
 /  Required Global Variable:                                /
 /      %XVARS -   Macro variable containing the names of    /
 /                 the X-variables in a string with space    /
 /                 delimiters.                               /
 ************************************************************/

 /***********************************************************
 /  Determine the largest label to be put on plot           /
 ***********************************************************/

%let name_len=1;

data _NULL_; set &ds;
  call symput('num_x',_N_);
run;

%do i=1 %to &num_x;
  %let temp=%scan(&xvars,&i,%str( ));
  %if %length(&temp)>&name_len %then %do;
    %let name_len=%length(&temp);
  %end;
%end;

 /***********************************************************
 /  Plot X-weights for each PLS component                   /
 ***********************************************************/

%do i=1 %to %eval(&max_lv-1);

  %let j=%eval(&i+1);

  data wt_anno;              *** Annotation Data Set for Plot ***;
    length text $ &name_len;
    retain function 'label' position '5' hsys '3' xsys '2' ysys '2' ;
    set &ds;
    text=%str(_name_); x=w&i; y=w&j;
  run;

  axis1 label=(angle=270 rotate=90 "X weight &j")
        major=(number=5) minor=none;

  axis2 label=("X-weight &i") minor=none;

  proc gplot data=&ds;
    plot w&j*w&i/anno=wt_anno vaxis=axis1 haxis=axis2 frame;
    symbol1 v=none i=none;
  run;

%end;

%mend;

%macro pltwtfrq(ds,
                plotyvar=W,
                plotxvar=f,
                max_lv=&lv,
                label=Weight);

 /************************************************************
 /  Plots X-Weights or X-Loadings versus the frequency for   /
 /  spectrometric calibration data sets.                     /
 /  Variables:                                               /
 /      DS -       Data set containing the weights/loadings  /
 /                 as variables with each observation        /
 /                 representing the weights for a particular /
 /                 X-variable, which in this case is a       /
 /                 frequency.                                /
 /      PLOTYVAR - The name (excluding the component number) /
 /                 of the weight/loading variables.  For     /
 /                 example, PLOTYVAR=w if the variables to   /
 /                 be plotted are w1, w2, w3,...             /
 /      PLOTXVAR - The variable name of the frequency        /
 /                 variable.                                 /
 /      MAX_LV -   Number of PLS components to be plotted    /
 /      LABEL -    The label for the vertical axis in the    /
 /                 plot.                                     /
 ************************************************************/

axis1 label=(angle=270 rotate=90 "&label")
      major=(number=5) minor=none;
axis2 label=("Frequency") minor=none;

%let plotvars=%str( );

%do i=1 %to &max_lv;
  %let plotvars=%str(&plotvars &plotyvar&i);
%end;

proc gplot data=&ds;
  plot (&plotvars)*&plotxvar/overlay legend vaxis=axis1
                             haxis=axis2 vref=0 lvref=2 frame;
  symbol1 v=none i=spline;
run;

%mend;

%macro getxload(dsoutmod,
                dsxload=xloads);

 /***********************************************************
 /  Gets X-loadings p from OUTMODEL data set:               /
 /    1.  Gets appropriate section of OUTMODEL data set.    /
 /    2.  Transposes it so the p's are column vectors.      /
 /    3.  Renames the columns to p1 - pA, where A is the    /
 /        number of PLS components in the final model.      /
 /  Variables:                                              /
 /      DSOUTMOD - Name of the OUTMODEL data set produced   /
 /                 by proc PLS.                             /
 /      DSXLOAD -  Name of the data set to contain the      /
 /                 X-loadings as variables.                 /
 /  Required Global Variable:                                /
 /      %XVARS -   Macro variable containing the names of    /
 /                 the X-variables in a string with space    /
 /                 delimiters.                               / 
 ***********************************************************/

data &dsxload; set &dsoutmod(keep=_TYPE_ &xvars);
  if _TYPE_='PQ' then output;

proc transpose data=&dsxload out=&dsxload; run;

data &dsxload; set &dsxload;
  n=_N_;
run;

%do i=1 %to &lv;

  data &dsxload; set &dsxload;
    rename col&i=p&i;
  run;

%end;

%mend;

%macro pltxload(ds,
                max_lv=&lv);

 /************************************************************
 /  Plots X-loadings for a given number of PLS components    /
 /  vs. those of the preceding PLS component.                /
 /  Variables:                                               /
 /      DS -       Name of the data set containing the       /
 /                 loadings as variables p1-pA, where A=LV,  /
 /                 the number of PLS components, and a       /
 /                 character variable _NAME_ containing the  /
 /                 X-variable names.                         /
 /      MAX_LV -   Number of the last PLS component to have  /
 /                 its loadings plotted.                     /
 /  Required Global Variable:                                /
 /      %XVARS -   Macro variable containing the names of    /
 /                 the X-variables in a string with space    /
 /                 delimiters.                               /
 ************************************************************/

 /***********************************************************
 /  Determine the largest label to be put on plot           /
 ***********************************************************/

%let name_len=1;

data _NULL_; set &ds;
  call symput('num_x',_N_);
run;

%do i=1 %to &num_x;
  %let temp=%scan(&xvars,&i,%str( ));
  %if %length(&temp)>&name_len %then %do;
    %let name_len=%length(&temp);
  %end;
%end;

 /***********************************************************
 /  Plot X-loadings for each PLS component                  /
 ***********************************************************/

%do i=1 %to %eval(&max_lv - 1);

  %let j=%eval(&i+1);

  data pltanno;             *** Annotation Data Set for Plot ***;
    length text $ &name_len;
    retain function 'label' position '5' hsys '3' xsys '2' ysys '2' ;
    set &ds;
    text=%str(_name_); x=p&i; y=p&j;
  run;

  axis1 label=(angle=270 rotate=90 "X loading &j")
        major=(number=5) minor=none;

  axis2 label=("X-loading &i") minor=none;

  proc gplot data=&ds;
    plot p&j*p&i/anno=pltanno vaxis=axis1 haxis=axis2 frame;
    symbol1 v=none i=none;
  run;

%end;

%mend;

%macro getyload(dsoutmod,
                dsyload=yloads);

 /***********************************************************
 / Gets Y-loadings q from OUTMODEL data set:                /
 /    1.  Gets appropriate section of OUTMODEL data set.    /
 /    2.  Transposes it so the q's are column vectors.      /
 /    3.  Renames the columns to q1 - qA, where A is the    /
 /        number of latent variables in the final model.    /
 /  Variables:                                              /
 /      DSOUTMOD - Name of the OUTMODEL data set produced   /
 /                 by proc PLS.                             /
 /      DSYLOAD -  Name of the data set to contain the      /
 /                 Y-loadings as variables.                 /
 /  Required Global Variable:                               /
 /      %YVARS -   Macro variable containing the names of   /
 /                 the Y-variables in a string with space   /
 /                 delimiters.                              /
 ***********************************************************/

data &dsyload; set &dsoutmod(keep=_TYPE_ _LV_ &yvars);
  if _TYPE_='PQ' then output;

proc transpose data=&dsyload out=&dsyload; run;

data &dsyload; set &dsyload;
  if _NAME_='_LV_' then delete;
run;

%do i = 1 %to &lv;

  data &dsyload; set &dsyload;
    rename col&i=q&i;
  run;

%end;

%mend;

%macro plt_y_lv(dsoutmod);

 /***********************************************************
 /  Plots Y-loadings for each Y-variable versus the PLS     /
 /  component.                                              /
 /  Variable:                                               /
 /      DSOUTMOD - The OUTMODEL data set from proc PLS.     /
 /  Required Global Variable:                               /
 /      %YVARS -   Macro variable containing the names of   /
 /                 the Y-variables in a string with space   /
 /                 delimiters.                              /
 ***********************************************************/

data dsyload; set &dsoutmod(keep=_TYPE_ _LV_ &yvars);
  if _TYPE_='PQ' then output;

goptions reset=symbol;

axis1 label=(angle=270 rotate=90 'Y loading')
      major=(number=5) minor=none;
axis2 label=('PLS Component') order=(1 to &lv by 1) minor=none;

proc gplot data=dsyload;
  plot (&yvars)*_LV_/overlay legend vaxis=axis1 haxis=axis2
                     vref=0 lvref=2 frame;
run;

%mend;

%macro pltyload(ds,
                max_lv=&lv);

 /************************************************************
 /  Plots Y-loadings for a given number of PLS components    /
 /  vs. those of the preceding PLS component.                /
 /  Variables:                                               /
 /      DS -       Name of the data set containing the       /
 /                 loadings as variables q1-qA, where A=LV,  /
 /                 the number of PLS components, and a       /
 /                 character variable _NAME_ containing the  /
 /                 Y-variable names.                         /
 /      MAX_LV -   Number of the last PLS component to have  /
 /                 its loadings plotted.                     /
 /  Required Global Variable:                                /
 /      %YVARS -   Macro variable containing the names of    /
 /                 the Y-variables in a string with space    /
 /                 delimiters.                               /
 ************************************************************/

 /***********************************************************
 /  Determine the largest label to be put on plot           /
 ***********************************************************/

data _NULL_; set &ds;
  call symput('num_y',_N_);
run;

%let name_len=1;

%do i=1 %to &num_y;
  %let temp=%scan(&yvars,&i,%str( ));
  %if %length(&temp)>&name_len %then %do;
    %let name_len=%length(&temp);
  %end;
%end;

 /***********************************************************
 /  Plot Y-loadings for each PLS component                  /
 ***********************************************************/

%do i=1 %to %eval(&max_lv-1);

  %let j=%eval(&i+1);

  data pltanno;                *** Annotation Data Set for Plot ***;
    length text $ &name_len;
    retain function 'label' position '5' hsys '3' xsys '2' ysys '2' ;
    set &ds;
    text=%str(_NAME_); x=q&i; y=q&j;
  run;

  axis1 label=(angle=270 rotate=90 "Y loading &j")
        major=(number=5) minor=none;

  axis2 label=("Y-loading &i") minor=none;

  proc gplot data=&ds;
    plot q&j*q&i/anno=pltanno vaxis=axis1 haxis=axis2;
    symbol1 v=none i=none;
  run;

%end;

%mend;

%macro pltxywts(dsxwts,
               dsyloads,
               norm=1,
               max_lv=&lv);

 /************************************************************
 /  Plots X-weights and Y-loadings on the same axes for a    /
 /  given number of PLS components vs. those of the          /
 /  preceding PLS component.                                 /
 /  Variables:                                               /
 /      DSXWTS -   Name of the data set containing the X-    /
 /                 weights as variables w1-wA, where A=LV,   /
 /                 the number of PLS components, and a       /
 /                 character variable _NAME_ containing the  /
 /                 X-variable names.                         /
 /      DSYLOADS - Name of the data set containing the Y-    /
 /                 loadings as variables q1-qA, where A=LV,  /
 /                 the number of PLS components, and a       /
 /                 character variable _NAME_ containing the  /
 /                 Y-variable names.                         /
 /      NORM -     Indicates whether to normalize the        /
 /                 X-weights and Y-loadings or not.          /
 /                 (1=Yes, 0=No)                             /
 /      MAX_LV -   Number of the last PLS component to have  /
 /                 its weights plotted.                      /
 /  Required Global Variables:                               /
 /      %XVARS -   Macro variable containing the names of    /
 /                 the X-variables in a string with space    /
 /                 delimiters.                               /
 /      %YVARS -   Macro variable containing the names of    /
 /                 the Y-variables in a string with space    /
 /                 delimiters.                               /
 ************************************************************/

 /***********************************************************
 /  Determine the largest label to be put on plot           /
 ***********************************************************/

data _NULL_; set &dsxwts;
  call symput('num_x',_N_);
run;

data _NULL_; set &dsyloads;
  call symput('num_y',_N_);
run;

%let name_len=1;

%do i=1 %to &num_x;
  %let temp=%scan(&xvars,&i,%str( ));
  %if %length(&temp)>&name_len %then %do;
    %let name_len=%length(&temp);
  %end;
%end;

%let nameleny=1;

%do i=1 %to &num_y;
  %let temp=%scan(&yvars,&i,%str( ));
  %if %length(&temp)>&nameleny %then %do;
    %let nameleny=%length(&temp);
  %end;
%end;

%if &name_len < &nameleny %then %let name_len = &nameleny;

 /***********************************************************
 /  Normalize weights if desired                            /
 ***********************************************************/

%if %eval(&norm) %then %do;

proc iml;
  use &dsxwts;
  read all var ("w1":"w&max_lv") into W;
  use &dsyloads;
  read all var ("q1":"q&max_lv") into Q;
  W=W#sqrt(1/W[##,]);  *** Normalize X-weights ***;
  Q=Q#sqrt(1/Q[##,]);  *** Normalize Y-loadings ***;
  w_col=("WQ1":"WQ&max_lv");
  _NAME_={&xvars};
  create dsxwts from W[colname=w_col rowname=_NAME_];
  append from W[rowname=_NAME_];
  q_col=("WQ1":"WQ&max_lv");
  _NAME_={&yvars};  
  create dsyloads from Q[colname=q_col rowname=_NAME_];
  append from Q[rowname=_NAME_];
quit;  
  
%end;

%else %do;

data dsxwts; set &dsxwts;

data dsyloads; set &dsyloads;

%end;


 /***********************************************************
 /  Plot X-weights and Y-loadings for each PLS component    /
 ***********************************************************/

%if &name_len>&nameleny %then %do;
  data ds; set dsxwts dsyloads;
%end;
%else %do;
  data ds; set dsyloads dsxwts;
%end;

%do i=1 %to %eval(&max_lv-1);

  %let j=%eval(&i+1);

  data wt_anno;              *** Annotation Data Set for Plot ***;
    length text $ &name_len;
    retain function 'label' position '5' hsys '3' xsys '2' ysys '2' ;
    set ds;
    text=%str(_name_); x=wq&i; y=wq&j;
  run;

  axis1 label=(angle=270 rotate=90 "Component &j Weight")
        major=(number=5) minor=none;

  axis2 label=("Component &i Weight") minor=none;

  proc gplot data=ds;
    plot wq&j*wq&i/anno=wt_anno vaxis=axis1 haxis=axis2 frame
                              vref=0 href=0;
    symbol1 v=none i=none;
  run;

%end;

%mend;

%macro get_bpls(dsoutmod,
                dsout=bpls);

 /************************************************************
 /  Gets B(PLS), the matrix of PLS regression coefficients   /
 /  of Y on X.  For each Y, the values represent the         /
 /  importance of each X-variable in the modeling of the     /
 /  corresponding Y-variable.                                /
 /  Variables:                                               /
 /      DSOUTMOD - Name of the OUTMODEL data set produced    /
 /                 by proc PLS.                              /
 /      DSOUT -    Name of the data set to contain the       /
 /                 regression coefficients, with the         /
 /                 variables representing columns in         /
 /                 B(PLS), and one variable naming the       /
 /                 X-variable for each row of B(PLS).        /
 /  Required Global Variables:                               /
 /      %XVARS -   Macro variable containing the names of    /
 /                 the X-variables in a string with space    /
 /                 delimiters.                               /
 /      %YVARS -   Macro variable containing the names of    /
 /                 the Y-variables in a string with space    /
 /                 delimiters.                               /
 /      %NUM_Y   - Number of response variables Y.           /
 ************************************************************/

data est_wb; set &dsoutmod;  if _TYPE_='WB' then output; run;
data est_pq; set &dsoutmod;  if _TYPE_='PQ' then output; run;

proc iml;
  use est_wb;
  read all var {&xvars} into w_prime;
  read all var {_Y_} into b;
  use est_pq;
  read all var {&xvars} into p_prime;
  read all var {&yvars} into q_prime;
  W=w_prime`;
  P=p_prime`;
  Q=q_prime`;
  B_PLS = W*inv(P`*W)*diag(b)*Q`;
  b_col=('B1':"B&num_y");
  x_var={&xvars};
  create &dsout from B_PLS[colname=b_col rowname=x_var];
  append from B_PLS[rowname=x_var];
quit;

%mend;

%macro plt_bpls(ds);

 /***********************************************************
 /  Plot the PLS predictor (regression) coefficients in     /
 /  B(PLS) vs. the frequency, for each response variable.   /
 /  Variables:                                              /
 /      DS -       Data set containing the columns of       /
 /                 B(PLS) as variables, as well as a        /
 /                 variable for the frequency.              /
 /  Required Global Variable:                               /
 /      %NUM_Y   - Number of response variables Y.          /
 ***********************************************************/

data &ds; set &ds;
  f=_n_;
run;

%let plotvars=%str( );

%do i=1 %to &num_y;
  %let plotvars=%str(&plotvars b&i);
%end;

axis1 label=(angle=270 rotate=90 'Coefficient')
      major=(number=5) minor=none;
axis2 label=('Frequency') minor=none;

proc gplot data=&ds;
  plot (&plotvars)*f / overlay legend vaxis=axis1 haxis=axis2
                       vref=0 lvref=2 frame;
  symbol1 v=none i=spline;
run;

%mend;

%macro get_vip(dsoutmod,
               dsvip=vip_data);

 /************************************************************
 /  Calculate VIP:  Variable Importance for the Projection.  /
 /  This represents the importance of each X-variable in     /
 /  the PLS modeling of both the X- and Y-variables.         /
 /  Variables:                                               /
 /      DSOUTMOD - Name of the OUTMODEL data set produced    /
 /                 by proc PLS.                              /
 /      DSVIP -    Name of the data set to contain the       /
 /                 variable named 'VIP' and the names of     /
 /                 X-variables.                              /
 /  Required Global Variables:                               /
 /      %XVARS -   Macro variable containing the names of    /
 /                 the X-variables in a string with space    /
 /                 delimiters.                               /
 /      %YVARS -   Macro variable containing the names of    /
 /                 the Y-variables in a string with space    /
 /                 delimiters.                               /
 ************************************************************/

data dsxwts; set &dsoutmod(keep=_TYPE_ _LV_ &xvars);
  if _TYPE_='WB' then output;

data y_rsq; set &dsoutmod(keep=_LV_ _TYPE_ &yvars _Y_);
  if _TYPE_='V' then output;
  drop _TYPE_;
run;

data y_rsq; merge y_rsq dsxwts; by _LV_;
  if _LV_=0 then delete;
run;

proc iml;
  use y_rsq;
  read all var {_Y_} into rsq_y;
  read all var {&xvars} into w_prime;
  A=nrow(rsq_y);
  K=ncol(w_prime);
  W=w_prime`;
  Wnorm=W#(1/sqrt(W[##,]));
  if A > 1 then do;
    part_rsq=rsq_y-(0//rsq_y[1:(A-1),]);
    tot_rsq=rsq_y[A,];
    vip_sq=((Wnorm##2)*part_rsq)#(K/tot_rsq);
    VIP=sqrt(vip_sq);
  end;
  else VIP=Wnorm#sqrt(K); 
  x_var={&xvars};
  create &dsvip from VIP[colname='VIP' rowname=x_var];
  append from VIP[rowname=x_var];
quit;

%mend;

%macro plot_vip(ds);

 /************************************************************
 /  Plot the VIP:  Variable Importance for the Projection.   /
 /  Variables:                                               /
 /      DS -       Data set containing the frequencies       /
 /                 the VIP for each frequency.               /
 ************************************************************/

data &ds; set &ds;
  f=_N_;
run;

axis1 label=(angle=270 rotate=90 'VIP')
      major=(number=10) minor=none;
axis2 label=('Frequency') minor=none;

proc gplot data=&ds;
  plot vip*f / overlay vaxis=axis1 haxis=axis2 vref=0.8 lvref=2
               frame;
  symbol1 v=none i=join;
run;

%mend;


%macro get_dmod(dsoutput,
                dsdmod=dmod,
                qresname=qres,
                id=n);

 /************************************************************
 /  Calculate the distance from each data point to the model /
 /  in both the X-space (DMODX) and in the Y-space (DMODY).  /
 /  Variables:                                               /
 /      DSOUTPUT - OUTPUT data set from proc PLS.            /
 /      DSDMOD -   Data set to contain the distances to      /
 /                 the model.                                /
 /      QRESNAME - Suffix of variable names for XQRES and    /
 /                 YQRES assigned by the user in the         /
 /                 proc PLS OUTPUT statement.                /
 /      ID -       Observation identification variable       /
 /                 in input data set.                        /
 /  Required Global Variable:                                /
 /      %LV -      Number of Latent Variables (PLS factors)  /
 ************************************************************/
  
data trn_out; set &dsoutput;
  if y&qresname ^= . then output;
run;

proc means data=trn_out noprint;
  var xqres;
  output out=outmeans n=n mean=xqres_mn;
run;

data _NULL_; set outmeans;
  call symput('num_trn',n);
  call symput('xqres_mn', xqres_mn);
run;

proc iml;
  use &dsoutput;
  read all var {x&qresname} into xqres;
  read all var {y&qresname} into yqres;
  read all var{&id} into id;
  dmodx=sqrt(xqres/&xqres_mn);
  
  do i=1 to nrow(xqres);
    if yqres[i]=. then 
       dmodx[i]=dmodx[i]/sqrt(&num_trn/(&num_trn-&lv-1));
  end;
  
  dmody=sqrt(yqres*(&num_trn/(&num_trn-&lv-1)));
  dmodboth=id||dmodx||dmody;
  col={&ID DMODX DMODY};
  create &dsdmod from dmodboth[colname=col];
  append from dmodboth;
quit;
  
%mend;

%macro cont_scr(est,
                out,
                dsout,
                obsnum,
                idvar=n,
                a=1);

  /************************************************************
  /  Calculates and plots the contribution of each X-variable /
  /  to an X-score for a particular observation.  This is     /
  /  useful in diagnosing the cause of an outlying X-score.   /
  /      EST -      Name of the OUTMODEL data set from        /
  /                 PROC PLS containing the weights w1-wA     /
  /                 where A=LV, the number of PLS components. /
  /      OUT -      Name of the OUTPUT data set containing    /
  /                 the X-scores as variables.                /
  /      DSOUT -    Name of the data set to be generated by   /
  /                 the macro containing the contributions    /
  /                 for each X-variable.                      /
  /      OBSNUM -   Number of the observation in              /
  /                 question.                                 /
  /      IDVAR -    Name of variable in the OUT data set      /
  /                 numbering the data points.  The default   /
  /                 is N.                                     /
  /      A -        PLS component associated with the         /
  /                 X-score of interest. The default is 1.    /
  /  Required Global Variables:                               /
  /      %XVARS -   Macro variable containing the names of    /
  /                 the X-variables in a string with space    /
  /                 delimiters.                               /
  /      %LV -      Number of Latent Variables (PLS factors)  /
  ************************************************************/

data est; set &est;
  if (_TYPE_='WB' or _TYPE_='PQ') then output;
run;

data out; set &out;
  if &idvar=&obsnum then output;
run;

proc iml;
  use est;
  read all var {&xvars} into WP;
  W=WP[1:&lv,];
  P=WP[(&lv+1):(2*&lv),];
  Wstar=W`*inv(P*W`);
  use &out;
  read all var {&xvars} into X;
  use out;
  read all var {&xvars} into x_i;
  contrib=(Wstar[,&a])` # (x_i - X[:,]);
  quantity=('contrib');
  xvar={&xvars};
  create &dsout from contrib[rowname=quantity colname=xvar];
  append from contrib[rowname=quantity];
quit;

proc transpose data=&dsout out=&dsout; run;

data &dsout; set &dsout;
  rename col1=contrib;
run;

axis1 label=(angle=270 rotate=90 'Contribution');

axis2 label=('X-variable');

proc gplot data=&dsout;
  plot contrib * _NAME_ / haxis=axis2 vaxis=axis1;
  symbol1 i=needles v=dot;
run;
quit;

%mend;

%macro cont2scr(est,
                out,
                dsout,
                obsnum,
                idvar=n,
                a1=1,
                a2=2);

  /************************************************************
  /  Calculates and plots the contribution of each X-variable /
  /  to a pair of X-scores for a particular observation.      /
  /  This is useful in diagnosing the cause of an outlying    /
  /  pair of X-scores.                                        /
  /      EST -      Name of the OUTMODEL data set from        /
  /                 PROC PLS containing the weights w1-wA     /
  /                 where A=LV, the number of PLS components. /
  /      OUT -      Name of the OUTPUT data set containing    /
  /                 the X-scores as variables.                /
  /      DSOUT -    Name of the data set to be generated by   /
  /                 the macro containing the contributions    /
  /                 for each X-variable.                      /
  /      OBSNUM -   Number of the observation in              /
  /                 question.                                 /
  /      IDVAR -    Name of variable in the OUT data set      /
  /                 numbering the data points.  The default   /
  /                 is N.                                     /
  /      A1 -       PLS component associated with the first   /
  /                 X-score of interest.  The default is 1.   /
  /      A2 -       PLS component associated with the second  /
  /                 X-score of interest.  The default is 2.   /
  /  Required Global Variables:                               /
  /      %XVARS -   Macro variable containing the names of    /
  /                 the X-variables in a string with space    /
  /                 delimiters.                               /
  /      %XSCRNAME - Prefix given to X-scores when OUT= data  /
  /                  set was defined.                         /
  /      %LV -      Number of Latent Variables (PLS factors)  /
  ************************************************************/

data est; set &est;
  if (_TYPE_='WB' or _TYPE_='PQ') then output;
run;
 
data out; set &out;
  if &idvar=&obsnum then output;
run;

proc iml;
  use est;
  read all var {&xvars} into WP;
  W=WP[1:&lv,];
  P=WP[(&lv+1):(2*&lv),];
  Wstar=W`*inv(P*W`);  
  use &out;
  read all var {&xscrname&a1 &xscrname&a2} into T;
  read all var {&xvars} into X;
  use out;
  read all var {&xscrname&a1 &xscrname&a2} into t_i;
  read all var {&xvars} into x_i;
  delta_t1=t_i[,1]-T[:,1];
  delta_t2=t_i[,2]-T[:,2];
  sd_t1=sqrt((T[##,1]-nrow(T)*T[:,1]**2)/(nrow(T)-1));
  sd_t2=sqrt((T[##,2]-nrow(T)*T[:,2]**2)/(nrow(T)-1)); 
  w1star=Wstar[,&a1];
  w2star=Wstar[,&a2];
  v_sq=(delta_t1/sd_t1)**2*(w1star)`##2+
       (delta_t2/sd_t2)**2*(w2star)`##2;
  v=sqrt(v_sq);
  delta_x=x_i-X[:,];
  contrib=(v#delta_x);
  quantity=('contrib');
  xvar={&xvars};
  create &dsout from contrib[rowname=quantity colname=xvar];
  append from contrib[rowname=quantity];
quit;

proc transpose data=&dsout out=&dsout; run;

data &dsout; set &dsout;
  rename col1=contrib;
run;

axis1 label=(angle=270 rotate=90 'Contribution');

axis2 label=('X-variable');

proc gplot data=&dsout;
  plot contrib * _NAME_ / haxis=axis2 vaxis=axis1;
  symbol1 i=needles v=dot;
run;
quit;

%mend;