/* Find best Box-Cox lambda parameter using algorithm 1 or 2. 
   Algorithm 1 (alg=1) is the optimized version of the original code which
     does a grid search. 
   Algorithm 2 first obtains an initial estimate of lambda using proc transreg, and
     then finds the optimal lambda by assuming the curve of the SSE values with 
	 respect to lambda is unimodal or monotonic from the initial estimate to
	 the optimal vlaue. 
	 
   The original %find_best_lambda function has been renamed to %find_best_lambda_orig.
*/
%macro find_best_lambda(input_data=,
												row_subset=,
												variable=,
												covariates=,
												weight=,
												lambda_start=,
												lambda_increment=,
												num_lambdas=,
												outlib=,
												outbc=,
												alg=2,
												DEBUG=0);
												
	%local wgtstr retdata;
	
	%let retdata = &outlib..&outbc.; 			 /* Output data set */
	%let wgtstr  = ;                       /* For a weight statement */
	%if &weight. ^= %str() %then %let wgtstr = weight &weight.;
	%let covariates = %left(%trim(&covariates));
												
	**Subset to non-zero values and apply other user-specified subsetting;
	data _subset;
		set &input_data.;
		
		%if &row_subset. ^= %str() %then %do;
			if (&row_subset.) and &variable. > 0 then _var_value = &variable.;
			else _var_value = .;
		%end;
		%else %do;
			if &variable. > 0 then _var_value = &variable.;
			else _var_value = .;
		%end;
		keep _var_value &weight &covariates;
	run;
	
	%if %eval(&alg = 1) %then %do;
		/* Original algorithm (grid search) optimized code */
		%find_best_lambda_alg1(&retdata, &variable, &covariates, &wgtstr,
							   &lambda_start, &lambda_increment, &num_lambdas);
	%end;
	%else %do;
		/* New algorithm */
		%find_best_lambda_alg2(&retdata, &variable, &covariates, &wgtstr,
						   	   &lambda_start, &lambda_increment, &num_lambdas, DEBUG=&DEBUG);
	%end;

%mend find_best_lambda;

