generated from CDCgov/template
-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathemkf_macro.sas
13903 lines (12178 loc) · 594 KB
/
emkf_macro.sas
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
/*
* Version 1.4 10-Aug-2024
*
* eMKF: Expansion of RAND's MKF macro to:
*
* - allow for unequally-spaced time points.
* - allow for quadratic and cubic trends in both MLE-based and Bayesian estimation settings.
* - allow for model averaging over (orthogonal) polynomial trend model sequences up to cubic in both MLE-based and Bayesian settings.
* - allow for common as well as group-specific AR parameters rho and tausq for the random effects in the Bayesian setting.
* - allow for random sampling variances in the Bayesian setting.
* - implement Gibbs sampling in PROC MCMC using user defined samplers (UDSs) that are precompiled with PROC FCMP, replacing .exe file (C code).
* - expand calculations of between-group disparities at the latest time point in the Bayesian setting.
*
* See README.md file for additional details on methodological differences between this eMKF version and RAND's MKF.
*
* Comments and code inserted for the eMKF version in this file are prefixed with *eMKF or otherwise indicated.
* All other comments and code are from the developers of the original MKF macro.
*
* Makram Talih, Ph.D. (NCHS contractor)
*
* Macros defined in this file:
*
* - MKF (main) .............................................................. line 41
* - bayesfit (Bayesian estimation workhorse) ................................ line 3645
* - bayesBMA (Bayesian model averaging workhorse, via mixture prior) ........ line 4925
* - htrp (ML-based estimation workhorse) .................................... line 6205
* - htrp2d (ML-based estimation workhorse for 2 outcomes) ................... line 7697
* - reformat (set up of dataset for analysis) ............................... line 9664
* - _counts_, zeros, thevacompr, etc. (various utility macros) .............. line 10082
* - gibbs_uds_compile_EP (Gibbs sampler for true state predictions) ......... line 10245
* - gibbs_uds_compile_RP (Gibbs sampler for random sampling variances) ...... line 10327
* - gibbs_uds_compile_MP (Gibbs samplers for mean hyper-parameters) ......... line 10364
* - gibbs_uds_compile_CP (Gibbs samplers for regression coefficients) ....... line 10540
* - gibbs_uds_compile_FP (Gibbs samplers for model flags) ................... line 12689
*
*/
data _null_;
run;
/* MKF macro -- eMKF expansion in 2023 Q1-Q2 by M. Talih
Macro defined 03-16-2007, Revised 05-08-09 D. McCaffrey. Revised 03-15-2010
It allows for the estimation of parameters from the non-linear mixed effect model.
In addition it allows for the estimation of model averaging as well as bayesian estimates
1- Printing (modified for eMKF)
modelprint : YES will print the proc nlmixed and/or proc mcmc results and NO will not. Default is NO.
finalprint : YES will print MKF estimates for last time point at the end of the estimations and NO will not. Default is YES.
2- data and variables (format 2 strongly advised for eMKF, as eMKF was not fully tested with format 1)
data : name of the dataset.
It can be a data where each year outcome is listed differently, i.e., Y1-Y8 (format 1)
or it can be an outcome where one year is stacked up on top of the other (format 2) and in that case
the time variable will need to be specified
outcome : the outcome of interest. It can be in two format:
a- multiple outcomes can be specified Y1-Y8 which means the data is in the format (1)
b- a single outcome with data format (2) and the time variable will need to be specified
se : the standard errors of the outcomes. It also needs to be in the a or b format to match the outcome format
neff : (eMKF) (Effective) sample sizes of the outcomes for the random sampling variance case in the Bayesian model
If empty, but random sampling variances are requested, an error is returned.
Ignored in the MLE-based models.
time : the time variable for data format (2). If missing, will assume format (1)
by : allows models to run for multiple strata at the same time and can be used for simulations
group : the different groups (e.g., race/ethnicity group) variable
3- Modeling option: The linear mixed model
slopes : the type of slope one needs: eMKF options are the following:
indep_cubic : The values of the parameters b1, b2, and b3 are computed for each group
indep_quad : b3=0. The values of the parameters b1 and b2 are computed for each group
indep_linear : b3=0 and b2=0. The value of the slope b1 is computed for each group
common_cubic : The values of each of the parameters b1, b2, and b3 are assumed to be the same across groups
common_quad : b3=0. The values of each of the parameters b1 and b2 are assumed to be the same across groups
common_linear : b3=0 and b2=0. The value of the slope b1 is assumed to be the same across groups
dropped : A model without time trend is computed (b3=0, b2=0, b1=0).
(empty) : DEFAULT when using one outcome -- results in Bayesian approach only.
eMKF: the number of time points will be checked and an error returned if:
* Cubic trend is requested when there are less than 5 time points
* Quadratic trend is requested when there are less than 4 time points
* Linear trend is requested when there are less than 3 time points
* There are less than 2 time points
If multiple options are selected, model averaging is computed using those options
* eMKF: When using two outcomes, DEFAULT is model average over all models up to cubic.
* eMKF: Combinations will be checked to ensure a common reference (shared descendent) model for the Bayes factor calculations.
** If both indep_quad and common_cubic are selected, then common_quad will be added if needed.
** If both indep_linear and common_cubic are selected, then common_linear will be added if needed.
** If both indep_linear and common_quad are selected, then common_linear will be added if needed.
4- Modeling option: The Bayesian model
Bayesmodel : This is the Bayesian model one wants to fit. eMKF options are the following:
bma_cubic : (DEFAULT) Bayesian model averaging implemented using a mixture of polynomial trend models up to unconstrained cubic.
bma_quad : Bayesian model averaging implemented using a mixture of polynomial trend models up to unconstrained quadratic.
bma_linear : Bayesian model averaging implemented using a mixture of polynomial trend models up to unconstrained linear.
full_cubic : Fully Bayesian cubic trends. The b1s, b2s, and b3s are independent draws from multilevel priors.
full_quad : Fully Bayesian quadratic trends. Here, b3s are 0, and b1s and b2s are independent draws from multilevel priors.
full_linear : Fully Bayesian linear trends. Here, b3s and b2s are 0, and b1s are indep. draws from a multilevel prior.
indep_cubic : Bayesian cubic trends. The b1s, b2s, and b3s are independent draws from uninformative (flat) priors.
indep_quad : Bayesian quadratic trends. Here, b3s are 0 and b1s and b2s are independent draws from uninformative (flat) priors.
indep_linear : Bayesian linear trends. Here, b3s and b2s are 0, and b1s are independent draws from an uninformative (flat) prior.
common_cubic : Bayesian model with common cubic trend. The common b1, b2, and b3 are each drawn from an uninformative (flat) prior.
common_quad : Bayesian model with common quadratic trend. b3=0 and the common b1 and b2 are each drawn from an uninformative (flat) prior.
common_linear : Bayesian model with common linear trend. b3=0, b2=0, and the common b1 is drawn from an uninformative (flat) prior.
dropped : Bayesian model with no time trend is computed (b3=0, b2=0, b1=0).
BayesmodelAvg : (eMKF) If multiple options are selected in Bayesmodel and BayesmodelAvg = NO or empty,
all selected models are included in the output datasets.
However, only the last model in the specified sequence is used for printing the MKF estimates for last time point.
If BayesmodelAvg = YES (DEFAULT) and more than one model is specified, then model averaging up to largest
unconstrained polynomial trend that was listed is conducted.
ARmodel : (eMKF) Default is common_ar.
If indep_ar, each group will have its own set of AR parameters rho and tausq, with common prior distribution.
randomVars : (eMKF) Default is YES.
Sampling variances will be treated as scaled chi-squared random variables with an inverse gamma prior.
5- Output prefix (modified for eMKF)
out : The name for the output prefix. Default is param.
Outputs (prefix + suffix) are saved. Here are some of the relevant suffixes:
_pred : Kalman prediction of the outcome of interest includes original values as well as parameters
_bayes : Kalman prediction of the outcome from bayesian modeling
(e.g., for OUT=param then PARAM_PRED will be the Kalman prediction data of the outcome of interest.)
xtrakeep: Any variable one wants to keep in the data while running models: weights, ... (eMKF: could be used to retain labels for multiyear data)
comparedata, comparedto : (eMKF) options from MKF allowing for estimating disparities/differences in Bayesian setting
6- Parameters that should only be changed by experienced users (modified for eMKF)
pdigit : Number of decimal digits for the printed outputs. Default is set to 4.
MLE-based model parameters (eMKF: passed on to proc nlmixed)
_rho_ : the value of the true rho that generated the data, if known. If not given, it will be estimated.
_tausq_ : the value of the true tau-square that generated the data, if known. If not given, it will be estimated.
DF : Non-linear model degrees of freedom. The default is set pretty high at 1,0000
Bayes model parameters (eMKF: passed on to proc mcmc)
chains : number of chains to run for the Bayesian estimation. Default is 4, as per Vehtari et al (2021; DOI 10.1214/20-BA1221)
GRthreshold : threshold to use for Gelman-Rubin convergence diagnostic, usually no larger than 1.1. Default is 1.01, as per Vehtari et al.
mcmcplot : if YES, within-chain trace/diagnostics plots from proc mcmc will be included (default is NO)
mcmclog : if YES, the full posterior samples will be retained in the work directory (default is NO, as dataset will be large)
seed : random number generating seed that will allow the user to reproduce the same results in the Bayesian model
maxtune : maximum number of proposal tuning loops (if empty, proc mcmc default of 24 is used; if 0, tuning will be skipped)
ntu : number of tuning iterations to use in each MCMC proposal tuning phase (if empty, proc mcmc default of 500 is used)
nbi : number of burn-in in MCMC iterations (if empty, proc mcmc default of 1000 is used; if 0, burn-in will be skipped)
nmc : number of post-burn-in MCMC iterations (if empty, proc mcmc default of 1000 is used)
thin : controls thinning rate (if empty, proc mcmc default of 1 is used)
accepttol : tolerance for target acceptance probabilities (targetaccept +|- accepttol). If empty, proc mcmc defaults are used.
targetaccept : target acceptance rate for random walk Metropolis. If empty, proc mcmc defaults are used (see proc mcmc documentation).
propcov : method used in constructing initial covariance matrix for the MH algorithm (see proc mcmc documentation).
If empty, proc mcmc default of IND will be used.
init : Option for generating initial values for the parameters. If empty, proc mcmc defaults to MODE (see proc mcmc documentation).
eMKF default is REINIT, which resets model parameters to the user-supplied initial values after tuning.
slicesampler : YES to use the slice sampler instead of MH algorithm for parameters that are not included in the Gibbs sampling steps.
Default is NO due to heavier computational load.
checkSampleSize : YES (default) to check sample size is large enough before proceeding. Set to NO to replicate earlier MKF from RAND.
orpoly : YES (default) for pre-transforming the design matrix using SAS IML orpol function. NO for "raw" polynomials.
If YES, regression coefficients will be reverse-transformed prior to macro end.
However, prior parameters below are assumed to be for the coefficients of the orthogonal polynomial regression if orpoly=YES.
malpha , palpha : prior mean and precision for alphas
mbeta1 : prior mean for mean linear coefficient across groups -- used for cubic, quadratic, or linear trends (full_, indep_, common_)
pbeta1 : prior precision for mean linear coefficient across groups -- used for SD hyperprior(s) in all 3 full_ models
mbeta2 : prior mean for mean quadratic coefficient across groups -- used for cubic or quadratic trends (full_, indep_, common_)
pbeta2 : prior precision for mean quadratic coefficient across groups -- used for SD hyperprior(s) in full_cubic or full_quad
mbeta3 : prior mean for mean cubic coefficient across groups -- used for cubic trends (full_, indep_, and common_)
pbeta3 : prior precision for mean cubic coefficient across groups -- used for SD hyperprior(s) in full_cubic
beta1l, beta1u : bounds for U(a,b) prior for SD of linear coefficients across groups -- used for hyperprior(s) in all 3 full_ models
beta2l, beta2u : bounds for U(a,b) prior for SD of quadratic coefficients across groups -- used for hyperprior(s) in full_cubic or full_quad
beta3l, beta3u : bounds for U(a,b) prior for SD of cubic coefficients across groups -- used for hyperprior(s) in full_cubic
mrho, prho : prior mean and precision for transformed rho -- ie., psi = ln[(1-rho)/(1+rho)]
taul, tauu : bounds for U(a,b) prior for tau (SD of innovation variance tausq)
vshape : Shape parameter for inverse gamma prior distribution of the variance when randomVars = YES
vscale : Scale parameter for inverse gamma prior distribution of the variance when randomVars = YES
wshape : Common shape parameter to use for Dirichlet prior on model indicators in mixture priors (default = 2)
cmploc : Desired location of CMP library (sasuser.funcs or work.funcs if sasuser is write-protected)
*/
%macro mkf(
data = ,
outcome = ,
se = ,
neff = ,
outcome2 = ,
se2 = ,
neff2 = ,
by = ,
group = ,
time = ,
slopes = ,
Bayesmodel = bma_cubic,
BayesmodelAvg = YES,
randomVars = YES,
ARmodel = common_ar,
xtrakeep = ,
out = param,
comparedto = ,
comparedata = ,
modelprint = NO,
finalprint = YES,
/* eMKF: The parameters from here down should not be modified unless the user really knows what s/he is doing*/
pdigit = 4,
_rho_ = ,
_tausq_ = ,
DF = 10000,
chains = 4 ,
GRthreshold = 1.01,
seed = 1235,
maxtune = 50,
ntu = 1000,
nbi = 10000,
nmc = 50000,
thin = 1,
accepttol = ,
targetaccept = ,
propcov = ,
init = reinit,
slicesampler = NO,
checkSampleSize= YES,
orpoly = YES,
/* eMKF: If parameters below are empty, the data will be used to inform the starting values */
malpha = ,
palpha = ,
mbeta1 = 0,
pbeta1 = ,
mbeta2 = 0,
pbeta2 = ,
mbeta3 = 0,
pbeta3 = ,
beta1l = 0,
beta1u = ,
beta2l = 0,
beta2u = ,
beta3l = 0,
beta3u = ,
mrho = 0,
prho = 1,
taul = 0.0001,
tauu = ,
vshape = ,
vscale = ,
wshape = 2,
mcmcplot = NO,
mcmclog = NO,
cmploc = work.funcs
);
%local _oo1_ _oo2_ ui uii uj uk un um uvar newuvar ug uloc
flag1 flag2 flag3 flag4 flag5 flag6 flag7 flag1f flag2f flag3f flag1a flag2a flag3a
crep run1 run2 run3 Bayesian _chainseed
xtrakeep22 toprint toprint2 _ssby _ssn _thekeeps _thekeepsb _thekeep1 _thekeep1b _thekeep2 _thekeep3 _thet
emkfkeep emkfrename _BMAmodel _slopes _rtimess _igrp_ _t_ _comp2;
/* eMKF: Added error check for the length of &out prefix */
%if &out = %str() or %length(&out) > 16 %then %do;
%put ERROR: Please specify a non-empty prefix (out=&out) that is no longer than 16 characters in length;
proc iml;
print " Error Note:";
print " Please specify a non-empty prefix (out=&out) that is no longer than 16 characters in length";
quit;
%return;
%end;
/* eMKF: Added error check for the length of &outcome variable name */
%if &outcome = %str() or %length(&outcome) > 18 %then %do;
%put ERROR: Please specify a non-empty variable name for outcome (&outcome) that is no longer than 18 characters in length;
proc iml;
print " Error Note:";
print " Please specify a non-empty variable name for outcome (&outcome) that is no longer than 18 characters in length";
quit;
%return;
%end;
/* eMKF: Added error check for the length of &outcome2 variable name (if applicable) */
%if &outcome2 ^= %str() and %length(&outcome2) > 18 %then %do;
%put ERROR: Please specify a variable name for outcome2 (&outcome2) that is no longer than 18 characters in length;
proc iml;
print " Error Note:";
print " Please specify a variable name for outcome2 (&outcome2) that is no longer than 18 characters in length";
quit;
%return;
%end;
/* eMKF: Added error check if both slopes and bayesmodel are empty */
%if &Bayesmodel = %str() and &slopes = %str() %then %do;
%put ERROR: There appears to be no estimation method specified: Please check!;
proc iml;
print " Error Note:";
print " There appears to be no estimation method specified. Please check. ";
quit;
%return;
%end;
/*eMKF: Check for any mispellings/misspecification of the model(s) */
%if &slopes ^= %str() %then %do;
%let uvar=; %let ui=0;
%do ui=1 %to %_counts_(&slopes);
%let uvar=;
%let uvar=%scan(&slopes, &ui);
%if %upcase(&uvar) ^=INDEP_CUBIC and %upcase(&uvar) ^=INDEP_QUAD and %upcase(&uvar) ^=INDEP_LINEAR and
%upcase(&uvar) ^=COMMON_CUBIC and %upcase(&uvar) ^=COMMON_QUAD and %upcase(&uvar) ^=COMMON_LINEAR and %upcase(&uvar) ^=DROPPED %then %do;
%put ERROR: &uvar is not a supported ML model specification;
%put ERROR- Model(s) must be one (or more) of indep_cubic, indep_quad, indep_linear, common_cubic, common_quad, common_linear, or dropped;
proc iml;
print " Error Note:";
print " Specified ML model(s) not supported.";
print " Model(s) must be one (or more) of indep_cubic, indep_quad, indep_linear, common_cubic, common_quad, common_linear, or dropped. ";
quit;
%return;
%end;
%end;
%end;
/*eMKF: Bayesian modeling flag and check for any mispellings/misspecification of the model(s) */
%let Bayesian=;
%if &Bayesmodel ^=%str() %then %do;
%let Bayesian= Yes;
%let uvar=; %let ui=0;
%do ui=1 %to %_counts_(&Bayesmodel);
%let uvar=;
%let uvar=%scan(&Bayesmodel, &ui);
%if %upcase(&uvar) ^=BMA_CUBIC and %upcase(&uvar) ^=BMA_QUAD and %upcase(&uvar) ^=BMA_LINEAR and
%upcase(&uvar) ^=FULL_CUBIC and %upcase(&uvar) ^=FULL_QUAD and %upcase(&uvar) ^=FULL_LINEAR and
%upcase(&uvar) ^=INDEP_CUBIC and %upcase(&uvar) ^=INDEP_QUAD and %upcase(&uvar) ^=INDEP_LINEAR and
%upcase(&uvar) ^=COMMON_CUBIC and %upcase(&uvar) ^=COMMON_QUAD and %upcase(&uvar) ^=COMMON_LINEAR and %upcase(&uvar) ^=DROPPED %then %do;
%put ERROR: &uvar is not a supported Bayesian model specification;
%put ERROR- Model(s) must be one (or more) of bma_cubic, bma_quad, bma_linear, full_cubic, full_quad, full_linear, indep_cubic, indep_quad, indep_linear, common_cubic, common_quad, common_linear, or dropped;
proc iml;
print " Error Note:";
print " Specified Bayesian model(s) not supported.";
print " Model(s) must be one (or more) of bma_cubic, bma_quad, bma_linear, full_cubic, full_quad, full_linear, indep_cubic, indep_quad, indep_linear, common_cubic, common_quad, common_linear, or dropped. ";
quit;
%return;
%end;
%end;
%end;
%else %let Bayesian = No; /* eMKF: Added to allow for Bayes approach to be explicitly disabled */
/*eMKF: Set up internal &by variable name _reps, e.g., for use in calls to HTRP and HTRP2D */
%let _ssby=;
%if &by ^=%str() %then %let _ssby=_reps;
/*eMKF: Set up dataset name for differences/disparities if requested */
%if &comparedata = %str() and &comparedto ^= %str() %then %let comparedata = &out._diff;
/*eMKF: Set up internal names for the outcome(s) variable(s) */
%let _oo1_=; %let _oo2_=;
%if &outcome ^=%str() %then %let _oo1_=%scan(&outcome, 1);
%if &outcome2 ^=%str() %then %let _oo2_=%scan(&outcome2, 1);
/*eMKF: Initialize macro variables for use in symbolic calculations/operations */
%let _thet=; %let _thekeeps=; %let _thekeepsb=; %let _thekeep1=; %let _thekeep1b=;
%let _thekeep2= &by &group &time;
%let _thekeep3=;
%if &time = %str() %then %let _thekeep2 = &_thekeep2 _time; /*eMKF: this was introduced in MKF to deal with data in format 1*/
%if %scan(&outcome2,1) =%str() or %scan(&se2,1) =%str() %then %do;
%if %scan(&outcome,2) ^=%str() %then %let _thekeep2 = &_thekeep2 _y _se;
%if %scan(&outcome,2) =%str() %then %let _thekeep2 = &_thekeep2 &outcome &se;
%if &by ^=%str() %then %let _thekeep3= &_thekeep3 &out.(keep= &by);
%let _thekeep3= &_thekeep3 &out.(keep= &group);
%if &time ^=%str() %then %let _thekeep3= &_thekeep3 &out.(keep= &time); /*eMKF: format 2 */
%if &time =%str() %then %let _thekeep3= &_thekeep3 &out.(keep= _time); /*eMKF: format 1 */
%if %scan(&outcome,2) ^=%str() %then %let _thekeep3= &_thekeep3 &out.(keep= _y) &out.(keep= _se);
%if %scan(&outcome,2) =%str() %then %let _thekeep3= &_thekeep3 &out.(keep= &outcome) &out.(keep= &se);
%end;
%if %scan(&outcome2,1) ^=%str() and %scan(&se2,1) ^=%str() %then %do;
%if %scan(&outcome,2) ^=%str() %then %let _thekeep2 = &_thekeep2 _y _se _y2 _se2;
%if %scan(&outcome,2) =%str() %then %let _thekeep2 = &_thekeep2 &outcome &se &outcome2 &se2;
%if &by ^=%str() %then %let _thekeep3= &_thekeep3 &out.(keep= &by);
%let _thekeep3= &_thekeep3 &out.(keep= &group);
%if &time ^=%str() %then %let _thekeep3= &_thekeep3 &out.(keep= &time);
%if &time =%str() %then %let _thekeep3= &_thekeep3 &out.(keep= _time);
%if %scan(&outcome,2) ^=%str() %then %let _thekeep3= &_thekeep3 &out.(keep= _y) &out.(keep= _se) &out.(keep= _y2) &out.(keep= _se2);
%if %scan(&outcome,2) =%str() %then %let _thekeep3= &_thekeep3 &out.(keep= &outcome) &out.(keep= &se) &out.(keep= &outcome2) &out.(keep= &se2);
%end;
/* eMKF: Bayesian estimation not yet implemented for two outcomes--use MLE-based estimation instead */
%if %scan(&outcome2,1) ^=%str() and %scan(&se2,1) ^=%str() %then %do;
%let Bayesian= No;
%let Bayesmodel= ;
/* eMKF: If slopes parameter is unspecified, default to MA up to cubic */
%if &slopes=%str() %then %let slopes= indep_cubic common_cubic indep_quad common_quad indep_linear common_linear dropped;
%end;
/******************************************************/
/* eMKF: Reformat the data only once at the top level */
/******************************************************/
data _bayesdata_ ;
run;
/*eMKF: macro variable &_ssn will be used in call to bayesfit to indicate generic variable name for neff */
%let _ssn = ;
/*eMKF: reformat data if needed and calculate variables needed for processing */
%if %upcase(&Bayesian) = YES and %upcase(&randomVars) = YES %then %do;
%if &neff = %str() %then %do; /*eMKF: if no effective sample sizes were provided, return an error */
%put ERROR: (Effective) sample sizes neff must be specified to fit random sampling variances;
proc iml;
print " Error Note:";
print " (Effective) sample sizes neff must be specified to fit random sampling variances. ";
quit;
%return;
%end;
%reformat(data=&data, outcome=&outcome, se=&se, neff=&neff, outcome2=&outcome2, se2=&se2, neff2=&neff2,
group=&group, time=&time, by=&by, randomVars=YES, outformat= _bayesdata_ );
%let _ssn = _n;
%end;
%else %do;
/*eMKF: Ignore effective samples sizes if random variances not requested in Bayesian setting, or for MLE-based estimation */
%reformat(data=&data, outcome=&outcome, se=&se, outcome2=&outcome2, se2=&se2,
group=&group, time=&time, by=&by, randomVars=NO, outformat= _bayesdata_ );
%let randomVars = NO;
%end;
proc sort data= _bayesdata_;
by &by _group_ _time ;
run;
/*eMKF: Note that _nlmixdata_ (formerly _bayesdata00_) has renamed variables */
/*eMKF: Added columns with new variables instead of re-naming, as original MKF code was re-creating generic variables _y, _time, etc. */
/*eMKF: Create copies of _group_ and _rep variables for use in proc iml matrix calculations */
%if %upcase(&Bayesian) ^= YES %then %do;
data _nlmixdata_;
set _bayesdata_;
_groupnum = _group_;
%if &by ^=%str() %then _reps = _rep;;
run;
%end;
/* eMKF: Datasets to check number of time points and groups are consistent within and across strata */
data _junk_ _freq1_ _freq2_ _freq3_;
run;
/* eMKF: Original RAND macro assumed < 10 groups and timepoints at this point in the code. Here, we do away with this assumption. */
data _junk_;
set _bayesdata_;
if _y=. then delete;
if _se=. then delete;
%if %scan(&outcome2,1) ^=%str() and %scan(&se2,1) ^=%str() %then if _y2=. then delete;;;
%if %scan(&outcome2,1) ^=%str() and %scan(&se2,1) ^=%str() %then if _se2=. then delete;;;
run;
proc freq data=_junk_ noprint;
tables _group_ /list out=_freq1_(rename=(count=ntime));
%if &by ^=%str() %then by &by;;
run;
proc freq data=_junk_ noprint;
tables _time /list out=_freq2_(rename=(count=ngroup));
%if &by ^=%str() %then by &by;;
run;
/*eMKF: added if-clause to avoid creating dataset _finalprint_ if it remains unused */
%if %upcase(&finalprint) = YES %then %do;
options formchar="|----|+|---+=|-/\<>*"; /*eMKF: system options for correct SAS monospace font in printed tables*/
data _finalprint_;
run;
data _finalprint_;
set _freq2_;
_rp_=1;
run;
data _finalprint_;
set _finalprint_;
keep=0;
%if &by ^=%str() %then by &by _time;;
%if &by =%str() %then by _rp_ _time;;
%if &by ^=%str() %then if last.&by and last._time then keep=1;;;
%if &by =%str() %then if last._rp_ and last._time then keep=1;;;
drop _rp_;
run;
data _finalprint_;
set _finalprint_;
if keep=1;
keep &by _time;
run;
%end;
proc sort data=_freq1_ nodupkey;
by &by ntime;
run;
proc sort data=_freq2_ nodupkey;
by &by ngroup;
run;
data _freq1_;
merge _freq1_ _freq2_;
%if &by ^=%str() %then by &by;;
drop _group_ _time percent;
run;
data _freq1_;
set _freq1_;
%if &by ^=%str() %then by &by;;
%if &by ^=%str() %then if first.&by then id =0;;;
stop=0;
id +1;
run;
data _freq2_;
run;
data _freq2_;
set _freq1_;
if id=1;
run;
data _freq1_;
set _freq1_;
stop=1;
if id=2;
run;
proc sort data= _freq1_ nodupkey;
by &by id;
run;
data _freq2_;
merge _freq2_ _freq1_(drop=id);
%if &by ^=%str() %then by &by;;
run;
%let run1=0; %let run2=0;
data _null_;
set _freq2_;
one=1;
if stop = 1 then call symput("run1" , stop);
if stop ne 1 then call symput("run2" , one);
run;
%let run2= %eval(&run2 +0);
%let run1= %eval(&run1 +0);
/* eMKF: Precaution to ensure groups and timepoints are consistent across &by strata */
%if &by ^= %str() and &run1 = 0 %then %do;
%let run3 = 0;
proc freq data=_freq2_ noprint;
tables ngroup /list out=_freq3_;
run;
data _freq3_;
set _freq3_;
stop = 0;
_ng_ + 1;
if _ng_ > 1 then stop = 1;
call symput("run3", stop);
run;
%let run3= %eval(&run3 +0);
%if &run3 = 0 %then %do;
proc freq data=_freq2_ noprint;
tables ntime /list out=_freq3_;
run;
data _freq3_;
set _freq3_;
stop = 0;
_nt_ + 1;
if _nt_ > 1 then stop = 1;
call symput("run3", stop);
run;
%end;
%let run1= %eval(&run3 +0);
%end;
proc datasets nolist;
delete _junk_ _freq1_ _freq2_ _freq3_;
run ;
%if &run1=0 and &run2=0 %then %do;
%put ERROR: There appears to be no data to work with: Please check!;
proc iml;
print " Error Note:";
print " There appears to be no data to work with. Please check. ";
quit;
%return; /*eMKF: added return functionality for easier debugging*/
%end;
%if &run1=1 and &run2=0 %then %do;
%put ERROR: Please check your data: The number of valid timepoints is inconsistent across groups;
proc iml;
print " Error Note:";
print " An error occurred with your data. ";
print " Check the data to make sure there are no missing values for means/SEs or only zero SEs, ";
print " and that all the groups have data for the same number of timepoints (&time). ";
print " You may need to combine some groups and/or timepoints to avoid this problem. ";
quit;
%return; /*eMKF: added return functionality for easier debugging*/
%end;
%if &run1=1 and &run2=1 %then %do;
%put ERROR: Please check your data: The number of valid timepoints is inconsistent across groups and strata;
%put ERROR- This problem occurred in at least one subgroup;
proc iml;
print " Error Note:";
print " An error occurred with your data. ";
print " Check the data to make sure there are no missing values for means/SEs or only zero SEs, ";
print " and that all the groups and strata have data for the same number of timepoints (&time). ";
print " You may need to combine some groups and/or timepoints to avoid this problem. ";
quit;
%return; /*eMKF: added return functionality for easier debugging*/
%end;
%if %scan(&outcome2,1) ^=%str() and %scan(&se2,1) ^=%str() %then %do;
%let ug=0;
data _null_;
set _bayesdata_;
one=1;
if _se2 ne . then call symput('ug',one);
run;
%let ug = %eval(&ug + 0);
%if &ug =0 %then %let run1=1;
%if &ug =0 %then %do;
%put ERROR: Please check your data: You are using two outcomes &outcome and &outcome2;
%put ERROR- Those outcomes or their standard errors should be of the same format: Check that and rerun the model!;
proc iml;
print " Error Note:";
print " An error occurred with your data. ";
print " Check the data to make sure that both outcomes and SEs are in the same format. ";
quit;
%return; /*eMKF: added return functionality for easier debugging*/
%end;
%end;
%if &run1=0 and &run2=1 %then %do;
/* Setting up the number of time points and the number of groups */
%let ug=0; /* eMKF: Number of groups */
data _freqg_;
run;
proc freq data=&data noprint;
tables &group /list out=_freqg_;
run;
data _freqg_;
set _freqg_;
_group_ +1;
call symput('ug',_group_);
keep _group_ &group;
run;
%let ug= %eval(0 + &ug);
%let un=0; /*eMKF: Macro variable for the number of time points */
data _freqn_;
run;
proc freq data=_bayesdata_ noprint;
tables _rtime /list out=_freqn_;
run;
data _freqn_;
set _freqn_;
_time+1;
call symput('un',_time);
keep _rtime _time;
run;
%let un= %eval(0 + &un);
%let _rtimess = ; /*eMKF: Macro variable for the real times to use in calculations */
data _freqn_;
set _freqn_;
retain _rts;
if _n_= 1 then _rts = cat(_rtime);
else _rts = catx(" ", _rts, _rtime);
call symput('_rtimess', _rts);
drop _rts;
run;
proc datasets nolist;
delete _freqn_ _freqg_ ;
run ;
%if &slopes ^=%str() %then %do;
/* eMKF: Guard against too many groups/timepoints: code for current implementation cannot handle more than 204 groups or more than 5508 data points */
%if (&ug > 204 or %eval(&ug*&un) > 5508) %then %do;
%put ERROR: Code for current implementation of this macro cannot handle more than 204 groups or more than 5508 data points (group by time) per stratum;
%put ERROR- Please reduce the number of groups and/or timepoints;
proc iml;
print " Error Note:";
print " Code for current implementation of this macro cannot handle more than 204 groups or more than 5508 data points (group by time) per stratum. ";
print " Please reduce the number of groups and/or timepoints. ";
quit;
%return;
%end;
/* eMKF: Make sure there are enough time points to fit the desired trends (up to cubic) */
%if (&un < 2) %then %do; /*eMKF: minimal check for 2+ points */
%put ERROR: There are not enough timepoints to work with! At least 2 timepoints are required (4+ recommended) to use this macro;
proc iml;
print " Error Note:";
print " There are not enough timepoints to work with! At least 2 timepoints are required (4+ recommended) to use this macro. ";
quit;
%return;
%end;
%if %upcase(&checkSampleSize) = YES %then %do;
%if (&un < 4) %then %do; /*eMKF: added check for 4+ points */
%put ERROR: There are not enough timepoints to work with! At least 4 timepoints are recommended to use this macro;
proc iml;
print " Error Note:";
print " There are not enough timepoints to work with! At least 4 timepoints are recommended to use this macro. ";
quit;
%return;
%end;
%let ui=0;
%do ui=1 %to %_counts_(&slopes);
%let uvar=;
%let uvar=%scan(&slopes, &ui);
%if (%upcase(&uvar) = INDEP_LINEAR or %upcase(&uvar) = COMMON_LINEAR) and (&un < 5) %then %do;
%put ERROR: At least 5 timepoints are recommended to fit linear trends: You may try dropping the time trends instead;
proc iml;
print " Error Note:";
print " At least 5 timepoints are recommended to fit linear trends. You may try dropping the time trends instead. ";
quit;
%return;
%end;
%if (%upcase(&uvar) = INDEP_QUAD or %upcase(&uvar) = COMMON_QUAD) and (&un < 6) %then %do;
%put ERROR: At least 6 timepoints are recommended to fit quadratic trends: You may try linear trends instead;
proc iml;
print " Error Note:";
print " At least 6 timepoints are recommended to fit quadratic trends. You may try linear trends instead. ";
quit;
%return;
%end;
%if (%upcase(&uvar) = INDEP_CUBIC or %upcase(&uvar) = COMMON_CUBIC) and (&un < 7) %then %do;
%put ERROR: At least 7 timepoints are recommended to fit cubic trends: You may try quadratic or linear trends instead;
proc iml;
print " Error Note:";
print " At least 7 timepoints are recommended to fit cubic trends. You may try quadratic or linear trends instead. ";
quit;
%return;
%end;
%end;
%end;
%let uvar =; %let toprint2=;
/* Compute the estimation for all different slopes estimation the user desires*/
%let uj=0;
%let uj= %eval(&uj + %_counts_(&slopes) );
%if &uj=1 %then %do;
%let uvar = &slopes;
/* eMKF: Strings used to account for various combinations of available parameter estimates (up to cubic only) */
%let emkfkeep = ; %let emkfrename = ;
%if %scan(&outcome2,1) =%str() or %scan(&se2,1) =%str() %then %do; /* one outcome */
%let emkfkeep = a; %let emkfrename = a=a_&uvar;
%if %upcase(&uvar) = INDEP_CUBIC or %upcase(&uvar) = COMMON_CUBIC %then %do;
%let emkfkeep = &emkfkeep b1 b2 b3;
%let emkfrename = &emkfrename b1=b1_&uvar b2=b2_&uvar b3=b3_&uvar;
%end;
%if %upcase(&uvar) = INDEP_QUAD or %upcase(&uvar) = COMMON_QUAD %then %do;
%let emkfkeep = &emkfkeep b1 b2;
%let emkfrename = &emkfrename b1=b1_&uvar b2=b2_&uvar;
%end;
%if %upcase(&uvar) = INDEP_LINEAR or %upcase(&uvar) = COMMON_LINEAR %then %do;
%let emkfkeep = &emkfkeep b1;
%let emkfrename = &emkfrename b1=b1_&uvar;
%end;
%end;
%else %do; /* two outcomes */
%let emkfkeep = o1a o2a; %let emkfrename = o1a=o1a_&uvar o2a=o2a_&uvar;
%if %upcase(&uvar) = INDEP_CUBIC or %upcase(&uvar) = COMMON_CUBIC %then %do;
%let emkfkeep = &emkfkeep o1b1 o2b1 o1b2 o2b2 o1b3 o2b3;
%let emkfrename = &emkfrename o1b1=o1b1_&uvar o2b1=o2b1_&uvar o1b2=o1b2_&uvar o2b2=o2b2_&uvar o1b3=o1b3_&uvar o2b3=o2b3_&uvar;
%end;
%if %upcase(&uvar) = INDEP_QUAD or %upcase(&uvar) = COMMON_QUAD %then %do;
%let emkfkeep = &emkfkeep o1b1 o2b1 o1b2 o2b2;
%let emkfrename = &emkfrename o1b1=o1b1_&uvar o2b1=o2b1_&uvar o1b2=o1b2_&uvar o2b2=o2b2_&uvar;
%end;
%if %upcase(&uvar) = INDEP_LINEAR or %upcase(&uvar) = COMMON_LINEAR %then %do;
%let emkfkeep = &emkfkeep o1b1 o2b1;
%let emkfrename = &emkfrename o1b1=o1b1_&uvar o2b1=o2b1_&uvar;
%end;
%end;
%let xtrakeep22=;
%let xtrakeep22= &xtrakeep &outcome &se &group &time impute inputorder &by;
%if &by ^=%str() %then %let xtrakeep22= &xtrakeep22 imputeb;
%let toprint2=&slopes;
/* eMKF: Simplified suffix assignment */
%let newuvar = ;
%if %upcase(&uvar) = INDEP_CUBIC %then %let newuvar = _GC;
%if %upcase(&uvar) = INDEP_QUAD %then %let newuvar = _GQ;
%if %upcase(&uvar) = INDEP_LINEAR %then %let newuvar = _GL;
%if %upcase(&uvar) = COMMON_CUBIC %then %let newuvar = _1C;
%if %upcase(&uvar) = COMMON_QUAD %then %let newuvar = _1Q;
%if %upcase(&uvar) = COMMON_LINEAR %then %let newuvar = _1L;
%if %upcase(&uvar) = DROPPED %then %let newuvar = _0;
%let _thekeep1 = &_thekeep1 pred_&uvar =pred&newuvar predse_&uvar=se&newuvar;
%let _thekeep2 = &_thekeep2 pred_&uvar predse_&uvar;
%let _thekeep3 = &_thekeep3 &out.(keep= pred_&uvar) &out.(keep= predse_&uvar);
%let _thekeeps = &_thekeeps pred&newuvar =&_oo1_._pred&newuvar se&newuvar =&_oo1_._se&newuvar;
/* eMKF: corrected issue in original MKF with re-labeling pred and se using outcome label when there are two outcomes */
%if %scan(&outcome2,1) ^=%str() and %scan(&se2,1) ^=%str() %then %do;
%let _thekeep1 = &_thekeep1 pred2_&uvar =pred2&newuvar pred2se_&uvar=se2&newuvar;
%let _thekeep2 = &_thekeep2 pred2_&uvar pred2se_&uvar;
%let _thekeep3 = &_thekeep3 &out.(keep= pred2_&uvar) &out.(keep= pred2se_&uvar);
%let _thekeeps = &_thekeeps pred2&newuvar =&_oo2_._pred&newuvar se2&newuvar= &_oo2_._se&newuvar;
%end;
%if %scan(&outcome2,1) =%str() or %scan(&se2,1) =%str() %then %do;
%htrp( data=_nlmixdata_, outcome=_y, se=_se, group=_groupnum, time=_time, by=&_ssby,
xtrakeep=&xtrakeep22, orpoly=&orpoly,
_rho_=&_rho_ , _tausq_=&_tausq_ , bvalue= &uvar , DF=&DF, out=&out , print=&modelprint );
data &out._pred;
set &out._pred(rename=(prediction=pred_&uvar predMSE=predMSE_&uvar predSE=predSE_&uvar )
drop=_2loglike _rho_ _tausq_
&emkfkeep /* eMKF: Additional quadratic and cubic terms to drop */
);;
run;
%end;
%else %do;
%htrp2d( data=_nlmixdata_, outcome=_y, se=_se, outcome2=_y2, se2=_se2, group=_groupnum, time=_time, by=&_ssby,
xtrakeep=&xtrakeep22, orpoly=&orpoly,
_rho_=&_rho_ , _tausq_=&_tausq_ , bvalue= &uvar , DF=&DF, out=&out , print=&modelprint );
data &out._pred;
set &out._pred(rename=(prediction=pred_&uvar predMSE=predMSE_&uvar predSE=predSE_&uvar
prediction2=pred2_&uvar predMSE2=pred2MSE_&uvar predSE2=pred2SE_&uvar )
drop=_2loglike _rho_ _tausq_
&emkfkeep /* eMKF: Additional quadratic and cubic terms to drop */
moddelta _tausq2_ _rho2_ err3 _se3 gamma gamma2 gammaeff predeff predeff2
gamma_blup gamma_se pred_blup pred_blup2
);;
run;
%end;
/* eMKF: clean-up. Deleting model-specific results but those could have useful information for the advanced user */
proc datasets nolist;
delete &out._covY &out._data &out._H &out._predvar ;
run ;
quit;
%end; /*End of condition if uj = 1 */
/* If there are more than one type of estimation methods requested then conduct a model averaging estimation */
/* eMKF: 1=indep_cubic, 2=indep_quad, 3=indep_linear, 4=common_cubic, 5=common_quad, 6=common_linear, 7=dropped */
%let flag1 =; %let flag2 =; %let flag3 =; %let flag4=; %let flag5=; %let flag6=; %let flag7=;
%if &uj > 1 %then %do;
/** eMKF: Check that a common reference/descendant model is included, and add to list if not:
** - If both indep_quad & common_cubic are selected, common_quad will be added to the list if not already specified.
** - If both indep_linear & common_cubic are selected, common_linear will be added to the list if not already specified.
** - If both indep_linear & common_quad are selected, common_linear will be added to the list if not already specified.
**/
%let flag1=0; %let flag2=0; %let flag3=0; %let flag4=0; %let flag5=0; %let flag6=0; %let flag7=0;
%let uvar=; %let ui = 0;
%do ui=1 %to %_counts_(&slopes); /* eMKF: first pass to set flags */
%let uvar=;
%let uvar=%scan(&slopes, &ui);
%if %upcase(&uvar) = INDEP_CUBIC %then %let flag1 = 1;
%if %upcase(&uvar) = INDEP_QUAD %then %let flag2 = 1;
%if %upcase(&uvar) = INDEP_LINEAR %then %let flag3 = 1;
%if %upcase(&uvar) = COMMON_CUBIC %then %let flag4 = 1;
%if %upcase(&uvar) = COMMON_QUAD %then %let flag5 = 1;
%if %upcase(&uvar) = COMMON_LINEAR %then %let flag6 = 1;
%if %upcase(&uvar) = DROPPED %then %let flag7 = 1;
%end;
%let uvar=;
%let _slopes = &slopes;
%if &flag5 = 1 and &flag3 = 1 and &flag6 ^= 1 and &flag7 ^= 1 %then %do;
%put ;
%put Both the common_quad and indep_linear models were specified with no shared descendant. Adding the common_linear model...;
%let _slopes = &_slopes common_linear;
%let flag6 = 1;
%put &slopes;
%put &_slopes;
%end;
%if &flag4 = 1 and &flag3 = 1 and &flag6 ^= 1 and &flag7 ^= 1 %then %do;
%put ;
%put Both the common_cubic and indep_linear models were specified with no shared descendant. Adding the common_linear model...;
%let _slopes = &_slopes common_linear;
%let flag6 = 1;
%put &slopes;
%put &_slopes;
%end;
%if &flag4 = 1 and &flag2 = 1 and &flag5 ^= 1 and &flag6 ^= 1 and &flag7 ^= 1 %then %do;
%put ;
%put Both the common_cubic and indep_quad models were specified with no shared descendant. Adding the common_quad model...;
%let _slopes = &_slopes common_quad;
%let flag5 = 1;
%put &slopes;
%put &_slopes;
%end;
/*eMFK: End check for nested models */
%let uvar=ModelAvg;
%let _thekeep1 = &_thekeep1 pred_&uvar =pred_MA predse_&uvar=se_MA;
%let _thekeep2 = &_thekeep2 pred_&uvar predse_&uvar;
%let _thekeep3 = &_thekeep3 &out.(keep= pred_&uvar) &out.(keep= predse_&uvar);
%let _thekeeps = &_thekeeps pred_MA=&_oo1_._pred_MA se_MA=&_oo1_._se_MA;
%if %scan(&outcome2,1) ^=%str() and %scan(&se2,1) ^=%str() %then %do;
%let _thekeep1 = &_thekeep1 pred2_&uvar =pred2_MA pred2se_&uvar=se2_MA;
%let _thekeep2 = &_thekeep2 pred2_&uvar pred2se_&uvar;
%let _thekeep3 = &_thekeep3 &out.(keep= pred2_&uvar) &out.(keep= pred2se_&uvar);
%let _thekeeps = &_thekeeps pred2_MA=&_oo2_._pred_MA se2_MA=&_oo2_._se_MA;
%end;
data &out._pred &out._ests _loglike_;
run;
%let uvar=; %let ui=0;
%do ui=1 %to %_counts_(&_slopes); /*eMKF: slopes replaced with _slopes */
%let uvar=;
%let uvar=%scan(&_slopes, &ui);
/* eMKF: Strings used to account for various combinations of available parameter estimates (up to cubic only) */
%let emkfkeep = ; %let emkfrename = ;
%if %scan(&outcome2,1) =%str() or %scan(&se2,1) =%str() %then %do; /* one outcome */
%let emkfkeep = a; %let emkfrename = a=a_&uvar;
%if %upcase(&uvar) = INDEP_CUBIC or %upcase(&uvar) = COMMON_CUBIC %then %do;
%let emkfkeep = &emkfkeep b1 b2 b3;
%let emkfrename = &emkfrename b1=b1_&uvar b2=b2_&uvar b3=b3_&uvar;
%end;
%if %upcase(&uvar) = INDEP_QUAD or %upcase(&uvar) = COMMON_QUAD %then %do;
%let emkfkeep = &emkfkeep b1 b2;
%let emkfrename = &emkfrename b1=b1_&uvar b2=b2_&uvar;
%end;
%if %upcase(&uvar) = INDEP_LINEAR or %upcase(&uvar) = COMMON_LINEAR %then %do;
%let emkfkeep = &emkfkeep b1;
%let emkfrename = &emkfrename b1=b1_&uvar;
%end;
%end;
%else %do; /* two outcomes */
%let emkfkeep = o1a o2a; %let emkfrename = o1a=o1a_&uvar o2a=o2a_&uvar;
%if %upcase(&uvar) = INDEP_CUBIC or %upcase(&uvar) = COMMON_CUBIC %then %do;
%let emkfkeep = &emkfkeep o1b1 o2b1 o1b2 o2b2 o1b3 o2b3;
%let emkfrename = &emkfrename o1b1=o1b1_&uvar o2b1=o2b1_&uvar o1b2=o1b2_&uvar o2b2=o2b2_&uvar o1b3=o1b3_&uvar o2b3=o2b3_&uvar;
%end;
%if %upcase(&uvar) = INDEP_QUAD or %upcase(&uvar) = COMMON_QUAD %then %do;
%let emkfkeep = &emkfkeep o1b1 o2b1 o1b2 o2b2;
%let emkfrename = &emkfrename o1b1=o1b1_&uvar o2b1=o2b1_&uvar o1b2=o1b2_&uvar o2b2=o2b2_&uvar;
%end;
%if %upcase(&uvar) = INDEP_LINEAR or %upcase(&uvar) = COMMON_LINEAR %then %do;
%let emkfkeep = &emkfkeep o1b1 o2b1;
%let emkfrename = &emkfrename o1b1=o1b1_&uvar o2b1=o2b1_&uvar;
%end;
%end;
/* eMKF: Simplified suffix assignment */
%let newuvar = ;
%if %upcase(&uvar) = INDEP_CUBIC %then %let newuvar = _GC;
%if %upcase(&uvar) = INDEP_QUAD %then %let newuvar = _GQ;
%if %upcase(&uvar) = INDEP_LINEAR %then %let newuvar = _GL;
%if %upcase(&uvar) = COMMON_CUBIC %then %let newuvar = _1C;
%if %upcase(&uvar) = COMMON_QUAD %then %let newuvar = _1Q;
%if %upcase(&uvar) = COMMON_LINEAR %then %let newuvar = _1L;
%if %upcase(&uvar) = DROPPED %then %let newuvar = _0;
%let _thekeep1 = &_thekeep1 pred_&uvar =pred&newuvar predse_&uvar=rmse&newuvar;
%let _thekeep2 = &_thekeep2 pred_&uvar predse_&uvar;
%let _thekeep3 = &_thekeep3 &out.(keep= pred_&uvar) &out.(keep= predse_&uvar);
%let _thekeeps = &_thekeeps pred&newuvar=&_oo1_._pred&newuvar rmse&newuvar=&_oo1_._se&newuvar;
%if %scan(&outcome2,1) ^=%str() and %scan(&se2,1) ^=%str() %then %do;
%let _thekeep1 = &_thekeep1 pred2_&uvar =pred2&newuvar pred2se_&uvar=rmse2&newuvar;
%let _thekeep2 = &_thekeep2 pred2_&uvar pred2se_&uvar;
%let _thekeep3 = &_thekeep3 &out.(keep= pred2_&uvar) &out.(keep= pred2se_&uvar);
%let _thekeeps = &_thekeeps pred2&newuvar=&_oo2_._pred&newuvar rmse2&newuvar=&_oo2_._se&newuvar;
%end;
%if %upcase(&uvar)=INDEP_CUBIC %then %let flag1=YES;
%if %upcase(&uvar)=INDEP_QUAD %then %let flag2=YES;
%if %upcase(&uvar)=INDEP_LINEAR %then %let flag3=YES;
%if %upcase(&uvar)=COMMON_CUBIC %then %let flag4=YES;
%if %upcase(&uvar)=COMMON_QUAD %then %let flag5=YES;
%if %upcase(&uvar)=COMMON_LINEAR %then %let flag6=YES;
%if %upcase(&uvar)=DROPPED %then %let flag7=YES;