/* 	Filename: 	LosMuertosCA.sas
	Purpose:	CA and MCD analysis of Type counts from Los Muertos Pueblo. 
				We also sort the data on the CA Dim 1 scores, so that they
				can be output for plotting.  Finally we do a PCA to compare 
				with the CA.
				Data from Andrew I. Duff 1996, Am Antiquity 61(1):89-101.
	Last Update: FDN 4.28.02
*/

/* import the counts */
PROC IMPORT OUT= WORK.losmuertos 
            DATAFILE= "f:\courses\anth588\LosMuertos.xls" 
            DBMS=EXCEL2000 REPLACE;
     GETNAMES=YES;
RUN;


data losmuertos; set losmuertos;
unitlevel=left(Unit||'_'||level);
drop unit level;

/* do the CA */
proc corresp data=LosMuertos out=ca dim=3;
var HeshotauthlaPoly  KwakinaPoly   PinedalePoly 
SpringervillePoly StJohnsPoly  TularosaBW;
id unitlevel;
run;

proc sort;
  by _type_ inertia;

data ca_obs;
  set ca;
  if _type_ eq 'OBS';
  rank=_n_;
proc plot;
  plot inertia*rank = '*' $ unitlevel;
  plot sqcos1*sqcos2 = '*' $  unitlevel;
  
data ca_vars(rename=(unitlevel=ware));
  set ca;
  if _type_ eq 'VAR' or _type_='SUPVAR';
  rank= _n_;

 proc plot;
   plot inertia*rank = '*' $  ware;
   plot sqcos1*sqcos2 = '*' $  ware;
  * plot sqcos3*sqcos4 = '*' $  ware;

proc plot data=ca vtoh=2;
  where _type_ eq 'VAR' or _type_='SUPVAR';
  plot dim2*dim1= '*' $  unitlevel ; * HAXIS=BY .1 VAXIS=BY .1;


proc plot data=ca;
  where _type_ eq 'OBS';
  plot dim2*dim1= '*' $  unitlevel ; *  HAXIS=BY .1 VAXIS=BY .1;
 
/* now compute the MCD's for the smoothed data and compare them to the CA Dim 1 scores */
/* start by transposing the data */

Proc sort data=losmuertos;
by unitlevel;
proc transpose data=losmuertos out=losmuertos1(rename=(col1=count)) name=type;
by unitlevel;
proc sort;
by type unitlevel;

/*Read in the median dates */
data mediandates; 
length type $ 20;
input Type $	MedianDate;
cards;
TularosaBW	1250
StJohnsPoly	1250
SpringervillePoly	1275
PinedalePoly	1300
KwakinaPoly	1337.5
HeshotauthlaPoly	1350
;

proc print;
run;

proc sort;
by type;

data both;      * merge the type median dates with the type counts;
merge losmuertos1 mediandates;
by type;

proc print;

proc summary data=both;
class unitlevel;
var mediandate;
weight count;
output out=meandates
	mean(mediandate)=MeanCeramicDate  ;
data meandates; set meandates;
if _type_ ne 0;
drop _type_ ;


proc sort data=ca_obs;
by unitlevel;
proc sort data=meandates;
by unitlevel;
data dates; merge ca_obs meandates;
by unitlevel;

proc plot;
plot meanceramicdate*dim1='*' $ unitlevel;
run;

proc corr spearman;
var meanceramicdate; with dim1;
run;
												  
/*merge the CA scores with the data and sort them for a seriation graph */
data out; merge ca_obs(keep=dim1 unitlevel) losmuertos;
by unitlevel; 
proc sort;
by dim1;
proc print; run;



/* convert counts to relative freqs for a PCA */
data LosMuertospct; set losMuertos;
array nums {*} HeshotauthlaPoly  KwakinaPoly   PinedalePoly 
			   SpringervillePoly StJohnsPoly  TularosaBW;
total=sum(of HeshotauthlaPoly  KwakinaPoly   PinedalePoly 
		SpringervillePoly StJohnsPoly  TularosaBW);
do i=1 to dim(nums);
	nums{i}= nums{i}/total;
end;
drop i total;


proc print;

/* do the PCA */
proc princomp data =losmuertospct out=pca outstat=eigenstructure covariance;
var HeshotauthlaPoly  KwakinaPoly   PinedalePoly 
SpringervillePoly StJohnsPoly TularosaBW;
run;


proc plot data=pca;
plot prin2*prin1 = '*' $ unitlevel;
run;


/* plot the eigenvectors */
data eigenvectors; set eigenstructure;
where _type_='SCORE';
proc transpose data = eigenvectors out=eigenvectors name=type; 
proc print data=eigenvectors;
proc plot data=eigenvectors;
plot prin2*prin1= '*' $ type;
run;