/* This code performs a two block partial least squares analysis, configured here for 
two blocks of shape coordinates (i.e., a singualar warps analysis). Unlike other routines 
in this series, this one is set up for a different dataset. The sample dataset was divided 
into two blocks of data and those were separately aligned by GPA. They can be read in 
separately in this routine from the new data files (block1.dta and block2.dta) provided. 
Results of this program include singular values, singular vectors and singular warps scores, 
which can be plotted separately in Interactive Data Analysis. */

data one; 
infile 'C:\...\block1.dta'
firstobs=2 obs=41; input cat $; run;

data block1; set one; 
infile 'C:\...\block1.dta'
firstobs=43; input x1-x30; run;

data two; 
infile 'C:\...\block2.dta'
firstobs=2 obs=41; input cat $; run;

data block2; set two; 
infile 'C:\...\block2.dta'
firstobs=43; input x1-x18; run;

proc iml; 
DF=7;
use block1; read all into X1; use block2; read all into X2;
N=nrow(X1); LM1=ncol(X1); LM2=ncol(X2);

if LM1<LM2 then do; 
	print '*****  The number of variables in block1 MUST be greater   *****',
	      '*****  than or equal to the number of variables in block2  *****',,
	      'Switch data blocks and re-run.'; 
	end;

else do;

lab1='B1SW1':'B1SW18'; lab2='B2SW1':'B2SW18'; 		
labSW=lab1||lab2;

meanX1=X1[+,]/N; X1bar=shape(meanX1,N,LM1); Dif1=X1-X1bar;
meanX2=X2[+,]/N; X2bar=shape(meanX2,N,LM2); Dif2=X2-X2bar;
CCV=(Dif1`*Dif2)/(N-1);

call SVD(u,q,v,CCV);

q=q[1:LM2-DF]; q2=q##2; var=q2[+,]; 
pcv=j(nrow(q),1); SW=j(nrow(q),1); 
do i=1 to nrow(q); SW[i,]=i; pcv[i,1]=q2[i,]*100/var; end;

swsX1 = Dif1*u;	swsX2 = Dif2*v;
sws=swsX1||swsX2; r=corr(sws);
r=r[1:LM2-DF,LM2+1:2#LM2-DF]; 
r=vecdiag(r); 

labstat={'SW' 'sing_val' '%_covar' 'corr. (r)'};
stat=SW||q||pcv||r;
reset autoname;
print stat[colname=labstat];

create stats from stat[colname=labstat]; append from stat;
create sws from sws[colname=labSW]; append from sws; 
create B1vectors from u; append from u;
create B2vectors from v; append from v;
end;
abort;

data singwarps; merge one sws; 
sps=substr(cat,1,2); sex=substr(cat,4,1); 
run;