I have n forests that I want to compute Moristia index. I have the code for one forest but I do not know how to repeat the program n times. The code for Morisita Index is the following:
data data;
infile "B:\sas\data.csv" delimiter=',' firstobs=2; * read the csv data with coordinates of the trees and multiple forests;
input forestid $ x y diameter;
if forestid ne 1 then delete; * select only one forest;
run;
*start Morosita's index computations;
proc iml;
/** Count the number of points in each cell of a 2D grid.
Input: p - a kx2 matrix of 2D points
XDiv - a subdivision of the X coordinate direction
YDiv - a subdivision of the Y coordinate direction
Output: The number of points in p that are contained in
each cell of the grid defined by XDiv and YDiv. **/
start Bin2D(p, XDiv, YDiv);
x = p[,1]; y = p[,2];
Nx = ncol(XDiv); Ny = ncol(YDiv);
/** Loop over the cells. Do not loop over observations! **/
count = j(Ny-1, Nx-1, 0); /** initialize counts for each cell **/
do i = 1 to Ny-1; /** for each row in grid **/
idx = loc((y >= YDiv) & (y < YDiv[i+1]));
if ncol(idx) > 0 then do;
t = x[idx]; /** extract point in this row (if any) **/
do j = 1 to Nx-1; /** count points in each column **/
count[i,j] = sum((t >= XDiv) & (t < XDiv[j+1]));
end;
end;
end;
return( count );
finish;
use data ; *bring external data into IML environment (SAS data but not data curently inside IML procedure);
read all var {x y} into p[colname=varNames]; *place your data in the matrix p that can be used by bining module;
/** define bins: 50 columns and 50 rows on [0,1000]x[0,1000] **/
cellsize=100;
XMin = 0; XMax = 1000; dx = cellsize;
YMin = 0; YMax = 1000; dy = cellsize;
XDiv = do(XMin, XMax, dx);
YDiv = do(YMin, YMax, dy);
count = Bin2D(p, XDiv, YDiv);
*create the dataset for MORISITA Index of Dispersion (1959);
create quad from count;
append from count;
close quad;
run;quit;
*MORISITA Index of Dispersion (1959);
data quad;
set quad;
index=catx("_", "Row",_n_); *create a new varaibale that concatenate the word Row with the row counting;
run;
*prepare data for computing Morisita index;
proc sort data=quad;
by index;
run;
*create one column from the matrix counting the points inside bins;
Proc transpose data=quad out=mor (rename=(col1=count));
by index;
run;
*start computing Morisita index;
data mor;
set mor;
c=count*(count-1);
run;
proc means data=mor noprint;
output out=morisita sum(c)=sumc sum(count)=points ;
run;quit;
data morisita;
set morisita;
MInd=sumc*_freq_/(points*(points-1));
pvalue=1-probchi((mind*(points-1)+_freq_-points)/(_freq_-1),_freq_-1);
run;
proc print data=morisita;
run; quit;
What a want is to have the Morisita's index computed for all forests and saved in one dataset. I tried macro, as they seems to be recommended to do repetitive tasks, but the tutorials and helps are designed for so simple tasks that I did not know why one should bother using them, while the complicated programs are so complicated that are useless. I was expecting some examples for this kind of situations, as it seems to be encountered very often in practice.
Any advice?
