-
Notifications
You must be signed in to change notification settings - Fork 0
/
relogitq.ado
654 lines (624 loc) · 29.6 KB
/
relogitq.ado
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
*! Version 1.1.1 October 28, 1999
* (C) Copyright 1999 Michael Tomz, Gary King, and Langche Zeng.
* Current version available at http://gking.harvard.edu
* RElogitQ: Quantities of Interest from Rare Events Logistic Regression
program define relogitq, rclass
version 6.0
if "`e(cmd)'" ~= "relogit" {
di in r _n "You must run relogit before relogitq."
exit 198
}
ck_xc /* is xc_mrt present? */
_strip `"`0'"' /* strip fd wrappers */
local fd `r(fd)' /* save arg in fd(arg) */
local 0 `r(rest)' /* rest of command line */
syntax [, PR RR(str) CHANGEx(str) LISTX MLE /*
*/ UNBIased SIMS(numlist integer max=1 >0) /*
*/ BAYES Level(numlist integer max=1 >0 <100)]
local analyt : word count `bayes'`unbiased'`mle' /* analytical methods? */
if `analyt' > 1 {
di in r _n "Error: Bayes, Unbiased & MLE options are mutually exclusive"
exit 198
}
local p `e(p)' /* proportion of 1's in pop */
gettoken pLO pHI : p /* upper and lower bound */
if `analyt' > 0 & "`pHI'" ~= "" {
di in r _n "Bayes, unbiased and MLE options are not available when " _c
di in r "the fraction of " _n "1's in the population is not known."
exit 198
}
if "`level'" == "" { local level $S_level }/* level for confid interval */
local lo = (100-`level')/2 /* lower bound of confid int */
local hi = `lo' + `level' /* upper bound of confid int */
if "`fd'" ~= "" {
if "`changex'" == "" {
di in r _n "If you request a first difference, you must specify " _c
di in r "the values of" _n "X by using the changex() option"
exit 198
}
ck_rr `fd'
}
else {
if "`changex'" ~= "" {
di in r _n "If you include the changex() option, you must request "
di in r "a first-difference by wrapping one of the options in fd()."
exit 198
}
if "`pr'`rr'" == "" { local pr pr } /* calculate pr by default*/
}
* PARSE CHANGEX AND RELATIVE RISK REQUESTS
_changex `"`changex'"' /* parse changex requests */
local nx `r(nx)' /* # sets of changed x's */
local i 1
while `i' <= `nx' {
local slist`i' `r(slist`i')' /* list of changes: start */
local elist`i' `r(elist`i')' /* list of changes: end */
local i = `i' + 1
}
_changex `"`rr'"' /* parse relrisk requests */
local nrr `r(nx)' /* # sets of rel risks */
local i 1
while `i' <= `nrr' {
local srr`i' `r(slist`i')' /* list of changes: start */
local err`i' `r(elist`i')' /* list of changes: end */
local i = `i' + 1
}
* RUN SETX TO PRODUCE CHANGEX AND RR SCENARIOS
local i 1
while `i' <= `nx' { /* FIRST DIFFERENCES */
if "`pHI'" == "" {
tempname xcS`i' xcE`i' /* xcStart, xcEnd vectors */
setx `slist`i'' $mrt_seto keepmrt /* set starting x's */
matrix `xcS`i'' = r(mrt_xc) /* save starting x's */
setx `elist`i'' $mrt_seto keepmrt /* set ending x's */
matrix `xcE`i'' = r(mrt_xc) /* save ending x's */
}
else {
tempname xcSL`i' xcSH`i' xcEL`i' xcEH`i'
setx `slist`i'' $mrt_seto keepmrt /* set starting x's */
matrix `xcSL`i'' = r(mrt_xcL)
matrix `xcSH`i'' = r(mrt_xcH)
setx `elist`i'' $mrt_seto keepmrt /* set ending x's */
matrix `xcEL`i'' = r(mrt_xcL)
matrix `xcEH`i'' = r(mrt_xcH)
}
local i = `i' + 1
}
local i 1
while `i' <= `nrr' { /* RELATIVE RISKS */
if "`pHI'" == "" {
tempname rrS`i' rrE`i'
setx `srr`i'' $mrt_seto keepmrt
matrix `rrS`i'' = r(mrt_xc)
setx `err`i'' $mrt_seto keepmrt
matrix `rrE`i'' = r(mrt_xc)
}
else {
tempname rrSL`i' rrSH`i' rrEL`i' rrEH`i'
setx `srr`i'' $mrt_seto keepmrt
matrix `rrSL`i'' = r(mrt_xcL)
matrix `rrSH`i'' = r(mrt_xcH)
setx `err`i'' $mrt_seto keepmrt
matrix `rrEL`i'' = r(mrt_xcL)
matrix `rrEH`i'' = r(mrt_xcH)
}
local i = `i' + 1
}
* SIMULATE PARAMETERS FROM MULTIVARIATE NORMAL
tempname b V A row
matrix `b' = e(b) /* parameter estimates */
matrix `V' = e(V) /* var-cov matrix of est*/
local nptosim = rowsof(`V') /* # params to simulate */
tempvar keepme /* mark the original ds */
mark `keepme' /* to be restored later */
if "`sims'" == "" { local sims 1000 } /* 1000 sims by default */
if `sims' > _N { qui set obs `sims' } /* expand ds to fit sims*/
local i 1
while `i' <= `nptosim' {
tempvar c`i'
qui gen `c`i''=invnorm(uniform()) in 1/`sims'
local cnames `cnames' `c`i'' /* collect names of vars*/
local i = `i' + 1
}
_chol `V' `nptosim' `A' /* Cholesky decomp of v */
matrix colnames `A' = `cnames' /* rename cols to */
local i 1
while `i' <= `nptosim' {
matrix `row' = `A'[`i',1...] /* get i^th row of A */
tempvar b`i' /* temporary variable */
local bsims `bsims' `b`i'' /* names of all b's */
matrix score `b`i'' = `row' /* c(NxK) * row(1xK)' */
qui replace `b`i'' = `b`i'' + `b'[1,`i'] /* add mean */
local i = `i' + 1
}
qui drop `cnames'
* CREATE LOW AND HIGH INTERCEPT
if "`pHI'" ~= "" {
summarize `e(depvar)' if e(sample), meanonly
local ybar = r(mean)
tempvar bHI bLO
qui gen `bLO'=`b`nptosim''-ln(1-`pLO')+ln(`pLO')-ln(`ybar')+ln(1-`ybar')
qui gen `bHI'=`b`nptosim''-ln(1-`pHI')+ln(`pHI')-ln(`ybar')+ln(1-`ybar')
local i 1
while `i' < `nptosim' {
local bsimsNC `bsimsNC' `b`i'' /* sims w/o constant */
local i = `i' + 1
}
local bsimsLO `bsimsNC' `bLO' /* sims w/ lo constant */
local bsimsHI `bsimsNC' `bHI' /* sims w/ hi constant */
}
if `nptosim' > `e(df_m)' { local k ,1 } /* append constant to xc? */
if "`listx'" ~= "" { setx } /* list chosen x-values */
* PROBABILITIES
tempname xc xcL xcH AnaXb AnaPr C
tempvar SimXb SimPr SimXbLO SimXbHI SimPrLO SimPrHI
if "`pr'" ~= "" {
if "`pHI'" == "" {
matrix `xc' = (mrt_xc`k') /* xchosen (vector) */
matrix colnames `xc' = `bsims' /* rename cols of xc */
qui matrix score `SimXb' = `xc' /* simulated xbetas */
qui gen `SimPr' = 1/(1+exp(-`SimXb')) /* simulated Pr(Y=1) */
if `analyt' == 0 {
_pctile `SimPr', p(50) /* simulated median */
return scalar Pr = r(r1)
}
else {
matrix `AnaXb' = (mrt_xc`k')*`b'' /* Analytical xbeta */
scalar `AnaPr' = 1/(1+exp(-`AnaXb'[1,1])) /* Analytical Pr(Y=1) */
if "`mle'" == "" {
mat `C' = (`AnaPr'-`AnaPr'^2)*(.5-`AnaPr') /* correction */
mat `C' = `C'*(mrt_xc`k')*`V'*(mrt_xc`k')' /* correction */
if "`bayes'"~="" {scalar `AnaPr'=`AnaPr'+`C'[1,1]} /* bayes */
else { scalar `AnaPr' = `AnaPr' - `C'[1,1] } /* unbias */
_recode `AnaPr' /* keep wi [0,1] bound */
}
return scalar Pr = `AnaPr' /* Analytical mean */
}
qui _pctile `SimPr', p(`lo',`hi') /* calc percentiles */
return scalar PrL = r(r1) /* lower bound of ci */
return scalar PrU = r(r2) /* upper bound of ci */
}
else {
matrix `xcL' = (mrt_xcL`k') /* xchosen|pLO (vector)*/
matrix colnames `xcL' = `bsimsLO' /* rename cols */
qui matrix score `SimXbLO' = `xcL' /* simd xbetas|pLO */
qui gen `SimPrLO' = 1/(1+exp(-`SimXbLO')) /* simd Pr(Y=1|pLO) */
matrix `xcH' = (mrt_xcH`k') /* xchosen|pHI (vector)*/
matrix colnames `xcH' = `bsimsHI' /* rename cols */
qui matrix score `SimXbHI' = `xcH' /* simd xbetas|pHI */
qui gen `SimPrHI' = 1/(1+exp(-`SimXbHI')) /* simd Pr(Y=1|pHI) */
qui _pctile `SimPrLO', p(`lo') /* calc percentiles */
return scalar PrL = r(r1) /* lower bound of ci */
qui _pctile `SimPrHI', p(`hi') /* calc percentiles */
return scalar PrU = r(r1) /* upper bound ci */
return scalar Pr = (`return(PrU)'+`return(PrL)')/2 /*mean of bounds */
}
di _n in g "Probabilities and confidence intervals based on initial X-values"
di " Pr(`e(depvar)'==1) = " %8.5f `return(Pr)'
di " `level'% Confidence Interval = " %8.5f `return(PrL)' " to " `return(PrU)'
}
* FIRST DIFFERENCES
tempname xcS xcE AnaXbS AnaXbE AnaPrS AnaPrE xcSL xcSH xcEL xcEH
tempname xcdiffL xcplus dPrL dPrH xcdiffH
tempvar SimXbS SimXbE diff SimXbSL SimXbSH SimXbEL SimXbEH diffL diffH
tempvar ORLO ARpiLO bsLO nonmon nonmonL nonmonU monL monU ARLLO ARULO
tempvar ORHI ARpiHI bsHI ARLHI ARUHI
local i 1
while `i' <= `nx' {
if "`pHI'" == "" {
matrix `xcS' = (`xcS`i''`k') /* starting values for X */
matrix `xcE' = (`xcE`i''`k') /* ending values for X */
matrix colnames `xcS' = `bsims' /* rename cols of startvals*/
qui matrix score `SimXbS' = `xcS' /* simulated starting Xb's */
matrix colnames `xcE' = `bsims' /* rename cols of endvals */
qui matrix score `SimXbE' = `xcE' /* simulated ending Xb's */
qui gen `diff' = 1/(1+exp(-`SimXbE')) - 1/(1+exp(-`SimXbS'))
if `analyt' == 0 {
_pctile `diff', p(50)
return scalar dPr_`i' = r(r1) /* simulated median 1stdiff*/
}
else {
matrix `AnaXbS' = (`xcS`i''`k')*`b'' /* Analyt XbS */
scalar `AnaPrS' = 1/(1+exp(-`AnaXbS'[1,1])) /* Analyt Pr(Y=1|S)*/
matrix `AnaXbE' = (`xcE`i''`k')*`b'' /* Analyt XbE */
scalar `AnaPrE' = 1/(1+exp(-`AnaXbE'[1,1])) /* Analyt Pr(Y=1|E)*/
if "`mle'" == "" {
mat `C' = (`AnaPrS'-`AnaPrS'^2)*(.5-`AnaPrS') /* correctn|S */
mat `C' = `C'*(`xcS`i''`k')*`V'*(`xcS`i''`k')' /* correctn|S */
if "`bayes'"~=""{scalar `AnaPrS'=`AnaPrS'+`C'[1,1]} /* bayes */
else { scalar `AnaPrS' = `AnaPrS' - `C'[1,1] } /* unbias*/
mat `C'= (`AnaPrE'-`AnaPrE'^2)*(.5-`AnaPrE') /* correctn|E */
mat `C' = `C'*(`xcE`i''`k')*`V'*(`xcE`i''`k')' /* correctn|E */
if "`bayes'"~=""{scalar `AnaPrE'=`AnaPrE'+`C'[1,1]} /* bayes */
else { scalar `AnaPrE' = `AnaPrE' - `C'[1,1] } /* unbias*/
_recode `AnaPrS' `AnaPrE' /* keep in [0,1] */
}
return scalar dPr_`i' = `AnaPrE' - `AnaPrS'
}
qui _pctile `diff', p(`lo',`hi') /* calculate percentiles */
return scalar dPrL_`i' = r(r1) /* lower bound of ci */
return scalar dPrU_`i' = r(r2) /* upper bound of ci */
drop `SimXbS' `SimXbE' `diff'
}
else {
matrix `xcSL' = (`xcSL`i''`k') /* starting values for X|pLO*/
matrix `xcEL' = (`xcEL`i''`k') /* ending values for X|pLO */
matrix `xcSH' = (`xcSH`i''`k') /* starting values for X|pHI*/
matrix `xcEH' = (`xcEH`i''`k') /* ending values for X|pHI */
* Simulate 1st differences under pLO and pHI
matrix colnames `xcSL' = `bsimsLO'
qui matrix score `SimXbSL' = `xcSL' /* sim XB|S,pLO */
matrix colnames `xcSH' = `bsimsHI'
qui matrix score `SimXbSH' = `xcSH' /* sim XB|S,pHI */
matrix colnames `xcEL' = `bsimsLO'
qui matrix score `SimXbEL' = `xcEL' /* sim XB|E,pLO */
matrix colnames `xcEH' = `bsimsHI'
qui matrix score `SimXbEH' = `xcEH' /* sim XB|E,pHI */
qui gen `diffL' = 1/(1+exp(-`SimXbEL')) - 1/(1+exp(-`SimXbSL'))
qui gen `diffH' = 1/(1+exp(-`SimXbEH')) - 1/(1+exp(-`SimXbSH'))
* CASE 1: pLO
* Calculate ARpi assuming pLO
matrix `xcdiffL' = `xcEL' - `xcSL'
matrix colnames `xcdiffL' = `bsims'
qui matrix score `ORLO' = `xcdiffL'
qui replace `ORLO' = exp(`ORLO') /* odds ratio */
qui gen `ARpiLO' = (sqrt(`ORLO')-1)/(sqrt(`ORLO')+1) /* AR_piLO */
* Intercept at which AR is maximized, given pLO
matrix `xcplus' = -.5*(`xcEL' + `xcSL') /* neg avg xcE, xcS */
matrix `xcplus'[1,`nptosim'] = 0 /* omit constant */
matrix colnames `xcplus' = `bsims'
qui matrix score `bsLO' = `xcplus' /* intercept value */
* bounds for AR in monotonic and non-monotinic regions
qui gen `nonmon'=(`bLO'<=`bsLO')&(`bsLO'<=`bHI') /* nonmonotonic region*/
qui egen `nonmonL' = rmin(`diffL' `diffH' `ARpiLO') /* nonmonotonic LO*/
qui egen `nonmonU' = rmax(`diffL' `diffH' `ARpiLO') /* nonmonotonic Hi*/
qui egen `monL' = rmin(`diffL' `diffH') /* monotonic LO */
qui egen `monU' = rmax(`diffL' `diffH') /* monotonic HI */
qui gen `ARLLO' = (`nonmon')*`nonmonL' + (1-`nonmon')*`monL'
qui gen `ARULO' = (`nonmon')*`nonmonU' + (1-`nonmon')*`monU'
* bounds
qui _pctile `ARLLO', p(`lo') /* calc percentiles */
scalar `dPrL' = r(r1) /* lower bound of ci */
qui _pctile `ARULO', p(`hi') /* calc percentiles */
scalar `dPrH' = r(r1) /* upper bound ci */
qui drop `ORLO' `ARpiLO' `bsLO' `nonmon' `nonmonL' `nonmonU' /*
*/ `monL' `monU' `ARLLO' `ARULO'
* CASE2: piHI
* Calculate ARpi assuming pHI
matrix `xcdiffH' = `xcEH' - `xcSH'
matrix colnames `xcdiffH' = `bsims'
qui matrix score `ORHI' = `xcdiffH'
qui replace `ORHI' = exp(`ORHI') /* odds ratio */
qui gen `ARpiHI' = (sqrt(`ORHI')-1)/(sqrt(`ORHI')+1) /* AR_piHI */
* Intercept at which AR is maximized, given pHI
matrix `xcplus' = -.5*(`xcEH' + `xcSH') /* neg avg xcE, xcS */
matrix `xcplus'[1,`nptosim'] = 0 /* omit constant */
matrix colnames `xcplus' = `bsims'
qui matrix score `bsHI' = `xcplus' /* intercept value */
* bounds for AR in monotonic and non-monotinic regions
qui gen `nonmon'=(`bLO'<=`bsHI')&(`bsHI'<=`bHI') /* nonmonotonic region*/
qui egen `nonmonL' = rmin(`diffL' `diffH' `ARpiHI') /* nonmonotonic LO*/
qui egen `nonmonU' = rmax(`diffL' `diffH' `ARpiHI') /* nonmonotonic Hi*/
qui egen `monL' = rmin(`diffL' `diffH') /* monotonic LO */
qui egen `monU' = rmax(`diffL' `diffH') /* monotonic HI */
qui gen `ARLHI' = (`nonmon')*`nonmonL' + (1-`nonmon')*`monL'
qui gen `ARUHI' = (`nonmon')*`nonmonU' + (1-`nonmon')*`monU'
* bounds
qui _pctile `ARLHI', p(`lo') /* calc percentiles */
if r(r1) < `dPrL' { scalar `dPrL' = r(r1) }
qui _pctile `ARUHI', p(`hi') /* calc percentiles */
if r(r1) > `dPrH' { scalar `dPrH' = r(r1) }
qui drop `ORHI' `ARpiHI' `bsHI' `nonmon' `nonmonL' `nonmonU' /*
*/ `monL' `monU' `ARLHI' `ARUHI' `SimXbSL' `SimXbSH' /*
*/ `SimXbEL' `SimXbEH' `diffL' `diffH'
* FINAL BOUNDS
return scalar dPrL_`i' = `dPrL' /* lower bound of ci */
return scalar dPrU_`i' = `dPrH' /* upper bound ci */
return scalar dPr_`i' = (`return(dPrL_`i')'+`return(dPrU_`i')')/2
}
gettoken change changex : changex, parse("&")
di _n in g "First difference `i': `change'"
di " dPr(`e(depvar)'==1) = " %8.5f `return(dPr_`i')'
di " `level'% Confidence Interval = " %8.5f `return(dPrL_`i')' " to " `return(dPrU_`i')'
gettoken ampers changex : changex, parse("&")
local i = `i' + 1
}
* RELATIVE RISKS
tempvar SimXbS SimXbE RR
local i 1
while `i' <= `nrr' {
if "`pHI'" == "" {
matrix `xcS' = (`rrS`i''`k')
matrix `xcE' = (`rrE`i''`k')
matrix colnames `xcS' = `bsims' /* rename cols of startvals*/
qui matrix score `SimXbS' = `xcS' /* simulated starting Xb's */
matrix colnames `xcE' = `bsims' /* rename cols of endvals */
qui matrix score `SimXbE' = `xcE' /* simulated ending Xb's */
qui gen `RR' = (1/(1+exp(-`SimXbE'))) / (1/(1+exp(-`SimXbS')))
if `analyt' == 0 {
_pctile `RR', p(50)
return scalar rr_`i' = r(r1) /* simulated median relrisk*/
}
else {
matrix `AnaXbS' = (`rrS`i''`k')*`b'' /* Analyt XbS */
scalar `AnaPrS' = 1/(1+exp(-`AnaXbS'[1,1])) /* Analyt Pr(Y=1|S)*/
matrix `AnaXbE' = (`rrE`i''`k')*`b'' /* Analyt XbE */
scalar `AnaPrE' = 1/(1+exp(-`AnaXbE'[1,1])) /* Analyt Pr(Y=1|E)*/
if "`mle'" == "" {
mat `C' = (`AnaPrS'-`AnaPrS'^2)*(.5-`AnaPrS') /* correctn|S */
mat `C' = `C'*(`rrS`i''`k')*`V'*(`rrS`i''`k')' /* correctn|S */
if "`bayes'"~=""{scalar `AnaPrS'=`AnaPrS'+`C'[1,1]} /* bayes */
else { scalar `AnaPrS' = `AnaPrS' - `C'[1,1] } /* unbias*/
mat `C'= (`AnaPrE'-`AnaPrE'^2)*(.5-`AnaPrE') /* correctn|E */
mat `C' = `C'*(`rrE`i''`k')*`V'*(`rrE`i''`k')' /* correctn|E */
if "`bayes'"~=""{scalar `AnaPrE'=`AnaPrE'+`C'[1,1]} /* bayes */
else { scalar `AnaPrE' = `AnaPrE' - `C'[1,1] } /* unbias*/
_recode `AnaPrS' `AnaPrE' /* keep in [0,1] */
}
return scalar rr_`i' = `AnaPrE' / `AnaPrS'
}
qui _pctile `RR', p(`lo',`hi') /* calculate percentiles */
return scalar rrL_`i' = r(r1) /* lower bound of ci */
return scalar rrU_`i' = r(r2) /* upper bound of ci */
drop `SimXbS' `SimXbE' `RR'
}
else {
tempvar SimXbSL SimXbSH SimXbEL SimXbEH rrL rrH RRL RRH
matrix `xcSL' = (`rrSL`i''`k')
matrix `xcSH' = (`rrSH`i''`k')
matrix `xcEL' = (`rrEL`i''`k')
matrix `xcEH' = (`rrEH`i''`k')
matrix colnames `xcSL' = `bsimsLO'
qui matrix score `SimXbSL' = `xcSL' /* sim XB|S,pLO */
matrix colnames `xcSH' = `bsimsHI'
qui matrix score `SimXbSH' = `xcSH' /* sim XB|S,pHI */
matrix colnames `xcEL' = `bsimsLO'
qui matrix score `SimXbEL' = `xcEL' /* sim XB|E,pLO */
matrix colnames `xcEH' = `bsimsHI'
qui matrix score `SimXbEH' = `xcEH' /* sim XB|E,pHI */
* simd rel risks pLO and pHI
qui gen `rrL' = (1+exp(-`SimXbSL')) / (1+exp(-`SimXbEL'))
qui gen `rrH' = (1+exp(-`SimXbSH')) / (1+exp(-`SimXbEH'))
qui egen `RRL' = rmin(`rrL' `rrH')
qui egen `RRH' = rmax(`rrL' `rrH')
* final bounds
qui _pctile `RRL', p(`lo') /* calc percentiles */
return scalar rrL_`i' = r(r1) /* lower bound of ci */
qui _pctile `RRH', p(`hi') /* calc percentiles */
return scalar rrU_`i' = r(r1) /* upper bound ci */
return scalar rr_`i' = (`return(rrL_`i')'+`return(rrU_`i')')/2
qui drop `SimXbSL' `SimXbSH' `SimXbEL' `SimXbEH' `rrL' `rrH' /*
*/ `RRL' `RRH'
}
gettoken change rr : rr, parse("&")
di _n in g "Relative Risk `i': `change'"
di " rr(`e(depvar)'==1) = " %8.5f `return(rr_`i')'
di " `level'% Confidence Interval = " %8.5f `return(rrL_`i')' " to " `return(rrU_`i')'
gettoken ampers changex : changex, parse("&")
local i = `i' + 1
}
qui keep if `keepme' == 1
end
****************************** SUB-PROCEDURES ********************************
program define ck_xc
* Checks for presence of xc matrix, simulated parameters, depvar
version 6.0
tokenize "`e(p)'"
if "`2'" == "" {
capture mat list mrt_xc
if _rc == 0 {
local xccols : colnames(mrt_xc)
if "`xccols'" == "`e(rhsvars)'" { local xcfound yes }
}
}
else {
capture mat list mrt_xcL /* does xcL mat exist? */
local rc = _rc
capture mat list mrt_xcH /* does xcH mat exist? */
local rc = `rc' + _rc
if `rc' == 0 { /* if xc does exist... */
local xccolsL : colnames(mrt_xcL) /* names of vars in xcL*/
local xccolsH : colnames(mrt_xcH) /* names of vars in xcH*/
if "`xccolsL'"=="`e(rhsvars)'" & /* if names match ivars
*/ & "`xccolsH'"=="`e(rhsvars)'" /*
*/ { local xcfound yes } /* then xc matrix is OK*/
}
}
if "`xcfound'" ~= "yes" {
di in r _n "X-values missing or invalid. Please rerun -setx-."
exit 198
}
end
program define _strip, rclass
* Strips fd and changex from command line
version 6.0
args cmdline
tokenize `"`cmdline'"', parse("@")
if "`1'" == "@" | "`2'" == "@" {
di in r _n "The -relogitq- command cannot contain the @ character."
exit 198
}
_parsep fd `cmdline' /* strip-off changex(arg) */
return local fd `r(Arg)'
return local rest `r(RestOpt)'
end
program define _parsep, rclass
* Parsing with embedded parentheses
* Based on Jeroen Weesie, "Parsing Options with Embedded Parentheses,"
* Stata Technical Bulletin 40 (Nov 1997), insert ip22.
* Limitation: input can't contain any @'s
version 6.0
if "`*'" == "" {
return local Arg /* pass empty argument */
return local RestOpt /* pass empty remainder */
exit
}
local optname `1' /* pluck-off option name */
mac shift
while "`1'" ~= "" { /* separate all tokens with @ */
local optamp `optamp'`1'@
mac shift
}
local H 0
local Mode0 None
local ProcArg 0
tokenize "`optamp'", parse("@()[]")
while "`1'" ~= "" {
if "`1'" == "`optname'" & `H' == 0 {
if "`2'" == "(" {
local H = `H' + 1
local Mode`H' p
local ProcArg 1
mac shift
}
}
else {
if "`1'" == "(" {
local H = `H' + 1
local Mode`H' p
}
else if "`1'" == "[" {
local H = `H' + 1
local Mode`H' b
}
else if "`1'" == ")" {
if "`Mode`H''" ~= "p" {
di in r "too many ')' or ']'"
exit 132
}
local H = `H' - 1
}
else if "`1'" == "]" {
if "`Mode`H''" ~= "b" {
di in r "too many ')' or ']'"
exit 132
}
local H = `H' - 1
}
else if "`1'" == "@" {
local 1 " "
}
if `ProcArg' == 1 {
if `H' > 0 { local Arg "`Arg'`1'" }
}
else { local RestOpt "`RestOpt'`1'" }
if `H' == 0 {
local ProcArg 0
local Arg "`Arg' "
}
}
mac shift
}
return local Arg `Arg'
return local RestOpt `RestOpt'
end
program define ck_fd
tokenize "`changex'", parse(",")
if "`1'" == "," | "`2'" == "," {
di in r _n "The changex option cannot contain a comma."
exit 198
}
end
program define _changex, rclass
* parses the changex commands
version 6.0
args changex
local nrhsvar : word count `e(rhsvars)' /* # rhs vars */
local nx 0 /* # sets of x */
while "`changex'" ~= "" {
local nx = `nx' + 1 /* next set */
gettoken changes changex : changex, parse("&") /* changes in */
while "`changes'" ~= "" { /* this set */
gettoken vars changes : changes, match(paren) /* variable(s) */
gettoken fun1 changes : changes, match(paren1) /* function 1 */
gettoken fun2 changes : changes, match(paren2) /* function 2 */
while "`vars'" ~= "" {
gettoken var vars : vars /* get varname */
if "`paren1'" == "" {local slist`nx' `slist`nx'' `var' `fun1'}
else {local slist`nx' `slist`nx'' `var' (`fun1')}
if "`paren2'" == "" {local elist`nx' `elist`nx'' `var' `fun2'}
else {local elist`nx' `elist`nx'' `var' (`fun2')}
}
}
return local slist`nx' `slist`nx'' /* return starting list #nx */
return local elist`nx' `elist`nx'' /* return ending list #nx */
gettoken amper changex : changex, parse("&") /* strip-off & */
}
return local nx `nx'
end
program define ck_pr, rclass
version 6.0
args pr prval
tempname catdi allcats
if "`pr'" == "no" { /* pr was not specified */
matrix `catdi' = J(1,2,0)
while "`prval'" ~= "" {
gettoken val prval : prval
if `val' == 0 { matrix `catdi'[1,1] = 1 }
else if `val' == 1 { matrix `catdi'[1,2] = 1 }
else {
di in r _n "Error: `val' is not a valid outcome for `e(depvar)'"
exit 198
}
}
}
else { /* pr was specified */
if "`prval'" ~= "" {
di in g _n "Note: You specified both the pr and the prval() " /*
*/ "options." _n "pr takes precedence over prval(), so -relogitq-"/*
*/ " will display " _n "Pr(Y=0) and Pr(Y=1), rather than the " /*
*/ "values listed in prval()."
}
matrix `catdi' = J(1,2,1)
}
return matrix catdi `catdi'
end
program define _recode
version 6.0
while "`0'" ~= "" {
gettoken p 0 : 0
if `p' < 0 { scalar `p' = 0 }
else if `p' > 1 { scalar `p' = 1 }
}
end
*! version 1.3 April 24, 1999 Michael Tomz
* Cholesky decomposition of an arbitrary matrix
* Input: V, the original matrix
* k, the dimension of the matrix
* A, the name of a new matrix to be created
* Output: A, the lower-triangle cholesky of V
program define _chol
version 6.0
args V k A
capt matrix `A' = cholesky(`V') /* square root of VC matrix */
if _rc == 506 { /* If VC ~ pos definite, ... */
tempname eye transf varianc vcd tmp
mat `eye' = I(`k') /* identity matrix */
mat `transf' = J(1,`k',0) /* initialize transf matrix */
local i 1
while `i' <= `k' {
scalar `varianc' = `V'[`i',`i'] /* variance of parameter `i' */
if `varianc' ~= 0 { /* if has variance, add row */
mat `tmp' = `eye'[`i',1...] /* of `eye' to transf matrix */
mat `transf' = `transf' \ `tmp'
}
local i = `i' + 1
}
mat `transf' = `transf'[2...,1...] /* drop 1st row of transf mat*/
mat `vcd' = `transf'*`V'*`transf'' /* decomposed VC (no 0 rows) */
mat `A' = cholesky(`vcd') /* square root of decomp VC */
mat `A' = `transf''*`A'*`transf' /* rebuild full sq root mat */
}
else if _rc { matrix `A' = cholesky(`V') } /* redisplay error message */
end
program define ck_rr
version 6.0
args fd
local 0 ,`fd'
syntax [, PR RR(string) ]
if "`rr'" ~= "" {
di in r _n "Error: relogitq does not allow first differences of " /*
*/ "relative risks."
exit 198
}
end