/* Find best Box-Cox lambda parameter using algorithm 2 */
%macro find_best_lambda_alg2(ret, variable, covars, wgtstr,
							lambda_start, lamda_increment, num_lambdas, DEBUG=0);
				
	/* The input data set will be called _subset containing variable _var_value */
				
	%local lambda0 lambda SSE0 SSE iter tmp lambda_stop probitVals probitSumSq tmpdat;
	%local contflag missval var plusMinus covflag left right step tstep maxiter;
												
	%let tmpdat  = _BC_alg2_tmp;      /* Temporary dat set */
	%let missval = 9.99e99;           /* Missing value for SSE */
	%let var     = _var_value;        /* Variable on data _subset */ 
	%let covflag = %length(&covars);
	%let tstep   = 0.1;               /* Step size for proc transreg */
	%let maxiter = &num_lambdas;      /* A maximum number of iterations to prevent an infinite loop */
	
    *Largest lambda value;
	%let lambda_stop = %sysevalf(&lambda_start + (&num_lambdas - 1) * &lambda_increment);
	
	* Initial estimate for lambda in case proc transreg fails;
	%let lambda0     = %sysevalf((&lambda_start + &lambda_stop)/2);     
	
	/* Get probit values that are used later */
	%let tmp         = %BC_getProbitVals();
	%let probitSumSq = %scan(%bquote(&tmp), 1, |); /* sum of squares  */
	%let probitVals  = %scan(%bquote(&tmp), 2, |); /* for array N later */
	
	%deleteData(&tmpdat);
	
	* Get initial estimate for lambda using proc transreg;
	%if &covflag %then %do;
		proc transreg data=_subset(where=(&var > 0)) 
			outtest=&tmpdat(keep=statistic value where=(statistic="Lambda"));
			model BoxCox(&var / lambda = &lambda_start to &lambda_stop by &tstep) = identity(&covars) / noprint; 
		run;
	%end;
	%else %do;
		*Add intercept and use the nozeroconstant option in proc transreg ;
		data _subset; set _subset;
			_intercept_ = 1;
		run;
		proc transreg data=_subset(where=(&var > 0)) 
			nozeroconstant
			outtest=&tmpdat(keep=statistic value where=(statistic="Lambda"));
			model BoxCox(&var / lambda = &lambda_start to &lambda_stop by &tstep) = identity(_intercept_) / noprint; 
		run;
	%end;
	
	* Store initial estimate for lambda as lambda0;
	%if %sysfunc(exist(&tmpdat)) %then %do;
		data _NULL_; set &tmpdat;
			call symputx("lambda0", value);
		run;
	%end;
	%if &DEBUG %then %put lambda0 = &lambda0;	
	
	*Compute SSE for initial lambda, store it in output data &tmpdat;
	%BC_mainLambdaCode(&tmpdat, &lambda0, &wgtstr, &covars, %bquote(&probitVals), &probitSumSq);
	
	* Store initial SSE in SSE0;
	data _NULL_; set &tmpdat;
		call symputx("SSE0", SSE);
	run;
	
	* Get the direction of decreasing SSE, try to the right first;
	%let lambda = %sysevalf(&lambda0 + &lambda_increment);
	%if %sysevalf(&lambda > &lambda_stop) %then %let lambda = &lambda_stop;
	%BC_mainLambdaCode(&tmpdat, &lambda, &wgtstr, &covars, %bquote(&probitVals), &probitSumSq);
	data _NULL_; set &tmpdat;
		call symputx("tmp", SSE);
	run;
	%if %sysevalf(&tmp < &SSE0) %then %do;
		* Go to the right, update lambda0 and SSE0 since the new lambda is better;
		%let lambda0   = &lambda;
		%let SSE0      = &tmp;
		%let plusMinus = +;
	%end;
	%else %do;
		*Go to the left, no need to update lambda0;
		%let plusMinus = -;
	%end;
	
	* Initialize variables for the while loop below;
	%let iter     = 0;
	%let contflag = 1;
	%let left     = 0;
	%let right    = 0;
	%let step     = %sysevalf(&tstep/2);  /* Initial step size in algorithm */
	
	* Find optimal lambda ;
	%do %while (&contflag);
		%let iter = %eval(&iter + 1);
		%let SSE  = &missval;
		%if "&plusMinus" = "+" %then %let right = %eval(&right + 1);
		%else %let left = %eval(&left + 1);
		
		* Try a new lambda;
		%let lambda = %sysevalf(&lambda0 &plusMinus &step); 
		
		* Check that lambda is inside the range of lambdas;
		%if %sysevalf(&lambda >= &lambda_start) or %sysevalf(&lambda <= &lambda_stop) %then %do;
			* Compute SSE for lambda;
			%BC_mainLambdaCode(&tmpdat, &lambda, &wgtstr, &covars, %bquote(&probitVals), &probitSumSq);
			data _NULL_; set &tmpdat;
				call symputx("SSE", SSE);
			run;
		%end;

		*See if SSE is smaller and check for missing value;
		%if %sysevalf(&SSE < &SSE0) and %sysevalf(&SSE < &missval) %then %do;
			*Smaller SSE, so update and continue;
			%let lambda0 = &lambda;
			%let SSE0    = &SSE;
			%if &DEBUG %then %put iter = &iter lambda = &lambda0;
		%end;
		%else %do;
			*Larger SSE, reduce step size and/or change direction;
			%let step = %sysevalf(&step/2);
			
			*If we have tried both directions, then reset;
			%if &right and &left %then %do;
				*Do not change direction;
			%end;
			%else %if &right %then %do;
				*Go to the left;
				%let plusMinus = -;
			%end;
			%else %do;
				*Go to the right;
				%let plusMinus = +;
			%end;
			*Reset since the step size has changed;
			%let right = 0;
			%let left  = 0;
		%end;
		
		*If the best lambda (lambda0) ios outside the grid values, then stop;
		%if %sysevalf(&lambda0 < &lambda_start) or %sysevalf(&lambda0 > &lambda_stop) %then %let contflag = 0;
		
		*Check the step size and stop if it is smaller then lambda_increment;
		%if %sysevalf(&step < &lambda_increment) %then %let contflag = 0;

		*Check max number of iterations to prevent infinite loop;
		%if %eval(&iter >= &maxiter) %then %let contflag = 0;
	%end;
	
	* Get lambda closest to lambda0 from above loop, and call it lambda0;
	%getClosestLambda(&tmpdat, &lambda0, &lambda_start, &lambda_increment, &num_lambdas);
	data _NULL_; set &tmpdat;
		call symputx("lambda0", lambda);
	run;
	%if &DEBUG %then %put closest lambda = &lambda0;
	
	*Final calculation for lambda0 and lambda0 +- increment;
	%BC_mainLambdaCode(&tmpdat, &lambda0, &wgtstr, &covars, %bquote(&probitVals), &probitSumSq);
	data _NULL_; set &tmpdat;
		call symputx("SSE0", SSE);
	run;
	%let lambda = %sysevalf(&lambda0 + &lambda_increment); 
	%if %sysevalf(&lambda > &lambda_stop) %then %let lambda = &lambda_stop;
	%BC_mainLambdaCode(&tmpdat, &lambda, &wgtstr, &covars, %bquote(&probitVals), &probitSumSq);
	data _NULL_; set &tmpdat;
		call symputx("SSE", SSE);
	run;
	%if %sysevalf(&SSE < &SSE0) %then %do;
		*Final lambda;
		%let lambda0 = &lambda;
		%let SSE0    = &SSE;
	%end;
	%else %do;
		*Try to the left;
		%let lambda = %sysevalf(&lambda0 - &lambda_increment); 
		%if %sysevalf(&lambda < &lambda_start) %then %let lambda = &lambda_start;
		%BC_mainLambdaCode(&tmpdat, &lambda, &wgtstr, &covars, %bquote(&probitVals), &probitSumSq);
		data _NULL_; set &tmpdat;
			call symputx("SSE", SSE);
		run;
		%if %sysevalf(&SSE < &SSE0) %then %do;
			%let lambda0 = &lambda;
			%let SSE0    = &SSE;
		%end;
	%end;
	
	%deleteData(&tmpdat);
	%if &DEBUG %then %put Number of iterations = &iter  Final lambda = &lambda0;	

	data &ret.;
		length variable $32.;
		variable    = "&variable.";
		tran_lambda = &lambda0;
	run;
	
