/****************************************************************/ /* S A S A U T O C A L L L I B R A R Y */ /* */ /* NAME: NLMeans */ /* TITLE: Differences, ratios or contrasts of means in */ /* generalized linear models */ /* PRODUCT: STAT */ /* SYSTEM: ALL */ /* KEYS: multiple comparisons contrasts */ /* PROCS: NLMIXED */ /* DATA: */ /* UPDATE: 04Oct2021 */ /* REF: */ /* MISC: The following SAS products are required to run */ /* this macro: BASE, STAT. The NLEST macro */ /* (Version 1.8 or later) is also required. */ /* */ /****************************************************************/ /*---------------------------------------------------------------------- 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. Both the NLMeans macro, and the NLEST macro that it calls, must be defined in the current SAS session before using the NLMeans macro. The NLMeans macro always attempts to check for a later version of itself. If it is unable to do this (such as if there is no active internet connection available), the macro will issue the following message in the log: NOTE: Unable to check for newer version The computations performed by the macro are not affected by the appearance of this message. The NLMeans macro can be used after fitting a generalized linear model (GLM) and estimating multiple linear combinations of model parameters, LBeta, which are response means when the inverse of the GLM link function, g, is applied. That is, mean = ginv(LBeta), where ginv is the inverse of g. GLMs can be fit by several SAS/STAT procedures including GENMOD, GLIMMIX, LOGISTIC, PROBIT, SURVEYLOGISTIC, and GEE. LBeta estimates are available using the LSMEANS, SLICE, and ESTIMATE statements in the modeling procedure. The ILINK option in these statements provides estimates of the means. If the STORE statement is used to save the fitted model, these statements are also available at any time using PROC PLM. While the DIFF option in the LSMEANS and SLICE statements provide pairwise differences on the link scale, LBeta_i-LBeta_j, differences on the mean scale, are not available. Similarly, in an ESTIMATE statement that defines a difference of LBeta the ILINK option applies the inverse of the link function to the difference, ginv(LBeta_i-LBeta_j) rather than computing the difference of the inverse linked estimates ginv(LBeta_i)-ginv(LBeta_j) = mean_i-mean_j. The quantity ginv(LBeta_i-LBeta_j) is generally only of interest in the case where the link function, g, is the log. In this case, the ILINK (or EXP) option estimates the ratio of means. Note that the NLMeans macro is not needed for models that use the identity link. This includes models fit by the REG, GLM, MIXED, or ORTHOREG procedures and others. For these models, differences of LBeta are differences of means. Consequently, the results of the DIFF option in the LSMEANS or SLICE statement, or the results of an ESTIMATE statement that defines a difference of LBeta directly provide the differences of means. You can also use the LSMESTIMATE statement to estimate or test differences or contrasts of means. The NLMeans macro can be used to provide estimates and tests of differences of means, ratios of means, or more complex contrasts of means. Standard errors are obtained using the delta method. To do this, you supply the saved model and a data set containing the coefficients used by the LSMEANS, SLICE, or ESTIMATE statements. You also indicate the link function used in the model. The model is best saved using the STORE statement in the modeling procedure. The coefficients can be saved by including the E option in any LSMEANS, SLICE, or ESTIMATE statement specified in the procedure, and by including an ODS OUTPUT statement to save the displayed table of coefficients in a data set. In most cases, the name of the coefficients table is Coef, so the following statement saves it in a data set. ods output coef=data-set-name; See the list of macro parameters below for details about how to provide the saved model and coefficients to the macro, about how to request means, ratios, or contrasts, and other options. The macro can process more than one set of estimates. Multiple sets of estimates occur when the modeling procedure includes: a SLICE statement, multiple LSMEANS or ESTIMATE statements, or a combination of these statements. The macro processes each set and a table of results is displayed for each set of estimates. By default, the results from each set are saved in separate data sets named EST1, EST2, EST3, ... . Specify options=append to create a single data set named EST_ALL of all results from all sets. If EST_ALL already exists, new results are appended to it. The last results set is also stored in data set EST. For ordinal multinomial models, the macro by default produces comparisons within each of the response functions. To do comparisons across all response functions, specify options=difall. For nominal multinomial models (link= glogit), the macro produces comparisons within each of the response level probabilities. When specifying a contrasts= data set, each contrast (row) of the data set must contain coefficients for all estimates in all response functions (ordinal) or in all response levels (nominal). Since the LSMEANS and SLICE statements require GLM parameterized CLASS variables (PARAM=GLM in the CLASS statement), the NLEST macro called by the NLMeans macro will typically display the following Warning message in this log. This Warning can be ignored. WARNING: The final Hessian matrix is not positive definite, and therefore the estimated covariance matrix is not full rank and may be unreliable. The variance of some parameter estimates is zero or some parameters are linearly related to other parameters. LIMITATION: The NLMeans macro is not intended for use with survival models (LIFEREG, PHREG, SURVEYPHREG), models from PROC CATMOD, or from SAS/ETS procedures. The macro cannot be used to compare means from zero-inflated models fit by PROC GENMOD. For models fit by the GENMOD or GLIMMIX procedures, use the LSMEANS, SLICE, or ESTIMATE statements in PROC PLM after saving the fitted model using the STORE statement. Do not use those statements in GENMOD or GLIMMIX. Note that the NLMeans macro cannot use coefficients produced by the LSMESTIMATE statement. The macro requires the coefficients that define the individual means. The use of inest= and incovb= should rarely, if ever, be necessary. If used, it is likely that modifications will be needed to those input data sets particularly when GLM parameterization of CLASS variables is used. In such cases, it might prove more convenient to use a full-rank parameterization and the NLEST macro directly. ------------------------------------------------------------------------ EXAMPLE: Estimate difference in group means in a log-linked gamma model Using data from the example titled "Gamma Distribution Applied to Life Data" in the GENMOD documentation, these statements fit the log-linked gamma model. The STORE statement saves the model in an item store for later use by the NLMeans macro. The LSMEANS statement with the ILINK option is used in PROC PLM to estimate the manufacturer means. proc genmod data = lifdat; class mfg / param=glm; model lifetime = mfg / dist=gamma link=log; store out=gammod; run; proc plm restore=gammod; lsmeans mfg / e ilink diff exp; ods output coef=coeffs; run; The difference of the manufacturer means is estimated by this call of the NLMeans macro. %NLMeans(instore=gammod, coef=coeffs, link=log, title=Difference of MFG means) The ratio of means can also be provided with this call of the NLMeans macro. %NLMeans(instore=gammod, coef=coeffs, link=log, options=ratio, title=Ratio of MFG means) ADDITIONAL EXAMPLES AVAILABLE AT support.sas.com: Estimating differences in probabilities with confidence interval http://support.sas.com/kb/37228.html Estimating rate differences (with confidence interval) using a Poisson model http://support.sas.com/kb/37344.html Estimating relative risks in a multinomial response model http://support.sas.com/kb/57798.html Estimating a relative risk (also called risk ratio, prevalence ratio) http://support.sas.com/kb/23003.html Estimating the difference in differences of means http://support.sas.com/kb/61830.html ------------------------------------------------------------------------ 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 Argument List Start =====================; %macro NLMeans( version, /* Any specified text as first argument displays the */ /* current version of the NLMeans macro in the log. */ instore=, /* Specifies the fitted model that was saved using the */ /* STORE statement in the modeling procedure. The OUT= */ /* option in the STORE statement saves the model in a */ /* file known as an item store. This is the preferred */ /* method for providing the required model information. */ /* However, if the modeling procedure does not offer the */ /* STORE statement, then you might be able to use the */ /* inest= and incovb= parameters. */ inest=, /* Specifies the data set of parameter estimates saved */ /* using an ODS OUTPUT statement in the modeling */ /* procedure. The parameter estimates of the model should */ /* be stored in a variable named ESTIMATE in this data */ /* set. An error is issued if the ESTIMATE variable is */ /* not found. */ incovb=, /* Specifies the data set containing the variance- */ /* covariance matrix of model parameters saved using an */ /* ODS OUTPUT statement in the modeling procedure. */ /* Typically, an option such as COVB is required in the */ /* modeling procedure to make this matrix available for */ /* saving. If the saved data set contains numeric */ /* variables other than those containing the covariance */ /* matrix, they should be removed before specifying the */ /* data set in this macro parameter in order to avoid a */ /* compatibility error. The incovb= data set should have */ /* the same number of observations (rows) and variables */ /* (columns) as the number of rows in the inest= data set */ /* in order to be compatible. */ coef=, /* Specifies the data set of estimate coefficients saved */ /* using an ODS OUTPUT statement in the modeling */ /* procedure. The required table must first be produced */ /* by including the E option in the LSMEANS, SLICE, or */ /* ESTIMATE statement in the modeling procedure. The */ /* table name is typically Coef and can be saved in a */ /* data set by specifying: ods output Coef=data-set-name; */ where=, /* Specifies a condition to subset the input data sets. */ link=, /* Specifies the link function used in the modeling */ /* procedure, typically in the LINK= option in the MODEL */ /* statement. Note that for models fit by PROC PROBIT */ /* which specifies the link with the DIST= option, */ /* specify link=probit or cumprobit when DIST=NORMAL; */ /* link=logit or cumlogit when DIST=LOGISTIC; or link=cll */ /* or cumcll when DIST=EXTREME or GOMPERTZ. */ diff=all, /* To estimate and test all pairwise differences of means,*/ /* specify diff=all. Sequential differences of means */ /* (mean1-mean2, mean2-mean3, ...) are provided when */ /* diff=seq. To obtain all differences with a control */ /* level, specify diff=number, where number is the */ /* position of the control level in the ordered list of */ /* levels shown in the results from the modeling */ /* procedure. For example, if level A is the control in a */ /* variable whose levels are shown in the order A, B, C, */ /* then specify diff=1. When there are multiple sets of */ /* estimates from the modeling procedure (presented in */ /* multiple tables), such as when the SLICE / E SLICEBY= */ /* statement is used or when multiple LSMEANS /E, */ /* SLICE /E, or ESTIMATE /E statements are used, then you */ /* can either specify a single type to apply to all */ /* estimate sets, or a list of types to apply different */ /* types to the sets. If diff= is omitted, all pairwise */ /* differences are done for all sets (diff=all). */ contrasts=, /* Specifies an optional data of coefficients defining */ /* contrasts among the estimates in each set of estimates.*/ /* The specified data set must contain a variable named */ /* SET and k variables named K1, K2, ... , Kk, where k is */ /* the number of Row variables in coef= data set. */ /* Optionally, a LABEL variable can be included to label */ /* each of the specified contrasts. For a given */ /* row in the data set, the value of SET indicates */ /* the set of estimates that the contrast in that row */ /* applies to. The SET value must match one of the values */ /* of the LMatrix variable in the coef= data set. The K */ /* variables correspond to the Row variables in the coef= */ /* data set which define each estimate. For multinomial */ /* models, k=mf where m is the number of estimates */ /* requested in each response function and f is the */ /* number of response functions. For link=glogit models, */ /* k=mf+m. The additional set of m variables correspond */ /* to the last response level. See the description of */ /* options=ratio for limitations on the contrast */ /* coefficients when mean ratios are desired. If */ /* contrasts= and diff= are both specified, diff= is */ /* ignored. */ df=, /* Specifies the degrees of freedom to be used in the */ /* tests and confidence intervals computed for the */ /* estimated functions. If omitted, large-sample Wald */ /* statistics are given. The degrees of freedom for */ /* testing a linear combination of parameters in a linear */ /* model would typically be the number of observations */ /* used in fitting the model minus the number of */ /* parameters estimated in the model - essentially, the */ /* error degrees of freedom. */ alpha=0.05, /* Specifies the alpha level to be used in computing */ /* confidence limits. If omitted, alpha=0.05. */ title=, /* Specifies a title for the table(s) of results. The */ /* title must not contain quotes (" or '), ampersands (&),*/ /* commas, or parentheses. If omitted, */ /* title=Nonlinear Function Estimate. */ null=0, /* Specifies c in the null hypothesis f(mu)=c, where f(mu)*/ /* is the function of means being tested. */ covdrop=, /* Specify any variables in the INCOVB= data set that */ /* do not contain columns of the covariance matrix. */ options=NOJOINT NONAMES NOREVERSE NORATIO DIFINFNS NOAPPEND PRINT /* */ /* Use options= to enable or disable any of several */ /* binary options: JOINT|NOJOINT, NAMES|NONAMES, */ /* REVERSE|NOREVERSE, RATIO|NORATIO, DIFINFNS|DIFALL. */ /* JOINT combines multiple sets of estimates into a */ /* single set before differencing or applying contrasts. */ /* NAMES displays the names of the model parameters used */ /* by the NLEST macro. REVERSE reverses the */ /* direction of differencing done by diff=. REVERSE is */ /* ignored when contrasts= is specified. RATIO requests */ /* estimation of ratios rather than differences of pairs */ /* of means. RATIO can be used with diff= or contrasts=, */ /* but with contrasts= each contrast must specify exactly */ /* one 1 coefficient to select the numerator, one -1 */ /* coefficient to select the denominator, and 0 */ /* coefficients otherwise. For multinomial models, DIFALL */ /* applies the differencing type specified in diff= */ /* across all estimates rather than within the response */ /* functions as done by DIFINFNS. This option is ignored */ /* for nonmultinomial models. APPEND merges results from */ /* all estimate sets into a single results data set named */ /* EST_ALL. NOAPPEND saves results sets in separate data */ /* sets (EST1, EST2, ...). The last results set is always */ /* saved in data set EST as well. If not specified, */ /* options=nojoint nonames noreverse noratio difinfns */ /* noappend. NOPRINT suppresses the table of estimates. */ ) / minoperator; ;*=========================== Macro Start =============================; %local notesopt; %let notesopt = %sysfunc(getoption(notes)); %let timenlm = %sysfunc(datetime()); %let _version=1.4; %if &version ne %then %put NOTE: &sysmacroname macro Version &_version..; %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; %if &version ne %then %let dbgopt=v,; %else %let dbgopt=; %end; /* Check for newer version */ %let _notfound=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 %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;; /* 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;