Macro for Multiple Comparisons with the Best Mean
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: MCB */
/* TITLE: Macro for Multiple Comparisons with the Best Mean */
/* PRODUCT: STAT */
/* SYSTEM: ALL */
/* KEYS: multiple comparisons, ranking and selection, */
/* KEYS: bioequivalence */
/* PROCS: MIXED, SORT, TRANSPOSE, PRINT, DATASETS */
/* DATA: */
/* */
/* SUPPORT: Randy Tobias UPDATE: 16Jan96 */
/* REF: Hsu, 1995 */
/* MISC: */
/****************************************************************/
/*-------------------------------------------------------------------
DISCLAIMER:
THIS INFORMATION IS PROVIDED BY SAS INSTITUTE INC. AS A SERVICE TO
ITS USERS. IT IS PROVIDED "AS IS". THERE ARE NO WARRANTIES,
EXPRESSED OR IMPLIED, AS TO MERCHANTABILITY OR FITNESS FOR A
PARTICULAR PURPOSE REGARDING THE ACCURACY OF THE MATERIALS OR CODE
CONTAINED HEREIN.
-------------------------------------------------------------------*/
/********************************************************************
Abstract:
Given a SAS data set with a grouping variable and a numeric
response variable, the macros defined here compute the within-group
means of the response, confidence intervals, and p-values for
multiple comparisons with the best mean (MCB). Confidence
intervals may be either constrained to contain 0 (%MCB and %MCW) or
unconstrained (%UMCB and %UMCW); comparisons may be with either the
maximum mean (%MCB and %UMCB) or the minimum mean (%MCW and %UMCW).
Introduction:
You are conducting an experiment on the effects of several
alternative drugs for treating a certain disease. The goal is to
determine which drugs are most effective. This is the task of
multiple comparisons---finding which groups are superior to which
other groups, and doing so in a manner which controls the over-all
probability of an incorrect inference. However, in this case not
all pairwise differences are of interest: you only want to compare
each drug with the true best drug.
This situation is called multiple comparisons with the best or MCB
(Hsu, 1996). It is related to several other multiple inference
techniques, such as bioequivalence testing and ranking and
selection (op cit.). MCB may assert certain treatments to be
inferior to the true best, and other treatments to be within such a
small distance of the true best that you may consider them to be
practically equivalent to the best.
By giving up the ability to say precisely how inferior the
not-the-best treatments are, MCB provides sharper inference than
can be achieved by evaluating all pairwise comparisons. On the
other hand, if you need to know how inferior the not-the-best
treatments are, unconstrained multiple comparisons with the best
(UMCB) provides this sort of analysis. MCB is executed by multiply
performing a one-sided Dunnett's test for comparisons with a
control, in turn treating each of the alternative drugs as the
control which is potentially the best; UMCB deduces from two-sided
Dunnett's tests (or alternatively from Tukey's all-pairwise test.)
You can use the %MCB macro to perform MCB analysis, and the %UMCB
macro to perform UMCB analysis, where in both cases the "best"
population mean is defined as the maximum one; use %MCW and %UMCW
if you want to treat the minimum population mean as the "best".
These macros use the MIXED procedure and the output manager to
perform Dunnett's and Tukey's tests and write the results to SAS
data sets, which are then processed to compute the standard form of
MCB and UMCB analysis, respectively.
Syntax:
The following arguments are required by each of the macros. They
must be the first three arguments and they must be in this order.
Do not use keywords for these arguments.
- the SAS data set containing the data to be analyzed
- the response variable
- the grouping variable
The following additional arguments may be listed in any order,
separated by commas:
MODEL= a linear model for the response, specified using the
effects syntax of GLM. The default is a one-way model
in the required grouping variable.
CLASS= classification variables involved in the linear model.
The default is the required grouping variable.
ALPHA= the level of significance for comparisons among the
means. The default is 0.05.
OUT= the name of the output dataset containing the MCB
analysis. The default is MCBOUT.
OPTIONS= a string containing either of the following options
NOPRINT - suppresses printed output of results
NOCLEAN - suppresses deletion of temporary datasets
Output:
The output dataset contains one observation for each group in the
dataset. For all four macros the output data set contains the
following variables:
LEVEL - formatted value of this group
LSMEAN - sample mean response within this group
SE - standard error of the sample mean for this group
CLLO - lower confidence limit for the difference between the
population mean of this group and the best population
mean
CLHI - upper confidence limit for the difference between the
population mean of this group and the best population
mean
To facilitate ranking and selection inference, the output dataset
for the %MCB and %MCW macros contains the following additional
variables:
RVAL - the smallest alpha level at which the population mean
of this group can be rejected as the best, for all
groups but the one with the best sample mean
SVAL - the smallest alpha level at which the population mean
of this group can be selected as the best treatment,
for the group with the best sample mean
Example: Comparison of filters
Hsu (1984) reports the results of a study undertaken to compare
seven different brands of water filter. For each brand, samples of
water were run through three filters and then the filters were
incubated; the response is the number of bacterial colonies grown
on filter. Some of the data are missing. The following data step
creates the FILTER dataset:
data filter;
do brand = 1 to 7;
do i = 1 to 3;
input ncolony @@;
output;
end;
end;
datalines;
69 122 95
118 154 102
171 132 182
122 119 .
204 225 190
140 130 127
170 165 .
;
A better filter is one that captures more bacteria and thus has a
higher colony count. Thus, the %MCB macro is appropriate:
%inc 'mcb.sas';
%mcb(filter,ncolony,brand);
This yields the following results:
+------------------------------------------------------------------+
| EFFECT BRAND LSMEAN SE CLLO CLHI RVAL SVAL |
| |
| BRAND 1 95.333 11.707 -153.94 0.0000 0.0001 . |
| BRAND 2 124.667 11.707 -124.61 0.0000 0.0009 . |
| BRAND 3 161.667 11.707 -87.608 0.0000 0.0418 . |
| BRAND 4 120.500 14.339 -133.84 0.0000 0.0013 . |
| BRAND 5 206.333 11.707 -7.950 86.8428 . 0.1006 |
| BRAND 6 132.333 11.707 -116.94 0.0000 0.0019 . |
| BRAND 7 167.500 14.339 -86.843 7.9498 0.1006 . |
+------------------------------------------------------------------+
The filter brand with the highest colony count was number 5, but
since the lower endpoint of the 95% confidence interval for the
difference between it and the best is negative, we cannot assert
that this particular brand is the best. However, we can say that
either brand 5 or 7 is the best, since these are the only two
brands for which the confidence interval properly contains 0.
References:
Hsu, Jason C. (1984) "Ranking and Selection and Multiple
Comparisons with the Best." In _Design_of_Experiments:_Ranking_
_and_Selection_, eds. Thomas J. Santner and Ajit C. Tamhane.
Marcel Dekker, NY.
Hsu, Jason C. (1996). _Multiple_Comparisons:_Theory_and_methods_,
Chapman & Hall, NY.
*********************************************************************/
%macro n2cclass(dsin,dsout,class,options = );
%let clean = 1;
%let iopt = 1;
%do %while(%length(%scan(&options,&iopt)));
%if (%upcase(%scan(&options,&iopt)) = NOCLEAN) %then
%let clean = 0;
%else
%put Warning: Unrecognized option %scan(&options,&iopt).;
%let iopt = %eval(&iopt + 1);
%end;
ods listing close;
ods output Variables=_Var;
proc contents data=&dsin;
run;
ods listing;
%let numc =;
%let nnumc = 0;
data _null_; set _Var;
_nc = 1; _cvar = trim(left(upcase(scan("&class",_nc))));
do while (_cvar ^= ' ');
if ( (trim(left(upcase(Variable)))=_cvar )
& (trim(left(Type ))='Num')) then do;
call symput('numc' ,
trim(left( symget('numc')||' '
||trim(left(Variable)))));
call symput('nnumc',trim(left(1+symget('nnumc'))));
end;
_nc = _nc + 1; _cvar = trim(left(upcase(scan("&class",_nc))));
end;
run;
data &dsout; set &dsin; run;
%do inumc = 1 %to &nnumc;
%let numcvar = %scan(&numc,&inumc);
data &dsout; set &dsout;
_tempc = trim(left(put(&numcvar,best16.)));
data &dsout; set &dsout;
drop &numcvar;
data &dsout; set &dsout;
rename _tempc=&numcvar;
run;
%end;
%let allv =;
%let nallv = 0;
proc sort data=_Var out=_Var; by Num;
data _null_; set _Var;
call symput('allv' ,
trim(left( symget('allv')||' '
||trim(left(Variable)))));
call symput('nallv',trim(left(1+symget('nallv'))));
run;
data _temp; set &dsout;
data &dsout; if (0); run;
%do ivar = 1 %to &nallv;
%let var = %scan(&allv,&ivar);
data &dsout; merge &dsout _temp(keep=&var);
run;
%end;
%if (&clean) %then %do;
proc datasets library=work nolist;
delete _Var _temp;
run;
%end;
%mend;
/*-------------------------------------------------------------------*/
/* Constrained MC with the best */
/*-------------------------------------------------------------------*/
%macro mcb(data ,
resp ,
mean ,
model = &mean ,
class = &mean ,
alpha = 0.05 ,
out = mcbout,
options = );
/*
/ Retrieve options.
/---------------------------------------------------------------------*/
%let print = 1;
%let clean = 1;
%let iopt = 1;
%do %while(%length(%scan(&options,&iopt)));
%if (%upcase(%scan(&options,&iopt)) = NOPRINT) %then
%let print = 0;
%else %if (%upcase(%scan(&options,&iopt)) = NOCLEAN) %then
%let clean = 0;
%else
%put Warning: Unrecognized option %scan(&options,&iopt).;
%let iopt = %eval(&iopt + 1);
%end;
/*
/ Count number of variables in grouping effect.
/---------------------------------------------------------------------*/
%let ivar = 1;
%do %while(%length(%scan(&mean,&ivar,*)));
%let var&ivar = %upcase(%scan(&mean,&ivar,*));
%let ivar = %eval(&ivar + 1);
%end;
%let nvar = %eval(&ivar - 1);
/*
/ Compute ANOVA and LSMEANS
/---------------------------------------------------------------------*/
ods listing close;
proc mixed data=&data;
class &class;
model &resp = &model;
lsmeans &mean;
ods output LSMeans=&out;
run;
%n2cclass(&out,&out,&class);
ods listing;
data &out; set &out; orig_n = _n_;
proc sort data=&out out=&out; by &mean;
run;
/*
/ Retrieve the levels of the classification variable.
/---------------------------------------------------------------------*/
data &out; set &out;
drop tvalue probt;
length level $ 20;
level = '';
%do ivar = 1 %to &nvar;
level = trim(left(level)) || ' ' || trim(left(&&var&ivar));
%end;
call symput('nlev',trim(left(_n_)));
call symput('lev'||trim(left(_n_)),level);
run;
/*
/ Now, perform Dunnett's comparison-with-control test with each
/ level as the control.
/---------------------------------------------------------------------*/
ods listing close;
proc mixed data=&data;
class &class;
model &resp = &model / dfm=sat;
%do ilev = 1 %to &nlev;
%let control =;
%do ivar = 1 %to &nvar;
%let control = &control "%scan(&&lev&ilev,&ivar)";
%end;
lsmeans &mean / diff=controlu(&control) cl alpha=&alpha
adjust=dunnett;
%end;
ods output Diffs=_mcb;
run;
%n2cclass(_mcb,_mcb,&class);
ods listing;
data _mcb; set _mcb;
length level1 $ 20 level2 $ 20;
level1 = '';
level2 = '';
%do ivar = 1 %to &nvar;
%let v1 = &&var&ivar;
%let v2 = _&&var&ivar;
%if (%length(&v2) > 8) %then
%let var2 = %substr(&v2,1,8);
level1 = trim(left(level1)) || ' ' || trim(left(&v1));
level2 = trim(left(level2)) || ' ' || trim(left(&v2));
%end;
run;
/*
/ Sort results by first and second level, respectively.
/---------------------------------------------------------------------*/
proc sort data=_mcb out=_tmcb1; by level1 level2;
proc transpose data=_tmcb1 out=_tmcb1 prefix=lo;
by level1; var AdjLower;
data _tmcb1; set _tmcb1; ilev = _n_;
proc sort data=_mcb out=_tmcb2; by level2 level1;
proc transpose data=_tmcb2 out=_tmcb2 prefix=lo;
by level2; var AdjLower;
data _tmcb2; set _tmcb2; ilev = _n_;
run;
/*
/ From Hsu (1996), p. 94:
/ Di+ = +( min_{j!=i} m_i - m_j + d^i*s*sqrt(1/n_i + 1/n_j))^+
/ = +(-max_{j!=i} m_j - m_i - d^i*s*sqrt(1/n_i + 1/n_j))^+
/ G = {i : min_{j!=i} m_i - m_j + d^i*s*sqrt(1/n_i + 1/n_j) > 0}
/ Di- = 0 if G = {i}
/ = min_{j!=i} m_i - m_j + d^j*s*sqrt(1/n_i + 1/n_j) otherwise
/---------------------------------------------------------------------*/
data clhi; set _tmcb2; keep level2 clhi ilev;
rename level2=level;
clhi = -max(of lo1-lo%eval(&nlev-1));
if (clhi < 0) then clhi = 0;
data _g; set clhi; if (clhi > 0);
run;
%let ng = 0;
%let g = 0;
data _null_; set _g;
call symput('ng',_n_ );
call symput('g' ,ilev);
run;
data cllo; set _tmcb1; keep level1 cllo ilev;
rename level1=level;
if ((&ng = 1) & (&g = ilev)) then cllo = 0;
else cllo = min(of lo1-lo%eval(&nlev-1));
run;
data cl; merge cllo clhi;
by level;
data &out; merge &out cl;
drop df ilev;
run;
/*
/ Compute RVAL and SVAL. RVAL is just the p-value for Dunnett's
/ test for all means except the best, and SVAL is the maximum RVAL.
/---------------------------------------------------------------------*/
data _slev; set &out; _i_ = _n_;
proc sort data=_slev out=_slev; by descending estimate;
%let ibest = 0;
data _null_; set _slev;
if (_n_ = 1) then call symput('ibest',_i_);
proc sort data=_mcb out=_pval; by level2 adjp;
proc transpose data=_pval out=_pval prefix=p; by level2; var adjp;
data _pval; set _pval; keep level2 rval;
rename level2=level;
if (_n_ = &ibest) then rval = .;
else rval = p1;
proc sort data=_pval out=_spval; by descending rval;
data _null_; set _spval; if (_n_ = 1) then call symput('sval',rval);
data _pval; set _pval;
if (_n_ = &ibest) then sval = &sval;
data &out; merge &out _pval; by level; drop level;
proc sort data=&out out=&out; by orig_n;
data &out; set &out; drop orig_n;
run;
/*
/ Print and clean up.
/---------------------------------------------------------------------*/
%if (&print) %then %do;
proc print uniform data=&out noobs;
run;
%end;
%if (&clean) %then %do;
proc datasets library=work nolist;
delete cllo clhi cl _slev _spval _pval _mcb _tmcb1 _tmcb2 _g;
run;
%end;
%mend;
/*-------------------------------------------------------------------*/
/* Constrained MC with the worst */
/*-------------------------------------------------------------------*/
%macro mcw(data,
resp ,
mean,
model = &mean,
class = &mean,
alpha = 0.05 ,
out = mcbout ,
options = );
/*
/ Retrieve options.
/---------------------------------------------------------------------*/
%let print = 1;
%let clean = 1;
%let iopt = 1;
%do %while(%length(%scan(&options,&iopt)));
%if (%upcase(%scan(&options,&iopt)) = NOPRINT) %then
%let print = 0;
%else %if (%upcase(%scan(&options,&iopt)) = NOCLEAN) %then
%let clean = 0;
%else
%put Warning: Unrecognized option %scan(&options,&iopt).;
%let iopt = %eval(&iopt + 1);
%end;
/*
/ Copy the dataset but reverse the sign of the response, so that
/ the best is the maximum response.
/---------------------------------------------------------------------*/
data _tmpds; set &data; &resp = -&resp; run;
%mcb(_tmpds,
&resp ,
&mean ,
model = &model ,
class = &class ,
alpha = &alpha ,
out = &out ,
options = &options);
/*
/ Reverse the sign of the results, so that the best is again the
/ minimum response.
/---------------------------------------------------------------------*/
data &out; set &out;
rename cllo=cllo;
rename clhi=clhi;
estimate = -estimate;
tvalue = -tvalue;
_temp = -cllo; cllo = -clhi; clhi = _temp; drop _temp;
run;
/*
/ Print and clean up.
/---------------------------------------------------------------------*/
%if (&print) %then %do;
proc print uniform data=&out noobs;
run;
%end;
%if (&clean) %then %do;
proc datasets library=work nolist;
delete _tmpds;
run;
%end;
%mend;
/*-------------------------------------------------------------------*/
/* Constrained MC with the best */
/*-------------------------------------------------------------------*/
%macro umcb(data,
resp ,
mean,
model = &mean,
class = &mean,
alpha = 0.05 ,
out = mcbout ,
method = EH ,
options = );
/*
/ Retrieve options.
/---------------------------------------------------------------------*/
%let print = 1;
%let clean = 1;
%let iopt = 1;
%do %while(%length(%scan(&options,&iopt)));
%if (%upcase(%scan(&options,&iopt)) = NOPRINT) %then
%let print = 0;
%else %if (%upcase(%scan(&options,&iopt)) = NOCLEAN) %then
%let clean = 0;
%else
%put Warning: Unrecognized option %scan(&options,&iopt).;
%let iopt = %eval(&iopt + 1);
%end;
/*
/ Count number of variables in grouping effect.
/---------------------------------------------------------------------*/
%let ivar = 1;
%do %while(%length(%scan(&mean,&ivar,*)));
%let var&ivar = %upcase(%scan(&mean,&ivar,*));
%let ivar = %eval(&ivar + 1);
%end;
%let nvar = %eval(&ivar - 1);
/*
/ Compute ANOVA and LSMEANS
/---------------------------------------------------------------------*/
ods listing close;
proc mixed data=&data;
class &class;
model &resp = &model;
lsmeans &mean;
ods output LSMeans=&out;
run;
%n2cclass(&out,&out,&class);
ods listing;
data &out; set &out; orig_n = _n_;
proc sort data=&out out=&out; by &mean;
run;
/*
/ Retrieve the levels of the classification variable.
/---------------------------------------------------------------------*/
data &out; set &out;
drop tvalue probt;
length level $ 20;
level = '';
%do ivar = 1 %to &nvar;
level = trim(left(level)) || ' ' || trim(left(&&var&ivar));
%end;
call symput('nlev',trim(left(_n_)));
call symput('lev'||trim(left(_n_)),level);
run;
%if (%upcase(&method) = TK) %then %do;
ods listing close;
proc mixed data=&data;
class &class;
model &resp = &model;
lsmeans &mean / diff=all cl alpha=&alpha adjust=tukey;
ods output Diffs=_mcb;
run;
%n2cclass(_mcb,_mcb,&class);
ods listing;
proc sort data=_mcb out=_mcb;
by &mean _&mean;
run;
/*
/ Add reverse differences.
/---------------------------------------------------------------------*/
data _mcb; set _mcb; keep level1 level2 AdjLower AdjUpper adjp;
length level1 $ 20 level2 $ 20;
level1 = '';
level2 = '';
%do ivar = 1 %to &nvar;
%let v1 = &&var&ivar;
%let v2 = _&&var&ivar;
%if (%length(&v2) > 8) %then
%let var2 = %substr(&v2,1,8);
level1 = trim(left(level1)) || ' ' || trim(left(&v1));
level2 = trim(left(level2)) || ' ' || trim(left(&v2));
%end;
output;
_tmplev = level1; level1 = level2; level2 = _tmplev;
_tmpcl = -AdjLower; AdjLower = -AdjUpper; AdjUpper = _tmpcl;
output;
run;
/*
/ Confidence limits are the minimum lower and upper CL's for each
/ level.
/---------------------------------------------------------------------*/
proc sort data=_mcb out=_mcb; by level1 level2;
proc transpose data=_mcb out=cllo prefix=lo;
by level1; var AdjLower;
proc transpose data=_mcb out=clhi prefix=hi;
by level1; var AdjUpper;
data cllo; set cllo;
rename level1=level;
cllo = min(of lo1-lo%eval(&nlev-1));
data clhi; set clhi;
rename level1=level;
clhi = min(of hi1-hi%eval(&nlev-1));
data cl; merge cllo(keep=level cllo) clhi(keep=level clhi);
run;
data &out; merge &out cl; drop level;
run;
%if (&clean) %then %do;
proc datasets library=work nolist;
delete _mcb cllo clhi cl;
run;
%end;
%end;
%else %do;
/*
/ Now, perform Dunnett's comparison-with-control test with each
/ level as the control.
/---------------------------------------------------------------------*/
ods listing close;
proc mixed data=&data;
class &class;
model &resp = &model / dfm=sat;
%do ilev = 1 %to &nlev;
%let control =;
%do ivar = 1 %to &nvar;
%let control = &control "%scan(&&lev&ilev,&ivar)";
%end;
lsmeans &mean / diff=control(&control) cl alpha=&alpha
adjust=dunnett;
%end;
ods output Diffs=_mcb;
run;
%n2cclass(_mcb,_mcb,&class);
ods listing;
data _mcb; set _mcb;
length level1 $ 20 level2 $ 20;
level1 = '';
level2 = '';
%do ivar = 1 %to &nvar;
%let v1 = &&var&ivar;
%let v2 = _&&var&ivar;
%if (%length(&v2) > 8) %then
%let var2 = %substr(&v2,1,8);
level1 = trim(left(level1)) || ' ' || trim(left(&v1));
level2 = trim(left(level2)) || ' ' || trim(left(&v2));
%end;
proc sort data=_mcb out=_mcb; by level2 level1;
data cl; keep cllo clhi;
array m{&nlev,&nlev}; /* m[i1]-m[i2] - |d|^i2*s[i1,i2] */
array p{&nlev,&nlev}; /* m[i1]-m[i2] + |d|^i2*s[i1,i2] */
array s{&nlev};
array l{&nlev};
array u{&nlev};
do i = 1 to &nlev; do j = 1 to &nlev;
m[i,j] = .; p[i,j] = .;
end; end;
do obs = 1 to %eval(&nlev*(&nlev-1));
set _mcb point=obs;
j = mod((obs-1),%eval(&nlev-1)) + 1;
i2 = int((obs-1)/%eval(&nlev-1)) + 1;
if (j < i2) then i1 = j;
else i1 = j + 1;
m[i1,i2] = AdjLower;
p[i1,i2] = AdjUpper;
end;
/*
/ From Hsu (1996), p. 120:
/ S = {i : min_{j!=i} m_i - m_j + |d|^i*s[i,j] > 0}
/ = {i : min_{j!=i} -(m_j - m_i - |d|^i*s[i,j]) > 0}
/ = {i : min_{j!=i} -m[j,i] > 0}
/---------------------------------------------------------------------*/
ns = 0;
do i = 1 to &nlev;
minmmji = 1e12;
do j = 1 to &nlev; if (j ^= i) then do;
if (-m[j,i] < minmmji) then minmmji = -m[j,i];
end; end;
s[i] = (minmmji > 0);
ns = ns + s[i];
end;
/*
/ From Hsu (1996), p. 115:
/ Lij = (i ^= j) * (m_i - m_j + |d|^j*s[i,j])
/ = (i ^= j) * p[i,j]
/ Li = min_{j in S} Lij
/
/ Uij = (i ^= j) * -(m_i - m_j + |d|^j*s[i,j])^-
/ = (i ^= j) * min(0,p[i,j])
/ Ui = max_{j in S} Uij
put "Edwards-Hsu intervals";
do i = 1 to &nlev;
li = 1e12;
do j = 1 to &nlev; if (s[j]) then do;
if (i = j) then lij = 0;
else lij = m[i,j];
if (lij < li) then li = lij;
end; end;
ui = -1e12;
do j = 1 to &nlev; if (s[j]) then do;
if (i = j) then uij = 0;
else uij = min(0,p[i,j]);
if (uij > ui) then ui = uij;
end; end;
put li 7.3 " < mu" i 1. " - max_j muj < " ui 7.3;
end;
/---------------------------------------------------------------------*/
/*
/ From Hsu (1996), p. 120:
/ If S = {i} then
/ Li* = (min_{j!=i} m_i - m_j - |d|^i*s[i,j] )^+
/ = (min_{j!=i} -(m_j - m_i + |d|^i*s[i,j]))^+
/ = (min_{j!=i} -p[j,i])^+
/ Otherwise
/ Li* = min_{j in S,j!=i} m_i - m_j - |d|^j*s[i,j]
/ = min_{j in S,j!=i} m[i,j]
/---------------------------------------------------------------------*/
do i = 1 to &nlev;
if ((ns = 1) & s[i]) then do;
minmpji = 1e12;
do j = 1 to &nlev; if (j ^= i) then do;
if (-p[j,i] < minmpji) then minmpji = -p[j,i];
end; end;
l[i] = max(0,minmpji);
end;
else do;
minpmij = 1e12;
do j = 1 to &nlev; if (s[j] & (j ^= i)) then do;
if (m[i,j] < minpmij) then minpmij = m[i,j];
end; end;
l[i] = minpmij;
end;
end;
/*
/ From Hsu (1996), p. 120:
/ If i in S then
/ Ui* = min_{j!=i} m_i - m_j + |d|^i*s[i,j]
/ = min_{j!=i} -(m_j - m_i - |d|^i*s[i,j])
/ = min_{j!=i} -m[j,i]
/ Otherwise
/ Ui* = -(max_{j in S,} m_i - m_j + |d|^j*s[i,j])^-
/ = -(max_{j in S,} p[i,j])^-
/---------------------------------------------------------------------*/
do i = 1 to &nlev;
if (s[i]) then do;
minmmji = 1e12;
do j = 1 to &nlev; if (j ^= i) then do;
if (-m[j,i] < minmmji) then minmmji = -m[j,i];
end; end;
u[i] = minmmji;
end;
else do;
minppij = -1e12;
do j = 1 to &nlev; if (s[j]) then do;
if (p[i,j] > minppij) then minppij = p[i,j];
end; end;
u[i] = minppij;
end;
end;
do i = 1 to &nlev;
cllo = l{i}; clhi = u{i};
output;
end;
stop;
data &out; merge &out cl; drop level;
run;
%if (&clean) %then %do;
proc datasets library=work nolist;
delete _mcb cl;
run;
%end;
%end;
proc sort data=&out out=&out; by orig_n;
data &out; set &out; drop orig_n;
run;
/*
/ Print and clean up.
/---------------------------------------------------------------------*/
%if (&print) %then %do;
proc print uniform data=&out noobs;
run;
%end;
%mend;
/*-------------------------------------------------------------------*/
/* Unconstrained MC with the worst */
/*-------------------------------------------------------------------*/
%macro umcw(data,
resp ,
mean,
model = &mean ,
class = &mean ,
alpha = 0.05 ,
out = mcbout,
method = EH ,
options = );
/*
/ Retrieve options.
/---------------------------------------------------------------------*/
%let print = 1;
%let clean = 1;
%let iopt = 1;
%do %while(%length(%scan(&options,&iopt)));
%if (%upcase(%scan(&options,&iopt)) = NOPRINT) %then
%let print = 0;
%else %if (%upcase(%scan(&options,&iopt)) = NOCLEAN) %then
%let clean = 0;
%else
%put Warning: Unrecognized option %scan(&options,&iopt).;
%let iopt = %eval(&iopt + 1);
%end;
/*
/ Copy the dataset but reverse the sign of the response, so that
/ the best is the maximum response.
/---------------------------------------------------------------------*/
data _tmpds; set &data; &resp = -&resp; run;
%umcb(_tmpds,
&resp ,
&mean ,
model = &model ,
class = &class ,
alpha = &alpha ,
out = &out ,
method = &method ,
options = &options);
/*
/ Reverse the sign of the results, so that the best is again the
/ minimum response.
/---------------------------------------------------------------------*/
data &out; set &out;
rename cllo=cllo;
rename clhi=clhi;
estimate = -estimate;
tvalue = -tvalue;
_temp = -cllo; cllo = -clhi; clhi = _temp; drop _temp;
run;
/*
/ Print and clean up.
/---------------------------------------------------------------------*/
%if (&print) %then %do;
proc print uniform data=&out noobs;
run;
%end;
%if (&clean) %then %do;
proc datasets library=work nolist;
delete _tmpds;
run;
%end;
%mend;