%mend find_best_lambda_alg2;

%macro getClosestLambda(ret, value, lambda_start, lambda_increment, num_lambdas);

	data &ret;
		keep lambda;
		retain min_lambda min_dist;
		value            = &value;
		lambda_start     = &lambda_start;
		lambda_increment = &lambda_increment;
		num_lambdas      = &num_lambdas;
		min_lambda       = lambda_start;
		lambda           = min_lambda;
		min_dist         = abs(min_lambda - value);
		do i = 1 to (num_lambdas - 1);
			lambda = lambda_start + i*lambda_increment;
			dist   = abs(lambda - value);
			if (dist <= min_dist) then do;
				min_dist   = dist;
				min_lambda = lambda;
			end;
			else do;
				lambda = min_lambda;				
				leave;
			end;				
		end;
		output;
		stop;
	run;
		
%mend getClosestLambda;

/* Get the probit sum of squares and the values for the array N. 
   This will return the objects probitVals and probitSumSq for macro 
   function BC_mainLambdaCode. */
%macro BC_getProbitVals(a=1, b=99, step=0.01);

  %local ret i tmp sumsq flag;
  
  %let ret   = ;
  %let sumsq = 0;
  %let flag  = 0;
  %do i = &a %to &b;
	%let tmp = %sysfunc(probit(&step * &i));
	%if &flag %then %let ret = &ret. , &tmp;
	%else %do;
		%let ret  = &tmp;
		%let flag = 1;
	%end;
	%let sumsq = %sysevalf(&sumsq + &tmp * &tmp);
  %end;
  %let ret = &sumsq | &ret;
  
  &ret
  
%mend BC_getProbitVals;

/* Main code for computing the SSE for a given lambda */
%macro BC_mainLambdaCode(out, lambda, wgtstr, covariates, probitVals, probitSumSq, missval=9.99e99);

	%local nocovflag;
	%if %length(&covariates) %then %let nocovflag = 0;
	%else %let nocovflag = 1;

	%if %sysevalf(&lambda = 0) %then %do;
		data _subset; set _subset;
			_lambda       = 0;
			_boxcox_value = log(_var_value);
		run;
	%end;
	%else %do;
		data _subset; set _subset;
			_lambda       = &lambda;
			_boxcox_value = (_var_value**_lambda - 1)/_lambda;
		run;
	%end;
		
	%if &nocovflag %then %do;
		proc univariate data=_subset noprint;
			var _boxcox_value;
			&wgtstr;
			output out=_percentiles pctlpts=1 to 99 by 1 pctlpre=P;
		run;
	%end;
	%else %do;
		ods select none;
		proc surveyreg data=_subset;
			model _boxcox_value = &covariates.;
			&wgtstr;
			output out=_reg_out r=residual;
		run;
		ods select all;
		proc univariate data=_reg_out noprint;
			var residual;
			&wgtstr;
			output out=_percentiles pctlpts=1 to 99 by 1 pctlpre=P;
		run;	
	%end;
	
	data &out (keep = SSE);
		set _percentiles;
		
		array N{99} (&probitVals); **Normal percentiles;
		array P{99} P1-P99;        **Box-Cox percentiles;
		
		P_mean        = mean(of P1-P99);
		bc_dev_sq     = 0;
		norm_bc_cross = 0;
		do i = 1 to 99;
			tmp           = P{i} - P_mean;
			bc_dev_sq     = bc_dev_sq + tmp*tmp;
			norm_bc_cross = norm_bc_cross + N{i}*tmp;
		end;
		
		SSE = &probitSumSq - (norm_bc_cross*norm_bc_cross)/bc_dev_sq;
		if (missing(SSE)) then SSE = &missval;
	run;

