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/
