Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

st: Stata Macro Issue

From   Kathryn.A.Sabadosa@Dartmouth.EDU (Kathryn A. Sabadosa)
Subject   st: Stata Macro Issue
Date   14 Oct 2002 15:28:13 EDT


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.


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


*! 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

* -----------
 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
		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

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:

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