-
Notifications
You must be signed in to change notification settings - Fork 2
/
Glimmix_R2_V3.sas
2018 lines (1768 loc) · 61.2 KB
/
Glimmix_R2_V3.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
*******************************************************************
* MACRO: Glimmix_R2
*
* DESCRIPTION: Calculates the R2 statistic for fixed effects
* in the Generalized linear mixed model
*
* REFERENCES: Edwards et al (2008)
* 'An R2 statistic for fixed effects in the linear mixed model.'
* Statistics in Medicine, 27, 6137-6157.
*
* Jaeger et al. (2015)
* 'An R^2 statistic for Fixed Effects in the GLMM Using Penalized Quasi Likelihood'
*
* PROGRAMMERS: Nicolas Ballarini / Byron Casey Jaeger
*
* Affiliation: UNC at Chapel Hill
*
* DATE: 08/01/2016
*
* HISTORY: Mixed_R2 -- Ballarini -- 12/04/2014 --
* Glimmix_R2 -- Ballarini / Jaeger -- 03/20/2015 --
* Glimmix_R2 -- Ballarini / Jaeger -- 06/08/2015 -- Added _ERROR option to run faster if Normal
* Glimmix_R2 -- Jaeger -- 08/01/2016 -- Added Partial R2 for Normal Errors
*
* LANGUAGE: SAS VERSION 9.3
*
* INPUT: SAS data set containing longitudinal data
*
* OUTPUT: R^2 for model and semi-partial R^2 for all
* fixed effects
*******************************************************************;
*******************************************************************
SUBMACRO: Dummy
A SMART GUIDE TO DUMMY VARIABLES: FOUR APPLICATIONS AND A MACRO
Susan Garavaglia and Asha Sharma Dun & Bradstreet
Murray Hill, New Jersey 07974
DESCRIPTION: Creates Dummy variables for ONE categorical variable
NOTE: The output dataset is the same as the input dataset
MACRO PARAMETERS:
dsn = input dataset name ,
var = variable to be categorized ,
prefix = categorical variable prefix ,
*******************************************************************;
%macro Dummy (dsn = , var =, prefix =);
proc summary data = &dsn nway noprint;
class &var ;
output out = x(keep=&var );
*proc print ;
*;
data _null_ ;
set x nobs=totx end=last;
if last then call symput ( 'tot', trim(left(put( totx, best. ) ) ) ) ;
call symput ( 'z' || trim ( left (put ( _n_,best. ) ) ),trim ( left ( &var ) ) ) ;
data _null_ ;
set x nobs=last;
if _n_=1 then call symput ( 'num',trim(left(put( last, best. ) ) ) ) ;
call symput ( 'c' || trim ( left (put ( _n_, best. ) ) ),trim ( left (&var)) );
run ;
%put #
%let list_&var.=;
%do k=2 %to #
%let list_&var.=&&list_&var &prefix&&c&k;
%end;
*%put &&list_&var.;
%let CLASS=&CLASS. &&list_&var.; ;
data &dsn.;
set &dsn. nobs=last;
array ct (&num) %do k=1 %to #
&prefix&&c&k
%end; ;
%do i = 1 %to #
select;
when (&var="&&c&i" ) ct(&i)=1;
otherwise ct(&i)=0; end;
%end; run ;
%mend Dummy ;
*******************************************************************
SUBMACRO: Domino
DESCRIPTION: Creates Dummy variables for categorical variables
using Dummy submacro. The categorical variables are listed in
a macro variable _CATEGORICAL, which is a parameter of the macro MIXED_R2
NOTE: The output dataset is the same as the input dataset
MACRO PARAMETERS:
dsn = input dataset name ,
lib = library of the dataset ,
*******************************************************************;
%macro Domino(dataset=,lib=WORK);
%let dataset=%upcase(&dataset);
%let lib=%upcase(&lib);
%put &_CATEGORICAL;
%global CLASS;
%let CLASS=;
%let j=1;
%do %until (%scan(&_CATEGORICAL,&j)= );
%let SelVar = %scan(&_CATEGORICAL,&j);
%Dummy(dsn=&lib..&dataset., var=&selvar., prefix=&selvar._ ) ;
%put &Selvar;
%let j = %eval(&j+1);
%put CLASS= &CLASS;
%end;
%mend Domino;
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: GLIMMIX */
/* TITLE: Macro for Generalized Linear Mixed Models */
/* PRODUCT: STAT */
/* SYSTEM: ALL */
/* KEYS: GENERALIZED MIXED LINEAR MODELS, */
/* PROCS: MIXED, MACRO */
/* DATA: */
/* */
/* SUPPORT: Russ Wolfinger */
/* REF: */
/* MISC: */
/* */
/****************************************************************/
/*-----------------------------------------------------------
TITLE
-----
GLIMMIX: a SAS macro for fitting generalized linear mixed
models using Proc Mixed and the Output Delivery System (ODS).
Requires SAS/STAT Version 8.
SUPPORT
-------
Russ Wolfinger, SAS Institute Inc. The original version was
written by Jason Brown, formerly of SAS Institute Inc.
Please send email to [email protected] with any suggestions
or corrections.
HISTORY
-------
initial coding 01Jun92 jbb
a few changes and additions 09Oct92 rdw
corrections from Dale McLerran, FHCRC 16Feb94 rdw
suggestions from David Murray, U. Minnesota 21Sep95 rdw
suggestions from Ken Goldberg, Wyeth-Ayerst 27Oct95 rdw
various minor updates 06Apr96 rdw
per suggestions from Ken Goldberg, INITIAL
option changed, INTERCEPT= option dropped,
and FITTING, NOPREV, and NOTEST options
added. 12Mar97 rdw
more Goldberg ideas: NOTES option added,
PARMS specification is only used in the
first iteration unless you also specify
NOPREV, some clean up 19May97 rdw
7.01 conversion 01Jul97 rdw
switched XBETA= and PRED= 14Nov97 rdw
save spatial coordinates as suggested by
Michael O'Kelly, Quintiles Dublin 01Dec97 rdw
eliminated LSMEANS / OM check 05Dec97 rdw
titling code from Dale McLerran, FHCRC 20Feb98 rdw
made PRINTLAST and FITTING the default 25Mar98 rdw
fixed problem with TYPE=SP(EXP) 30Apr98 rdw
allowed METHOD=MIVQUE0 to persist as
suggested by Svetlana Rudnaya, Ford 14Aug98 rdw
changed output data set as suggested
by Carol Gotway-Crawford, CDC 25Sep98 rdw
DESCRIPTION
-----------
The macro uses iteratively reweighted likelihoods to fit the
model; refer to Wolfinger, R. and O'Connell, M., 1993,
``Generalized Linear Mixed Models: A Pseudo-Likelihood
Approach,'' Journal of Statistical Computation and
Simulation, 48.
By default, GLIMMIX uses restricted/residual psuedo
likelihood (REPL) to find the parameter estimates of the
generalized linear mixed model you specify. The macro calls
Proc Mixed iteratively until convergence, which is decided
using the relative deviation of the variance/covariance
parameter estimates. An extra-dispersion scale parameter is
estimated by default.
There are a few macros at the beginning; all are used in the
main macro, GLIMMIX. This macro will work on any type of
model with the error distributions and link functions given
in the ERRLINK macro. In addition, you can specify your
own error and/or link functions. In order to do this, you
must specify error=user and/or link=user in conjunction with
the errvar=, errdev=, linku=, linkud=, linkui=, and linkuid=
options.
The relevant information is saved using the MAKE statement
of Proc Mixed, which is a part of ODS.
The following are reserved variable names and should not be
used in your input SAS data set:
col, deta, dmu, eta, lowereta, lowermu, mu, pred, resraw,
reschi, stderreta, stderrmu, uppereta, uppermu, var, _offset,
_orig, _w, _wght, _y, _z
The following data sets are created by the macro and exist
after completion unless certain options exclude them:
_class, _con, _cov, _diff, _dim, _ds, _est, _fitstats, _lsm,
_model, _pred, _predm, _slice, _soln, _solnr, _tests3
To see how each of these data sets are created, search the macro
code below. If the data set you want is not one of these, add
an appropriate MAKE statement to your STMTS= specification.
CAUTION: This macro can produce biased results for repeated
binary data with few repeats on each subject. Refer to
Breslow and Clayton (1993, JASA, 9-25).
SYNTAX
------
Syntax for the macro is similar to that of Proc Mixed.
There are other options that are macro-specific, however.
%glimmix(data=,
procopt=,
stmts=,
weight=,
freq=,
error=,
errvar=,
errdev=,
link=,
linkn=,
linknd=,
linkni=,
linku=,
linkud=,
linkui=,
linkuid=,
numder=,
cf=,
converge=,
maxit=,
offset=,
out=,
outalpha=,
options=
)
where
data specifies the data set you are using. It can either
be a regular input data set or the _DS data set
from a previous call to GLIMMIX. The latter is used
to specify starting values for GLIMMIX and should be
accompanied by the INITIAL option described below.
procopt specifies options appropriate for a PROC
MIXED statement. Refer to the Proc Mixed
documentation for more information.
stmts specifies Proc Mixed statements for the analysis,
separated by semicolons and listed as a single
argument to the %str() macro function. Statements
may include any of the following: CLASS, MODEL,
RANDOM, REPEATED, PARMS, ID, CONTRAST, ESTIMATE,
and LSMEANS. Syntax and options for each
statement are exactly as in the Proc Mixed
documentation. If you wish to use the OM option
with the LSMEANS statement, you should specify
OM=dataset to avoid conflicts with weights.
weight specifies a weighting variable for the analysis
This allows you to construct your own weights
which can modify or replace the ones constructed
by GLIMMIX.
freq specifies a frequency variable for the analysis.
It replicates observations with the number of
replicates being equal to the value of the FREQ
variable.
error specifies the error distribution. Valid types are:
binomial|b, normal|n, poisson|p, gamma|g,
invgaussian|ig, and user|u
When you specify error=user, you must also provide
the errvar= and errdev= options. The default
error distribution is binomial.
errvar specifies the user-defined variance function. It
must be expressed as a function the argument "mu"
(see examples).
errdev specifies the user-defined deviance function. It
must be expressed as a function the arguments
"_y", which is the response variable, and "mu",
which is the mean. You are allowed to use "_wght"
also, which corresponds to the denominator of a
binomial response. Typical deviance functions are
as follows:
normal (_y-mu)**2
poisson 2*_y*log(_y/mu);
binomial 2*_wght*(_y*log(_y/mu)+
(1-_y)*log((1-_y)/(1-mu)))
gamma -2*log(_y/mu)
invgaussian (((_y-mu)**2)/(_y*mu*mu))
The default deviance is binomial.
link specifies the link function. Valid types are
logit, probit, cloglog, loglog, identity,
power(), log, exp, reciprocal, nlin, and user.
(warning: nlin has not been tested, and it currently
uses an MQL-type estimation scheme.)
When you specify link=nlin, you must also provide
the linkn=, linknd=, and linkni= options. When
you specify link=user, you must also provide the
ulink=, dulink=, and iulink= options. The default
link is different for each error distribution and
is as follows:
Distribution Default Link
------------ ------------
Binomial Logit
Poisson Log
Normal Identity
Gamma Reciprocal
Invgaussian Power(-2)
linkn specifies a nonlinear link function. It must be
enclosed in %str() and assign a value to "mu" by
using parameters "b1" - "bk".
linknd specifies the derivative of the nonlinear link
function.
linkni specifies the initial values for the nonlinear
link function.
linku specifies a user-defined link function. It must
be expressed as a function with the argument "mu".
linkud specifies the derivative of the user-defined link
function with respect to mu. It must be expressed
as a function with argument "mu". For an
approximation, use the formula
(u(mu+h)-u(mu-h))/(2*h)
where u() is the link and h is a small number.
linkui specifies the inverse of the user-defined link.
It must be expressed as a function with argument
"eta".
linkuid specifies the derivative of the inverse of the
user-defined link. It must be expressed as a function
with argument "eta".
numder specifies the tolerance used to numerically differentiate
certain link functions (e.g. probit and power). It has
a default value of 1e-5.
cf specifies the correction factor added to the data
in order to avoid singularities in the initial
iteration. It has a default value of 0.5.
converge sets the convergence criterion for the GLIMMIX
macro. This is not the convergence criteria used
for each internal Proc Mixed call, but rather the
criterion used to assess convergence of the entire
macro algorithm. It has a default value of 1e-8.
maxit specifies the maximum number of iterations for the
GLIMMIX macro to converge. It has a default value of
20.
offset specifies the offset variable. By default no offset
is used.
out specifies a name for an output data set. This data
set is the predicted value data set from Proc Mixed with
the following additional variables:
eta = linear predictor (xbeta) + offset
stderreta = approximate std err of eta
lowereta = lower confidence limit for eta
uppereta = upper confidence limit for eta
mu = inverse link transform of eta
dmu = derivative of mu with respect to eta
stderrmu = approx std err of mu via delta method
lowermu = lower cl for mu, inv link transform of lowereta
uppermu = upper cl for mu, inv link transform of uppereta
var = variance
resraw = raw residual, y - mu
reschi = scaled residual, (y-mu)/sqrt(phi*var)
deta = derivative of eta with respect to mu
_w = weight used in final Proc Mixed call
_z = dependent variable used in final Proc Mixed call
If none is given, then a default name
of _OUTFILE is used.
outalpha specifies an alpha level for the confidence limits
in the out= data set.
options specifies GLIMMIX macro options separated by
spaces:
INITIAL specifes that the input data set is actually
the _DS data set from a previous call to
GLIMMIX. This allows you to restart a
problem that stopped or to specify starting
values.
MQL computes MQL estimates (see Breslow and
Clayton, 1993, JASA, 9-25). The default
is PQL with an extra-dispersion
parameter.
NOPREV prevents use of previous covariance parameter
estimates as starting values for the next
iteration.
NOPRINT suppresses all printing.
NOITPRINT suppresses printing of the iteration
history.
NOTES requests printing of SAS notes, date, and page
numbers during macro execution. By default,
the notes, date, and numbers are turned off
during macro execution and turned back on after
completion.
PRINTALL prints all Proc Mixed runs.
PRINTDATA prints the pseudo data after each
iteration.
OUTPUT
------
The output from this macro is a printout of selected tables
from the final iteration of Proc Mixed. All of these tables
are stored in data sets whose names begin with an
underscore; you can scan the macro code to find the name of
the data set you wish to use.
EXAMPLE SYNTAX
--------------
1) Both of the following examples specifiy the same
analysis: logistic regression, no random effects
%glimmix(data=ingots,
stmts=%str(
class soak;
model nready/ntotal=soak heat;
)
)
%glimmix(data=ingots,
stmts=%str(
class soak;
model nready/ntotal=soak heat;
),
error=user,errvar=mu*(1-mu),
errdev=2*_wght*(_y*log(_y/mu) +
(1-_y)*log((1-_y)/(1-mu))),
link=user,
linku=log(mu/(1-mu)),
linkud=1/(mu*(1-mu)),
linkui=exp(eta)/(1+exp(eta)),
linkuid=-exp(eta)/(1+exp(eta))**2;
)
Here _wght corresponds to ntotal and _y to nready/ntotal.
2) This example uses the random, lsmeans, and options
arguments:
%glimmix(data=salaman1,
stmts=%str(
class fpop fnum mpop mnum;
model y = fpop|mpop;
random fpop*fnum mpop*mnum;
lsmeans fpop|mpop / cl;
)
options=noitprint
)
3) This example uses the procopt, random, and offset
arguments:
%glimmix(data=ship,
procopt=order=data,
stmts=%str(
class type year period;
model y=type;
random year period;
),
error=poisson,link=log,offset=service
)
4) This example uses the repeated argument:
%glimmix(data=salaman1,
stmts=%str(
class fpop fnum mpop mnum;
model y = fpop|mpop;
repeated / type=ar(1) sub=fpop*fnum;
lsmeans fpop|mpop / cl;
)
)
DISCLAIMER
-----------
THIS INFORMATION IS PROVIDED BY SAS INSTITUTE INC. AS A
SERVICE TO ITS USERS. IT IS PROVIDED "AS IS". THERE ARE NO
WARRANTIES, EXPRESSED OR IMPLIED, AS TO MERCHANTABILITY OR
FITNESS FOR A PARTICULAR PURPOSE REGARDING THE ACCURACY OF
THE MATERIALS OR CODE CONTAINED HEREIN.
-------------------------------------------------------------*/
/*--------------------------------------------------------------*/
/* */
/* %mvarlst */
/* Make a variable list from the class list, model */
/* specification, and random specification. */
/* */
/*--------------------------------------------------------------*/
%macro mvarlst;
%let varlst =;
%let mdllst = &mdlspec;
/*---get response variable---*/
%if %index(&response,/) %then
%let varlst = %scan(&response,1,/) %scan(&response,2,/) &varlst;
%else %let varlst = &response &varlst;
/*---get fixed effects---*/
%if %index(&mdllst,@) %then %do;
%let j = 1;
%let mdl = &mdllst;
%let mdllst=;
%do %while(%length(%scan(&mdl,&j,' ')));
%let var=%scan(&mdl,&j,' ');
%if %index(&var,@) %then %do;
%let b = %eval(%index(&var,@)-1);
%let mdllst = &mdllst %substr(%quote(&var),1,&b);
%end;
%else %let mdllst = &mdllst &var;
%let j = %eval(&j+1);
%end;
%end;
%let iv = 1;
%do %while (%length(%scan(&mdllst,&iv)));
%let varlst = &varlst %scan(&mdllst,&iv);
%let iv = %eval(&iv + 1);
%end;
/*---get random effects---*/
%let iv = 1;
%do %while (%length(%scan(&rndlst,&iv)));
%let temp = %scan(&rndlst,&iv);
%if &temp ne INT and &temp ne INTERCEPT %then
%let varlst = &varlst &temp;
%let iv = %eval(&iv + 1);
%end;
/*---get repeated effects---*/
%let iv = 1;
%do %while (%length(%scan(&replst,&iv)));
%let temp = %scan(&replst,&iv);
%if &temp ne DIAG %then %let varlst = &varlst &temp;
%let iv = %eval(&iv + 1);
%end;
%let varlst = &varlst &class &id &freq &weight;
%mend mvarlst;
/*--------------------------------------------------------------*/
/* */
/* %trimlst */
/* Get rid of repetitions in a list */
/* */
/*--------------------------------------------------------------*/
%macro trimlst(name,lst);
%let i1 = 1;
%let tname =;
%do %while (%length(%scan(&lst,&i1,%str( ))));
%let first = %scan(&lst,&i1,%str( ));
%let i2 = %eval(&i1 + 1);
%do %while (%length(%scan(&lst,&i2,%str( ))));
%let next = %scan(&lst,&i2,%str( ));
%if %quote(&first) = %quote(&next) %then %let i2=10000;
%else %let i2 = %eval(&i2 + 1);
%end;
%if (&i2<10000) %then %let tname = &tname &first;
%let i1 = %eval(&i1 + 1);
%end;
%let &name = &tname;
%mend trimlst;
/*--------------------------------------------------------------*/
/* */
/* %remove */
/* Remove a word from a string */
/* */
/*--------------------------------------------------------------*/
%macro remove(sntnc,wrd);
%let sentence=%str( )%nrbquote(&sntnc);
%if &sentence^=%str( ) %then %do;
%let word=%str( )%nrbquote(&wrd);
%let answer=;
%let i=%index(&sentence,&word);
%if &i and &word^=%str( ) %then %do;
%if &i>1 %then %let answer=%qsubstr(&sentence,1,&i-1);
%let j=%eval(&i+%index(%qsubstr(&sentence,&i+1),%str( )));
%if &j>&i %then
%let answer=&answer%qsubstr(&sentence,&j);
%end;
%else %let answer=&sentence;
%unquote(&answer)
%end;
%mend remove;
/*--------------------------------------------------------------*/
/* */
/* %errlink */
/* Create macro variables that contain the link, invlink, */
/* derivative, and the variance funtion for the following */
/* error distributions and link functions: */
/* */
/* error distn: normal, binomial, poisson, gamma, and */
/* inverse gaussian */
/* link func: logit, probit, complementary log-log, */
/* log-log, identity, power(), log, exp, and */
/* reciprocal. */
/* */
/* The user-defined specification is given by leaving the */
/* error distribution field blank and then giving the link, */
/* the derivative of the link, the inverse link, and the */
/* variance function. The parameters for each are: */
/* */
/* mu: variance function, link, and the derivative */
/* of the link */
/* eta: inverse link; */
/* */
/*--------------------------------------------------------------*/
%macro errlink;
/*---error distributions: set variance and deviance functions---*/
%let exiterr = 0;
%if %length(&error)=0 %then %do;
%if %length(&errvar) and %length(&errdev) %then %let error=USER;
%else %let error=BINOMIAL;
%end;
%if %length(&linkn) and %length(&link)=0 %then %let link=NLIN;
%if %length(&linku) and %length(&link)=0 %then %let link=USER;
%if &error=BINOMIAL or &error=B %then %do;
%let errorfn=BINOMIAL;
%let varform=mu*(1-mu);
%let devform=_wght*2*(_y*log(_y/mu) + (1-_y)*log((1-_y)/(1-mu)));
%if %length(&link)=0 %then %let link=LOGIT;
%end;
%else %if &error=POISSON or &error=P %then %do;
%let errorfn=POISSON;
%let varform=mu;
%let devform=_wght*2*_y*log(_y/mu);
%if %length(&link)=0 %then %let link=LOG;
%end;
%else %if &error=NORMAL or &error=N %then %do;
%let errorfn=NORMAL;
%let varform=1;
%let devform=_wght*(_y-mu)**2;
%if %length(&link)=0 %then %let link=IDENTITY;
%end;
%else %if &error=GAMMA or &error=G %then %do;
%let errorfn=GAMMA;
%let varform=mu**2;
%let devform=_wght*-2*log(_y/mu);
%if %length(&link)=0 %then %let link=RECIPROCAL;
%end;
%else %if &error=INVGAUSSIAN or &error=IG %then %do;
%let errorfn=INVGAUSSIAN;
%let varform=mu**3;
%let devform=_wght*(((_y-mu)**2)/(_y*mu*mu));
%if %length(&link)=0 %then %let link=INVGAUSSIAN;
%end;
%else %if &error=USER or &error=U %then %do;
%let errorfn=USER;
%if %length(&errvar) %then %let varform=&errvar;
%else %let exiterr = 1;
%if %length(&errdev) %then %let devform=&errdev;
%else %let exiterr = 1;
%if %length(&link)=0 %then %let link=LOGIT;
%end;
/*---truncate link function, so we can match if a power link---*/
%if %length(&link)>5 %then %let trlink=%substr(&link,1,5);
%else %let trlink=&link;
/*---link functions; set eta, mu, and derivative formulas---*/
%if &trlink=LOGIT %then %do;
%let linkfn=LOGIT;
%let etaform=log(mu/(1-mu));
%let detaform=1/(mu*(1-mu));
%let muform=exp(eta)/(1+exp(eta));
%let dmuform=exp(eta)/(1+exp(eta))**2;
%end;
%else %if &trlink=PROBI %then %do;
%let linkfn=PROBIT;
%let etaform=probit(mu);
%let detaform=(probit(mu+&numder)-probit(mu-&numder))/
(2*&numder);
/*
%let detaform=(probit(min(1,mu+&numder*(1+abs(mu)))) -
probit(max(1,mu-&numder*(1+abs(mu))))) /
(2*&numder*(1+abs(mu)));
*/
%let muform=probnorm(eta);
%let dmuform=(probnorm(eta+&numder)-probnorm(eta-&numder))/
(2*&numder);
%end;
%else %if &trlink=CLOGL %then %do;
%let linkfn=COMPLEMENTARY LOG LOG;
%let etaform=log(-log(1-mu));
%let detaform=-1/((1-mu)*log(1-mu));
%let muform=1-exp(-exp(eta));
%let dmuform=exp(-exp(eta))*exp(eta);
%end;
%else %if &trlink=LOGLO %then %do;
%let linkfn=LOG LOG;
%let etaform=-log(-log(mu));
%let detaform=-1/(mu*log(mu));
%let muform=exp(-exp(-eta));
%let dmuform=exp(-exp(-eta))*exp(-eta);
%end;
%else %if &trlink=IDENT %then %do;
%let linkfn=IDENTITY;
%let etaform=mu;
%let detaform=1;
%let muform=eta;
%let dmuform=1;
%end;
%else %if &trlink=LOG %then %do;
%let linkfn=LOG;
%let etaform=log(mu);
%let detaform=1/mu;
%let muform=exp(eta);
%let dmuform=exp(eta);
%end;
%else %if &trlink=EXP %then %do;
%let linkfn=EXPONENTIAL;
%let etaform=exp(mu);
%let detaform=exp(mu);
%let muform=log(eta);
%let dmuform=1/eta;
%end;
%else %if &trlink=RECIP %then %do;
%let linkfn=INVERSE;
%let etaform=1/mu;
%let detaform=-1/mu**2;
%let muform=1/eta;
%let dmuform=-1/eta**2;
%end;
%else %if &trlink=POWER %then %do;
%let linklen = %eval(%length(&link)-7);
%let expon=%substr(&link,7,&linklen);
%let linkfn=POWER(&expon);
%let etaform=mu**(&expon);
%let detaform=(&expon)*mu**(&expon-1);
/*
%let detaform=((mu+&numder)**(&expon)-(mu-&numder)**(&expon))/
(2*&numder);
%let detaform=((mu+&numder*(1+abs(mu)))**(&expon) -
(mu-&numder*(1+abs(mu)))**(&expon))/
(2*&numder*(1+abs(mu)));
*/
%let muform=eta**(1/(&expon));
%let dmuform=(1/(&expon))*eta**(1/(&expon)-1);
%end;
%else %if &trlink=INVGA %then %do;
%let linkfn=POWER(-2);
%let etaform=mu**(-2);
%let detaform=-2*mu**(-3);
%let muform=eta**(-1/2);
%let dmuform=(-1/2)*eta**(-3/2);
%end;
%else %if &trlink=BOXCO %then %do;
%let linkfn=BOX-COX;
%let linklen = %eval(%length(&link)-8);
%let expon=%substr(&link,8,&linklen);
%let etaform=(mu**(&expon)-1)/(&expon);
%let detaform=mu**((&expon)-1);
%let muform=((&expon)*eta + 1)**(1/(&expon));
%let dmuform=((&expon)*eta + 1)**(1/(&expon)-1);
%end;
%else %if &trlink=USER %then %do;
%let linkfn=USER;
%if %length(&linku) %then %let etaform=&linku;
%else %let exiterr = 1;
%if %length(&linkud) %then %let detaform=&linkud;
%else %let exiterr = 1;
%if %length(&linkui) %then %let muform=&linkui;
%else %let exiterr = 1;
%if %length(&linkuid) %then %let dmuform=&linkuid;
%else %let exiterr = 1;
%end;
%else %if &trlink=NLIN %then %do;
%let linkfn=NONLINEAR;
%if %length(&linkn) %then %let nlinform=&linkn;
%else %let exiterr = 1;
%if %length(&linknd) %then %let nlinder=&linknd;
%else %let exiterr = 1;
%end;
%if %index(&options,DEBUG) %then %do;
%put options = &options;
%put intopt = &intopt;
%put varlst = &varlst;
%put error = &errorfn;
%put variance = &varform;
%put deviance = &devform;
%put link: eta = &etaform;
%put dlink: deta = &detaform;
%put invlink: mu = &muform;
%end;
%mend errlink;
/*--------------------------------------------------------------*/
/* */
/* %init */
/* Sets the initial values for the iterations. */
/* */
/*--------------------------------------------------------------*/
%macro init;
%let off = &offset;
%if %index(&intopt,NLIN) %then %do;
/*---determine number of parameters---*/
%let nb = 0;
%let i = 1;
%do %while(%index(&linkni,B&i));
%let nb = %eval(&nb + 1);
%let i = %eval(&i + 1);
%end;
%let nu = 0;
%let ns = 0;
%let nus = 0;
%let varlst = &varlst one mu;
%end;
data _ds;
set %unquote(&data);
%if not %index(&options,INITIAL) %then %do;
/*---move away from parameter space boundary for the binomial
error situation---*/
%if %index(&response,/) %then %do;
mu = (%scan(&response,1,/) + &cf)/(%scan(&response,2,/) +
2*&cf);
_wght = %scan(&response,2,/) ;
%end;
%else %if &errorfn=BINOMIAL %then %do;
mu = (&response + &cf)/(1 + 2*&cf);
_wght = 1;
%end;
%else %do;
mu = &response + &cf;
_wght = 1;
%end;
%if %length(&weight) %then %do;
_wght = &weight * _wght;
%end;
_y = &response;
var = &varform;
_offset = &off;
%if %index(&intopt,NLIN) %then %do;
array b{&nb} b1-b&nb;
array db{&nb} db1-db&nb;
one = 1;
%do i = 1 %to &nb;
%let idx = %index(&linkni,B&i);
b&i = %scan(%substr(&linkni,&idx),2,'=' ' ');
%end;
&nlinform
_z = _y - mu;
&nlinder
do i = 1 to &nb;
_z = _z + db{i}*b{i};
end;
_w = _wght / var;
%end;
%else %do;
eta = &etaform ;
deta = &detaform ;
_w = _wght / ((deta**2)*(var));
_z = (_y-mu)*deta + eta - _offset;
%end;
%if %length(&freq) %then %do;
do i = 1 to &freq;
if i=1 then _orig='y';
else _orig='n';
output;
end;
%end;
%else %do;
_orig='y';
%end;
if (_w = .) then _w = 1;
/*
keep _y _z _w _offset _wght _orig &varlst;
*/
%end;
run;
%if %index(&options,PRINTDATA) %then %do;
proc print;
run;
%end;
%let iter = 0;
%mend init;
/*--------------------------------------------------------------*/
/* */
/* %newdata */
/* Create the new data set with the updated values */
/* */
/*--------------------------------------------------------------*/
%macro newdata;
/*---save previous parameter estimates---*/
data _oldsoln;
set _soln;
run;
/*---save previous estimates of covariance matrix---*/
data _oldcov;
set _cov;
%let covsaved = 1;
run;
%if %index(&options,DEBUG) %then %put Creating new pseudo data.;
/*---create new data set---*/
data _ds;
%if %index(&intopt,NLIN) %then %do;
set _ds;
array b{&nb} b1-b&nb;
array db{&nb} db1-db&nb;
/*
array u{&nu,&ns} u1-u&nus;
array du{&nu,&ns} du1-du&nus;
array _r{&nu} _r1-_rν
*/
&nlinform
_z = _y - mu;
&nlinder
do i = 1 to &nb;
_z = _z + db{i}*b{i};
end;
/*
do i = 1 to ν
_r{i} = du{i,subject};
_z = _z + _r{i}*u{i,subject};
end;
end;