%macro boxcox(phimin=,phimax=, steps=, model=, class=, stmts=, data=, response=);

/***********************************************************/
/* This macro performs the Box-Cox-Power transformation    */
/* Box, G. E. P. and Cox, D. R. (1964) An analysis of      */
/* transformations (with discussion).                      */
/* J. R. Statist. Soc. B., 26, 211-246.                    */
/*                                                         */
/*   phimin   = smallest phi                               */
/*   phimax   = largest phi                                */
/*   steps    = number of steps in grid search             */
/*   class    = list of variables for class statment       */
/*   model    = list of terms in model statement           */
/*   stmts    = %str(other mixed statements, incl. random) */
/*   data     = dataset with original data                 */
/*   response = name of response variable in dataset       */
/*                                                         */
/***********************************************************/

data t;
set &data;
phi=&phimin;
size=(&phimax-&phimin)/&steps;
do i=0 to &steps;
  if phi ne 0 then _y_=(&response**phi-1)/phi;
  if abs(phi)<1e-10 then _y_=log(&response);
  jac=(phi-1)*log(&response); 
  output;
  phi=phi+size;
end;

proc sort data=t out=t;
by i phi;
run;

ods output fitstatistics=fit;
proc mixed data=t method=ml;
class &class;
model _y_=&model;
&stmts;
by i phi;
run;

data fit2;
set fit;
if descr='-2 Log Likelihood';
loglik=-value/2;

proc means data=t noprint;
var jac;
output out=jac sum=;
by i phi;
run;

data loglik;
merge jac fit2;
loglik=loglik+jac;

proc print data=loglik;
run;

symbol value=dot i=join;
proc gplot data=loglik;
plot loglik*phi;
run;

proc means data=loglik noprint;
var loglik;
id phi;
output out=max max=max;
run;

data r;
if _n_=1 then set max;
set loglik;
if loglik=max;
deviance=-2*loglik;

title 'Optimal value of phi';
proc print data=r;
var phi loglik deviance;
run;

title ' ';

%mend;


/*Example: Number of weed plants per plot in a completely randomized 
  experiment comparing three herbicides (trt=1, 2, 3) and a control (trt=4).
*/

data s;
do rep=1 to 6;
  do trt=1 to 4;
    input y@@; output;
  end;
end;
datalines;
	4	8	25	33
	5	11	28	21
	2	9	20	48
	5	12	15	18
	4	7	14	53
	1	7	30	31
;

%boxcox(phimin=0,phimax=1, steps=10, model=trt, class=trt, data=s, response=y);

/*
  Turnip data. Taken from:
  Lee, C.J., M. O'Donnell, and M. O'Neill. 2008. 
  Statistical analysis of field trials with changing treatment variance. 
  Agronomy Journal 100:484-489.
*/
data w;
do variety=1 to 2;
do date=1 to 2;
do density=1 to 4;
do block=1 to 4;
input y@@; output;
end;
end;
end;
end;
datalines;
2.7 1.4 1.2 3.8
7.3 3.8 3.0 1.2
6.5 4.6 4.7 0.8
8.2 4.0 6.0 2.5
4.4 0.4 6.5 3.1
2.6 7.1 7.0 3.2
24.0 14.9 14.6 2.6
12.2 18.9 15.6 9.9
1.2 1.3 1.5 1.0
2.2 2.0 2.1 2.5
2.2 6.2 5.7 0.6
4.0 2.8 10.8 3.1
2.5 1.6 1.3 0.3
5.5 1.2 2.0 0.9
4.7 13.2 9.0 2.9
14.9 13.3 9.3 3.6
;

%boxcox(phimin=0.0,phimax=1.0, steps=10, model=variety|date|density, stmts = %str(random block;), class=variety date density block, data=w, response=y);
/*reduce range of grid from 0.1 to 0.3*/
%boxcox(phimin=0.1,phimax=0.3, steps=20, model=variety|date|density, stmts = %str(random block;), class=variety date density block, data=w, response=y);
/*reduce range of grid from 0.22 to 0.24*/
%boxcox(phimin=0.22,phimax=0.24, steps=20, model=variety|date|density, stmts = %str(random block;), class=variety date density block, data=w, response=y);