%mend BC_mainLambdaCode;		
		
/* Find best Box-Cox lambda parameter by grid search. This function is the optimized
   version of the original %find_best_lambda. The original %find_best_lambda function 
   is below and has been renamed to %find_best_lambda_orig.
   The data set used in _subset. */
%macro find_best_lambda_alg1(retdata, variable, covariates, wgtstr,
							 lambda_start, lambda_increment, num_lambdas);
												
	%local lambda i nocovflag;
												
	%if %length(&covariates) %then %let nocovflag = 0;
	%else %let nocovflag = 1;
	
	**Calculate 1st-99th percentiles of Box-Cox transformations for the lambda grid;
	%do i = 0 %to %eval(&num_lambdas. - 1);
		%let lambda = %sysevalf(&lambda_start. + &i.*&lambda_increment.);

		%if %sysevalf(&lambda = 0) %then %do;
		  data _subset; set _subset;
			_lambda       = 0;
			_boxcox_value = log(_var_value);
		  run;
		%end;
		%else %do;
		  data _subset; set _subset;
			_lambda       = &lambda;
			_boxcox_value = (_var_value**_lambda - 1)/_lambda;
		  run;
		%end;
		
		%if &nocovflag %then %do;
			proc univariate data=_subset noprint;
				var _boxcox_value;
				&wgtstr;
				output out=_percentiles pctlpts=1 to 99 by 1 pctlpre=P;
			run;
		%end;
		%else %do;
			ods select none;
			proc surveyreg data=_subset;
				model _boxcox_value = &covariates.;
				&wgtstr;
				output out=_reg_out r=residual;
			run;
			ods select all;
			proc univariate data=_reg_out noprint;
				var residual;
				&wgtstr;
				output out=_percentiles pctlpts=1 to 99 by 1 pctlpre=P;
			run;	
		%end;
		/* Add lambda */
		data _percentiles;
			set _percentiles;
			lambda = &lambda.;
		run;
		
		%if &i. = 0 %then %do;
			data _lambda_percentiles; set _percentiles; run;
		%end;
		%else %do;
			proc append base=_lambda_percentiles
                data=_percentiles
                force;
			run;
		%end;
	%end;
	
	**Calculate SSE for each lambda, save for min SSE;
	data &retdata (keep = variable tran_lambda);
		set _lambda_percentiles end=last;
		
		retain N1-N99 SSE_min lambda_min norm_dev_sq;
		
		array N{99} N1-N99; **Normal percentiles;
		array P{99} P1-P99; **Box-Cox percentiles;
		
		if _N_ = 1 then do;
			SSE_min     = .;
			lambda_min  = .;
			norm_dev_sq = 0;
			do i = 1 to 99;
				tmp         = probit(0.01*i);
				N{i}        = tmp;
				norm_dev_sq = norm_dev_sq + tmp*tmp;
			end;
		end;
		
		P_mean        = mean(of P1-P99);
		bc_dev_sq     = 0;
		norm_bc_cross = 0;
		do i = 1 to 99;
			tmp           = P{i} - P_mean;
			bc_dev_sq     = bc_dev_sq + tmp*tmp;
			norm_bc_cross = norm_bc_cross + N{i}*tmp;
		end;
	
		SSE = norm_dev_sq - (norm_bc_cross*norm_bc_cross)/bc_dev_sq;
		if missing(SSE_min) or SSE < SSE_min then do;
			SSE_min    = SSE;
			lambda_min = lambda;
		end;

		if last then do;
			length variable $32.;
			variable    = "&variable.";
			tran_lambda = lambda_min;
			output;
		end;
	run;
	
%mend find_best_lambda_alg1;

