/** Latest version */ /** Update: This program doesn't work for byvar yet, delete this parameter */ /** coped from Ellen /udd/stleh/ehmac/relrss8.sas, March 8, 2004, modified by Ruifeng and saved in /udd/strui/relrisk **/ /** This program only works with repeated measurment with empirical (robust) variance**/ %macro relrisk8(data=, depend=, independ=, notes=nonotes, id=, empcal=F, where=, extravar=, reptype=ind, class=, initmeth=poisson, fintercept=, finitval=, withinvar=, noint=F, cov=F, corr=F, contrast=, estimate=, titlenum=4, outdat=, gmprint=F, poisprint=F, showsump=F); %let _fdl = %sysfunc(getoption(formdlim)); %let _nt = %sysfunc(getoption(notes)); options ¬es validvarname=upcase replace formdlim='='; options nocenter ps=78 ls=80; /*********************** parameters explanation ********************* Macro relrisk obtains relative risk estimates using proc GENMOD with (usually) initial values taken from some other source. The initial values are needed because the method for obtaining relative risk - running GENMOD with the log link and binomial distribution - can have problems converging. By default, initializations for GENMOD are obtained by running PROC GENMOD with the log link and the Poisson distribution. Used with the empirical variance, this method produces estimates fairly close to those that would be produced by PROC GENMOD with log binomial. This program was originally written using PROC NLIN, which is still available. For NLIN, the initialization of the intercept is adjusted if any predicted values from NLIN are greater than 1. You can also try to get the initial values using GENMOD run as a logistic or you can force the initial values. If you are adding one variable to an analysis, you can use the estimates for the original variables to provide initial values for the next analysis. This works if the added variable does not affect the others too much. If you only want to run RELRISK8 to get the relative risk printout, you can set the method to NONE. Error messages are printed if parameters or data do not conform to requirements. All error messages begin with "ERROR". If any error occurs, the program terminates and no output is produced. ; I have added several things that you can take or leave. (Sally Skinner, sally.skinner@tufts.edu). 7/21/03. Along with PARAM, which contains the estimate and p-value information, the ODS output files for the model fitting tests (geefit), the model infor (mod) and the convergence statust (conv) are available for additional ODS work. There are also two global macro variables (&nlinstatus and &nlinreason) which give convergence status information for any analysis used to produce initial values. These can be used to shut down an automatic process if NLIN does not converge. ; %* Parameters: data: The dataset. depend: The dependent variable. independ: The independent variable. notes: (default=nonotes) Set equal to NOTES if you want the SAS notes to print. id: The ID number for each person used for SUBJECT= in the REPEATED statement. empcal: (default=F) Set empcal to T to get the empirical (robust) variance. NOTE: for now, it can only be T byvar: (default=mergvbl_) If you are using a b-variable to subset the analysis, put it here. The default, MERGVBL_, is set to the same value for everyone. where: Add a WHERE statement to the data step, i.e. where=%str(zb_femal=0 & zb_hart1=1). Anytime there is punctuation or operators in the value of a parameter, you must put the macro function %str() around it. You may get extremely strange results if you do not. extravar: Additional variables used in the process. initmeth: (default=Poisson) The method for finding the initial values for each parameter. INITMETH has five choices: POISSON - which uses PROC GENMOD with the POISSON distribution. This is the default NLIN - which uses PROC NLIN to determine the initial values. PREV uses the previous estimates to initialize genmod. It can only be used if the previous run of RELRISK has converged and one or two variables have been replaced or added. FORCE allows the user to force either the intercept, the initial estimates, or both, using FINTERCEP and FINITVAL. Here is an example of the call: INITMETH=force, FINTERCEPT=-5, FINITVAL=-1 -1 2 This sets the initial value of the intercept to -5 and the initial estimates for the first three variables to -1, -1 and 2, respectively. NONE - which does not produce initial values for the log binomial analysis. This is useful if you just want to have the relative risks calculated automatically, or if a particular model runs fine without any initial values. fintercept: Initial value of the intercept when INITMETH = FORCE. It can be left blank for no-intercept models finitval: Initial value of the betas for the independent variables. reptype: (default=ind) Type of covariance matrix used (IND, EXCH, UN, TOEP (not TOEP(2)!! ). withinvar: The index variable for within-group records. This must have a value here if REPTYPE=TOEP or UN. noint: If NOINT=T, a no-intercept model is run. This cannot be run with INITMETH=NLIN. cov: If COV=T then the covariance matrix is printed out and sent to ODS. corr: If CORR=T, then the correlation matrix is printed out and sent to ODS. contrast: This is where you put the contrast statements. You MUST use %str(), i.e contrast=%str(contrast zb_Femal 1...). estimate: This is where you put the estimate statments. They also must have %str() around them. titlenum: (default=4) More titles have been added to make it easier to navigate through the output. In order to avoid having these titles knock out your titles, use TITLENUM to state what line the MACRO titles should be in. The default is line 4. If you are using title1 and title2, TITLENUM should be set to 3. ; **************************** end of parameters explanation *****************/ /* upcasing parameters */ %let cov=%upcase(&cov); %let corr=%upcase(&corr); %let empcal=%upcase(&empcal); %let initmeth=%upcase(&initmeth); %let reptype=%upcase(&reptype); %let gmprint=%upcase(&gmprint); %let poisprint=%upcase(&poisprint); /* changed macro so that poisson always runs with empirical estimates, even if empcal is F */ /* clearing titles from previous runs */ title&titlenum ' '; /* if one obs per subject, we want the results of the proc genmod without the repeated statement */ %if "&reptype" eq "IND" %then %let empcal=F; /* localizing new macro vbls */ %local callerr typerr predmax predmin interdif letters numivar modnlin i dimen indep nlinstatus nlinreason initgen initintr noconv _maxpred_ _pct_ _nsubj daterr diffpred _sumpred_ _nevents_ poisconvst ; %let nlinstatus = 0; * Initialize macro variables to avoid carryover when there NLIN does not run properly (SCS 7/15/03); %let typerr=0; %let callerr=0; %let daterr=0; %let predmax=; %let predmin=; %let interdif=; %let letters=; %let numivar=; %let modnlin=; options ¬es; /* check that the required variables are input */ /* check for all errors before aborting--eh 5/12/04 */ %if %length(&data)=0 %then %do; %let callerr=1; data _null_; put "ERROR in macro call: The name of the input dataset was omitted."; put 'The RELRISK8 macro stopped.'; file print; put "ERROR in macro call: The name of the input dataset was omitted."; put 'The RELRISK8 macro stopped.'; %end; %if %length(&depend)=0 %then %do; %let callerr=1; data _null_; put "ERROR in macro call: The name of the dependent variable was omitted."; put 'The RELRISK8 macro stopped.'; file print; put "ERROR in macro call: The name of the dependent variable was omitted."; put 'The RELRISK8 macro stopped.'; %end; %if %length(&independ)=0 %then %do; %let callerr=1; data _null_; put "ERROR in macro call: "; put " The name(s) of the independent variable(s) was(were) omitted."; put 'The RELRISK8 macro stopped.'; file print; put "ERROR in macro call: "; put " The name(s) of the independent variable(s) was(were) omitted."; put 'The RELRISK8 macro stopped.'; run; %end; %if &initmeth ^= NLIN & &initmeth ^= PREV & &initmeth ^= NONE & &initmeth ^= FORCE & &initmeth ^=LOGISTIC & &initmeth ^=POISSON %then %do; %let callerr=1; data _null_; put "ERROR in macro call: "; put " INITMETH must be either NONE, NLIN, PREV, FORCE, LOGISTIC or POISSON."; put " You have INITMETH=&initmeth ."; file print; put "ERROR in macro call: "; put " INITMETH must be either NONE, NLIN, PREV, FORCE, LOGISTIC or POISSON."; put " You have INITMETH=&initmeth ."; %end; %if &empcal eq T and "&id" eq "" %then %do; %let callerr=1; data _null_; put "ERROR in macro call: You asked for a GEE model (empcal=T),"; put " but did not specify a variable to serve as the subject (ID)."; file print; put "ERROR in macro call: You asked for a GEE model (empcal=T),"; put " but did not specify a variable to serve as the subject (ID)."; run; %end; %if &initmeth = FORCE & %length(&fintercept)=0 & %length(&finitval)=0 & &noint^=T %then %do; %let callerr=1; data _null_; put "ERROR in macro call: INITMETH=FORCE "; put " but there are no values for FINTERCEPT or FINITVAL"; put 'The RELRISK8 macro stopped.'; file print; put "ERROR in macro call: INITMETH=FORCE "; put " but there are no values for FINTERCEPT or FINITVAL"; put 'The RELRISK8 macro stopped.'; %end; %if &initmeth = NLIN & &noint=T %then %do; %let callerr=1; data _null_; put "ERROR in macro call: PROC NLIN CANNOT BE RUN WITH THE NO-INTERCEPT MODEL"; put " PLEASE USE ANOTHER METHOD OR SHUT OFF THE PARAMETER NOINT"; put 'The RELRISK8 macro stopped.'; file print; put "ERROR in macro call: PROC NLIN CANNOT BE RUN WITH THE NO-INTERCEPT MODEL"; put " PLEASE USE ANOTHER METHOD OR SHUT OFF THE PARAMETER NOINT"; put 'The RELRISK8 macro stopped.'; %end; %if &initmeth = FORCE & %length(&fintercept)>0 & &noint=T %then %do; %let callerr=1; data _null_; put "ERROR in macro call: Intercept value is forced for a no-intercept analysis"; put 'The RELRISK8 macro stopped.'; file print; put "ERROR in macro call: Intercept value is forced for a no-intercept analysis"; put 'The RELRISK8 macro stopped.'; %end; %if "&reptype" ne "EXCH" and "&reptype" ne "IND" and "&reptype" ne "CS" and "&reptype" ne "SIMPLE" and %length(&withinvar)=0 %then %do; %let callerr=1; data _null_; put "ERROR in macro call: REPTYPE=&reptype, but WITHINVAR has no value."; put " For unstructured covariance, there must be a within-id counter"; put 'The RELRISK8 macro stopped.'; file print; put "ERROR in macro call: REPTYPE=&reptype, but WITHINVAR has no value."; put " For unstructured covariance, there must be a within-id counter"; put 'The RELRISK8 macro stopped.'; %end; /* check that data set &data exists */ %str( data _null_; if 0 then set &data; stop; run; ); %if &syserr ne 0 %then %do; %let daterr=1; data _null_; put "ERROR in macro call: The dataset you specified does not exist."; put 'The RELRISK8 macro stopped.'; file print; put "ERROR in macro call: The dataset you specified does not exist."; put 'The RELRISK8 macro stopped.'; %end; %if &daterr eq 1 %then %goto errend; /* determine the number of independent variables - "numivar" */ %let i=1; %let numivar=%numargs(&independ); /* RL: what is this for ? indep becomes ^A^S^B&independent^A^T^B */ %let indep=%str(%()&independ%str(%)); %let dimen=%eval(&numivar+1); data _dat_998a ; set &data; keep &depend &independ &extravar &id &withinvar ; %if %length(&where) > 0 %then %do; where &where; %end; run; %if &syserr ne 0 %then %do; %let daterr=1; data _null_; put "ERROR in macro run: "; put " An error occurred during the creation of the working dataset,"; put " while keeping DEPEND and INDEPEND from &data"; put " The names of your variables may be incorrect."; put " dependent vbl is &depend . "; put " predictors are &independ ."; %if "&empcal" eq "T" %then %do; put " subject is &id ."; if "&reptype" ^in ("IND", "EXCH", "SIMPLE", "CS") then do; put " WITHINVAR is &withinvar ."; end; %end; put 'The RELRISK8 macro stopped.'; file print; put "ERROR in macro run: "; put " An error occurred during the creation of the working dataset,"; put " while keeping DEPEND and INDEPEND from &data"; put " The names of your variables may be incorrect."; put " dependent vbl is &depend . "; put " predictors are &independ ."; %if "&empcal" eq "T" %then %do; put " subject is &id ."; if "&reptype" ^in ("IND", "EXCH", "SIMPLE", "CS") then do; put " WITHINVAR is &withinvar ."; end; %end; put 'The RELRISK8 macro stopped.'; %end; %if &daterr eq 1 %then %goto errend; data _dat_998; set _dat_998a; /* keep only observations with known values for the dependent and independent variables */ if nmiss(of &independ &depend) gt 0 then delete; %if %length(&id) gt 0 and "&empcal" eq "T" %then %do; if nmiss(&id ) eq 0; %end; %if %length(&withinvar) gt 0 and "&empcal" eq "T" %then %do; if nmiss(&withinvar) eq 0 ; %end; mergvbl_=1; _id_=-_n_; run; %if %length(&id) eq 0 %then %let id=_id_; %nobsnvar(data=_dat_998); %if &nobs eq 0 %then %do; %let daterr=1; data _null_; put 'ERROR in macro run: '; put ' The working data set has no observations with non-missing data.'; put " The names of your variables may be incorrect."; %if "&empcal" eq "T" %then %do; if "&reptype" in ("EXCH", "IND", "CS", "SIMPLE") then do; put " The subject variable is &id . "; end; else do; put " subject vbl is &id . vbl determining order of observations is &withinvar ."; end; %end; put " dependent variable is &depend ."; put " predictors are &independ ."; put 'The RELRISK8 macro stopped.'; file print; put 'ERROR in macro run: '; put ' The working data set has no observations with non-missing data.'; put " The names of your variables may be incorrect."; %if "&empcal" eq "T" %then %do; if "&reptype" in ("EXCH", "IND", "CS", "SIMPLE") then do; put " The subject variable is &id . "; end; else do; put " subject vbl is &id . vbl determining order of observations is &withinvar ."; end; %end; put " dependent variable is &depend ."; put " predictors are &independ ."; put 'The RELRISK8 macro stopped.'; run; proc print data=_dat_998a (obs=5); var &id &depend &independ &extravar ; run; %end; %if &daterr eq 1 %then %goto errend; proc contents noprint data= _dat_998 out= _dat_999 ; run; /* check that the variables are all numeric */ data _null_; set _dat_999 end=_end_; typerr=0; if name ne %upcase("&id") and type ne 1 then do; typerr=1; put 'ERROR in macro run: The type is not numeric for variable ' name '.'; put 'The RELRISK8 macro stopped.'; file print; put 'ERROR in macro run: The type is not numeric for variable ' name '.'; put 'The RELRISK8 macro stopped.'; end; call symput('typerr', typerr); run; %if &callerr eq 1 or &typerr=1 %then %goto errend; /**************************************/ /* initializing by NLIN */ /**************************************/ %if &initmeth=NLIN %then %do; %let initnlin = intercept = .1 ; %let letters = intercept; %do i=1 %to &numivar.; %let initnlin= &initnlin a&i =.1 ; %let letters = &letters a&i ; %end; %let modnlin=exp(intercept ; %do i=1 %to &numivar.; %let modnlin= &modnlin + a&i*%scan(&independ,&i) ; %end ; %let modnlin=&modnlin); /*gets nlin estimates*/ ods listing close; proc nlin data=&data outest=parms; model &depend=&modnlin.; parameters &initnlin.; output out=predict p=predict; ods output convergencestatus=nlconv; title&titlenum 'NLIN TO FIND INITIAL VALUES'; run; ods output close; data _null_; set nlconv; call symput ("nlinstatus",put (status,1.)); call symput ("nlinreason",put(reason,$100.)); run; %if &nlinstatus ne 0 %then %goto cucun; proc univariate data=predict noprint; var predict; output out=predstat max=predmax min=predmin; run; data _null_; set predstat; call symput('predmax',predmax); call symput('predmin',predmin); run; %if &predmax>1 %then %do; /* create macro variable "interdif", the amount by which the intercept estimate from NLIN must be adjusted in order to force all fitted values to be less than 1. */ data _null_; interdif=log(1/(&predmax+.0001)); call symput('interdif',interdif); run; %end; %else %let interdif=0; /* create macro variables "initgen1" and "initgen2" and "initintr" */ data _null_; set parms; file print; array vars[&dimen.] &letters.; length initgen $ %eval((&dimen+1)*30); if _type_='FINAL' then do; /* the adjustment factor, exp(&interdif)=1/(&predmax+.0001), insures that predicted values from NLIN are less than 1 */ initintr = "intercept= " ||left(trim(put(intercept+&interdif,14.9))); initgen = "initial= "; do ii=1 to &numivar; initgen=trim(initgen)||' '||left(trim(put(vars[ii+1],14.9))); end; call symput("initintr",initintr); call symput("initgen",initgen); end; run; cucun: %if &nlinstatus ne 0 %then %do; %let initmeth = POISSON ; ods listing; data _null_; put 'ERROR in macro run: Proc NLIN did not converge.'; put ' Try another initialization method.'; file print; put 'ERROR in macro run: Proc NLIN did not converge.'; put ' Try another initialization method.'; run; %end; %end; /*End of getting initial values using PROC NLIN */ /***************************************/ /* initmeth=FORCE */ /***************************************/ %if &initmeth=FORCE %then %do; ods listing close; data _null_; length fint $15 fval $ %eval((&dimen+1)*30); fint = "&fintercept"; fval = "&finitval"; if fint>'' then fint="intercept= " || fint; call symput("initintr",fint); if fval>'' then fval = "initial= " || fval; call symput("initgen",fval); run; /* %put "Initialization Method: User values" ; %put "User Intercept: " &initintr; %put "User initial values: " &initgen; */ title&titlenum 'User values (FORCE) used for initialization'; %end; /* End of process for forced initial values */ /******************************************/ /* initmeth=PREV */ /******************************************/ %if &initmeth=PREV %then %do; ods listing close; data _null_; set prevval end=eof; length intcept $15 initlist $ %eval((&dimen+1)*30); retain intcept '' initlist 'initial='; if upcase(parameter)="INTERCEPT" then intcept = 'intercept= ' ||left(trim(put (estimate, 14.9))); else if upcase(parameter)^="SCALE" then initlist=trim(initlist)||' '||left(trim(put(estimate,14.9))); call symput("initintr",intcept); call symput("initgen", initlist); run; title&titlenum 'Initialization from results of a previous run (PREV)'; %end; /* End of process for using previous betas */ /***************************************/ /* initmeth=LOGISTIC */ /***************************************/ %if &initmeth = LOGISTIC %then %do; ods listing close; proc genmod data=_dat_998 descending; class &id &withinvar &class ; model &depend=&independ /dist=bin link=logit %if &noint=T %then noint; ; ods output ConvergenceStatus=logitconv; %if &empcal eq T %then %do; ods output GEEEmpPEST=logitpar; repeated subject=&id / type=&reptype %if %length(&withinvar) > 0 %then %do; withinsubject=&withinvar %end; ; %end; %else %do; ods output ParameterEstimates=logitpar ; %end; title&titlenum 'Initial values from PROC GENMOD with logit link (initmeth=LOGISTIC)'; run; ods output close; data logitinit; length intcept $15 initlist $ %eval((&dimen+1)*30); set logitconv (in=c) logitpar (in=p) end=eof; retain stat 0 intcept '' initlist 'initial='; if c & status>stat then status=stat; if p then do; %if &noint=T %then %do; if upcase(parameter)="INTERCEPT" then intcept='';%end; %else %do; if upcase(parameter)="INTERCEPT" then intcept = 'intercept= ' ||left(trim(put (estimate,14.9))); else %end; if upcase(parameter)^="SCALE" then initlist=trim(initlist)||' '||left(trim(put(estimate,14.9))); end; call symput("initintr",intcept); call symput("initgen",initlist); %end; /* End of process using betas from a logistic regression*/ /******************************************/ /* initmeth=POISSON */ /******************************************/ %if &initmeth = POISSON %then %do; proc genmod data=_dat_998 ; ods listing close; class &id &withinvar &class ; model &depend=&independ / dist=poisson link=log %if &noint=T %then noint; ; ods output ConvergenceStatus=poisconv; ods output GEEEmpPEST=poisspar; %if "&cov" eq "T" %then %do; ods output geercov=poiscovemp; %end; %if "&corr" eq "T" %then %do; ods output geercorr=poiscorremp; %end; repeated subject=&id / type=&reptype %if %length(&withinvar) > 0 %then %do; withinsubject=&withinvar %end; %if "&cov" eq "T" %then %do; ecovb %end; %if "&corr" eq "T" %then %do; ecorrb %end; ; run; ods output close; ods listing; options notes; data poisconv; set poisconv; call symput('poisconvst', status); if status ne 0 then do; put 'ERROR in macro run: Proc genmod with Poisson variance'; put ' and log link and empirical variance did not converge'; file print; put 'ERROR in macro run: Proc genmod with Poisson variance'; put ' and log link and empirical variance did not converge'; end; run; ods listing close; %if &poisconvst ne 0 %then %let initmeth=NONE; %if &poisconvst ne 0 %then %goto noinit; %if &poisprint eq T %then %do; %pfrq; data ppar; set poisspar; length varname $20 ; if parameter ne ' ' then varname=parameter; else if parm ne ' ' then varname=parm; pval=max(probz, probchisq); format pval pvalue6.4; units=1; if varname ne 'Scale' and StdErr ne 0 then do; rrest=exp(units*estimate); rrcilo95=exp(units*estimate-1.96*units*stderr); rrcihi95=exp(units*estimate+1.96*units*stderr); end; label rrest = 'relative risk' rrcilo95 = '95% CI Low' rrcihi95 = '95% CI high' varname = 'Parameter' pval = 'P value' level1 = ' ' ; if level1 eq . then level1=.; keep level1 rrest rrcilo95 rrcihi95 varname pval estimate units stderr ; run; ods listing; title&titlenum "Results of running RELRISK8 on data &data"; title%eval(&titlenum + 1) 'Results of Poisson regression with empirical variance'; title%eval(&titlenum + 2) "Outcome is &depend . &_nevents_ events (&_pct_ %) in &nobs trials"; proc print noobs label data=ppar; %if %length(&class) eq 0 %then %do; var varname Estimate Stderr pval units rrest rrcilo95 rrcihi95; %end; %else %do; var varname level1 estimate stderr pval units rrest rrcilo95 rrcihi95; %end; run; title&titlenum ' '; ods listing close; %end; data poissinit; length intcept $15 initlist $ %eval((&dimen+1)*30); set poisconv (in=c) poisspar (in=p) end=eof; retain stat 0 intcept '' initlist 'initial='; if c & status>stat then status=stat; if p then do; %if &noint=T %then %do; if upcase(parameter)="INTERCEPT" then intcept='';%end; %else %do; if upcase(parameter)="INTERCEPT" then intcept = 'intercept= ' ||left(trim(put (estimate,14.9))); else %end; if upcase(parameter)^="SCALE" then initlist=trim(initlist)||' '||left(trim(put(estimate,14.9))); end; call symput("initintr",intcept); call symput("initgen",initlist); %end; /* End of process using betas from a poisson regression*/ %noinit: /***************************** initmeth=NONE *****************************/ %if &initmeth = NONE %then %do; %let initintr= ; %let initgen = ; %end; /* End of no method for initial values */ /****************************************/ /* finally to the main genmod */ /****************************************/ %if &gmprint eq T %then %do; ods listing; %end; %else %do; ods listing close; %end; proc genmod data=_dat_998 descending; class &id &withinvar &class ; model &depend=&independ / dist=bin link=log &initintr &initgen %if &cov = T %then covb ; %if &corr = T %then corrb ; %if &noint =T %then noint ;; output out= outstats pred=predmean u=upper l=lower hesswgt=hessian reschi=chiresid resdev=devresid reslik=lkhoodresid stdxbeta=se_xbeta stdreschi=standchires stdresdev=standdevres xbeta = x_beta ; %if &empcal eq T %then %do; ods output GEEEmpPEST=param ; ods output GEEModInfo=mod; %if "&corr" eq "T" %then %do; ods output geencorr=corrbinmb; %end; %if "&cov" eq "T" %then %do; ods output geencov=covbinmb; %end; %end; %else %do; ods output ParameterEstimates=param; ods output modelInfo=mod; %if "&corr" eq "T" %then %do; ods output covb=covbin; %end; %if "&cov" eq "T" %then %do; ods output corrb=corrbin; %end; %end; ods output ModelFit (persist=run) = geefit; ods output ConvergenceStatus=conv; ods output responseprofile = respprof; %if %length(&estimate) >0 %then %do; ods output estimates=ests; %end; %if %length(&contrast)> 0 %then %do; ods output contrasts=contr; %end; %if &empcal eq T %then %do; repeated subject=&id / type=&reptype %if %length(&withinvar) > 0 %then %do; withinsubject=&withinvar %end; %if "&cov" eq "T" %then %do; mcovb %end; %if "&corr" eq "T" %then %do; mcorr %end; ; %end; %if %length(&estimate) > 0 %then %do; &estimate; %end; %if %length(&contrast) > 0 %then %do; &contrast;%end; run; ods listing; data _null_; set conv; if Status ne 0 then nonconv=1; else nonconv=0; call symput('noconv', nonconv); run; %if &noconv ne 0 %then %do; title&titlenum ' '; ods listing ; data _null_; put "WARNING in macro run: The log-binomial model did not converge"; put ' The macro will show the results of the poisson model with'; put ' empirical variance below.'; put "The initial values were"; put " &initintr"; put " &initgen"; file print; put "WARNING in macro run: The log-binomial model did not converge"; put ' The macro will show the results of the poisson model with'; put ' empirical variance below.'; put "The initial values were"; put " &initintr"; put " &initgen"; run; proc print data=conv; run; options obs=max; %pfrq; /***** use the results of the poisson regression, if already done *****/ %if &poisconvst eq 0 %then %do; %if "&poisprint" ne "T" %then %do; data ppar; set poisspar; length varname $20 ; if parameter ne ' ' then varname=parameter; else if parm ne ' ' then varname=parm; pval=max(probz, probchisq); format pval pvalue6.4; units=1; if varname ne 'Scale' and StdErr ne 0 then do; rrest=exp(units*estimate); rrcilo95=exp(units*estimate-1.96*units*stderr); rrcihi95=exp(units*estimate+1.96*units*stderr); end; label rrest = 'relative risk' rrcilo95 = '95% CI Low' rrcihi95 = '95% CI high' varname = 'Parameter' pval = 'P value' level1 = ' ' ; if level1 eq . then level1=.; keep level1 rrest rrcilo95 rrcihi95 varname pval estimate units stderr ; run; ods listing; title&titlenum "Results of running RELRISK8 on data &data"; title%eval(&titlenum + 1) 'Results of Poisson regression with empirical variance'; title%eval(&titlenum + 2) "Outcome is &depend . &_nevents_ (&_pct_ %) in &nobs observations."; proc print noobs label data=ppar; %if %length(&class) eq 0 %then %do; var varname Estimate Stderr pval units rrest rrcilo95 rrcihi95; %end; %else %do; var varname level1 estimate stderr pval units rrest rrcilo95 rrcihi95; %end; run; %end; %end; %else %if %length(&poisconvst) eq 0 %then %do; ods listing close; %pfrq; proc genmod data=_dat_998; class &id &withinvar &class; model &depend = &independ / dist=poisson link=log %if &noint eq T %then noint; ; ods output convergencestatus=poisconv geeemppest=poisspar %if "&cov" eq "T" %then geercov=poiscovemp; %if "&corr" eq "T" %then geercorr=poiscorremp; ; repeated subject=&id / type=&reptype %if %length(&withinvar) gt 0 %then %do; withinsubject=&withinvar %end; %if "&cov" eq "T" %then %do; ecovb %end; %if "&corr" eq "T" %then %do; ecorrb %end; ; run; ods output close; ods listing; options notes; data poisconv; set poisconv; call symput ('poisconvst', status); if status ne 0 then do; put 'ERROR in macro run: Poisson regression '; put ' with empirical variance did not converge'; put 'ERROR in macro run: Poisson regression '; put ' with empirical variance did not converge'; file print; end; run; %if &poisconvst ne 0 %then %goto errend; options ¬es; data ppar; set poisspar; length varname $20 ; if parameter ne ' ' then varname=parameter; else if parm ne ' ' then varname=parm; pval=max(probz, probchisq); format pval pvalue6.4; units=1; if varname ne 'Scale' and StdErr ne 0 then do; rrest=exp(units*estimate); rrcilo95=exp(units*estimate-1.96*units*stderr); rrcihi95=exp(units*estimate+1.96*units*stderr); end; label rrest = 'relative risk' rrcilo95 = '95% CI Low' rrcihi95 = '95% CI high' varname = 'Parameter' pval = 'P value' level1 = ' ' ; if level1 eq . then level1=.; keep level1 rrest rrcilo95 rrcihi95 varname pval estimate units stderr ; run; ods listing; options notes; title&titlenum "Results of running RELRISK8 on data &data ."; title%eval(&titlenum + 1) 'Results of Poisson regression with empirical variance'; title%eval(&titlenum + 2) "Outcome is &depend . &_nevents_ events (&_pct_ %) in &nobs trials."; proc print noobs label data=ppar; %if %length(&class) eq 0 %then %do; var varname Estimate Stderr pval units rrest rrcilo95 rrcihi95; %end; %else %do; var varname level1 estimate stderr pval units rrest rrcilo95 rrcihi95; %end; run; %end; %if &poisconvst ne 0 %then %goto errend; %end; %if &noconv eq 0 %then %do; proc means noprint data=outstats; var predmean; output out=_maxpred_ max=maxpred sum=sumpred; run; data _maxpred_; set _maxpred_; call symput ('_maxpred_', maxpred); call symput ('_sumpred_', trim(left(round(sumpred)))); if maxpred gt 1 then do; put 'The maximum predicted probability was ' (maxpred) (5.2) ' .'; put ' This is NOT a legal value.'; put ' You need to rerun RELRISK8 with different initialization.'; put ' One possibility is to use the output of this procedure for all parameters'; put ' except the intercept. Then Force the intercept to be much lower.'; put ' Another option is to use the results of the Poisson regression with the empirical variance.'; file print; put 'The maximum predicted probability was ' (maxpred) (5.2) ' .'; put ' This is NOT a legal value.'; put ' You need to rerun RELRISK8 with different initialization.'; put ' One possibility is to use the output of this procedure for all parameters'; put ' except the intercept. Then Force the intercept to be much lower.'; put ' Another option is to use the results of the Poisson regression with the empirical variance.'; end; run; /* even though get an illegal value, will print out results, so can use for force in next run */ data respprof; set respprof; if orderedvalue eq 1 then nevents=count; run; proc means noprint data=respprof; var nevents; output out=_sumev_ mean=nevents; run; data _sumev_; set _sumev_; nobs=&nobs; pct=round(100*nevents/nobs, .1); call symput('_pct_', trim(left(pct))); call symput('_nevents_', trim(left(nevents))); sumpred=&_sumpred_; lf=log(sumpred/nevents); diffpred=0; if lf lt -.15 or lf gt .15 then diffpred=1; run; %if &empcal eq T %then %do; data mod; set mod; where label1 eq: 'N' and index(label1, 'uster') gt 0; run; data mod; set mod end=_end_; if _end_; call symput('_nsubj', trim(left(cvalue1))); run; %end; title%eval(&titlenum + 1) "Results of running RELRISK8 on data &data"; %if &empcal eq T %then %do; title%eval(&titlenum + 2) " GEE (generalized estimating equations) model with TYPE=&reptype "; title%eval(&titlenum + 3) "Outcome is &depend . &_nevents_ events (&_pct_ %) in &nobs trials of &_nsubj subjects."; %end; %else %do; title%eval(&titlenum + 3) "Outcome is &depend . &_nevents_ events (&_pct_ %) in &nobs trials."; %end; %if &showsump eq T %then %do; %if &diffpred eq 1 %then %do; title%eval(&titlenum + 4) "But model predicts &_sumpred_ events."; %end; %end; run; ods output close; * Save the parameters in case you want to use the previous betas for the initial values; data prevval; set param; run; /* calculate relative risk and 95% confidence bounds */ data param; set param; length varname $20 ; if parameter ne ' ' then varname=parameter; else if parm ne ' ' then varname=parm; pval=max(probz, probchisq); format pval pvalue6.4; units=1; if varname ne 'Scale' and StdErr ne 0 then do; rrest=exp(units*estimate); rrcilo95=exp(units*estimate-1.96*units*stderr); rrcihi95=exp(units*estimate+1.96*units*stderr); end; label rrest = 'relative risk' rrcilo95 = '95% CI Low' rrcihi95 = '95% CI high' varname = 'Parameter' pval = 'P value' level1 = ' ' ; if level1 eq . then level1=.; keep level1 rrest rrcilo95 rrcihi95 varname pval estimate units stderr ; run; proc print label data=param noobs; %if %length(&class) eq 0 %then %do; var varname Estimate Stderr pval units rrest rrcilo95 rrcihi95; %end; %else %do; var varname level1 estimate stderr pval units rrest rrcilo95 rrcihi95; %end; run; %if %length(&outdat) gt 0 %then %do; data &outdat; set param; %end; %end; proc datasets nolist; delete _dat_998 _dat_999 parms predstat predict poisspar ppar poissinit _maxpred_ respprof _sumev_ mod prevval param poisconv frqs ; run; quit; %goto realend; %errend: options notes; %realend: options notes obs=max formdlim="&_fdl" ; %mend; %macro numargs(arg, delimit); %if %quote(&arg)= %then %do; 0 %end; %else %do; %let n=1; %do %until (%qscan(%quote(&arg), %eval(&n), %str( ))=%str()); %let n=%eval(&n+1); %end; %eval(&n-1) %end; %mend numargs; %macro nobsnvar (data=, nvarsp=, nobsp=); %global dset nvars nobs; %let dset=&data; %let dsid=%sysfunc(open(&dset)); %if &dsid %then %do; %let nobs=%sysfunc(attrn(&dsid, NOBS)); %let nvars=%sysfunc(attrn(&dsid, NVARS)); %let rc=%sysfunc(close(&dsid)); %end; %else %put Could not open &data--%sysfunc(sysmsg()); %mend; %macro pfrq; proc freq noprint data=_dat_998; tables &depend / out=frqs; run; data frqs; set frqs end=_end_; retain nobs 0; nobs=nobs + count; if &depend eq 1 then do; call symput ('nevents', trim(left(count))); nevents=count; end; if _end_ then do; call symput ('nobs', trim(left(nobs))); pct=round(100* nevents / nobs, .1); call symput('_pct_', trim(left(pct))); end; run; %mend;