/*
Filename: 	Site7CA.sas
Purpose:	SAS code to do a CA on ceramic type counts from plowzone quadrats
			at Site 7, Monticello. First we smooth the data. Then do the CA on the 
			smoothed counts. Then compute the meanceramic dates for the smoothed 
			data and compare them against the CA dimension 1 scores. 
Last Update: FDN 4.29.02
*/

title 'CA of Site 7 ceramic assemblages';

data TypeCounts; set Site7PlowZone;
drop northing easting ;


data coords; set Site7PlowZone;
keep  prov northing easting;


proc plot data=coords;
plot northing * easting = '*' $ prov;

run;

proc iml;
**********************************************************************************
Function: 	 Smoothit
Description: Uses a moving average flter to smooth artifact counts  
Arguments:	 c:			an n X 2 matrix of northing and easting coordinates 
						for the centers of the sample locations.
			 x:			a sample location x type matrix of artifact counts. 
			 range:		The distance over which the counts shuld be smoothed
			 method:	"WMA" for a weighted moving average filter with bisquare weights
						"MA" for a simple moving average filter. 
			 treatment:	"Counts" to smooth the counts.
						"Percent" to smooth percentages computed from the counts.
Returns:	 SMOOTHX: 	A matrix ofsmoothed counts.
***********************************************************************************;
start Smoothit(c,x,range,method,treatment);
n=nrow(x);
cols= ncol(x);
wtmat= j(n,n,{0}); 				*dimension weights matrix;
if treatment="Percent" then  do;	*compute the percentages;
      do i={1} to n;
	     if x[i,+]  >  0 then total=x[i,+];
		 	else total=1;
         x[i,]= x[i,] / j(1,cols,total);
         end;
      end;
if method="WMA" then  do;
    do j={1} to n by{ 1};  *get distances between data points;       			
       do k={1} to n by{ 1};
          dist= ((( c[k,]- c[j,])*( c[k,]- c[j,])`)##{.5});
          if ( dist < (range+{.001})) then
               wtmat[j,k]= (1-(dist/range)##2)##2;
       end;
     end;
end; 
if method="MA" then  do;
     do j={1} to n by{ 1}; *get  distances between data points;      
         do k={1} to n by{ 1};
            dist= ((( c[k,]- c[j,])*( c[k,]- c[j,])`)##{.5});
            if ( dist < (range+{.001})) then  wtmat[j,k]={1};
         end;
      end;
end;
smoothx=wtmat*x;	 *apply the weights;  
*print c x dist wtmat smoothx ;	  *uncomment this line to see intemediate results';
return(SMOOTHX);
finish Smoothit;
*********************************************************************************;

/* OK let's read the data and do the business */
use TypeCounts;                  
read all into x [colname=vars rowname=prov];
use coords;                    
read all into c;	   * coords contains the quadrat northings and eastings;

smoothx=Smoothit(c,x,40,"WMA","Counts");
 
create smoothx  from smoothx [rowname=prov colname=vars ];
append from smoothx [rowname=prov];
close smoothx;

proc print;


proc corresp data =smoothx out=ca short  dimens=3;
var _numeric_ ;
id prov;

*plot the results;

proc sort;
  by _type_ inertia;

data ca_obs;
  set ca;
  if _type_ eq 'OBS';
  rank=_n_;
proc plot;
  plot inertia*rank = '*' $ prov;
  plot sqcos1*sqcos2 = '*' $  prov;
  plot sqcos3*sqcos4 = '*' $  prov;


data ca_vars(rename=(prov=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= '*' $  prov ; * HAXIS=BY .1 VAXIS=BY .1;
 * plot dim3*dim4= '*' $  prov ; * HAXIS=BY .1 VAXIS=BY .1;
  *plot dim2*dim3= '*' $  prov;

proc plot data=ca;
  where _type_ eq 'OBS';
  plot dim2*dim1= '*' $  prov ; *  HAXIS=BY .1 VAXIS=BY .1;
 * plot dim3*dim4= '*' $  prov ; *  HAXIS=BY .1 VAXIS=BY .1;
  *plot dim2*dim3= '*' $  prov;


run;

/* 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=smoothx;
by prov;
proc transpose data=smoothx out=smoothx1(rename=(col1=count)) name=type;
by prov;
proc sort;
by type prov;

/*Read in the median dates */
data mediandates; 
length type $ 20;
input Type $	MedianDate;
cards;
PearlwarePolyPainted	1812.5
WhiteSaltGlaze			1765
Delft					1762.5
Whieldon				1757.5
Cauliflower				1767
Creamware				1791
Pearlware				1804.5
PearlwareAnnular		1796
PearlwareBluePainted	1804.5
PearlwareTP				1812.5
ChinesePorcelain		.
Astbury					1737.5
Jackfield				1770
NMidSlip				1745
OtherSlip				.
BlackGlazed				.
CoarseEW				.
WesterwalSW				1757.5
BritBrownSW				1757.5
NottSW					1775
PearlwareShell			1807.5
;

proc print;

run;

proc sort;
by type;

data both; 
merge smoothx1 mediandates;
by type;

proc print;

proc summary data=both;
class prov;
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 prov;
proc sort data=meandates;
by prov;
data dates; merge ca_obs meandates;
by prov;

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

proc corr spearman;
var meanceramicdate; with dim1;
run;