*Original code for %find_best_lambda by grid search;
%macro find_best_lambda_orig(input_data=,
												row_subset=,
												variable=,
												covariates=,
												weight=,
												lambda_start=,
												lambda_increment=,
												num_lambdas=,
												outlib=,
												outbc=);
												
	%local lambda i;
												
	**Subset to non-zero values and apply other user-specified subsetting;
	data _subset;
		set &input_data.;
		
		%if &row_subset. ^= %str() %then %do;
		
			if (&row_subset.) and &variable. > 0 then _var_value = &variable.;
			else _var_value = .;
		%end;
		%else %do;
		
			if &variable. > 0 then _var_value = &variable.;
			else _var_value = .;
		%end;
	run;
	
	**Calculate 1st-99th percentiles of Box-Cox transformations for the lambda grid;
	%do i = 0 %to %eval(&num_lambdas. - 1);
	
		data _NULL_;
		
			lambda = &lambda_start. + &i.*&lambda_increment.;
			call symputx("lambda", lambda);
		run;
	
		data _subset;
			set _subset;
			
			_lambda = &lambda.;
			if _lambda = 0 then _boxcox_value = log(_var_value);
			else _boxcox_value = (_var_value**_lambda - 1)/_lambda;
		run;
		
		%if &covariates. = %str() %then %do;
		
			proc univariate data=_subset noprint;
			
				var _boxcox_value;
				%if &weight. ^= %str() %then %do;
					weight &weight.;
				%end;
				
				output out=_percentiles pctlpts=1 to 99 by 1 pctlpre=P;
			run;
		%end;
		%else %do;
		
			ods select none;
			
			proc surveyreg data=_subset;
			
				model _boxcox_value = &covariates.;
				%if &weight. ^= %str() %then %do;
					weight &weight.;
				%end;
				
				output out=_reg_out r=residual;
			run;
			
			ods select all;
			
			proc univariate data=_reg_out noprint;
			
				var residual;
				%if &weight. ^= %str() %then %do;
					weight &weight.;
				%end;
				
				output out=_percentiles pctlpts=1 to 99 by 1 pctlpre=P;
			run;
		%end;
		
		data _percentiles;
			set _percentiles;
			
			lambda = &lambda.;
		run;
		
		data _lambda_percentiles;
			set %if &i. ^= 0 %then _lambda_percentiles;
					_percentiles;
		run;
	%end;
	
	**Calculate SSE for each lambda;
	data _lambda_percentiles;
		set _lambda_percentiles;
		
		retain N1-N99;
		
		array N{99} N1-N99; **Normal percentiles;
		array P{99} P1-P99; **Box-Cox percentiles;
		
		if _N_ = 1 then do;
			do i = 1 to 99;
			
				N{i} = probit(0.01*i);
			end;
		end;
		
		P_mean = mean(of P1-P99);
		
		norm_dev_sq = 0;
		bc_dev_sq = 0;
		norm_bc_cross = 0;
		do i = 1 to 99;
		
			norm_dev_sq = norm_dev_sq + N{i}**2;
			bc_dev_sq = bc_dev_sq + (P{i} - P_mean)**2;
			norm_bc_cross = norm_bc_cross + N{i}*(P{i} - P_mean);
		end;
		
		SSE = norm_dev_sq - (norm_bc_cross**2)/bc_dev_sq;
	run;
	
	**Select lambda with lowest SSE;
	proc sort data=_lambda_percentiles; by SSE; run;
	data &outlib..&outbc. (keep = variable tran_lambda);
		set _lambda_percentiles;
		
		length variable $32.;
		
		if _N_ = 1 then do;
		
			variable = "&variable.";
			tran_lambda = lambda;
			output;
		end;
	run;
%mend find_best_lambda_orig;

