%macro glmpi(version, alpha=0.05, type=all, options=, /* inputs */ data=_last_, response=, trials=, dist=, scale=, power=, pred=, lower=, upper=, /* outputs */ out=GLM_PI, lclq=lclq, uclq=uclq, lclqu=lclqu, uclqu=uclqu ) / minoperator; %let time = %sysfunc(datetime()); %macro existchk(data=, var=, dmsg=e, vmsg=e); %global status; %let status=ok; %if &dmsg=e %then %let dmsg=ERROR; %else %if &dmsg=w %then %let dmsg=WARNING; %else %let dmsg=NOTE; %if &vmsg=e %then %let vmsg=ERROR; %else %if &vmsg=w %then %let vmsg=WARNING; %else %let vmsg=NOTE; %if %quote(&data) ne %then %do; %if %sysfunc(exist(&data)) ne 1 %then %do; %put &dmsg: Data set %upcase(&data) not found.; %let status=nodata; %end; %else %if &var ne %then %do; %let dsid=%sysfunc(open(&data)); %if &dsid %then %do; %let i=1; %do %while (%scan(&var,&i) ne %str() ); %let var&i=%scan(&var,&i); %if %sysfunc(varnum(&dsid,&&var&i))=0 %then %do; %put &vmsg: Variable %upcase(&&var&i) not found in data %upcase(&data).; %let status=novar; %end; %let i=%eval(&i+1); %end; %let rc=%sysfunc(close(&dsid)); %end; %else %put ERROR: Could not open data set &data.; %end; %end; %else %do; %put &dmsg: Data set not specified.; %let status=nodata; %end; %mend; %let _version=2.0; %if &version ne %then %put NOTE: &sysmacroname macro Version &_version; %if &data=_last_ %then %let data=&syslast; %let notesopt = %sysfunc(getoption(notes)); %let newchk=1; %let version=%upcase(&version); %if %index(&version,DEBUG) %then %do; options notes mprint %if %index(&version,DEBUG2) %then mlogic symbolgen; ; ods select all; %put _user_; %end; %else %do; %if %index(&version,NOCHK) %then %let newchk=0; options nonotes nomprint nomlogic nosymbolgen; ods exclude all; %end; /* Check for newer version */ %if &newchk %then %do; %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 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(&version,DEBUG)=0 %then options nonotes;; %end; %let type=%upcase(&type); %let proc=%upcase(%scan(&dist,2)); %if &proc= %then %let proc=null; %let dist=%upcase(%scan(&dist,1)); /* Process OPTIONS= */ %let validopts=COVER; %let cover=0; %let i=1; %do %while (%scan(&options,&i) ne %str() ); %let option&i=%upcase(%scan(&options,&i)); %if &&option&i=COVER %then %let cover=1; %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; /* Errors */ %if &dist= or &pred= or &alpha= or &out= or &lclq= or &uclq= %then %do; %put ERROR: DIST=, PRED=, ALPHA=, OUT=, LCLQ=, and UCLQ= are required.; %goto exit; %end; %if &type ne ALL and &type ne Q and &type ne QU %then %do; %put ERROR: Valid values of TYPE= are ALL, Q, or QU.; %goto exit; %end; %let validprocs=ADAPTIVEREG FMM GAMPL GENMOD GLIMMIX HPFMM HPGENSELECT HPLOGISTIC NLMIXED PROBIT SURVEYLOGISTIC; %if not (&proc in null COUNTREG GLM REG LIFEREG &validprocs) %then %do; %put ERROR: Valid procedures in DIST= are &validprocs..; %goto exit; %end; %let reqvar=&pred; %if &cover and &response ne %then %let reqvar=&reqvar &response; %if &dist=BINOMIAL and &trials ne %then %let reqvar=&reqvar &trials; %if (&type in QU ALL) and &lower ne %then %let reqvar=&reqvar &lower; %if (&type in QU ALL) and &upper ne %then %let reqvar=&reqvar &upper; %existchk(data=&data, var=&reqvar); %if (&status in nodata novar) %then %goto exit; /* Log notes */ %if %index(&version,DEBUG)=0 %then options notes;; %if &cover and &response= %then %do; %put NOTE: RESPONSE= is required with OPTIONS=COVER.; %let cover=0; %end; %if (&type in QU ALL) and (&lower= or &upper=) %then %do; %put NOTE: LOWER= and UPPER= are required for TYPE=QU.; %put NOTE- Intervals with mean uncertainty (TYPE=QU) will not be computed.; %if &type=ALL %then %let type=Q; %if &type=QU %then %goto exit; %end; %if %index(&version,DEBUG)=0 %then options nonotes;; /* Distribution parameters */ %let reqscale=0; %if &dist=GAMMA %then %do; %if (&proc in GLIMMIX ADAPTIVEREG) %then %do; %let parmsq=1/&scale, &pred*&scale; %let parmsquu=1/&scale, &upper*&scale; %let parmsqul=1/&scale, &lower*&scale; %end; %else %do; %let parmsq=&scale, &pred/&scale; %let parmsquu=&scale, &upper/&scale; %let parmsqul=&scale, &lower/&scale; %end; %let qdist=GAMMA; %let reqscale=1; %end; %else %if &dist=POISSON %then %do; %let parmsq=&pred; %let parmsquu=&upper; %let parmsqul=&lower; %let qdist=POISSON; %end; %else %if &dist=NEGBIN %then %do; %let parmsq=1/(1+&pred*&scale), 1/&scale; %let parmsquu=1/(1+&upper*&scale), 1/&scale; %let parmsqul=1/(1+&lower*&scale), 1/&scale; %let qdist=NEGBIN; %let reqscale=1; %end; %else %if &dist=NORMAL %then %do; %if (&proc in GLIMMIX FMM HPGENSELECT GAMPL ADAPTIVEREG HPFMM NLMIXED) %then %let scale=%quote(sqrt(&scale)); %let parmsq=&pred, &scale; %let parmsquu=&upper, &scale; %let parmsqul=&lower, &scale; %let qdist=NORMAL; %let reqscale=1; %end; %else %if &dist=BINOMIAL %then %do; %if &trials= %then %do; %put ERROR: TRIALS= is required for DIST=&dist..; %goto exit; %end; %let parmsq=&pred, &trials; %let parmsquu=&upper, &trials; %let parmsqul=&lower, &trials; %let qdist=BINOMIAL; %end; %else %if &dist=IGAUSS %then %do; %if (&proc in GLIMMIX FMM HPGENSELECT GAMPL ADAPTIVEREG HPFMM) %then %let scale=%quote(sqrt(&scale)); %let parmsq=1/&scale**2, &pred; %let parmsquu=1/&scale**2, &upper; %let parmsqul=1/&scale**2, &lower; %let qdist=IGAUSS; %let reqscale=1; %end; %else %if &dist=GENPOISSON %then %do; %let parmsq=&pred*exp(-&scale), 1-exp(-&scale); %let parmsquu=&upper*exp(-&scale), 1-exp(-&scale); %let parmsqul=&lower*exp(-&scale), 1-exp(-&scale); %let qdist=GENPOISSON; %let reqscale=1; %end; %else %if &dist=TWEEDIE %then %do; %if &power= %then %do; %put ERROR: POWER= is required for DIST=&dist..; %goto exit; %end; %let parmsq=&power, &pred, &scale; %let parmsquu=&power, &upper, &scale; %let parmsqul=&power, &lower, &scale; %let qdist=TWEEDIE; %let reqscale=1; %end; %else %if &dist=WEIBULL %then %do; %if (&proc in LIFEREG) %then %let parmsq=1/&scale, exp(&pred); %else %do; %let parmsq=1/&scale, &pred/gamma(1+&scale); %let parmsquu=1/&scale, &upper/gamma(1+&scale); %let parmsqul=1/&scale, &lower/gamma(1+&scale); %end; %let qdist=WEIBULL; %let reqscale=1; %end; %else %if &dist=BETA %then %do; %let parmsq=&pred*&scale, (1-&pred)*&scale; %let parmsquu=&upper*&scale, (1-&upper)*&scale; %let parmsqul=&lower*&scale, (1-&lower)*&scale; %let qdist=BETA; %let reqscale=1; %end; %else %if &dist=LOGNORMAL %then %do; %if (&proc in FMM HPFMM) %then %let pred=log(&pred)-&scale/2; %if (&proc in LIFEREG) %then %let scale=%quote(&scale**2); %let parmsq=&pred, sqrt(&scale); %let parmsquu=&upper, sqrt(&scale); %let parmsqul=&lower, sqrt(&scale); %let qdist=LOGNORMAL; %let reqscale=1; %end; %else %do; %put ERROR: DIST= must be one of POISSON, NEGBIN, GAMMA, NORMAL, BINOMIAL,; %put ERROR- IGAUSS, GENPOISSON, TWEEDIE, WEIBULL, BETA, LOGNORMAL.; %goto exit; %end; %if &reqscale and &scale= %then %do; %put ERROR: SCALE= is required for DIST=&dist..; %goto exit; %end; /* Compute prediction intervals: TYPE=Q: quantile-based TYPE=QU: quantile-based including mean uncertainty */ data &out; set &data; /* Quantile-based interval */ %if &type=Q or &type=ALL %then %do; &lclq = quantile("&qdist",&alpha/2,&parmsq); &uclq = quantile("&qdist",1-&alpha/2,&parmsq); %if &cover %then %do; _CoverQ=(&lclq<=&response<=&uclq); label _CoverQ="Quantile-based"; %end; %end; /* Quantile-based interval including mean uncertainty */ %if &type=QU or &type=ALL %then %do; &lclqu = quantile("&qdist",&alpha/2,&parmsqul); &uclqu = quantile("&qdist",1-&alpha/2,&parmsquu); %if &cover %then %do; _CoverQU=(&lclqu<=&response<=&uclqu); label _CoverQU="Quantile+Uncertainty*"; %end; %end; run; %if %sysfunc(exist(&out)) %then %do; %if %index(&version,DEBUG)=0 %then options notes;; %let dsid=%sysfunc(open(&out)); %if &dsid %then %do; %let nobs=%sysfunc(attrn(&dsid,NOBS)); %let nvar=%sysfunc(attrn(&dsid,NVARS)); %let rc=%sysfunc(close(&dsid)); %put NOTE: The data set &out has &nobs observations and &nvar variables.; %end; %else %do; %put WARNING: Could not check size of OUT= data set.; %end; %if %index(&version,DEBUG)=0 %then options nonotes;; %end; %if &cover %then %do; %let nominal=%sysevalf(1-&alpha); ods select all; proc means data=&out n mean; var _cov:; output out=Cover mean=CovQ CovQU; title "Interval coverage"; title2 "Nominal: &nominal"; %if &type=QU or &type=ALL %then %do; footnote "*Quantile+Uncertainty intervals can exceed the nominal coverage level."; %end; run; data Cover; set Cover; Alpha=α N=_freq_; drop _:; run; %end; /* Clean up */ %exit: %if %index(&version,DEBUG) %then %do; options nomprint nomlogic nosymbolgen; %put _user_; %end; ods select all; options ¬esopt; title; footnote; %let time = %sysfunc(round(%sysevalf(%sysfunc(datetime()) - &time), 0.01)); %put NOTE: The &sysmacroname macro used &time seconds.; %mend;