/*
Filename:        BraunDIV.sas
Purpose:         SAS code to read Braun's data on stylistic element frequencies by thickness class
                 and region.  Thickness classes are grouped to lessen sample size variation.
				 Then compute diversity stats;
Last Update: FDN 4.8.02
*/

/* use IML to compute the diversity stats */
proc iml;
**************************************************************************
Function:       theta_solve
Description:iterative estimation of theta -- t_e (Ewens 1972, Neiman 1995)
Arguments:      n:              Sample size
                        s:      Number of types
Returns:        est:    a single theta estimate -- t_e;
start theta_solve(s,n);
        maxiter=1000 ;                          *maximum number of iterations;
        res= .0001   ;                          *convergence criterion;
        valuevec=j(maxiter,1,0);
        target=s;
        est=1;
        inc=.1;
        valuevec[1,]= sum((est#j(1,n,1))/(est+(0:(n-1))));   /*equation to solve*/
        i=2;
        do until ((i>maxiter) | abs((valuevec[i-1,]-target))<res);
        valuevec[i,]= sum((est#j(1,n,1))/(est+(0:(n-1))));
        if abs(valuevec[i-1,]-target)< abs(valuevec[i,]-target)  then
                 inc=-inc;       /*are we going in the right direction?*/
        if ((valuevec[i-1,]-target)>0)^= ((valuevec[i,]-target)>0)  then
                 inc= -(inc/2);      /*has an over shoot occurred?*/
        est=est+inc;
        i=i+1;
end;
if abs((valuevec[i-1,]-target))> res then est=.;
return(est);
finish theta_solve;
********************************************************************************;
********************************************************************************
Function:       theta_vec
Description: Turns function THETA_SOLVE loose on vectors of aguments and returns
                        a vector of estimates
Arguments:      Total:     a column vector of Sample sizes
                        n_types:   a column vector of Number of types (richness) in each sample
Requires:       Function theta_solve
Returns:        Theta_EST_VEC: a Vector of theta estimates;

start theta_vec(n_types,total);
        theta_est_vec=j(nrow(n_types),1,0);
do i=1 to nrow(n_types) by 1;
    theta_est_vec[i,]= theta_solve(n_types[i,],total[i,]);
end;
return (theta_est_vec);
finish theta_vec;
*********************************************************************************;
**********************************************************************************
Function:       diversity
Description: Computes diversity statistics from an input assmblage*type matrix
Arguments:
Returns:        h_shannon: a column vector of shannon diversity
                        e_shannon: a column vector of shannon evenness
                        d_simpson: a column vector of simpson diversity
                        e_simpson: a column vector of simpson evenness
                        t_f:       a column vector of t_f estimates (Neiman 1995)
                        total:     a column vector of Sample sizes
                        n_types:   a column vector of Number of types (richness) in each sample;

start diversity(INMAT,h_shannon,e_shannon,d_simpson,e_simpson,t_f,total,n_types);
nobs=nrow(inmat); nvars=ncol(inmat);
sums=j(nobs,nvars,0);                      * dimension  ;
h_shannon=j(nobs,1,0);             * dimension a vector for shannon ;
total=inmat[,+];                   * sample sizes;
n_types=(inmat>0)[,+];             * number of types (species, alleles);
  do i=1 to nvars;
    sums[,i]=total;
  end;
prop=inmat/sums;          *proportions;
f=(prop#prop)[,+];        *isonymy;
d_simpson=1/f  ;          *effective number of alleles, Simpson's D;
e_simpson = d_simpson#(1/n_types); * evenness for simpson's D;
t_f= d_simpson-1;         * estimate of theta_f (Neiman 1995) ;
do i=1 to nobs;
  x= prop[i,loc(prop[i,])];             *error trap for division by 0;
  h_shannon[i,]=-(sum(x#log(x)));       *shannon's H;
end;
e_shannon= h_shannon/log(total);        *evenness for shannon's H;
finish diversity;
************************************************************************;
/* OK let's read the data and do the business */

use braun1;
read all into inmat ;          *read in the numeric data;
read all var {group tclass};   *read in the character variable for labels;

call diversity(INMAT,h_shannon,e_shannon,d_simpson,e_simpson,t_f,total,n_types);
t_e=theta_vec(n_types,total);

create div var {group tclass h_shannon e_shannon d_simpson e_simpson
t_f t_e total n_types };
append;
close div ;

proc print;

proc plot;
plot h_shannon*tclass= '*' $ group; 
run;