**Suggested Winsorization report;
%macro find_suggested_winsorization(input_data=,
																		row_subset=,
																		lambda=,
																		variable=,
																		is_episodic=,
																		covariates=,
																		weight=,
																		id=,
																		repeat_obs=,
																		iqr_multiple=,
																		print_report=,
																		outlib=,
																		outwin=,
																		outlim=);
	
	**Subset of data used to find thresholds;
	data _subset;
		set &input_data.;
		
		%if &row_subset. ^= %str() %then %do;
		
			if (&row_subset.) and &variable. > 0 then _var_value = &variable.;
			else _var_value = .;
		%end;
		%else %do;
		
			if &variable. > 0 then _var_value = &variable.;
			else _var_value = .;
		%end;
	run;
	
	**Box-Cox transform variable;
	data _subset;
		set _subset;
		
		_lambda = &lambda.;
		if _lambda = 0 then _boxcox_value = log(_var_value);
		else _boxcox_value = (_var_value**_lambda - 1)/_lambda;
	run;
	
	**Find outlier thresholds on the Box-Cox transformed scale;
	%if &covariates. = %str() %then %do;
	
		proc univariate data=_subset noprint;
		
			var _boxcox_value;
			%if &weight. ^= %str() %then %do;
				weight &weight.;
			%end;
			
			output out=_percentiles pctlpts=25 75 pctlpre=P;
		run;
		
		data _thresholds;
			set _percentiles;
			
			iqr = P75 - P25;
			
			bc_lower = P25 - &iqr_multiple.*iqr;
			bc_upper = P75 + &iqr_multiple.*iqr;
		run;
	%end;
	%else %do;
	
		ods select none;
	
		proc surveyreg data=_subset;
		
			model _boxcox_value = &covariates.;
			%if &weight. ^= %str() %then %do;
				weight &weight.;
			%end;
			
			output out=_reg_out r=residual p=prediction;
		run;
		
		ods select all;
		
		proc univariate data=_reg_out noprint;
			
			var residual;
			%if &weight. ^= %str() %then %do;
				weight &weight.;
			%end;
			
			output out=_percentiles pctlpts=25 75 pctlpre=P;
		run;
		
		data _thresholds;
			set _reg_out;
			if _N_ = 1 then set _percentiles;
			
			retain resid_lower resid_upper;
			
			if _N_ = 1 then do;
			
				iqr = P75 - P25;
			
				resid_lower = P25 - &iqr_multiple.*iqr;
				resid_upper = P75 + &iqr_multiple.*iqr;
			end;
			
			bc_lower = resid_lower + prediction;
			bc_upper = resid_upper + prediction;
		run;
	%end;
	
	**Backtransform thresholds to original scale;
	data _thresholds;
		set _thresholds;
		
		lambda = &lambda.;
		if lambda = 0 then do;
		
			_lower = exp(bc_lower);
			_upper = exp(bc_upper);
		end;
		else do;
		
			if bc_lower < -1/lambda then _lower = 0;
			else _lower = (lambda*bc_lower + 1)**(1/lambda);
		
			if bc_upper < -1/lambda then _upper = 0;
			else _upper = (lambda*bc_upper + 1)**(1/lambda);
		end;
	run;
	
	**Dataset of Winsorization limits;
	%if &covariates. = %str() %then %do;
		
		**Dataset of lower and upper limits for whole dataset;
		data &outlib..&outlim. (keep = lower upper);
			set _thresholds (rename=(_lower=lower _upper=upper));
		run;
	%end;
	%else %do;
	
		**Dataset of lower and upper limits for every observation;
		data &outlib..&outlim. (keep = &id. &repeat_obs. lower upper);
			set _thresholds (keep = &id. &repeat_obs. _lower _upper
											 rename=(_lower=lower _upper=upper));
		run;
	%end;
	
	**Report of outlier values and suggested Winsorization;
	data &outlib..&outwin. (keep = &id. &repeat_obs. &variable. &variable._win);
		set &input_data.;
		%if &covariates. = %str() %then %do;
			if _N_ = 1 then set _thresholds (keep = _lower _upper);
		%end;
		%else %do;
			set _thresholds (keep = _lower _upper);
		%end;
		
		&variable._win = &variable.;
		
		%if %upcase("&is_episodic.") = "Y" %then %do;
		
			**for episodic variables, Winsorize low outliers to zero;
			if &variable. < _lower and &variable. > 0 then &variable._win = 0; 
			else if &variable. > _upper and &variable. > 0 then &variable._win = _upper;
		%end;
		%else %do;
		
			if &variable. < _lower then &variable._win = _lower;
			else if &variable. > _upper then &variable._win = _upper;
		%end;
		
		if &variable. ^= &variable._win then output;
	run;
	
	%if %upcase("&print_report.") = "Y" %then %do;
	
		proc print data=&outlib..&outwin.; 
		
			title "Suggested Winsorization for &variable.";
		run;
		title;
	%end;
%mend find_suggested_winsorization;

