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

/***********************************************************/
/* This macro performs the folded exponential transfor-    */
/* mation by                                               */
/*                                                         */
/* Piepho, H.P. (2003): The folded exponential             */
/* transformation for proportions. The Statistician        */
/*                                                         */
/*   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       */
/*   opt      =%str(options to mixed call                  */
/*                                                         */
/*  WARNING: When model contains continous covariates, all */
/*  records with missing values in at least one of the     */
/*  covarates must be deleted prior to submission to       */
/*  this macro                                             */
/*                                                         */
/***********************************************************/

data t;
phi=&phimin;
size=(&phimax-&phimin)/&steps;
do ii=0 to &steps;
  do jj=1 to filesize;
    set &data point=jj nobs=filesize;
    if abs(phi)> 1e-10 then _y_=(exp(phi*&response)-1)/phi/2 - (exp(phi*(1-&response))-1)/phi/2;
    if abs(phi)<1e-10 then _y_=&response;
    jac=( exp(phi*&response) + exp(phi*(1-&response)) )/2; 
    jac=log(jac); 
    output;
  end;
  phi=phi+size;
end;
stop;  
/*
proc sort data=t out=t;
by phi;
run;
*/

*proc print data=t; run;

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

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

proc means data=t;
var jac;
output out=jac sum=;
by 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;
var loglik;
id phi;
output out=max max=max;
run;

data r;
if _n_=1 then set max;
set loglik;
if loglik=max;

title 'optimal value for phi';
proc print data=r;
var phi loglik;
run;
title ' ';

%mend;



data d;
input
vers pos dg sorte dorma leba;
x=dorma/leba;
datalines;
1	.	1	1	21	93
1	.	1	1	25	95
1	.	1	1	16	94
1	.	1	1	35	98
2	5	4	1	8	99
2	1	4	1	4	96
2	2	4	1	13	96
2	4	4	1	15	94
3	1	5	1	10	99
3	2	5	1	7	98
3	4	5	1	9	95
3	5	5	1	8	96
1	.	1	2	64	88
1	.	1	2	68	93
1	.	1	2	51	86
1	.	1	2	55	85
2	1	4	2	56	94
2	2	4	2	45	99
2	4	4	2	53	96
2	5	4	2	45	95
3	1	5	2	56	95
3	2	5	2	46	87
3	4	5	2	56	94
3	5	5	2	43	89
1	.	1	3	3	98
1	.	1	3	0	99
1	.	1	3	0	96
1	.	1	3	2	100
2	3	4	3	1	99
2	1	4	3	1	100
2	6	4	3	0	94
2	4	4	3	1	96
3	1	5	3	0	100
3	2	5	3	0	98
3	5	5	3	0	98
3	6	5	3	0	100
1	.	1	4	12	96
1	.	1	4	9	98
1	.	1	4	11	99
1	.	1	4	8	98
2	3	4	4	4	99
2	3	4	4	2	97
2	2	4	4	3	99
2	6	4	4	5	94
3	2	5	4	2	100
3	3	5	4	8	100
3	5	5	4	1	99
3	6	5	4	1	97
1	.	1	5	5	95
1	.	1	5	7	99
1	.	1	5	14	98
1	.	1	5	21	98
2	2	4	5	1	100
2	6	4	5	3	97
2	5	4	5	5	97
2	3	4	5	0	96
3	1	5	5	7	98
3	3	5	5	2	93
3	4	5	5	1	96
3	6	5	5	4	97
1	.	1	6	10	97
1	.	1	6	10	98
1	.	1	6	5	93
1	.	1	6	9	96
2	6	4	6	5	97
2	4	4	6	4	95
2	5	4	6	3	96
2	1	4	6	2	94
3	3	5	6	5	96
3	3	5	6	0	90
3	4	5	6	2	95
3	6	5	6	2	93
1	1	2	7	52	95
1	2	2	7	36	93
1	3	2	7	50	92
1	4	2	7	40	93
2	1	3	7	15	94
2	2	3	7	27	90
2	3	3	7	35	95
2	4	3	7	21	91
3	1	6	7	20	92
3	2	6	7	13	90
3	3	6	7	20	92
3	4	6	7	20	94
1	1	2	8	36	86
1	2	2	8	44	94
1	3	2	8	46	90
1	4	2	8	37	87
2	1	3	8	30	90
2	2	3	8	26	91
2	3	3	8	37	96
2	4	3	8	37	94
3	1	6	8	23	88
3	2	6	8	23	81
3	3	6	8	15	86
3	4	6	8	26	93
1	1	2	9	35	92
1	2	2	9	36	96
1	3	2	9	39	97
1	4	2	9	44	99
2	1	3	9	7	99
2	2	3	9	14	96
2	3	3	9	12	99
2	4	3	9	7	97
3	1	6	9	14	96
3	2	6	9	11	100
3	3	6	9	16	97
3	4	6	9	15	98
1	1	2	10	63	91
1	2	2	10	67	89
1	3	2	10	65	88
1	4	2	10	71	88
2	1	3	10	17	92
2	2	3	10	.	44
2	3	3	10	28	94
2	4	3	10	39	97
3	1	6	10	15	95
3	2	6	10	22	90
3	3	6	10	22	95
3	4	6	10	33	93
;

%foldedexp(phimin=4.7,phimax=4.9, steps=100, class=vers dg sorte, model=sorte, stmts=%str(random vers vers*dg vers*sorte;), data=d, response=x);