options nodate nonumber ps=55 ls=80;

title 'Kuehl -- Example 4.1 & Levene (Med) Test';

data crab;
do site=1 to 6;
do rep=1 to 25;
input count @@;
output;
end; end;
datalines;
0 0 22 3 17 0 0 7 11 11 73 33 0 65 13 44 20 27 48 104 233 81 22 9 2
415 466 6 14 12 0 3 1 16 55 142 10 2 145 6 4 5 124 24 204 0 0 56 0 8
0 0 4 13 5 1 1 4 4 36 407 0 0 18 4 14 0 24 52 314 245 107 5 6 2
0 0 0 4 2 2 5 4 2 1 0 12 1 30 0 3 28 2 21 8 82 12 10 2 0
0 1 1 2 2 1 2 29 2 2 0 13 0 19 1 3 26 30 5 4 94 1 9 3 0
0 0 0 2 3 0 0 4 0 5 4 22 0 64 4 4 43 3 16 19 95 6 22 0 0
;
run;

proc means mean std min max;
class site;
var count;
run;

proc glm;
class site;
model count=site;
output out=err_out r=ehat p=ybar;
run;

proc sort;
by ehat;
run;

data two;
set err_out;
i=_n_;
f_i=(i-0.5)/150;
nq_i=probit(f_i);
run;


proc plot;
plot ehat*nq_i;
title2 'Normal Probability Plot';
run;

proc plot;
plot ehat*ybar;
title2 'Plot of Residuals vs Site Means';
run;

proc sort data=two;
by site ehat;
run;

data median;
set two;
by site ehat;
retain median ord;
if first.site then do;
ord=0;
median=0;
end;
ord=ord+1;
if ord < 13 then median=median+0;
else if ord=13 then median=median+ehat;
else if ord > 13 then median=median;
if last.site then do;
output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output;
end;
keep site median;
run;
quit;

data three;
merge two median;
z=abs(ehat-median);
run;

proc sort;
by site;
run;

data zibar;
set three;
by site;
retain zibar;
if first.site then do;
zibar=0;
end;
zibar=zibar+z;
if last.site then do;
zibar=zibar/25;
output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output;
end;
keep site zibar;
run;


data zbar;
set three;
retain zbar 0;
zbar=zbar+z;
if _n_=150 then do;
zbar=zbar/150;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
output; output; output; output; output; output; output; output; output; output;
end;
keep site zbar;
run;
quit;

data four;
merge three zibar zbar;
mst=((zibar-zbar)**2)/5;
mse=((z-zibar)**2)/144;
run;


proc print;
var site ehat median z zibar zbar mst mse;
sum mst mse;
title2 'Sum of MST and MSE used in Levene test (F0=MST/MSE)';
run;

quit;