**Influential subject report;
%macro find_influential_subjects(input_data=,
																 row_subset=,
																 lambda=,
																 variable=,
																 weight=,
																 id=,
																 repeat_obs=,
																 alpha=,
																 multiple_test=,
																 print_report=,
																 outlib=,
																 outinf=);
																 
	%local num_recalls;
																 
	**Remove non-consumption days;
	data _subset;
		set &input_data.;
		if &variable. > 0;
	run;
	
	**Assign numbers to recalls and calculate total number of recalls;
	proc sort data=_subset; by &repeat_obs.; run;
	data _subset;
		set _subset end=last_obs;
		by &repeat_obs.;
		
		retain _recallnum;
		
		if _N_ = 1 then _recallnum = 0;
		
		if first.&repeat_obs then _recallnum = _recallnum + 1;
		
		if last_obs = 1 then call symputx("num_recalls", _recallnum);
	run;
	
	**Transpose to one record per subject;
	proc sort data=_subset; by &id.; run;
	data _subject_data (keep = &id. &weight. &variable.1-&variable.&num_recalls.);
		set _subset;
		by &id.;
		
		retain &variable.1-&variable.&num_recalls.;
		
		array obs{&num_recalls.} &variable.1-&variable.&num_recalls.;
		
		if first.&id. then do;
		
			do _i = 1 to &num_recalls.;
			
				obs{_i} = .;
			end;
		end;
		
		obs{_recallnum} = &variable.;
	run;
	
	**Box-Cox transform the variable values;
	data _subject_data;
		set _subject_data;
		
		array obs{&num_recalls.} &variable.1-&variable.&num_recalls.;
		array bc{&num_recalls.} bc1-bc&num_recalls.;
		
		lambda = &lambda.;
		
		do i = 1 to &num_recalls.;
		
			if lambda = 0 then bc{i} = log(obs{i});
			else bc{i} = (obs{i}**lambda - 1)/lambda;
		end;
	run;
	
	**Calculate the variance of observations on the Box-Cox transformed scale for each subject;
	data _subject_data;
		set _subject_data;
		
		array bc{&num_recalls.} bc1-bc&num_recalls.;
		
		_subject_variance = var(of bc{*});
		_num_subjects = n(of bc{*});
		_subject_df = _num_subjects - 1;
	run;
	
	**Subset to subjects with at least 2 observations;
	data _subject_data;
		set _subject_data;
		if _num_subjects >= 2;
	run;
	
	**Calculate weighted sum of variances and total degrees of freedom;
	data _sums (keep = _sum_variance _sum_weights _sum_df);
		set _subject_data end=last_obs;
		
		retain _sum_variance _sum_weights _sum_df _num_subjects;
		
		if _N_ = 1 then do;
		
			_sum_variance = 0;
			_sum_weights = 0;
			_sum_df = 0;
			_num_subjects = 0;
		end;
		
		_sum_variance = _sum_variance + _subject_variance*_subject_df*&weight.;
		_sum_weights = _sum_weights + _subject_df*&weight.;
		_sum_df = _sum_df + _subject_df;
		_num_subjects = _num_subjects + 1;
		
		if last_obs = 1 then output;
	run;
	
	**Merge sums into dataset;
	data _subject_data;
		set _subject_data;
		
		_merge = 1;
	run;
	
	data _sums;
		set _sums;
		
		_merge = 1;
	run;
	
	data _subject_data;
		merge _subject_data
					_sums;
		by _merge;
	run;
	
	**F-test on variance for each subject against mean variance of other subjects;
	data _subject_data;
		set _subject_data;
		
		**Sum of variance and degrees of freedom for all subjects EXCEPT the current one;
		_sum_variance = _sum_variance - _subject_variance*_subject_df*&weight.;
		_sum_weights = _sum_weights - _subject_df*&weight.;
		_sum_df = _sum_df - _subject_df;
		
		**F-test;
		_mean_variance = _sum_variance/_sum_weights;
		_f_statistic = _subject_variance/_mean_variance;
		_p_value = 1 - probf(_f_statistic, _subject_df, _sum_df);
	run;
	
	proc sort data=_subject_data; by _p_value; run;
	
	**Determine if the F-test p-value is significant based on the provided alpha;
	%if &multiple_test. = %str() %then %do;
		
		data _subject_data;
			set _subject_data;
			
			if _p_value < &alpha. then _influential = 1;
			else _influential = 0;
		run;
	%end;
	%else %if %upcase("&multiple_test.") = "BONFERRONI" %then %do;
	
		data _subject_data;
			set _subject_data;
				
			if _p_value < &alpha./_num_subjects then _influential = 1;
			else _influential = 0;
		run;
	%end;
	%else %if %upcase("&multiple_test.") = "BENJAMINI-HOCHBERG" %then %do;
	
		%local last_significant;
	
		data _NULL_;
			set _subject_data end=last_obs;
			
			retain _last_significant;
			
			if _N_ = 1 then _last_significant = 0;
			
			if _p_value < (&alpha./_num_subjects)*_N_ then _last_significant = _N_;
			
			if last_obs = 1 then call symputx("last_significant", _last_significant);
		run;
		
		data _subject_data;
			set _subject_data;
			
			if _N_ <= &last_significant. then _influential = 1;
			else _influential = 0;
		run;
	%end;
	%else %do;
	
		%put Specified type of multiple testing correction is not available.;
	%end;
	
	**Create report of subjects found to be influential;
	data &outlib..&outinf. (keep = &id. _p_value &variable.1-&variable.&num_recalls.);
		set _subject_data;
		if _influential = 1;
	run;
	
	%if %upcase("&print_report.") = "Y" %then %do;
	
		proc print data=&outlib..&outinf.;
		
			title "Subjects with Influential Variance of &variable.";
		run;
		title;
	%end;
