/* SURV19str: 01 SEP 04 Miguel Hernan, Roger Logan, James Robins Copyright, 2005, President and Fellows of Harvard College SURV.SAS is a set of SAS macros to run either a marginal structural Cox model or a nested structural accelerated failure time model. References (introductory level): Hernán MA, Brumback B, Robins JM. Marginal structural models to estimate the causal effect of zidovudine on the survival of HIV-positive men. Epidemiology 2000;11:561-570. Hernán MA, Cole SR, Margolick JB, Cohen MH, Robins JM. Structural accelerated failure time models for survival analysis in studies with time-varying treatments. Pharmacoepidemiology and Drug Safety 2005 (in press). */ %macro surv( /* input data */ data=, /* data set name */ id=, /* subject identifier (variable name) */ time=, /* time of follow-up (variable name) first observation for each subject must be time 0 same time scale as 'T' variable below */ outcMSM=, /* outcome variable for MS Cox (pooled logistic) model 1: event, 0: no event */ covMSM=, /* list of variables included in MS Cox model (baseline variables and time-varying intercept) */ classMSM=, /* list of categorical variables in covMSM */ AMSM=, /* exposure variable in MS Cox model */ A=, /* dichotomous exposure variable in SN AFT model and outcome in models for treatment */ covAd=, /* list of variables included in treatment model (denominator) */ classAd=, /* list of categorical variables in covAd */ covAn=, /* list of variables included in treatment model (numerator) */ classAn=, /* list of categorical variables in covAn */ Cense=, /* dichotomous variable indicating administrative (by eof) censoring 1: yes, 0: no (MS Cox model only, optional) */ covCed=, /* list of variables included in adm. censoring model (denominator) */ classCed=, /* list of categorical variables in covCed */ covCen=, /* list of variables included in adm. censoring model (numerator) */ classCen=, /* list of categorical variables in covCen */ Cens=, /* dichotomous variable indicating censoring 1: yes, 0: no MSM Cox: 1) if using 'cense' for censoring by administrative eof, then cens should be 1 only when censored before eof 2) if not using cense, then cens is 1 when censored due to any cause (including administr. eof) SN AFT 1) for igncens=0, 'cens' is censoring before eof 2) for igncens=1, 'cens' is any censoring */ covCd=, /* list of variables included in censoring model (denominator) */ classCd=, /* list of categorical variables in covCd */ covCn=, /* list of variables included in treatment model (numerator) */ classCn=, /* list of categorical variables in covCn */ initA=, /* time of treatment initiation (for monotonic exposures only) */ T=, /* time of event, same time scale as 'time' variable */ maxC=, /* maximum time when death would have been observed for each individual */ LSM=, /* time-dependent variable in product term with exposure (SN AFT model only, optional, one-parameter model if no LSM specified */ eligib=, /* eligibility variable 1: yes, 0: no */ weight=, /* user-defined weights for SN AFT model */ /* data analysis specifications */ MSMCox= 0, /* 0 for SN AFT model, 1 for MS Cox model*/ bootstrap= 0, /* bootstrap estimate of variance */ boot_one_a_time= 1, /* 1: one bootstrap sample at a time, 0 by processing */ nboot= 200, /* number of bootstrap samples */ bseed= 12345, /* random number seed for boostrap */ /* data analysis specifications for MSMCox=0 */ itt= 0, /* "intention to treat" approach: 1 yes, 0 no */ igncens= 0, /* ignorable adm. censoring: 1 yes, 0 no (see above for details) */ maxCall=, /* maximum allowed time when death would have been observed for any subject */ fixedbaseline= 1, /* 1 indicates that maxCall is measured from start of f-u 0 indicates that maxCall is measured from each time NOTE: T, maxC, time are always measured from start of follow-up */ hfunc= 0, /* Function of H used in estimating equation 0: delta 1: delta*H 2: delta*log(H) NOTE: when dr=1 and hfunc = 2, then E[delta*log(H)] estimated by least squares */ conf_interval= 1, /* invert test for finding confidence interval (with numpsi = 1) */ minpA= 0, /* exclude observations w/ probability of A maxpA only with msmcox=0 */ rr= 0, /* computes causal rate ratio from AFT parameter only with msmcox=0 and numpsi= 1 rr=1 increases CPU time substantially. it is better no to use it simultaneously with conf_interval=1 */ /* computational specifications */ MSM_posdef= 0, /* use positive definite estimate in calculation of msm variance */ gridsearch= 0, /* 1: grid search 0: no grid search */ use_grid_start= 0, /* use grid search results to determine starting point of optimization method. Only with gridsearch= 1. */ Psi1upper= 3,Psi1lower= -3, /* bounds of optimization */ Psi1step= 0.1, Psi2upper= 1,Psi2lower= -1, Psi2step= 0.1, est_method= 1, /* 0 = No numerical optimization 1 = Nelder-Mead Simplex 2 = Simulated Annealing 3 = Find root of estimating equation Only with numpsi= 1 4 = Smoothes test statistic (spline) and finds minimum Only with gridsearch= 1 */ Psi1init= 0, Psi2init= 0, /* initial values for N-M Simplex */ conf_min= 0.001, /* lower bound for sumeef2 used to calculated confidence region. Also used as a potential stopping criteria for the simulated annealing method Also used for a way of determining if convergence was achieved in the Nelder-Mead simplex method for including a given bootstrap sample in the analysis */ countmax= 100, /* maximum number of evaluations for finding bounds of confidence interval */ /* parameters for Simulated Annealing */ t0_eval= 15, /* number of iterations used in finding an initial temperature */ sim_eval= 20, /* number of iterations or steps of decreasing the temperature */ ddwell= 25, /* number of evaluations used in the equilbrium module */ rho= 0.1, /* used in calculating new point x = rho * t * tan(theta) where theta is a random point in the inerval [-rrange,rrange] **/ rrange= 1.57, /* range of points used in evaluating tangent for finding new test point. default value is pi/2 */ tscale= 1.0, /* used in the calculation of the next temperatur t = t0/(1+i*tscale) */ tmin= 1.0 , /* a possible lower bound for annealing temperature */ K= 1.0, /* used in the probability of accepting a new point */ /* output options */ print_opt= 1, /* print out the results for the optimization method when est_mthod = 1 */ optn_level= 0, /* output level for Nelder Mead Simplex method : est_method = 1 0 = none 5 = all */ print_boot= 0, /* print out the results for all bootstrap samples */ plot_graph= 1, save_graph= 0, /* if plot_graph = 1 and save_graph = 1 then output graph in some format given by variables gdevice and ftype. default type: pdf */ gdevice= pdf, ftype= pdf, save_results= 1, /* save grid search/bootstrap results into sas data sets */ save_weights= 0, /* save the weights used in the MSM and SM models for censoring and/or sw */ libsurv= '.' , /* libname for saved data sets */ savename= results, /* name of file for saving results */ estprob= 1, rnseed= , debug= 0 /* 1 = keep intermediary data sets, 0 = delete intermediary data sets */ ); /* The program is structured into 9 main modules: %SURV, the main driver, inputs data, analysis and computing specifications %SETDATAFT prepares the data from SURV or GENDATA to be analyzed by the macros %MSMCOX and %SNMAFT %SNMAFT estimates the parameters of a structural nested accelerated failure time model with a dichotomous exposure. It implements non doubly-robust and doubly robust estimators, the minmin and the maxmin (beta version) methods for dealing with administrative censoring, and IPCW to deal with censoring by loss to follow-up/competing risks. It also allows to consider administrative censoring as ignorable. %MSMCOX estimates the parameters of a marginal structural Cox model as described in Hernan, Brumback, Robins (JASA 2001). To add: analytic variance, doubly-robust estimation %BIGBOOT, boostrap estimates of variance %FIND_ROOT, bisection method to find root of estimating equation of one-dimensional parameters in SNMAFT %SIM_ANN, simulated annealing to find minimum of objective function in SNMAFT %NUMARGS, auxillary macro for finding number of entries in a macro list %OUTPUT_RESULTS, macro for displaying final results %REMOVE_CLASS, remove entries in one list from a second list %CONVERT_CLASS, converts class variables to numeric variables so they can be used in the models for E[delta * f(H)| covDelta] in %SNMAFT %REMOVE_EXTRANEOUS, removes entries in a class list if they are not contained in the corresponding covariate list, e.g. classDelta and covDelta NOTES regarding censoring igncens=0 no assumptions regardind censoring by eof (non ignorable) the random potential censoring time maxC is always added to models for A and Cens because, even if maxC is not a risk factor for H(psi), it is a certainly risk factor for C(psi) and therefore predicts delta(psi). the option maxCall allows to fix a non random eof (shorter than maxC for the subject with the largest potential censoring time). that is, the length of the study is shortened. still need to adjust for maxC as some subjects' delta(psi) is predicted by it. because the macro does not use efficient estimators, this shortening (with the corresponding loss of events) may paradoxically increase efficiency */ %local individuals chilevel cov_betas n_betas num_events events_used ; /* consistency checks and input error management (to be further developed) */ %let sample = 1; %let get_cov = 0; %if %eval(&MSMCox) = 1 %then %do; %let gridsearch = 0; %let est_method = 0; %let msm_var = 1; %let conf_interval = 0 ; %let results = msmsw; %if %bquote(&Cens)= %then %do; %let covCd = ; %let classCd = ; %let covCn = ; %let classCn = ; %let covCed = ; %let classCed = ; %let covCen = ; %let classCen = ; %end; %if %bquote(&cense) = %then %do; %let covCed = ; %let classCed = ; %let covCen = ; %let classCen = ; %end; %end; %else %do; %let msm_var = 0; %let cov_betas = &A ; %let numpsi = 1 ; %if %bquote(&LSM)^= %then %do; %let cov_betas = &A &A.&LSM ; %let numpsi = 2; %end; %if %eval(&numpsi) ne 1 %then %let conf_interval = 0; %if %eval(&est_method) = 3 %then %let conf_interval = 0; %let n_betas = %eval(&numpsi); %if %eval(&numpsi)^= 1 %then %let rr = 0; %let results = snm ; %end; %if %eval(&gridsearch)= 0 %then %do; %let use_grid_start = 0; %if %eval(&est_method)^= 3 %then %let plot_graph= 0; %end; %if %eval(&gridsearch)= 1 & %eval(&use_grid_start)= 0 %then %let boostrap = 0; %if %eval(&bootstrap) = 1 %then %do ; %let conf_interval = 0; %let MSM_var = 0; %let get_cov = 1; %end; %if %eval(&bootstrap) = 0 %then %do ; %let nboot = 0; %let boot_one_a_time = 0; %end; %if %eval(&msmcox) = 1 & &msm_var = 1 %then %do; data varbetas ;run; %end; %if %eval(&est_method) = 2 %then %sim_ann;; %let bsample = 0 ; /* index bootstrap samples (setdata) */ %let _loop_start = 0; %let _loop_end = %eval(&nboot); %if &bootstrap = 1 & &boot_one_a_time = 0 %then %let _loop_end = 0; /* results data sets */ %if %eval(&MSMCox) = 1 %then data msmsw; %else data snm;; run; %setdataft ; /* process class statements and possible bootstrap samples */ %if %eval(&MSMCox) = 0 %then %snmaft; %if %eval(&MSMCox) = 1 %then %msmcox; /* output results */ %if (%eval(&est_method)^= 0 | %eval(&msmcox) = 1) %then %do; data raw_results ; set &results ; run; data name_order; length name $20; %do i = 1 %to &n_betas; %let vartmp = %scan(&cov_betas,%eval(&i)); name = "&vartmp"; order = %eval(&i); output; %end; run; proc sort data=name_order; by name; run; data _results; length name $20; set &results; %if &bootstrap = 1 & &msmcox = 0 %then if objfunc < &conf_min or _bsample_ = 0;; for_merge = 1; run; %if &bootstrap = 1 & &msmcox = 0 %then %do; /* number of bootstrap samples that converged to find an estimate of psi */ proc freq data=_results noprint ; table name/ out=_boot_freq; where name = "&A"; run; data _boot_freq; set _boot_freq; if count ne %eval(&nboot + 1) then do; true_count = count - 1; put '***************** WARNING ******************'; put 'Sample number ' sample ': only ' true_count " out of &nboot bootstrap samples converge to less than &conf_min" ; end; run; %end; %if %eval(&print_boot) = 1 %then %do; proc print data=_results noobs; var _bsample_ name estimate %if &msmcox = 0 %then objfunc; ; run; %end; proc sort data = _results ; by _bsample_; run; proc transpose data = _results out=_results_tr (drop = _name_ ); var estimate ; by _bsample_; run; %if &get_cov = 1 %then %do; proc corr data = _results_tr cov nocorr nosimple outp = _myout noprint; var col1-col%eval(&n_betas); %if &bootstrap = 1 %then %do; where _bsample_ ne 0 ; %end; run; data _myout2; set _myout (drop = _name_ ); if _type_ = 'MEAN' or _type_ = 'STD' ; rename %do i = 1 %to &n_betas; %let vartmp = %scan(&cov_betas,%eval(&i)); col&i = &vartmp %end;; for_merge = 1; run; proc transpose data = _myout2 out = _myout2t (rename = ( _name_ = name col1 = mean col2 = std)); %if &bootstrap = 0 %then by for_merge;; run; proc sort data = _myout2t; by name ; run; %end; /* get_cov = 1 */ proc sort data = _results; by name; run; data &results ; %if &get_cov = 1 %then merge _results _myout2t; %else set _results;; by name ; if _bsample_ = 0 ; %if &get_cov = 1 & &msm_var = 0 %then %do; label p_value = 'Pr > |Z|'; format p_value PVALUE8.4 ; lb = estimate - 1.96 * std; ub = estimate + 1.96 * std; z = estimate / std; p_value = 2*(1-probnorm(abs(z))); %if &bootstrap = 1 %then bbias = estimate - mean; ; %end; drop for_merge; run; %if %sysevalf(&MSMCox) = 0 %then %do; title "ITT= %eval(&itt), Ign.cens= %eval(&igncens), Max EOF= &maxCall."; %end; /* print out the results */ data &results; set &results ; %if &msmcox =1 & &msm_var = 1 %then %do; label p_value = 'Pr > |Z|' ; format p_value PVALUE8.4 ; %end; run; proc sort data = &results; by name; run; %if &get_cov = 1 %then %do; data _cov; length name $20; set _myout (rename = ( _NAME_ = name )) ; if _type_ = 'COV'; %do i = 1 %to &n_betas; %let vartmp = %scan(&cov_betas,%eval(&i)); rename col&i = &vartmp ; if name = "COL&i" then name = "&vartmp"; %end; drop _type_ ; run; proc sort data = _cov; by name; run; %end; data &results; merge &results name_order %if &get_cov = 1 %then _cov ;; by name ; %if &get_cov = 0 | &msm_var = 1 %then %do; rename %if &msm_var = 1 %then msm_std = std msm_lb = lb msm_ub = ub; %if &conf_interval = 1 %then ub1 = ub lb1 = lb;; %end; drop _bsample_ ; run; proc sort data = &results out = &results (drop = order ); by order; run; %if &get_cov = 1 & &msm_var = 0 %then %do; data cov_betas; merge _cov name_order; by name; run; proc sort data = cov_betas out = cov_betas (drop = order) ; by order; run; %end; %end; /* est_method ne 0 */ %if %eval(&gridsearch) = 1 or %eval(&est_method) ^= 0 or %eval(&msmcox) = 1 %then %output_results; %if %eval(&save_results) = 1 %then %do; libname libsurv &libsurv; data libsurv.&savename; set &results ; run; %let cov_name = cov_betas; %if &msm_var = 1 %then %let cov_name = varbetas; %if &bootstrap = 1 %then %do; data libsurv.boot_results ; set _results; run; %end; %if &msm_var = 1 | &get_cov = 1 %then %do; data libsurv.cov_betas ; set &cov_name ; run; %end; %if &save_weights = 1 %then %do; data libsurv.weights; set weights; run; %end; %end; %if %eval(&debug) = 0 %then %do; proc datasets library = work nolist; delete _unoboot_ _cov _myout _myout2 _myout2t _results_tr name_order _results all_results raw_results %if %eval(&rr) = 1 %then _parms_ ; %if %eval(&bootstrap) = 1 %then _justid ; %if %eval(&save_results) = 0 %then %do; %if %eval(&bootstrap) = 1 %then _boottmp_ _boot_freq; %if %eval(&est_method) = 2 %then sa_graph; %end;; run; quit; %end; title; %mend surv; %macro bigboot(datain=,_nboot=1,_id=,sortvar=,dataout=); %do i = 1 %to &_nboot; %if %eval(&boot_one_a_time) = 0 %then %let stmp = %sysevalf(&bseed * &i); %else %let stmp = %sysevalf(&bseed * &bsample) ; proc surveyselect data=_justid noprint method = urs n= &individuals seed = &stmp out = _idsamp; run; data _tmp; merge _idsamp (in = _a) &datain ; by &_id; if _a; _bsample_ = %if %eval(&boot_one_a_time) = 0 %then &i; %else &bsample; ; run; data _tmp; set _tmp ; previd = lag(&_id); prevhits = lag(numberhits); if _n_ = 1 then _tmpid = 0; else do; if previd ^= &_id then _tmpid = _tmpid + prevhits; end; do i = 1 to numberhits; _newid_ = _tmpid + i; output; end; retain _tmpid ; drop _tmpid previd prevhits i numberhits; run; proc sort data= _tmp; by _newid_ &sortvar; run; %if %eval(&i) = 1 %then %do; data &dataout; set _tmp; run; %end; %else %do; proc append base = &dataout data=_tmp; run; %end; %end; %if %eval(&debug) = 0 %then %do; proc datasets library = work nolist; delete _idsamp _tmp %if &boot_one_a_time = 0 %then _justid ; ; run; quit; %end; %mend bigboot; %macro setdataft; proc sort data= &data out=_cero_; by &id &time; run; %let nclass = 7 ; %let class_var = %nrstr( &classMSM &classAd &classAn &classCd &classCn &classCed &classCen); %let cov_var = %nrstr( &covMSM &covAd &covAn &covCd &covCn &covCed &covCen); %let aaaa = %unquote(&class_var) ; %if %bquote(&aaaa) ^= %then %do; %remove_extraneous; %convert_class; %end; %let initial_name = _cero_; %let logistic_data = _ceroboot_; %if &msmcox = 1 %then %do; %if &bootstrap = 0 | ( &boot_one_a_time = 1 & &bsample = 0 ) %then %let logistic_data = _cero_; %end; %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time) = 0 %then %let final_name = _unoboot_; %else %let final_name = _uno_; /* one dataset per iteration */ data _cero_ (index = (_bsample_)); set _cero_ end = _end_; by &id ; if first.&id then do; _subjects_ + 1; %if &msmcox = 0 %then if &T ^= . then _nevents + 1; ; end; if last.&id then do; %if &msmcox = 1 %then if &outcMSM = 1 then _nevents + 1 ;; end; _newid_ = &id; _bsample_ = 0; if _end_ then do; call symput('individuals',_subjects_); call symput('num_events',left(_nevents)); end; drop _subjects_ ; run; %if %eval(&bootstrap) = 1 %then %do; /* create the data set _justid once only */ %if %eval(&boot_one_a_time) = 0 OR %eval(&bsample) = 0 %then %do; data _justid ; set _cero_; by _newid_; if first._newid_; keep _newid_; run; %end; %if %eval(&boot_one_a_time) = 0 %then %do; /* create all bootstrap samples at once */ %bigboot(datain=_cero_ , _nboot=&nboot , _id=_newid_, sortvar=&time, dataout=_boottmp_ ); data _ceroboot_ ; set _cero_ _boottmp_ ; run; %let initial_name = _ceroboot_; %end; %if %eval(&boot_one_a_time) = 1 AND %eval(&bsample) > 0 %then %do; %bigboot( datain=_cero_ , _id=_newid_, sortvar=&time, dataout=_ceroboot_ ); %let initial_name = _ceroboot_ ; %end; %let id = _newid_; /* sequential id starting a 1 */ %end; %mend setdataft; %macro msmcox; title 'Marginal Structural Cox Model, Stabilized Weights'; %let cov_betas = Intercept &Amsm &covMSM ; %let n_betas = %numargs(&cov_betas) ; %do _loop = %eval(&_loop_start) %to %eval(&_loop_end); %let bsample = &_loop; /* name of results dataset */ %let logistic_data = _ceroboot_; %if &bootstrap = 0 | ( &boot_one_a_time = 1 & &bsample = 0 ) %then %let logistic_data = _cero_; %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time) = 0 %then %let final_name = _unoboot_; %else %let final_name = _uno_; /* one dataset per iteration */ %if %eval(&boot_one_a_time) = 1 AND %eval(&bsample) > 0 %then %do; %bigboot(datain=_cero_ ,_id=_newid_,sortvar=&time,dataout=_ceroboot_ ); %let initial_name = _ceroboot_ ; %let id = _newid_; /* sequential id starting a 1 */ %end; %let modelAd = ; %let modelAn = ; %let modelEd = ; %let modelEn = ; %let modelCd = ; %let modelCn = ; %if %eval(&estprob) = 1 %then %do; proc logistic data = &logistic_data descending %if (&bootstrap = 1 & (&boot_one_a_time = 0 OR &bsample ne 0) ) %then noprint; ; ods exclude ClassLevelInfo TypeIII Association FitStatistics GlobalTests Oddsratios; %if %bquote(&initA)^= %then where &time <= &initA | &initA = .;; model &A = &covAd ; %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time) = 0 %then by _bsample_ ;; output out = _modelAd_ (keep= _bsample_ &id &time pA_d) p = pA_d; title "Model for denominator of stabilized weights for &A"; run; %let modelAd = _modelAd_; %if %bquote(&Cense)^= and %eval(&MSMCox) = 1 %then %do; proc logistic data = &logistic_data %if (&bootstrap = 1 & (&boot_one_a_time = 0 OR &bsample ne 0) ) %then noprint; ; ods exclude ClassLevelInfo TypeIII Association FitStatistics GlobalTests Oddsratios; where &time>0 ; model &cense = &covCed; %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time) = 0 %then by _bsample_ ;; output out = _modelEd_ (keep= _bsample_ &id &time pE_d) p = pE_d; title "Model for denominator of stabilized weights for &cense "; run; %let modelEd = _modelEd_; %end; %if %bquote(&Cens)^= %then %do; proc logistic data = &logistic_data %if (&bootstrap = 1 & (&boot_one_a_time = 0 OR &bsample ne 0) ) %then noprint; ; where &cens ne .; ods exclude ClassLevelInfo TypeIII Association FitStatistics GlobalTests Oddsratios; model &Cens = &covCd; %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time) = 0 %then by _bsample_ ;; output out = _modelCd_ (keep= _bsample_ &id &time pC_d) p = pC_d; title "Model for denominator of stabilized weights for &cens"; run; %let modelCd = _modelCd_; %end; %end; /* estprob */ proc logistic data = &logistic_data descending %if (&bootstrap = 1 & (&boot_one_a_time = 0 OR &bsample ne 0) ) %then noprint;; ods exclude ClassLevelInfo TypeIII Association FitStatistics GlobalTests Oddsratios; %if %bquote(&initA)^= %then where &time <= &initA | &initA = . ;; model &A = &covAn; %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time) = 0 %then by _bsample_ ;; output out = _modelAn_ (keep= _bsample_ &id &time pA_n) p = pA_n; title "Model for numerator of stabilized weights for &A"; run; %let modelAn = _modelAn_; %if %bquote(&Cense)^= %then %do; proc logistic data = &logistic_data %if (&bootstrap = 1 & (&boot_one_a_time = 0 OR &bsample ne 0) ) %then noprint;; ods exclude ClassLevelInfo TypeIII Association FitStatistics GlobalTests Oddsratios; where &time>0 ; model &cense = &covCen; %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time) = 0 %then by _bsample_ ;; output out = _modelEn_ (keep= _bsample_ &id &time pE_n) p = pE_n; title "Model for numerator of stabilized weights for &cense"; run; %let modelEn = _modelEn_; %let numE = numE; %end; %if %bquote(&Cens)^= %then %do; proc logistic data = &logistic_data %if (&bootstrap = 1 & (&boot_one_a_time = 0 OR &bsample ne 0) ) %then noprint; ; where &cens ne .; ods exclude ClassLevelInfo TypeIII Association FitStatistics GlobalTests Oddsratios; model &Cens = &covCn; %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time) = 0 %then by _bsample_ ;; output out = _modelCn_ (keep= _bsample_ &id &time pC_n) p = pC_n; title "Model for numerator of stabilized weights for &cens"; run; %let modelCn = _modelCn_; %let numC = numC; %end; title; data &final_name (index = (_bsample_) ); merge &logistic_data &modelAd &modelEd &modelCd &modelAn &modelEn &modelCn end = end; by %if %eval(&bootstrap) = 1 %then _bsample_; &id &time; if first.&id then do; numA=1; denA=1; numE=1; denE=1; numC=1; denC=1; /* pE_n=1; pE_d=1; */ %if &msm_var = 1 %then %do; den = 1; num = 1; %end; _occasion_ = -1; end; retain numA denA numE denE numC denC ; %if &msm_var = 0 %then if &outcMSM ne .;; %if &msm_var = 1 %then %do; if pc_n=. then do; pc_n=1; pc_d=1; end; %end; %if %bquote(&Cens)^= %then %do; denC = denC* pC_d; numC = numC* pC_n; %end; %if %bquote(&Cense)^= %then %do; if pE_d = . then pE_d = 1; if pE_n = . then pE_n = 1; numE = numE* pE_n; denE = denE* pE_d; %end; %if %bquote(&initA)= %then %do; if &A = 1 then do; numA = numA*pA_n; denA = denA*pA_d; end; else if &A = 0 then do; numA = numA*(1-pA_n); denA = denA*(1-pA_d); end; %end; %else %do; if pA_d = . then pA_d = 1; if pA_n = . then pA_n = 1 ; if &time<&initA or &initA=. then do; numA = numA*(1-pA_n); denA = denA*(1-pA_d); end; else if &time = &initA then do; numA = numA*pA_n; denA = denA*pA_d; end; else do; numA = numA; denA = denA; end; %end; num = numA; den = denA; %if %bquote(&Cens)^= %then %do; num = num * numC ; den = den * denC ; %end; %if %bquote(&Cense)^= %then %do; num = num * numE; den = den * denE; %end; sw = num / den; %if &msm_var = 1 %then %do; dies = 0; if &T ne . then dies = 1 ; %end; _occasion_ + 1; if &outcMSM = 1 then _events_used + 1 ; if end then call symput('events_used',_events_used) ; keep &id &classMSM &outcMSM &A &AMSM &covMSM sw _bsample_ _occasion_ %if (&save_weights = 1) | (&msm_var = 1) %then &time ; %if &msm_var = 1 %then %do; dies &covAd &covCd pA_d pC_d &cens %if %bquote(&cense)^= %then &cense pE_d ; %end; ; run; %if &save_weights = 1 %then %do; data weights; set &final_name (keep = &id &time sw ) ; run; %end; %if (&bootstrap = 0 | (&boot_one_a_time = 1 AND &bsample = 0) ) %then %do; proc univariate data = &final_name; title 'Distribution of stabilized weights'; var sw; run; %end; %if &msm_var = 1 %then %do; %let n_covAd = %numargs(&covAd) ; /* number of variables in model for pA_d : covAd */ %let scoreAd = scoreA_int ; %let n_covMSM = %numargs(&covMSM); %do i = 1 %to &n_covAd; %let vartmp = %scan(&covAd, %eval(&i)); %let scoreAd = &scoreAd scoreA_&vartmp ; %end; %let scores = &scoreAd ; %let n_scores = %eval(&n_covAd + 1); %if %bquote(&Cens)^= %then %do; %let n_covCd = %numargs(&covCd); /* number of variables in model for pC_d: covCd */ %let scoreCd = scoreC_int; %do i = 1 %to &n_covCd; %let vartmp = %scan(&covCd, %eval(&i)); %let scoreCd = &scoreCd scoreC_&vartmp ; %end; %let scores = &scores &scoreCd ; %let n_scores = %eval(&n_scores + &n_covCd + 1 ); %if %bquote(&cense)^= %then %do; %let n_covCed = %numargs(&covCed); /* number of variables in model for pE_d : covCed */ %let scoreCed = scoreCe_int; %do i = 1 %to &n_covCed; %let vartmp = %scan(&covCed, %eval(&i)); %let scoreCed = &scoreCed scoreCe_&vartmp ; %end; %let scores = &scores &scoreCed ; %let n_scores = %eval(&n_scores + &n_covCed + 1 ); %end; %end; %else %do; %let n_covCd = 0; %let scoreCd = ; %end; data _for_follow; set &final_name /* _cero_ */ end = _end_; by &id &time; _keep = 0; if first.&id then _subjects + 1; if last.&id then do; if &T = . then followup = _occasion_ ; if &T ne . then followup = floor(&T); _keep = 1; end; if _end_ then call symput('individuals',_subjects); if _keep = 0 then delete; keep &id followup ; run; data &final_name; merge &final_name _for_follow; by &id; derien = 1; run; data nuisance; set &final_name; by &id &time; retain %do i = 1 %to &n_scores; %let vartmp = %scan(&scores,%eval(&i)); &vartmp 0 %end;; if &time=0 then do; %do _i = 1 %to &n_scores; %let vartmp = %scan(&scores, %eval(&_i)); &vartmp = 0; %end; end; _tmp = &A - pA_d ; %if %bquote(&Cens)^= %then %do; if &cens = . then _tmp = 0; %end; scoreA_int = _tmp + scoreA_int; %do _i = 1 %to %eval(&n_covAd); %let vartmp = %scan(&covAd, %eval(&_i)); _tmp = &vartmp * (&A - pA_d) ; %if %bquote(&Cens)^= %then %do; if &cens = . then _tmp = 0; %end; scoreA_&vartmp = _tmp + scoreA_&vartmp ; %end; %if %bquote(&Cens)^= %then %do; if &cens > . then do; scoreC_int = (&Cens - 1 + pC_d)+scoreC_int; %do _i = 1 %to &n_covcd; %let vartmp = %scan(&covcd ,%eval(&_i)); scoreC_&vartmp = &vartmp * (&Cens - 1 + pC_d ) + scoreC_&vartmp; %end; %if %bquote(&cense)^= %then %do; scoreCe_int = (&Cense - 1 + pE_d)+scoreCe_int; %do _i = 1 %to &n_covCed; %let vartmp = %scan(&covCed ,%eval(&_i)); scoreCe_&vartmp = &vartmp * (&Cense - 1 + pE_d ) + scoreCe_&vartmp; %end; %end; end; %end; if _occasion_ = followup then do; keep &id &scores ; output; end; run; data &final_name; merge &final_name nuisance; by &id ; if &outcMSM = . %if %bquote(&cens)^= %then or &cens = 1 or &cens = . ; then delete; rien = 1; run; %end; /* msm_var = 1*/ %if &debug = 0 %then %do; proc datasets library=work nolist; delete _modelAd_ &modelAn &modelEd &modelEn &modelCd &modelCn %if &boot_one_a_time = 0 %then _cero_ ; %if %eval(&debug) = 0 %then _ceroboot_ ;; run; quit; %end; proc logistic data= &final_name descending outest= emppest (keep = Intercept &AMSM &covMSM %if %eval(&bootstrap) = 1 AND %eval(&boot_one_a_time) = 0 %then _bsample_ ; ) %if (&bootstrap = 1 AND (&boot_one_a_time = 0 OR &bsample ne 0) ) %then noprint;; ods exclude ClassLevelInfo TypeIII Association FitStatistics GlobalTests; model &outcMSM = &AMSM &covMSM; weight sw; %if %eval(&bootstrap) = 1 AND %eval(&boot_one_a_time) = 0 %then by _bsample_ ;; title "Weighted model for outcome &outcMSM"; run; proc transpose data = emppest out = emppest (drop = _LABEL_ rename = (col1 = estimate _NAME_ = name)); %if %eval(&bootstrap) = 1 AND %eval(&boot_one_a_time) = 0 %then by _bsample_ ;; run; %if &msm_var = 1 %then %do; proc sort data= &final_name out= &final_name; by _occasion_ followup dies descending &id; run; data &final_name; set &final_name; by _occasion_ followup dies descending &id ; retain numeriq 0; if first.followup then numeriq=1; followup=followup+numeriq/1000000; numeriq=numeriq+1; keep &id &time &outcMSM sw &Amsm &covMSM &scores followup _occasion_; run; proc sort data= &final_name ; by descending _occasion_ descending followup; run; proc freq data=&final_name noprint; table _occasion_ / out = index_set; run; data index_set; set index_set end = _end_; retain countsum 0 ; tot = round( count /(0.01* percent)) ; endindex = tot - countsum ; countsum = countsum + count; startindex = endindex - count + 1 ; if _end_ then call symput('_tlevels_',_n_); rename count = _freq_ ; run; %if &debug = 0 %then %do; proc datasets library = work nolist; delete censtemp c_coeffs _for_follow hazpar drbeta ; quit; %end; proc iml; use index_set; read all var {_FREQ_} into indexes; read all var {startindex} into startindex; read all var {endindex} into endindex; use nuisance; read all var{ &scores } into AC_scoresII ; use &final_name; read all var{ &scores } into AC_scores where ( _occasion_ = 0); read all var{ &Amsm &covMSM } into PRED; read all var {&outcMSM} into DEAD; read all var {sw} into weight; final_indiv = nrow(pred); ones = J(final_indiv,1,1); pred = ones || pred ; use emppest; read all var {estimate } into betas; regressor = pred*betas; hazard = 1/(1+exp(-regressor)); U = weight#(dead-hazard); new_u = U # pred; individuals = nrow(AC_scores); n_beta = ncol(new_u); NUvec=J(individuals,n_beta,0); k=&_tlevels_-1; do while (k>=0); startk=startindex[k+1]; endk=endindex[k+1]; indexk=indexes[k+1]; NUvec[1:indexk,] = NUvec[1:indexk,] + new_u[startk:endk,]; k=k-1; end; lambda0 = t(NUvec) * NUvec; gamma0 = t(NUvec) * AC_scores; omega0 = t(AC_scoresII) * AC_scoresII; omega0= ginv(omega0); %if &msm_posdef = 1 %then %do; Unew = NUvec; Eus = gamma0; Ess_inv = omega0; part2 = Eus * Ess_inv ; part2 = part2 * t(AC_scores ) ; tmp = Unew - t(part2) ; meat = t(tmp) * tmp ; %end; %else %do; meat=(lambda0)-gamma0*omega0*t(gamma0); %end; /* derivative of U with respect to beta */ dE_dbeta = hazard # (1-hazard) ; dE_dbeta = -1*(weight # dE_dbeta) # pred ; n_beta = ncol(pred); new_deriv = J(n_beta,n_beta,0); do i = 1 to n_beta ; tmp = pred[,i]; tmp2 = dE_dbeta # tmp ; deriv_tmp = tmp2[+,]; new_deriv[i,] = deriv_tmp; end; deriv = new_deriv; varbetas= ( ginv(deriv)*meat*t( ginv(deriv))); create varbeta_kk from varbetas ; append from varbetas; sterrB=varbetas[2,2]##0.5; name = "12345678901234567890"; /* needed to define a character variable */ create _kk_ var{name msm_std msm_lb msm_ub z p_value }; %do i = 1 %to &n_betas; sterrb = varbetas[&i,&i]##0.5 ; %let vartmp = %scan(&cov_betas, %eval(&i)); name = "&vartmp" ; betasw = betas[&i]; msm_lb = betas[&i] - 1.96 * sterrb; msm_ub = betas[&i] + 1.96 * sterrb; msm_std = sterrb; z = betasw / sterrb; p_value = 2*(1.0 - probnorm(abs(z))); append; %end; quit; data _kk_; merge emppest _kk_ ; _bsample_ = &bsample; run; %end; %if &msm_var = 1 %then %do; data varbeta_kk; length name $20; set varbeta_kk; rename %do i = 1 %to %eval(&n_betas ) ; %let vartmp = %scan(&cov_betas,%eval(&i)); col&i = &vartmp %end;; array names[&n_betas] $20 names1-names&n_betas ; %do i = 1 %to %eval(&n_betas ); %let vartmp = %scan(&cov_betas, %eval(&i)); names&i = "&vartmp " ; %end; name = names[_n_ ]; drop names1-names&n_betas; run; data varbetas; set varbetas varbeta_kk; if &A ne . ; run; %end; /* msm_var = 1 */ %else %do; data _kk_; set emppest ; if _bsample_ = . then _bsample_ = &bsample ; run; %end; data msmsw; set msmsw _kk_; if estimate ne .; run; %if &debug = 0 %then %do; proc datasets library=work nolist; delete &final_name _kk_ emppest %if &msm_var = 1 %then varbeta_kk index_set nuisance; %if &bootstrap = 1 %then _boottmp_ ;; run; quit; %end; %end; /* end _loop */ title; %mend msmcox; %macro snmaft; %do _loop = %eval(&_loop_start) %to %eval(&_loop_end); %let bsample = %eval(&_loop); %let initial_name = _cero_; %if &bootstrap = 1 & &boot_one_a_time = 0 %then %let initial_name = _ceroboot_ ; %let logistic_data = _ceroboot_; %if &bootstrap = 0 | ( &boot_one_a_time = 1 & &bsample = 0 ) %then %let logistic_data = _cero_; %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time) = 0 %then %let final_name = _unoboot_; %else %let final_name = _uno_; /* one dataset per iteration */ %if %eval(&boot_one_a_time) = 1 AND %eval(&bsample) > 0 %then %do; %bigboot( datain=_cero_,_id=_newid_,sortvar=&time,dataout=_ceroboot_ ); %let initial_name = _ceroboot_ ; %let id = _newid_; /* sequential id starting a 1 */ %end; %let logistic_data = _ceroboot_; %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time) = 0 %then %let final_name = _unoboot_; %else %let final_name = _uno_; /* one dataset per iteration */ data _ceroboot_; set &initial_name end = _end_; by %if %eval(&bootstrap) = 1 %then _bsample_; &id ; if _end_ then do; _chilev_ = cinv(.95,&numpsi ); call symput('chilevel',_chilev_); end; ulti = _end_; if not ulti then set &initial_name (firstobs = 2 keep = &time rename = (&time = _nexttime_) ); %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time)= 0 %then if last._bsample_ then _nexttime_ = . ;; %if %eval(&igncens) = 0 %then %do; %if %bquote(&maxCall)^= %then %do; %if %eval(&fixedbaseline) = 1 %then &maxC = min(&maxCall,&maxC);; %if %eval(&fixedbaseline) = 0 %then &maxC = min(&maxCall+&time,&maxC);; /* q censoring in interval (t,t+1] not defined if censored by eof in same interval or at t */ %if %bquote(&cens)^= %then if &time <= &maxC < _nexttime_ or &time = &maxC then &cens = .;; %end; /* maxCall */ /* add maxC as baseline covariate in models */ %let maxc2= _maxc2_; %let covAd = &covAd &maxC &maxc2; %let covCd = &covCd &maxC &maxc2; %end; /* igncens = 0 */ %else %if %eval(&igncens) = 1 %then %do; %if %eval(&fixedbaseline) = 1 %then &maxC = &maxCall;; %if %eval(&fixedbaseline) = 0 %then &maxC = &maxCall+&time;; %if %bquote(&cens)= %then %do; %let cens = _censeof_; &cens = 0; %end; if &time <= &maxC < _nexttime_ or &time = &maxC then &cens = 1; %end; /* igncens = 1 */ if &time <= &maxC; if &T>&maxC then &T=.; %if %eval(&igncens) = 0 %then _maxc2_= &maxC * &maxC;; drop ulti _nexttime_; run; %let modelAd = ; %let modelCd = ; %if %eval(&estprob) = 1 %then %do; proc logistic data = &logistic_data descending %if (&bootstrap = 1 & (&boot_one_a_time = 0 OR &bsample ne 0) ) %then noprint; ; ods exclude ClassLevelInfo TypeIII Association FitStatistics GlobalTests Oddsratios; %if %bquote(&initA)^= | %bquote(&eligib)^= %then %do; %if %bquote(&eligib)^= & %bquote(&initA)^= %then %do ; where &eligib = 1 AND ( &time <= &initA | &initA = . ) ; %end; %else %if %bquote(&eligib)^= %then %do; where &eligib = 1 ; %end; %else %if %bquote(&initA)^= %then %do; where &time <= &initA | &initA = . ; %end; %end; %if %bquote(&classAd)^= %then class &classAd /param=ref descending;; model &A = &covAd ; %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time) = 0 %then by _bsample_ ;; output out = _modelAd_ (keep= _bsample_ &id &time pA_d) p = pA_d; title "Estimated probabilites for &A : E[&A | covA ]"; run; %let modelAd = _modelAd_; %if %bquote(&Cens)^= %then %do; proc logistic data = &logistic_data %if (&bootstrap = 1 & (&boot_one_a_time = 0 OR &bsample ne 0) ) %then noprint; ; where &cens ne .; ods exclude ClassLevelInfo TypeIII Association FitStatistics GlobalTests Oddsratios; %if %bquote(&classCd)^= %then class &classCd /param=ref descending;; model &Cens = &covCd; %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time) = 0 %then by _bsample_ ;; output out = _modelCd_ (keep= _bsample_ &id &time pC_d) p = pC_d; title "Model for denominator of weights for &Cens"; run; title ; %let modelCd = _modelCd_; * exclude censored by loss to follow-up for AFT model; %if %eval(&rr) = 1 %then %do; data _lastobs_; set &logistic_data; by _bsample_ &id; if last.&id; censlfu= &cens; keep _bsample_ &id censlfu; run; data &logistic_data; merge &logistic_data _lastobs_; by _bsample_ &id; /* if censlfu ne 1; */ run; proc datasets library=work nolist; delete _lastobs_ ; run; quit; %end; %end; /* cens */ %end; /* estprob */ data &final_name (index = (_bsample_) ); merge &logistic_data &modelAd &modelCd end= _end_; by %if %eval(&bootstrap) = 1 %then _bsample_; &id &time; if first.&id then do; denA=1; denC=1; _occasion_ = -1; if _bsample_ = 0 & . < &T < &maxC then _events_used + 1 ; end; retain denA denC; if pA_d = . then pA_d=1; if pC_d = . then pC_d=1; _ipc_= 1 %if %bquote(&Cens)^= %then / pC_d; ; %if %eval(&rr) = 1 %then if censlfu = 1 then _ipc_ = 0 ;; _occasion_+1; _last_ = last.&id; if _end_ then do; %if %eval(&boot_one_a_time) = 1 %then if _bsample_ = 0 then ; call symput('events_used',left(_events_used)); end; %if %eval(&rr) = 0 %then if .< &T <= &maxC ;; data &final_name; set &final_name end = ulti; _durat_ = &maxC - &time; %if %eval(&itt) = 1 %then _durat2_ = &T - &time;; ultim = ulti; if not ultim then set &final_name (firstobs = 2 keep = &time rename = (&time = _nexttime_)); %if %eval(&rr) = 0 %then %do; _interval_ = sum(_nexttime_*(1-_last_),&T*_last_) - &time; %end; %else %do; if .< &T <= &maxC then _interval_ = sum(_nexttime_*(1-_last_),&T*_last_) - &time; else _interval_ = -1; %end; run; %if &save_weights = 1 %then %do; data weights; set &final_name (keep = &id &time _ipc_ ) ; run; proc sort data = weights ; by &id descending &time ; run; data weights; set weights; by &id ; if first.&id then holder = 1; ipcw = _ipc_*holder; holder = ipcw; retain holder; drop holder; run; proc sort data = weights; by &id &time; run; %end; %else %do; %if %eval(&rr) = 1 %then %do; data &final_name; set &final_name; if censlfu ne 1 ; run; %end; %end; proc datasets library=work nolist; delete _modelAd_ &modelCd %if &boot_one_a_time = 0 %then _cero_ ; %if %eval(&debug) = 0 %then _ceroboot_ ;; run; quit; /* need to differentiate between bootstrap and non-bootstrap 1) non-bootstrap will need only one call of original snmaft macro. 2) for bootstrap: a) when boot_one_a_time = 1 will call original snmaft macro once. b) when boot_one_a_time = 0 will call original snmaft macro nboot + 1 times. */ %let _sloop_start = 0; %let _sloop_end = %eval(&nboot) ; %if %eval(&bootstrap) = 1 & %eval(&boot_one_a_time) = 1 %then %do; %let _sloop_start = %eval(&bsample); %let _sloop_end = %eval(&bsample); %end; %do _sloop = %eval(&_sloop_start) %to %eval(&_sloop_end); %if %eval(&bootstrap) = 1 %then %put Evaluating bootstrap sample &bsample;; %if %eval(&bootstrap) = 1 AND %eval(&boot_one_a_time) = 0 %then %do; %let bsample = &_sloop ; data _uno_; set _unoboot_ (idxname = _bsample_) ; if _bsample_ = &bsample; drop _bsample_; run; %end; proc iml; /* GAMMA Function */ start gammaf (Psigamma); load Amin Amax %if %eval(&numpsi)= 2 %then Lmin Lmax; ; Rmin = Psigamma[1] %if %eval(&numpsi) = 2 %then + Lmin*Psigamma[2]; ; Rmax = Psigamma[1] %if %eval(&numpsi) = 2 %then + Lmax*Psigamma[2]; ; %if %eval(&numpsi) = 1 %then %do; if Rmax < 0 then gammahold = exp(Amax*Rmax); else gammahold = exp(Amin*Rmax); %end; %else %do; if Psigamma[2] > 0 then do; if Rmin > 0 then gammahold = exp(Amin*Rmin); else gammahold = exp(Amax*Rmin); end; else do; if Rmax < 0 then gammahold = exp(Amax*Rmax); else gammahold = exp(Amin*Rmax); end; %end; return (gammahold); finish gammaf; /* SUMEE Function */ start sumeef(Psi); load _occasion_ &A &LSM pA_d ipcw _durat_ _interval_ _last_ %if %bquote(&eligib)^= %then &eligib; %if %eval(&itt) = 1 %then _durat2_ ; %if %bquote(&weight)^= %then weight ; last_for_rr; rows = nrow(_occasion_); newy = j(rows,1,0); holder = 0; numberpeople = 0; Smat = J(&individuals,&numpsi,0); %if %eval(&itt) = 1 %then %do; Hstep = _durat2_#exp(&A*Psi[1]); %end; %else %if %eval(&itt)^= 1 %then %do; model = _interval_#(exp(&A#(Psi[1] %if %eval(&numpsi) = 2 %then + Psi[2]*&LSM ;))); Hstep = J(rows,1,0);; %end; modelinc = exp(Psi[1] %if %eval(&numpsi) = 2 %then + Psi[2]*&LSM ; ); gamma = gammaf(Psi); yhold = _durat_ * gamma; yhold_g = yhold - gamma; yhold0 = yhold_g + 1; yhold1 = yhold_g + modelinc; numberpeople = 0; do i = rows to 1 by -1; if _last_[i] then do; holder = 0; delta = 1; numberpeople = numberpeople + 1; contrib = J(_occasion_[i]+1,&numpsi,0); cens_eof = 0; %if %eval(&rr)= 1 %then if _interval_[i] = -1 then cens_eof = 1;; end; if cens_eof then do; delta = 0; step = 0; end; else do; %if %eval(&itt)^= 1 %then %do; Hstep[i] = holder + model[i]; holder = Hstep[i]; %end; Y = Hstep[i] <= yhold[i]; delta = delta*Y; %if %sysevalf(&hfunc) = 0 %then step = 1; %else %if %sysevalf(&hfunc) = 1 %then step = Hstep[i]; %else step = log(Hstep[i]);; end; newY[i] = delta * step ; user_wt = 1; %if %bquote(&weight)^= %then user_wt = weight[i] ;; EA = pA_d[i]; %if %bquote(&eligib)^= %then if &eligib[i] then do;; filacontrib = _occasion_[i]+1; if pA_d[i]>&minpA & pA_d[i]<&maxpA then contrib[filacontrib,1] = user_wt * ipcw[i] * newY[i]* (&A[i] - EA) ; else contrib[filacontrib,1] = 0; %if %sysevalf(&numpsi) = 2 %then %do; if pA_d[i]>&minpA & pA_d[i]<&maxpA then contrib[filacontrib,2] = contrib[filacontrib,1]*&LSM[i]; else contrib[filacontrib,2] = 0; %end; %if %bquote(&eligib)^= %then end; ; if _occasion_[i] = 0 then Smat[numberpeople,] = contrib[+,]; end; Sval = Smat[+,]; %if %sysevalf(&est_method) ^= 3 %then %do; Save = Sval/&individuals; Smat = Smat - repeat(Save,&individuals); Sigma = J(&numpsi,&numpsi,0); /* covariance */ do j = 1 to &numpsi; do k = j to &numpsi; do l = 1 to nrow(Smat); Sigma[j,k] = Sigma[j,k] + Smat[l,j]*Smat[l,k]; end; Sigma[k,j] = Sigma[j,k]; end; end; store sval; if sigma = 0 then do; sigma = 1.0e-16; end; estimeq = (Sval)*Inv(Sigma)*Sval`; %end; %else estimeq = sval[1];; if last_for_rr= 1 then store Hstep Yhold;; return(estimeq); finish sumeef; /* SUMEE Function 2 */ start sumeef2(Psi); estimeq2 = ((sumeef(Psi) - &chilevel )##2); return(estimeq2); finish sumeef2; use _uno_; read all var{_occasion_ &A &LSM pA_d _ipc_ _durat_ _interval_ _last_ %if %bquote(&eligib)^= %then &eligib; %if %eval(&itt) = 1 %then _durat2_ ;}; %if %bquote(&weight)^= %then read all var {&weight} into weight ;; close _uno_; rows = nrow(_occasion_); last_for_rr = 0; Amin = min(&A); Amax = max(&A); %if %eval(&numpsi) = 2 %then %do; Lmin = min(&LSM); Lmax = max(&LSM); %end; store Amin Amax %if %eval(&numpsi) = 2 %then Lmin Lmax; ; /* computing ipc weights */ ipcw = J(rows,1,1); do i = rows to 1 by -1; if _last_[i] then holder=1; ipcw[i] = _ipc_[i]*holder; holder = ipcw[i]; end; store _occasion_ &A &LSM pA_d ipcw _durat_ _interval_ _last_ %if %bquote(&eligib)^= %then &eligib; %if %eval(&itt) = 1 %then _durat2_ ; %if %bquote(&weight)^= %then weight ; last_for_rr; %if %eval(&gridsearch) = 1 %then %do; %if %eval(&print_opt) ^= 0 %then %do; print "Grid search : psi1 : [&psi1Lower , &psi1upper] step size = &psi1step "; %if &numpsi = 2 %then print " psi2 : [&psi2Lower , &psi2upper] step size = &psi2step "; ; %end; ymin = 1e10; create _graph_ var {test psi_1 psi_2}; %let gridsize1 = %sysevalf( (&Psi1upper - &Psi1lower)/&Psi1step,ceil); %let gridsize2 = %sysevalf( (&Psi2upper - &Psi2lower)/&Psi2step,ceil); %if &numpsi = 1 %then %let upperlimit = 0; %else %if &numpsi = 2 %then %let upperlimit = &gridsize2; %do i=0 %to &gridsize1 %by 1; %let j = %sysevalf(&Psi1lower + &i*&Psi1step); %do k=0 %to &upperlimit %by 1; %if &numpsi = 1 %then %let l = 0; %else %let l = %sysevalf(&Psi2lower + &k*&Psi2step); Psi0 = {&j &l}; test = sumeef(Psi0); psi_1= &j; psi_2= &l; append var {test psi_1 psi_2}; %if %eval(&est_method) ^= 3 %then %do; if test < ymin then do; ymin = test; xmin = Psi0; end; %end; %else %do; /* using estimating equation, not test statistic */ if abs(test) < ymin then do; ymin = abs(test); xmin = Psi0; end; %end; %end; %end; close _graph_; objfunc = ymin; psi_1 = xmin[1]; psi_2 = xmin[2]; psi_1_low = 0; psi_1_high = 0; xr = xmin; %end; %if %sysevalf(&est_method) = 1 %then %do; %let outplev = &optn_level; %if %eval(&print_opt) = 1 and &bsample = 0 %then %let outplev = 3; %if %eval(&numpsi) = 1 %then con = { &psi1lower , &psi1upper } ; %else con = { &psi1lower &psi2lower , &psi1upper &psi2upper};; optn = {0 &outplev}; /* minimization print (all: 5, nothing: 0 */ %if %eval(&use_grid_start) = 1 %then %do; Psi0 = psi_1 ; %if %eval(&numpsi) = 2 %then Psi0 = psi_1|| psi_2 ; ; %end; %else %if %eval(&use_grid_start)^= 1 %then %do; Psi0 = {&Psi1init }; %if %eval(&numpsi) = 2 %then Psi0 = {&Psi1init &Psi2init }; ; %end; call nlpnms(rc,xr,"SUMEEF",psi0,optn,con); psi_1= xr[1]; %if %eval(&numpsi) = 2 %then psi_2= xr[2]; %else psi_2 = 0; ; objfunc = sumeef(xr); /* 95% CI for one dimensional psi */ %if %sysevalf(&numpsi) = 1 & %eval(&conf_interval) = 1 %then %do; %if &print_opt = 1 %then print 'Confidence interval'; ; /* first find rough estimate of confidence interval by stepping abs(xr) in both directions */ /* lower estimate */ start for_conf(x); tmp = sumeef(x) -&chilevel; return (tmp); finish for_conf; delta = 1; testlow = objfunc; psilow = xr; countlow = 0; do while ( testlow < &chilevel & countlow < &countmax ); psilow = psilow - delta; testlow = sumeef(psilow); countlow = countlow + 1; end; %if &print_opt = 1 %then print psilow testlow countlow; ; /* upper estimate */ psihigh = xr; testhigh = objfunc; counthigh = 0; do while (testhigh < &chilevel & counthigh < &countmax ); psihigh = psihigh + delta; testhigh = sumeef(psihigh); counthigh = counthigh + 1; end; %if &print_opt = 1 %then print psihigh testhigh counthigh ;; /* better estimate */ if (testhigh > &chilevel) & (testlow > &chilevel) then do; /* bisection method */ left = xr; fleft = objfunc - &chilevel; right = psihigh; fright = testhigh - &chilevel ; middle = (left + right)/2.0; fmiddle = for_conf(middle); count = 0; diff = right - left; do until(abs(fmiddle) < &conf_min | diff < 0.001 | count > &countmax ); test = fmiddle * fleft ; if test < 0 then do; right = middle; fright = fmiddle; end; else do; left = middle; fleft = fmiddle; end; middle = (left + right)/2.0; fmiddle = for_conf(middle); count = count + 1; diff = right - left; end; psi_high = middle; objfunc_high = fmiddle + &chilevel; psi_1_high = psi_high; %if &print_opt = 1 %then print psi_1_high objfunc_high count;; /* lower bound */ left = psilow; fleft = testlow - &chilevel; right = xr; fright = objfunc - &chilevel ; middle = (left + right)/2.0; fmiddle = for_conf(middle); count = 0; diff = right - left ; do until ( abs(fmiddle) < &conf_min | diff < 0.001 | count > &countmax ); test = fmiddle * fleft ; if test < 0 then do; right = middle; fright = fmiddle; end; else do; left = middle; fleft = fmiddle; end; middle = (left + right)/2.0; fmiddle = for_conf(middle); diff = right - left ; count = count + 1; end; psi_low = middle; objfunc_low = fmiddle + &chilevel; psi_1_low = psi_low; %if &print_opt = 1 %then print psi_1_low objfunc_low count ;; %if &print_opt = 1 %then print xr psi_1_low psi_1_high;; end; else do; print "Can not find estimate of confidence interval using less than &countmax function evaluations."; end; %end; %else %do; psi_1_low = 0; psi_1_high = 0; %end; %end; %else %if %sysevalf(&est_method) = 2 %then %do; t0 = 10; lbd = J(&numpsi,1,0); ubd = J(&numpsi,1,0); lbd[1] = &psi1lower; ubd[1] = &psi1upper; %if %eval(&numpsi) = 2 %then %do; lbd[2] = &psi2lower; ubd[2] = &psi2upper; %end; call find_t0("sumeef",&numpsi,&t0_eval,t0,lbd,ubd,xmin_sa,ymin_sa); %if %eval(&use_grid_start) = 1 %then %do; x = xmin; if ymin_sa < ymin then x = xmin_sa; %end; %else %if %eval(&use_grid_start)^= 1 %then x = xmin_sa;; t= anneal("sumeef",&numpsi,&sim_eval,x,t0,lbd,ubd); objfunc = sumeef(x); psi_1 = x[1]; psi_2 = 0; %if %eval(&numpsi) = 2 %then psi_2 = x[2];; %if &print_opt = 1 %then %do; print 'Simulated nnealing : Optimization results '; print " psi1 : [&psi1Lower , &psi1upper] "; %if &numpsi = 2 %then print " psi2 : [&psi2Lower , &psi2upper] "; ; print t0 t xmin_sa ymin_sa psi_1 %if &numpsi = 2 %then psi_2 ; objfunc ; %end; /* 95% CI for one dimensional psi */ %if %sysevalf(&numpsi) = 1 & %eval(&conf_interval) = 1 %then %do; %if &print_opt = 1 %then print 'Confidence interval '; ; lbd_lower = lbd; ubd_upper = ubd; lbd_upper = x[1]; ubd_lower = x[1]; call find_t0("sumeef2",&numpsi,&t0_eval,t0,lbd_lower,lbd_upper,xmin_sa,ymin_sa); xrlow = xmin_sa; t= anneal("sumeef2",&numpsi,&sim_eval,xrlow,t0,lbd_lower,lbd_upper); objfunc_low = sumeef2(xrlow); psi_1_low = xrlow[1]; %if &print_opt = 1 %then print psi_1_low objfunc_low; ; call find_t0("sumeef2",&numpsi,&t0_eval,t0,ubd_lower,ubd_upper,xmin_sa,ymin_sa); xrhigh = xmin_sa; t= anneal("sumeef2",&numpsi,&sim_eval,xrhigh,t0,ubd_lower,ubd_upper); objfunc_high = sumeef2(xrhigh); psi_1_high = xrhigh[1]; %if &print_opt = 1 %then print psi_1_high objfunc_high ;; %end; %else %do; psi_1_low = 0; psi_1_high = 0; %end; %end; %else %if %sysevalf(&est_method) = 3 %then %find_root; %else %if %sysevalf(&est_method) = 4 %then %do; use _graph_; read all var {test psi_1 psi_2} ; %if %eval(&numpsi) = 1 %then %do; for_spline = psi_1 || test ; call splinec(fitted,coeff,endval,for_spline); nf = nrow(fitted); ymin_s = ymin ; do i=1 to nf; yhat = fitted[i,2]; if yhat < ymin_s then do; ymin_s = yhat; xmin_s = fitted[i,1]; end; end; print xmin_s ymin_s; psi_1 = xmin_s; psi_2 = 0; psi = psi_1 // {0}; objfunc = sumeef(psi); /* value of sumeef(psi) lower than min from grid search? note that spline value may be negative */ if objfunc > ymin then do; psi_1 = xmin[1]; objfunc = ymin; end; psi_1_low = 0; psi_1_high = 0; %end; %else %if %eval(&numpsi) = 2 %then %do; x_s = psi_1 || psi_2 ; call tpspline(fitted,coeff,adiag,gcv,x_s,test); store x_s test coeff; psi = {0 0}; start eval_tspline(psi); load x_s coeff; call tpsplnev(pred,psi,x_s,coeff); tmp = pred[1]; return (tmp) ; finish eval_tspline ; optn = {0 1}; con = {&psi1lower &psi2lower, &psi1upper &psi2upper }; Psi0 = xmin; call nlpnms(rc,xr,"eval_tspline",psi0,optn,con); psi_1= xr[1]; psi_2= xr[2]; objfunc = sumeef(xr); if objfunc > ymin then do; print 'No improvement using the Nelder-Mead Simplex method with spline over grid search'; psi_1 = xmin[1]; psi_2 = xmin[2]; objfunc = ymin; end; psi_1_low = 0; psi_1_high = 0; %end; /* numpsi = 2 */ %end; /* est_method = 4 */ create _kk_ var {objfunc psi_1 psi_2 %if &numpsi = 1 and &conf_interval = 1 %then psi_1_low objfunc_low psi_1_high objfunc_high ; }; append var {objfunc psi_1 psi_2 %if &numpsi = 1 and &conf_interval = 1 %then psi_1_low objfunc_low psi_1_high objfunc_high ; }; close _kk_; name = "12345678901234567890"; create _drbeta var{name estimate objfunc %if %eval(&gridsearch) = 1 %then %do; grid_min grid_psi1 %if %eval(&numpsi) = 2 %then grid_psi2; %end; %if &conf_interval = 1 & &numpsi = 1 %then lb1 ub1 objfunc_low objfunc_high ; }; %do iii = 1 %to &n_betas ; %let vartmp = %scan(&cov_betas, %eval(&iii)); name = "&vartmp" ; estimate = psi_&iii; %if %eval(&gridsearch) = 1 %then %do; %do jjj = 1 %to %eval(&numpsi); grid_psi&jjj = xmin[&jjj]; %end; grid_min = ymin ; %end; %if &conf_interval = 1 & &numpsi = 1 %then %do; lb1 = psi_1_low ; ub1 = psi_1_high; %end; append; %end; close _drbeta; /* causal rate ratio */ %if %eval(&rr)= 1 %then %do; last_for_rr= 1; store last_for_rr; psiw = psi_1 // {0}; callkk = sumeef(psiw); load Hstep Yhold; create _weibull_ var {Hstep Yhold ipcw}; append var {Hstep Yhold ipcw}; %end; quit; /* causal rate ratio */ %if %eval(&rr)= 1 %then %do; data _weibull_; set _weibull_; if Hstep> Yhold or Hstep= 0 then do; H= Yhold; censweibull= 1; end; else do; H= Hstep; censweibull= 0; end; run; proc lifereg data= _weibull_ outest= _parms_ %if &print_opt = 0 %then noprint ; ; weight ipcw; model H*censweibull(1)= /d=weibull; run; proc datasets library=work nolist; delete _weibull_ ; run; quit; data _drbeta ; merge _parms_ (keep = _scale_ ) _drbeta; weibull_shape= 1/_scale_; crr= exp(weibull_shape*estimate /* psi_1 */); label crr='Causal Rate Ratio'; drop weibull_shape _scale_ ; %end; %if %eval(&plot_graph)= 1 %then %do; %if %eval(&save_graph) ^= 0 %then filename graphout "mygraph.&ftype";; goptions reset=global %if %eval(&save_graph) ^= 0 %then %do; device = &gdevice gsfname = graphout gsfmode = replace %end; gunit=pct border cback=white colors=(black blue green red) ftitle=swissb ftext=swiss htitle=4 htext=2; symbol1 color=black interpol= spline /* none step needle join r */ value= %if %sysevalf(&gridsize1) < 20 %then dot; %else none;; %if %bquote(&maxCall)= %then %let maxCall_tmp = Random; %else %let maxCall_tmp = &maxCall ; %if %eval(&est_method)^= 3 %then title 'Test statistic as a function of Psi'; %else title 'Estimating equation as a function of Psi';; footnote h=2 j=r "ITT= %sysevalf(&itt), Ign.cens= &igncens., Max EOF= &maxCall_tmp."; %if &numpsi = 1 %then %do; proc gplot data = _graph_; plot test*psi_1/ haxis = &Psi1lower to &Psi1upper %if %sysevalf(&gridsize1) < 20 %then by &Psi1step; hminor = 9 vminor = 4 grid %if %eval(&est_method) ^= 3 %then vref= &chilevel ; %else vref = 0 ; cvref= red; %end; /* numpsi = 1*/ %else %if &numpsi = 2 %then %do; proc g3d data= _graph_; plot psi_2*psi_1= test/ grid tilt= 70 rotate=0 45 90 xticknum=%sysevalf(((&Psi1upper- &Psi1lower)/&Psi1step)+1) yticknum=%sysevalf(((&Psi2upper- &Psi2lower)/&Psi2step)+1) zticknum=6 zmin=0 zmax= &chilevel ctop=green cbottom=red; proc gcontour data= _graph_; plot psi_2*psi_1= test /grid levels= &chilevel; %end; run; %end; /* do_graph = 1 */ data snm; set snm _drbeta; if estimate ne .; if _bsample_ = . then _bsample_ = &bsample ; run; %if &debug = 0 %then %do; proc datasets library= work nolist; delete _uno_ _kk_ _drbeta ; run; quit; %end; title; footnote; %end; /* end of _sloop loop for when bootstrap = 1 and boot_one_a_time = 0 */ %end; /* end of _loop loop */ %mend snmaft; %macro find_root; /* for one dimensional psi */ leftbnd = %sysevalf(&psi1lower); rightbnd = %sysevalf(&psi1upper); lf = sumeef(leftbnd); rf = sumeef(rightbnd); chk = lf * rf; if chk > 0 then do; /* print leftbnd lf rightbnd rf chk;*/ newleftbnd = leftbnd; newrightbnd = newleftbnd + %sysevalf(&psi1step); rf = sumeef(newrightbnd); chk = lf * rf; do while ( chk > 0 & newleftbnd < rightbnd ); newleftbnd = newrightbnd; lf = rf; newrightbnd = newrightbnd + %sysevalf(&psi1step); rf = sumeef(newrightbnd); chk = lf * rf; /* print newleftbnd lf newrightbnd rf chk; */ end; if newrightbnd <= rightbnd & chk < 0 then do; leftbnd = newleftbnd; rightbnd = newrightbnd; delta = rightbnd - leftbnd; end; else print "Bad starting points. Can not find end points having opposite signs "; end; delta = rightbnd - leftbnd; eval = 0; if chk < 0 then do; /* the end points have different signs */ middle = (leftbnd + rightbnd)/2.0; mf = sumeef(middle); do while ( abs(mf) > 0.0001 & delta > 1.e-6 ); delta = rightbnd - leftbnd; middle = (leftbnd + rightbnd)/2.0; mf = sumeef(middle); /* print eval leftbnd rightbnd lf rf mf; */ if lf*mf < 0 then do; rightbnd = middle; rf = mf; end; else do; leftbnd = middle; lf = mf; end; eval = eval + 1; end; psi_1 = middle; objfunc = mf; psi_2 = 0; psi_1_low = 0; psi_1_high = 0; end; /* chk < 0 */ else do; psi_1 = .; objfunc = .; psi_2 = .; covered = .; esteq1 = . ; psi_1_low = .; psi_1_high = .; end; /* could not find any good starting points */ %mend find_root; %macro sim_ann; proc iml; start find_t0(usrfunc,dimen,miter,t0,lbd,ubd,xmin,ymin); sdof = 0.0; range = ubd - lbd; newrho = ρ init = (lbd + ubd)/2.0; x = init; xnew = J(1,dimen,0); call execute("y0 =" + usrfunc + "(x);"); ymin = y0; yold = y0; xmin = x; inc = 0; do i=1 to miter; do j=1 to dimen; t1 = ranuni(300); /* theta = 2.0 * range*theta - range; */ theta = lbd[j] + range[j] * t1; xnew[j] = theta; end; call execute("ynew =" + usrfunc + "( xnew );"); if ynew < ymin then do; ymin = ynew; xmin = xnew; end; if ynew > yold then do; inc = inc + 1; dof = ynew - yold; sdof = sdof + dof; end; x = xnew ; yold = ynew; end; dfave = sdof/inc; t0 = 0.22314*dfave; return; finish find_t0; store module=find_t0; /* iterate a few times at the present temperature to reach thermal equilibrium */ start equilibrate( usrfunc,dimen,t, n,x,xnew,xbest,y,ynew,ybest,lbd,ubd); range = &rrange; newrho = ρ equil = 0; xtmp = J(1,dimen,0); delta = 1.0; do i = 1 to n ; do while (equil < 10); do j = 1 to dimen; theta = ranuni(300); theta = 2.0 * range*theta - range; xc = newrho * t * tan ( theta ); xnew[j] = x[j] + xc; if xnew[j] > ubd[j] then xnew[j] = ubd[j]; if xnew[j] < lbd[j] then xnew[j] = lbd[j]; end; /* "energy" */ call execute("ynew =" + usrfunc + "( xnew );"); if dimen = 1 then do; p1 = xnew[1]; end; else if dimen = 2 then do; p1 = xnew[1]; p2 = xnew[2]; end; c = ynew - y; if c < 0.0 then do; /* keep xnew if energy is reduced */ xtmp = x; x = xnew; xnew = xtmp; y = ynew; if y < ybest then do; do j = 1 to dimen; xbest[j] = x[j]; end; ybest = y; end; delta = abs( c ); if y ^= 0.0 then delta = delta/y; ; else if ynew ^= 0.0 then delta = delta/ynew; ; /* equilibrium is defined as a 10% or smaller change in 10 iterations */ if delta < 0.10 then equil = equil + 1; ; else equil = 0; ; end; else do; equil= equil + 1; end; /* end of c < 0 */ end; if equil > 9 then return (i); end; return (i); finish equilibrate; store module = equilibrate; /* cool the system with annealing */ start anneal(usrModule,dimen, iter,x,t0,lbd,ubd); xtmp = J(1,dimen,0); xnew = J(1,dimen,0); y = J(1,1,0); ynew = J(1,1,0); ybest = J(1,1,0); call execute("y =" + usrModule + "( x );"); xbest = x; ynew = y; ybest = y; newrho = ρ range = &rrange; newtscale = &tscale; n = iter; ntmp = %eval(10*&ddwell); tmp = equilibrate( usrModule,dimen,t0, ntmp ,x,xnew,xbest,y,ynew,ybest,lbd,ubd); told = t0; t = t0; i = 0; do while ( ( i <= n | t > &tmin ) & ybest > &conf_min); i = i + 1; t = t0 /(1.0 + i * newtscale ); told = t; tmp = equilibrate(usrModule,dimen,t, &ddwell,x,xnew,xbest,y,ynew,ybest,lbd,ubd); do j = 1 to dimen; theta1 = ranuni(100); theta = 2.0*range*theta1 - range; xc = (newrho * t * tan(theta)); xnew[j] = x[j] + xc; if xnew[j] > ubd[j] then xnew[j] = ubd[j]; if xnew[j] < lbd[j] then xnew[j] = lbd[j]; end; /* "energy" */ call execute("ynew =" + usrModule + "( xnew );"); if dimen = 1 then p1 = xnew[1]; else if dimen = 2 then do; p1 = xnew[1]; p2 = xnew[2]; end; c = ynew - y; if ynew<=y then do; /* keep xnew if energy is reduced */ xtmp = x; x = xnew; xnew = xtmp; y = ynew; if y theta ) then do; xtmp = x; x = xnew; xnew = xtmp; y = ynew; end; end; end; ntmp = %eval(10 * &ddwell); tmp = equilibrate(usrModule,dimen,t, ntmp ,x,xnew,xbest,y,ynew,ybest,lbd,ubd); x = xbest; return (t); finish anneal; store module= anneal; quit; %mend sim_ann; %macro numargs(arg); %let n = 1; %if &arg^= %then %do; %do %until (%scan(&arg,%eval(&n),%str( ))=%str()); %let word = %scan(&arg,&n); %let n = %eval(&n+1); %end; %end; %eval(&n-1) /* there is no ; here since it will be used as %let a = %numargs(&b) ; and the ; is included at the end of this line */ %mend numargs; %macro output_results; %local for_crr for_crr2 ; %if &get_cov = 1 %then %do; data _cov ; merge _cov name_order; by name; proc sort data = _cov out = _cov (drop = order ); by order; run; data _cov; set _cov; run; %end; data all_results; set &results %if &msm_var = 1 %then varbetas ; %else %if &get_cov = 1 %then _cov;; %if %eval(&msmcox) = 0 %then %do; /* at most two psi for aft */ est = left(estimate); if _n_ = 1 then call symput('psiopt1',est); if _n_ = 2 then call symput('psiopt2',est); %end; run; %if %eval(&gridsearch) = 1 %then %do; %let linegrid0 =Grid search: ; %let forgrid = grid_psi1 8.4; %if %eval(&numpsi) = 2 %then %let forgrid = '(' grid_psi1 8.4 ' , ' grid_psi2 8.4 ')' ; %if %eval(&est_method) = 0 %then %let indent = @32 ; %end; %if %eval(&est_method)^= 0 %then %do; %let indent = @35 ; %if %eval(&est_method) = 1 %then %let method = Nelder-Mead Simplex: ; %else %if %eval(&est_method) = 2 %then %let method = Simulated Annealing: ; %else %if %eval(&est_method) = 4 %then %do; %let method = Spline, using values from grid search: ; %let indent = @52 ; %end; %end; %if %eval(&msmcox) = 1 %then %let firstline = @35; %else %if %eval(&msmcox) = 0 %then %let firstline = @16; %if &msmcox = 0 %then %do; %let line0 = Nested Structural Accelerated Failure Time Model; %let line01 = ; %end; %else %do; %let line0 = Marginal Structural Cox Model; %let line01 = Variance not corrected for weights; %end; %if &n_betas =1 %then %let vartype = v; %else %let vartype = cov ; %if &bootstrap = 0 %then %do; %if &msm_var = 1 %then %do; %let line01 = Analytic &vartype.ariance for stabilized weights; %let line1 = @30 'Standard' @40 ' 95% Confidence' ; %let line2 = @5 'Parameter' @20 'Estimate' @30 ' Error' @40 ' Limits' @60 ' Z' @70 'Pr > |Z|' ; %let line3 = @5 name @20 estimate 8.4 @30 std 8.4 @40 lb 8.4 @50 ub 8.4 @60 Z 8.4 @70 p_value ; %end; %else %do; /* msm_var = 0 */ %if &conf_interval = 1 & &get_cov = 0 %then %do; %if %eval(&rr) = 1 %then %do; %let for_crr = @60 'Causal rate ratio' ; %let for_crr2 = @65 crr 8.4; %end; %let line01 = ; %let line02 = Inverting test to find confidence interval ; %let line1 = @35 ' 95% Confidence' ; %let line2 = @10 'Parameter' @25 'Estimate' @35 ' Limits' /* &for_crr */; %let line3 = @10 name @25 estimate 8.4 @35 lb 8.4 @45 ub 8.4 /* &for_crr2 */ ; %end; /* conf_interval = 1 */ %else %if &get_cov = 1 %then %do; %let line01 = ; %let line1 = @40 'Standard' @50 ' 95% Confidence' ; ; %let line2 = @5 'Parameter' @20 'Estimate' @30 ' Mean' @40 ' Error' @50 ' Limits' @70 ' Z' @80 'Pr > |Z|' ; ; %let line3 = @5 name @20 estimate 8.4 @30 mean 8.4 @40 std 8.4 @50 lb 8.4 @60 ub 8.4 @70 Z 8.4 @80 p_value ; ; %end; %else %do; %if %eval(&rr) = 1 %then %do; %let for_crr = @45 'Causal rate ratio' ; %let for_crr2 = @50 crr 8.4; %end; %let line1 = ; %let line2 = @10 'Parameter' @25 'Estimate' ; %let line3 = @10 name @25 estimate 8.4 ; %end; %end; /* msm_var = 0 */ %end; %else %do; %let line01 = Bootstrap Mean and &vartype.ariance using &nboot bootstrap samples ; %let line1 = @35 'Standard' @45 ' 95% Confidence' ; %let line2 = @1 'Parameter' @15 'Estimate' @25 ' Mean' @35 ' Error' @45 ' Limits' @65 ' Z' @75 'Pr > |Z|' ; %let line3 = @1 name @15 estimate 8.4 @25 mean 8.4 @35 std 8.4 @45 lb 8.4 @55 ub 8.4 @65 Z 8.4 @75 p_value ; %end; data _null_; set all_results; file print notitles; if _n_ = 1 then put &firstline "&line0" // @5 "Number of events : &num_events" / %if %bquote(&maxCall)^= %then @5 "Number of events used : &events_used " / ; ; %if %eval(&msmcox) = 0 %then %do; %if %eval(&est_method) ne 3 %then %do; if _n_ = 1 then put @5 'Minimum value of estimating function' ; %end; %if %eval(&gridsearch) = 1 & %eval(&est_method) ne 3 %then %do; if _n_ = 1 then put @10 'Grid search: ' &indent grid_min 8.4 ' at psi = ' &forgrid ; %end; %if %eval(&est_method)^= 0 %then %do; %if %eval(&est_method)^= 3 %then %do; _psi1 = &psiopt1 ; %if %eval(&numpsi) = 2 %then _psi2 = &psiopt2 ;; if _n_ = 1 then put @10 "&method" &indent objfunc e8.4 ' at psi = ' %if %eval(&numpsi) = 1 %then _psi1 8.4 ; %else %if %eval(&numpsi) = 2 %then '(' _psi1 8.4 ' , ' _psi2 8.4 ')' ; // %if %eval(&conf_interval) = 1 %then @5 "&line02" // ; ; /* end of put */ %end; %else %do; if _n_ = 1 then put @10 'Bisection method yields estimating equation = ' objfunc e8.4 ' at psi = ' estimate 8.4; if _n_ = 1 then put @10 'Confidence intervals can not be calculated using this method.'; %end; %end; %end; %if %bquote(&line01)^= %then if _n_ = 1 then put @5 "&line01" / ;; if _n_ = 1 then put // &line1 / &line2 /; if _n_ <= %eval(&n_betas) then put &line3; %if %eval(&rr) = 1 %then %do; lcrr = log(crr); put // @10 'Log causal rate ratio : ' lcrr 8.4 // ; %end; if _n_ > %eval(&n_betas) then do; %if %eval(&n_betas) < 5 & (&get_cov = 1 | &msm_var = 1) %then %do; if _n_ = %eval(&n_betas+1) then put // @40 "&vartype.ariance matrix " // @10 'Parameter' %let start = 30; %do i = 1 %to &n_betas; %let vartmp = %scan(&cov_betas,%eval(&i)); %let varlength&i = %length(&vartmp); %if %eval(&&varlength&i) < 11 %then %do ; %let spaces = %eval(10 - &&varlength&i); %if %eval(&spaces) > 1 %then %let spaces = %eval(&spaces / 2) ; %let position = %eval(&start + &spaces) ; %end; %else %do; %let position = %eval(&start) ; %let vartmp = %substr(&vartmp,1,9); %end; %let start = %eval(&start + 10) ; @%eval(&position) "&vartmp" %end; // ; if _n_ >= %eval(&n_betas + 1) then put @10 name %let start = 30; %do i = 1 %to &n_betas; %let vartmp = %scan(&cov_betas,%eval(&i)); @%eval(&start) &vartmp 8.4 %let start = %eval(&start + 10); %end;; /* put */ %end; /* n_betas < 5 */ %else %do; if _n_ = %eval(&n_betas + 2) then do; %if &msm_var = 1 %then put // @10 'Covariance matrix is contained in data set varbetas' ;; %if &get_cov =1 & &msm_var = 0 %then put // @10 'Covariance matrix is contained in data set cov_betas' ;; end; %end; end; run; %mend output_results; %macro remove_class; /* remove entries in list2 from list1 */ %local i j test word ; %let i = 1; %let _nonclass= ; %do %until(%scan(&_allcov,&i,%str( )) = %str() ); %let test = %scan(&_allcov,&i); %let j = 1; %let add = 1; %do %until(%scan(&_class,&j,%str( ))=%str()); %let word = %scan(&_class,&j); %if &test = &word %then %let add = 0 ; %let j = %eval(&j + 1); %end; %if %eval(&add) = 1 %then %let _nonclass = &_nonclass &test ;; %let i = %eval(&i + 1); %end; %mend remove_class; %macro convert_class ; %local cnt uniquecnt newlist _nonclass; %let keyfld = %unquote(&class_var) ; /* create list of unique elements of all the class variables */ %let i = 1; %do %until(%scan(&keyfld,&i,%str( ) ) = %str() ); %let var&i = %scan(&keyfld,&i,%str( )); %let i = %eval(&i + 1); %end; %let cnt = %eval(&i - 1); %let newlist = &var1 ; %let in = 1; %let test = 2; %do %until(%eval( &test) > %eval(&cnt)); %let testword = &&var&test; %let add = 1; %do i = 1 %to ∈ %let tmp = %scan(&newlist,&i,%str( )); %if &testword = &tmp %then %do; %let add = 0; %let i = %eval(&in + 1); %end; %end; %if &add = 1 %then %do; %let newlist = &newlist &testword ; %let in = %eval(&in + 1); %end; %let test = %eval(&test + 1); %end; %let uniquecnt = ∈ /* find number of levels for each class variable. store the results in a macro variable */ %do varnum = 1 %to &uniquecnt; %let vartmp = %scan(&newlist,&varnum); proc freq data= _cero_ noprint; tables &vartmp / out = &vartmp; run; data &vartmp; set &vartmp (drop = count percent) end = _end_; by &vartmp; if first.&vartmp then _levels_ + 1; call symput("_&vartmp"||left(put(_levels_,2.)),&vartmp); if _end_ then call symput("_&vartmp"||'levels',_levels_); retain _levels_; run; proc datasets library = work nolist; delete &vartmp; run; quit; %end; /* create the new class variables in the original data set */ data _cero_; set _cero_; %do varnum = 1 %to &uniquecnt; %let vartmp = %scan(&newlist, &varnum); %let new_levels = %eval(&&_&vartmp.levels - 1); %do i = 1 %to &new_levels; _&&vartmp.&i = 0; %let holder = &&_&vartmp.&i; if &vartmp = "&holder" then _&&vartmp.&i = 1; %end; %end; run; /* create new list of covariates with re-coded class variables */ %do i = 1 %to &nclass; %let var = %qscan(&class_var, &i,%str( )); %let cvar = %qscan(&cov_var,&i,%str( )); %let _class = %unquote(&var); %let _allcov = %unquote(&cvar); %if %bquote(&_class)^= %then %do; %remove_class; %let j = 1; %let in = 1; %let tmp = ; %let ctmp= ; %do %until(%scan(%unquote(&var),&j,%str( ) ) = %str() ); %let vartmp = %scan(%unquote(&var),&j,%str( )); %let j = %eval(&j + 1); %let new_levels = %eval(&&_&vartmp.levels - 1); %do k = 1 %to &new_levels; %let tmp = &tmp _&&vartmp.&k ; %end; %end; %let %substr(&cvar,2) = &_nonclass &tmp; %end; %end; %mend convert_class; %macro remove_extraneous; %do i = 1 %to &nclass; %let _var = %qscan(&class_var, &i,%str( )); %let _cvar = %qscan(&cov_var,&i,%str( )); %let _class = %unquote(&_var); %let _allcov = %unquote(&_cvar); %if %bquote(&_class)^= %then %do; /* remove entries in list1 that are not in list2 */ %local test word ; %let ii = 1; %let _class_used= ; %do %until(%scan(&_class,&ii,%str( )) = %str() ); %let test = %scan(&_class,&ii); %let j = 1; %let remove = 1; %do %until(%scan(&_allcov,&j,%str( ))=%str()); %let word = %scan(&_allcov,&j); %if &test = &word %then %let remove = 0 ; %let j = %eval(&j + 1); %end; %if %eval(&remove) = 0 %then %let _class_used = &_class_used &test ;; %let ii = %eval(&ii + 1); %end; %let %substr(&_var,2) = &_class_used; %end; %end; %mend remove_extraneous; /*************************************************************************/