[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

From |
Kathryn.A.Sabadosa@Dartmouth.EDU (Kathryn A. Sabadosa) |

To |
statalist@hsphsun2.harvard.edu |

Subject |
st: Stata Macro Issue |

Date |
14 Oct 2002 15:28:13 EDT |

Greetings, I am writing to ask if you would take a look at the enclosed do file and stata macro. I am trying to show an adjusted survival. The macro, adjbobo.ado, essentially generates 6 survival curves, the crude survival curves, the mean survival curves and the adjusted survival curves, however, the graph is not as a typical sts plot. If I am reading the macro code correctly, the y axis is set by the data. I would like to control the graphic display. I have tried to insert code to set the max and min, however, the macro failed. To run the macro I have used the following steps: run a do file to rename variables type adjbobo Revproc Sex at the command line Any help you could give me, would be greatly appreciated. If you would like a copy of the data set that I am working with to see the actual graph that is generated, please let me know and I will forward it to you. DO FILE CODE: rename sex Sex rename rf Rf rename pyrs Pyrs rename age2cat Age2cat rename revproc Revproc /* loc tim Pyrs loc dth death quietly { gen byte death= 1-dead su `tim' if `tim'>0 } noisily stset `tim', fail(`dth') */ noisily stset Pyrs dead THE MACRO ADO FILE CODE: *! program diradj version 2.0.0 - 03jan2000 - GvM *! *! usage: diradj varlist *! *! first var in varlist is <main> regressor (sought effect) *! followed by other regressors (that are to be adjusted for) *! *! data must be -stset- first; all regressors must be binary *! *! performs "direct adjustment" as described in Maruch (1981) in order to *! obtain better adjustment than via classical "mean subject" (see Ghali etal) *! *! computed as total probability= sum( Pr(T>t|pattern i) * Pr(pattern i) ) *! *! where the survival probabilities are obtained via Cox regression *! and Pr(pattern i) is estimated as (N in i / total N) *! thus producing weighted averages suW0 (for main reg=no) and suW1 (yes) *! *! program also computes Crude survival: Cru0 and Cru1 *! as well as "mean patient" adjustment: suM0 and suM1 *! *! --> final dataset contains: _t (time) Cru0 Cru1 suM0 suM1 suW0 suW1 *! prog def adjbobo * version 6.0 st_is 2 analysis syntax varlist unab vL: `varlist' gettoken z vL:vL /* main regressor : e.g. Diab */ loc regs `vL' /* other regressors: e.g. creat CHF male */ loc nreg: word count `regs' loc xreg=`nreg'+1 tempname m A hrCru Bz Bm hr qui { * ---- patterns sca npat= 2^`nreg' /* possible number of patterns */ g np=0 loc j 0 sca `m'=2^`nreg' while `j'<`nreg' { sca `m'= int(`m'/2+.01) loc j=`j'+1 loc x: word `j' of `regs' qui replace np= np+`m' if `x' } mat `A'=J(npat,1,0) /* will contain the Ni's */ loc NN 0 loc p 0 while `p'<npat { count if np==`p' & _st loc NN=`NN'+r(N) loc p =`p'+1 mat `A'[`p',1]=r(N) } * ---- cox reg stcox `z' , basesurv(Cru0) /* Cru0= crude baseline surv */ sca `hrCru'=exp(_b[`z']) stcox `z' `regs' , basesurv(s0) /* s0= full reg baseline surv */ sca `Bz'=_b[`z'] sca `Bm'= 0 * -- adjust: mean patient --> haz h(t)= h0(t) * exp(Xbar'Beta) loc j 0 while `j'<`nreg' { loc j=`j'+1 loc x: word `j' of `regs' su `x' if _st, meanonly sca `Bm'=`Bm' + _b[`x']*r(mean) } sort _t qui { by _t: keep if _n==_N keep _t Cru0 s0 g double Cru1= Cru0^`hrCru' g double suM0= s0^exp(`Bm') g double suM1= s0^exp(`Bm'+`Bz') order _t C* su* g double suW0=0 g double suW1=0 } * -- adjust: = direct adjustment loc p 1 while `p'<=npat { if `A'[`p',1] { /* some patterns may be empty */ loc j =`nreg' sca `hr'=0 loc i=`p'-1 while `i' { loc x: word `j' of `regs' sca `hr'=`hr'+ mod(`i',2)*_b[`x'] loc j=`j'-1 loc i= int(`i'/2+.01) } replace suW0= suW0+`A'[`p',1] * s0^exp(`hr') replace suW1= suW1+`A'[`p',1] * s0^exp(`hr'+`Bz') } loc p=`p'+1 } drop s0 replace suW0= suW0/`NN' replace suW1= suW1/`NN' #delimit ; drop if Cru0==Cru0[_n-1] & suM0==suM0[_n-1] & suW0==suW0[_n-1] & Cru1==Cru1[_n-1] & suM1==suM1[_n-1] & suW1==suW1[_n-1] ; #delimit cr } show end * ----------- prog def show /* shows name of surv curve next to curve stata won't do that alone so it must be programmed low-level instead of calling this <show> could just say gr su* C* _t, xla yla s(..oo..) c(llllll) to produce the 6 curves instead of that 1 line need 46 lines! to obtain same plus the name and axes */ loc su suM0 suM1 suW0 suW1 Cru0 Cru1 loc wc: word count `su' tempvar x y qui g `x'=. qui g str8 `y'= "" loc grmin= 1 loc j 0 qui while `j'<`wc' { loc j=`j'+1 loc v: word `j' of `su' loc va=`v'[_N] loc grmin= min(`grmin',`va') replace `x'= `va' in `j' replace `y'="`v'" in `j' } loc grmin=int(50*`grmin')/50 loc a `su' _t, s(..dd..) c(llllll) xla yla(`grmin'(.02)1.0) t1(" ") sort `x' tempname X mkmat `x' in 1/`wc', mat(`X') loc tx loc ht loc k=`wc'+1 while `k'>1 { loc k= `k'-1 loc u= `x'[`k'] loc ht "`ht' `u'" loc u= `y'[`k'] loc tx "`tx' `u'" } sort _t loc bb bbox(0,0,23000,28900,700,300,0) gr `a' `bb' border gph open gr gph pen 1 loc k 0 while `k'<`wc' { loc k=`k'+1 loc w: word `k' of `tx' loc h: word `k' of `ht' loc h=18000*(1-`h')/(1-`grmin')+2200 gph text `h' 29000 0 -1 `w' } gph close end * exit Kathy Sabadosa, MPH Dartmouth-Hitchcock Medical Center Clinical Research Section, Department of Medicine 562E, Borwell One Medical Center Drive Lebanon, NH 03756 Phone (603)-650-7751 Phone (603)-643-3548 Fax (603)-650-8972 * * For searches and help try: * http://www.stata.com/support/faqs/res/findit.html * http://www.stata.com/support/statalist/faq * http://www.ats.ucla.edu/stat/stata/

- Prev by Date:
**st: standard errors after nonlinear function of regression coefficients** - Next by Date:
**Re: st: standard errors after nonlinear function ofregressioncoefficients** - Previous by thread:
**st: standard errors after nonlinear function of regression coefficients** - Next by thread:
**Re: st: standard errors after nonlinear function ofregressioncoefficients** - Index(es):

© Copyright 1996–2015 StataCorp LP | Terms of use | Privacy | Contact us | What's new | Site index |