%mend find_influential_subjects;

**Generates Q-Q plots against standard normal distribution before and after Box-Cox transformation;
%macro make_qq_plots(input_data=,
										 row_subset=,
										 variable=,
										 lambda=,
										 weight=);
										 
	**Subset data;
	data _input (keep = _orig_value _boxcox_value &weight.);
		set &input_data.;
		
		%if &row_subset. ^= %str() %then %do;
		
			if &row_subset.;
		%end;
		
		**Remove non-consumption days;
		if &variable. > 0 then _orig_value = &variable.;
		else _orig_value = .;
		
		**Box-Cox transformation;
		%if &lambda. = 0 %then %do;
		
			_boxcox_value = log(_orig_value);
		%end;
		%else %do;
		
			_boxcox_value = (_orig_value**&lambda. - 1)/&lambda.;
		%end;
	run;
	
	**Calculate quantiles of original scale and Box-Cox transformed variable;
	proc univariate data=_input vardef=wdf noprint;
	
		var _orig_value _boxcox_value;
		
		%if &weight. ^= %str() %then %do;
		
			weight &weight.;
		%end;
		
		output out=_quantiles pctlpre=orig_P boxcox_P
													pctlpts=1 to 99 by 1
													
													mean=orig_mean boxcox_mean
													std=orig_std boxcox_std;
	run;
	
	**Transpose quantiles for plotting;
	data _qq_plot_data (keep = std_norm_quantile orig_quantile_actual orig_quantile_normal boxcox_quantile_actual boxcox_quantile_normal);
		set _quantiles;
		
		array orig_quantiles{99} orig_P1-orig_P99;
		array boxcox_quantiles{99} boxcox_P1-boxcox_P99;
		
		do i = 1 to 99;
		
			q = 0.01*i;
			
			**Standard normal quantiles;
			std_norm_quantile = quantile("Normal", q, 0, 1);
			
			**Original scale quantiles;
			orig_quantile_actual = orig_quantiles{i};
			orig_quantile_normal = quantile("Normal", q, orig_mean, orig_std);
			
			**Box-Cox scale quantiles;
			boxcox_quantile_actual = boxcox_quantiles{i};
			boxcox_quantile_normal = quantile("Normal", q, boxcox_mean, boxcox_std);
			
			output;
		end;
	run;
	
	**Q-Q plot of original scale variable vs. standard normal;
	proc sgplot data=_qq_plot_data;
	
		scatter x=orig_quantile_actual y=std_norm_quantile / name="Actual" legendlabel="Actual" markerattrs=(symbol=circle color=black);
		series x=orig_quantile_normal y=std_norm_quantile / name="Normal" legendlabel="Normal" lineattrs=(pattern=solid color=black thickness=2);
		
		xaxis label="&variable.";
		yaxis label="Standard Normal";
		
		keylegend "Actual" "Normal" / across=1 position=bottomright;
		
		title "Quantiles of &variable. vs. Standard Normal";
		title2 "Original Scale";
	run;
	
	**Q-Q plot of Box-Cox scale variable vs. standard normal;
	proc sgplot data=_qq_plot_data;
	
		scatter x=boxcox_quantile_actual y=std_norm_quantile / name="Actual" legendlabel="Actual" markerattrs=(symbol=circle color=black);
		series x=boxcox_quantile_normal y=std_norm_quantile / name="Normal" legendlabel="Normal" lineattrs=(pattern=solid color=black thickness=2);
		
		xaxis label="&variable.";
		yaxis label="Standard Normal";
		
		keylegend "Actual" "Normal" / across=1 position=bottomright;
		
		title "Quantiles of &variable. vs. Standard Normal";
		title2 "Box-Cox Scale, Lambda = &lambda.";
	run;
%mend make_qq_plots;