/*---------------------------------------------------------------------- The NLMeans macro provides multiple comparisons (pairwise, sequential, or against a control) on the mean scale among the levels of an effect in a generalized linear model. More complex contrasts of means and ratios of pairs of means can also be estimated and tested. In addition to Base SASŪ, SAS/STATŪ and the NLEST macro (Version 1.8 or later) are required to run this macro. Both the NLMeans macro and the NLEST macro that it calls must be defined in the current SAS session before using the NLMeans macro. For full documentation of the NLMeans macro, including descriptions of all available options and examples, see http://support.sas.com/kb/62362. ------------------------------------------------------------------------ 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. ----------------------------------------------------------------------*/ %macro NLMeans(version, instore=, inest=, incovb=, coef=, where=, link=, diff=all, contrasts=, df=, alpha=0.05, title=, null=0, covdrop=, options=NOJOINT NONAMES NOREVERSE NORATIO DIFINFNS NOAPPEND PRINT) / minoperator mindelimiter=' '; %local notesopt; %let notesopt = %sysfunc(getoption(notes)); %let timenlm = %sysfunc(datetime()); %let _version=2.0; %if &version ne %then %put NOTE: &sysmacroname macro Version &_version..; %let newchk=1; %if %index(%upcase(&version),DEBUG) %then %do; options notes mprint %if %index(%upcase(&version),DEBUG2) %then mlogic symbolgen; ; ods select all; %put _user_; %let dbgopt=debug,; %if %index(%upcase(&version),DEBUG2) %then %let dbgopt=debug2,; %end; %else %do; options nonotes nomprint nomlogic nosymbolgen; ods exclude all; %let dbgopt=; %if &version ne %then %let dbgopt=v,; %if %index(%upcase(&version),NOCHK) %then %do; %let newchk=0; %let dbgopt=NOCHK,; %end; %end; /* Check for newer version */ %if &newchk %then %do; %let _notfound=0; %let _newver=0; filename _ver url 'http://ftp.sas.com/techsup/download/stat/versions.dat' termstr=crlf; data _null_; infile _ver end=_eof; input name:$15. ver; if upcase(name)="&sysmacroname" then do; call symput("_newver",ver); stop; end; if _eof then call symput("_notfound",1); run; options notes; %if &syserr ne 0 or &_notfound=1 or &_newver=0 %then %put NOTE: Unable to check for newer version of &sysmacroname macro.; %else %if %sysevalf(&_newver > &_version) %then %do; %put NOTE: A newer version of the &sysmacroname macro is available at; %put NOTE- this location: http://support.sas.com/ ; %end; %if %index(%upcase(&version),DEBUG)=0 %then options nonotes;; %end; /* Check inputs */ %if (&instore= and &inest= and &incovb=) or (&instore= and (&inest= or &incovb=)) %then %do; %put ERROR: Either INSTORE= or both INEST= and INCOVB= must be specified.; %goto exit; %end; %if &coef= or %sysfunc(exist(&coef)) ne 1 %then %do; %put ERROR: COEF= data set not specified or not found.; %goto exit; %end; %if &contrasts ne and %sysfunc(exist(&contrasts)) ne 1 %then %do; %put ERROR: CONTRASTS= data set not found.; %goto exit; %end; %if %quote(&null)= %then %let null=0; %else %do; %if (%sysfunc(verify(%sysfunc(trim(%sysfunc(left(&null)))),'0123456789.-'))>0) + (%sysfunc(findc(%sysfunc(trim(%sysfunc(left(&null)))),'-'))>1) + (%sysfunc(count(&null,-))>1) + (%sysfunc(count(&null,.))>1) %then %do; %put ERROR: The NULL= value must be a numeric value.; %goto exit; %end; %if %quote(&null)=. or %quote(&null)=%quote(-) %then %do; %put ERROR: The NULL= value must be a numeric value.; %goto exit; %end; %end; %if %quote(&df) ne %then %do; %if (%sysfunc(verify(%sysfunc(trim(%sysfunc(left(&df)))),'0123456789.'))>0) + (%sysfunc(count(&df,.))>1) %then %do; %put ERROR: The DF= value must be a non-zero, positive value.; %goto exit; %end; %if %sysevalf(&df=0 or &df=.) %then %do; %put ERROR: The DF= value must be a non-zero, positive value.; %goto exit; %end; %end; %if %sysevalf(&alpha<=0 or &alpha>=1) %then %do; %put ERROR: The ALPHA= value must be between 0 and 1.; %goto exit; %end; /* Verify COEF= data set has necessary variables */ %let status=ok; %let dsid=%sysfunc(open(&coef)); %if &dsid %then %do; %if %sysfunc(varnum(&dsid,LMATRIX))=0 or %sysfunc(varnum(&dsid,ROW1))=0 or (%sysfunc(varnum(&dsid,PARAMETER))=0 and %sysfunc(varnum(&dsid,EFFECT))=0) %then %do; %put ERROR: Invalid COEF= data set. At least one variable Parameter (or; %put ERROR- Effect), LMatrix, or Row1 not found.; %let status=input_err; %end; %let rc=%sysfunc(close(&dsid)); %if &status=input_err %then %goto exit; %end; %else %do; %put ERROR: Could not open COEF= data set.; %goto exit; %end; %let validopts=JOINT NOJOINT NAMES NONAMES REVERSE NOREVERSE RATIO NORATIO DIFINFNS DIFALL APPEND NOAPPEND PRINT NOPRINT; %let joint=0; %let names=0; %let rdiff=0; %let ratio=0; %let difall=0; %let append=0; %let prt=1; %let i=1; %do %while (%scan(&options,&i) ne %str() ); %let option&i=%upcase(%scan(&options,&i)); %if &&option&i=JOINT %then %let joint=1; %else %if &&option&i=NAMES %then %let names=1; %else %if &&option&i=REVERSE %then %let rdiff=1; %else %if &&option&i=RATIO %then %let ratio=1; %else %if &&option&i=DIFALL %then %let difall=1; %else %if &&option&i=APPEND %then %let append=1; %else %if &&option&i=NOPRINT %then %let prt=0; %else %do; %let chk=%eval(&&option&i in &validopts); %if not &chk %then %do; %put ERROR: Valid values of OPTIONS= are &validopts..; %goto exit; %end; %end; %let i=%eval(&i+1); %end; %if %index(%quote(&title),%str(%")) or %index(%quote(&title),%str(%')) %then %do; %put ERROR: Do not use quotes (%str(%") or %str(%')) in TITLE=.; %goto exit; %end; %if &df= %then %let dfopt=; %else %let dfopt=, df=&df; %if %quote(&null)= %then %let nullopt=; %else %let nullopt=, null=&null; %if &alpha= %then %let alphaopt=; %else %let alphaopt=, alpha=α %if %quote(&title)= %then %let ttlopt=; %else %let ttlopt=, title=%quote(&title); %if &names=1 %then %let namopt=; %else %let namopt=, listnames=no; %if &prt=1 %then %let prtopt=; %else %let prtopt=, print=no; %if &rdiff %then %do; %let locode=-1; %let hicode=1; %end; %else %do; %let locode=1; %let hicode=-1; %end; /* Determine inverse link function */ %let link=%quote(%upcase(&link)); %if &link= %then %do; %put ERROR: The LINK= option is required.; %goto exit; %end; %else %if &link=LOGIT or &link=CUMLOGIT or &link=CLOGIT %then %let ilink="logistic("||trim(left(exp))||")"; %else %if &link=PROBIT or &link=CUMPROBIT or &link=CPROBIT %then %let ilink="probnorm("||trim(left(exp))||")"; %else %if &link=CLL or &link=CLOGLOG or &link=CUMCLL or &link=CCLL or &link=CCLOGLOG or &link=CUMCLOGLOG %then %let ilink="(1-exp(-exp("||trim(left(exp))||")))"; %else %if &link=LOGLOG or &link=CUMLOGLOG %then %let ilink="exp(-exp(-("||trim(left(exp))||")))"; %else %if &link=LOG or &link=GLOGIT %then %let ilink="exp("||trim(left(exp))||")"; %else %if %substr(&link,1,3)=POW %then %do; %let power=%quote(%substr(&link,6)); %if &power=0 %then %let ilink="exp("||trim(left(exp))||")"; %else %let ilink="exp(log("||trim(left(exp))||")/&power)"; %end; %else %if %upcase(&link)=IDENTITY %then %do; options notes; %put NOTE: This macro is not needed when LINK=IDENTITY.; %put NOTE- Use results of LSMEANS, SLICE or ESTIMATE.; options nonotes; %goto exit; %end; %else %do; %put ERROR: Valid LINK= values are LOG, LOGIT, CUMLOGIT, PROBIT, ; %put ERROR- CUMPROBIT, CLL, CUMCLL, LOGLOG, CUMLOGLOG, GLOGIT, ; %put ERROR- or POWERx where x is a number.; %goto exit; %end; %if &difall and &link=GLOGIT %then %do; %put ERROR: OPTIONS=DIFALL cannot be used with LINK=GLOGIT.; %goto exit; %end; /* Determine number of response functions for multinomial models */ %let nfn=1; %let cumlinks=CUMLOGIT CLOGIT CUMPROBIT CPROBIT CUMCLL CCLL CCLOGLOG CUMCLOGLOG CUMLOGLOG GLOGIT; %let cumchk=%eval(&link in &cumlinks); /* Subset the coef= data set if a where= condition is specified */ data _coefsub; set &coef; %if %quote(&where) ne %then where %str(&where);; run; %if &syserr %then %do; %put ERROR: Subsetting the COEF= data set failed.; %goto exit; %end; %if &cumchk %then %do; data _null_; set _coefsub end=eof; if index(Parameter,"Intercept") or index(Effect,"Intercept") then nfn+1; else do; call symput('nfn',cats(nfn)); stop; end; run; %end; /* Get number of sets of estimates */ %if &joint %then %do; %let nsets=1; %let set1=1; proc transpose data=_coefsub out=_coef(drop=_name_); by LMatrix; var row:; run; proc transpose data=_coef out=_coef(drop=_name_) prefix=Row; var col:; run; data _coef; set _coef; LMatrix=1; run; %end; %else %do; proc freq data=_coefsub nlevels; table LMatrix / out=_setnums noprint; ods output nlevels=_nlvl; run; data _null_; set _nlvl; call symput('nsets',cats(nlevels)); run; data _null_; set _setnums; call symput(cats('set',_n_),cats(lmatrix)); run; data _coef; set _coefsub; run; %end; %if &contrasts= %then %do; %let diff=%upcase(&diff); %let ndifvals=%sysfunc(countw(&diff)); %if &ndifvals ne 1 and &ndifvals ne &nsets %then %do; %put ERROR: There are &nsets sets of estimates. Specify a single value; %put ERROR- in DIFF= to apply to all sets, or one value for each set.; %put ERROR- Valid values are ALL, SEQ, or a positive integer.; %goto exit; %end; %let i=1; %do %while (%scan(&diff,&i) ne %str() ); %let diferr=0; %let setdif&i=%scan(&diff,&i); %if &&setdif&i ne ALL and &&setdif&i ne SEQ %then %do; %let diferr=1; %if %sysevalf(&&setdif&i>0 and &&setdif&i<1e8) %then %do; %if %sysevalf(%sysfunc(mod(&&setdif&i,1))=0) %then %let diferr=0; %end; %if &diferr %then %do; %put ERROR: Invalid DIFF= value found. Valid values are ALL, SEQ or; %put ERROR- a positive integer.; %goto exit; %end; %end; %let i=%eval(&i+1); %end; %end; /* Process each set of estimates in COEF= data set */ %do s=1 %to &nsets; /* Construct mean expressions for NLEST macro */ proc transpose data=_coef %if &joint=0 %then (where=(LMatrix=&&set&s)); out=_ct(keep=col: where=(col1 ne .)); by lmatrix; var row:; run; data _glognum(keep=fn1) _glogden(keep=sum); set _ct nobs=nlsm end=last; array lb (*) col:; array fn (&nfn) $32767; length exp sum sum1 $32767; do i=1 to &nfn; *initialize &nfn fns for each row in _ct; if lb(1) ne 0 then fn(i)=catx("*",lb(1),cats("B_p",i)); end; do i=2 to dim(lb); *cols; do j=1 to &nfn; *fns; if lb(i) ne 0 then fn(j)=catx("+",fn(j),catx("*",lb(i),cats("B_p",i+j-1))); end; end; do i=1 to &nfn; exp=fn(i); fn(i)=&ilink; if i=1 then sum1=fn(i); else sum1=catx("+",sum1,fn(i)); end; sum="(1+"||trim(left(sum1))||")"; call symput("nlsm",cats(nlsm)); output _glognum; if _n_<=nlsm/&nfn then output _glogden; %if &link=GLOGIT %then %do; if last and &nfn>1 then do; fn1="1"; do i=1 to nlsm/&nfn; output _glognum; end; end; %end; run; data _prexp; length exp $32767; %if &link=GLOGIT %then %do; set %do i=1 %to &nfn+1; _glogden %end; ; merge _glognum; exp=catx("/",fn1,sum); %if &ratio %then exp="("||trim(left(exp))||")";; %end; %else %do; set _glognum; exp=fn1; %end; keep exp; run; proc transpose data=_prexp out=_lbt(keep=col:); var exp; run; /* With DIFALL: don't restrict differencing to be within response functions */ %if &difall %then %let nfn=1; /* Number of estimates per response function in multinomial models */ %let nperfn=%sysevalf(&nlsm/&nfn); /* Comparisons: pairwise diffs, diff with control, or input contrasts */ %if &link=GLOGIT and &nfn>1 %then %do; %let nk=%sysevalf((&nfn+1)*&nperfn); %let badglgcoef=0; %if %sysevalf(%sysfunc(mod(&nperfn,1)) ne 0) %then %let badglgcoef=1; %else %do; data %do f=1 %to &nfn; %let first=%eval((&f-1)*&nperfn+1); %let last=%eval(&f*&nperfn); _glg&f(keep=row&first - row&last) %end; ; set _coef %if &joint=0 %then (where=(LMatrix=&&set&s)); ; %do i=1 %to &nfn-1; if mod(_n_,&nfn)=&i then output _glg&i; %end; if mod(_n_,&nfn)=0 then output _glg&nfn; run; %do f=2 %to &nfn; proc compare base=_glg1 compare=_glg&f noprint out=_glgcomp outnoequal; var row1-row&nperfn; %let first=%eval((&f-1)*&nperfn+1); %let last=%eval(&f*&nperfn); with row&first-row&last; run; %let uneq=0; data _null_; dsid=open("_glgcomp"); nobs=attrn(dsid, "nobs"); if nobs>0 then call symput ("uneq",1); rc=close(dsid); run; %if &uneq %then %let badglgcoef=1; %end; %end; %if &badglgcoef %then %do; %put ERROR: For LINK=GLOGIT, the requested estimates must be estimated in; %put ERROR- each of the &nfn logits. If the ESTIMATE statement was used,; %put ERROR- the CATEGORY=JOINT option is required.; %put ERROR: The above message is for estimate set &&set&s...; %goto endset; %end; %end; %else %let nk=&nlsm; %if &nperfn=1 %then %do; %put ERROR: Only one estimate - no differencing possible; %put ERROR: The above message is for estimate set &&set&s...; %goto endset; %end; %let label=0; %if &contrasts= %then %do; %if &ndifvals=1 %then %let setdif&s=&setdif1; %if %sysevalf(&&setdif&s>0 and &&setdif&s<1e8) %then %do; %if &&setdif&s>&nperfn %then %do; %put ERROR: The DIFF= value, &&setdif&s, for estimate set &&set&s must be an; %put ERROR- integer between 1 to the the number of estimates, &nperfn..; %goto exit; %end; %end; data _cntrs(keep=k:); array k (&nk) (&nk*0); do f=1 to &nfn %if &link=GLOGIT and &nfn>1 %then +1; ; first=(f-1)*&nperfn+1; last=f*&nperfn; %if &&setdif&s=ALL %then %do; do i=first to last-1; do j=i+1 to last; do h=1 to last; if h=i then k(h)=&locode; else if h=j then k(h)=&hicode; else k(h)=0; end; output; end; end; %end; %else %if &&setdif&s=SEQ %then %do; do i=first to last-1; do h=1 to last; if h=i then k(h)=&locode; else if h=i+1 then k(h)=&hicode; else k(h)=0; end; output; end; %end; %else %do; do j=first to last; do h=1 to last; if h=(first-1)+&&setdif&s then k(h)=&locode; else if h=j then k(h)=&hicode; else k(h)=0; end; if j ne (first-1)+&&setdif&s then output; end; %end; end; run; %let noobs=0; data _null_; dsid=open("_cntrs"); nobs=attrn(dsid, "nobs"); if nobs=0 then call symput ("noobs",1); rc=close(dsid); run; %if &noobs %then %do; %put ERROR: Could not produce differencing coefficients. If a multinomial; %put ERROR- model was used you might need to specify OPTIONS=DIFALL.; %put ERROR: The above message is for estimate set &&set&s...; %goto endset; %end; %end; %else %do; %let misscoef=0; %let miss_set=0; data _cntrs; set &contrasts; array k (*) k:; nk=dim(k); if set=&&set&s; if (nk-nmiss(of k:))<&nk then do; call symput("misscoef",1); stop; end; run; data _null_; dsid=open("_cntrs"); nobs=attrn(dsid, "nobs"); if nobs=0 then call symput ("miss_set",1); if varnum(dsid,"LABEL") then call symput("label",1); rc=close(dsid); run; %if &misscoef %then %do; %put ERROR: Too few contrast coefficients found in CONTRASTS= data set; %put ERROR- for estimate set &&set&s (LMatrix=&&set&s) in COEF= data set.; %put ERROR- No analysis done for estimate set &&set&s...; %goto endset; %end; %if &miss_set %then %do; %put ERROR: No contrast coefficients found in CONTRASTS= data found for; %put ERROR- estimate set &&set&s (LMatrix=&&set&s) in COEF= data set.; %goto endset; %end; %if &ratio %then %do; %let badcoef=0; data _null_; set _cntrs; badcoef=0; cnt1=0; cntm1=0; array k (*) k:; do i=1 to dim(k); if k(i) not in(-1,0,1,.) then badcoef=1; if k(i)=1 then cnt1+1; if k(i)=-1 then cntm1+1; end; if cnt1 ne 1 or cntm1 ne 1 then badcoef=1; if badcoef then call symput("badcoef","1"); run; %if &badcoef %then %do; %put ERROR: OPTIONS=RATIO requires all contrasts to contain a single 1; %put ERROR- for the numerator, a single -1 for the denominator, and; %put ERROR- zeros otherwise.; %put ERROR: The above message is for estimate set &&set&s...; %goto endset; %end; %end; %end; /* Construct fdata= data set for NLEST macro */ data _fd; if _n_=1 then set _lbt; set _cntrs; array lb (*) col:; array k (*) k:; length f $32767; %if &label=0 %then length label $32767;; do i=1 to dim(lb); %if &label=0 %then label=catx(" ",label,k(i));; %if &ratio %then %do; if k(i)=1 then num=lb(i); if k(i)=-1 then den=lb(i); f=catx("/",num,den); %end; %else if k(i) ne 0 then f=catx("+",f,catx("*",k(i),lb(i))); ; end; %if &ratio and &label=0 %then %do; negpos=index(label,"-"); substr(label,negpos,1)="/"; %end; keep label f; run; /* Call NLEST macro to estimate/test the contrasts */ %if &instore ne %then %NLEST(&dbgopt instore=&instore, where=&where, covdrop=&covdrop, fdata=_fd &dfopt &alphaopt &namopt &ttlopt &prtopt &nullopt); %else %if &inest ne and &incovb ne %then %NLEST(&dbgopt inest=&inest, incovb=&incovb, where=&where, covdrop=&covdrop, fdata=_fd &dfopt &alphaopt &namopt &ttlopt &prtopt &nullopt); /* Append results sets if requested or save as separate files */ %if &append %then %do; %if %sysfunc(exist(est_all)) and &s=1 %then %do; options notes; %put NOTE: Data set EST_ALL already exists. Appending to it.; %if %index(%upcase(&version),DEBUG)=0 %then options nonotes;; %end; proc append base=est_all data=est; run; %end; %else %if &nsets > 1 %then %do; data est&s; set est; run; %end; /* End of processing set of estimates */ %endset: %end; %exit: %if %index(%upcase(&version),DEBUG) %then %do; options nomprint nomlogic nosymbolgen; %put _user_; %end; %else %do; proc datasets nolist nowarn; delete _coef: _ct _setnums _nlvl _glognum _glogden _prexp _cntrs _fd _lbt _rows _nrows _glg:; run; quit; %end; ods select all; options ¬esopt; %let timenlm=%sysfunc(round(%sysevalf(%sysfunc(datetime())-&timenlm),0.01)); %put NOTE: The &sysmacroname macro used &timenlm seconds.; %mend;