********************************************************************** Attached below are the following ado files: stbtcalc.ado stbtplot.ado stgtcalc.ado stgtplot.ado mkname1.ado running.ado And help files: stbtcalc.hlp stbtplot.hlp stgtcalc.hlp stgtplot.hlp running.hlp ********************************************************************** stbtcalc.ado Instructions for use: Cut out the following program and put it in a file called "stbtcalc.ado". Be sure that the file is saved as a plain text (ASCII) file and that it has a hard return at the end of the last line in the file. Put stbtcalc.ado in your C:\ADO directory (you may have to create this directory) or your Stata working directory (i.e., your current directory). Note: Do not put it in the C:\STATA\ADO directory. Only official Stata ado and hlp files should be placed in this directory. For more information, see [2] ado in Vol. 1 of the Reference manuals. ---------------------------------cut here---start-of-stbtcalc.ado------ *! version 1.0.2 PR/PS 08May98. program define stbtcalc version 5.0 st_is local varlist "opt ex none" #delimit ; local options "Window(int 0) Increment(int 0) CLuster(string) noHR Level(integer $S_level) Robust noSHow SMooth COEF(str) SE(str) *" ; #delimit cr local if "opt" local in "opt" parse "`*'" if "`coef'"=="" {local coef BT} if "`se'" =="" {local se BS} local id : char _dta[st_id] local t : char _dta[st_t] local t0 : char _dta[st_t0] local d : char _dta[st_d] local id : char _dta[st_id] local w : char _dta[st_w] local wt : char _dta[st_wt] quietly { tempvar touse dead mark `touse' `if' `in' `w' markout `touse' `t' `t0' `d' `varlist' markout `touse' `cluster' `id', strok if "`wt'"=="pweight" { local robust "robust" } if "`robust'"!="" & "`cluster'"=="" & "`id'"!="" { local cluster "`id'" } if "`cluster'"!="" { local cluster "cluster(`cluster')" } if "`t0'" != "" { local t0 "t0(`t0')" } if "`d'" != "" { gen byte `dead'=`d'!=0 if `d'!=. } else { gen byte `dead'=1 } tempvar wdead wevents gen byte `wdead'=`dead' local D "dead(`d')" local d "dead(`wdead')" st_show `show' if "`id'"!="" & "`t0'"==""{ di in red "id not currently supported without t0" exit 198 } if "`w'`wt'"!="" { di in red "weights not currently supported" exit 198 } cox `t' `varlist' `w' if `touse', `robust' `cluster' `t0' nocoef /* */ `options' `D' tempname beta matrix `beta' = get(_b) local vars : colnames (`beta') if "`vars'"!="`varlist'" { noi di "Some variables not used by cox regression " /* */ "(presumably due to collinearity)" noi di "The following will be used by stbtcalc: " noi di "`vars'" } /* Set up _b() and _se() variables */ parse "`vars'", parse(" ") local nv 0 while "`1'"!="" { local nv = `nv'+1 local var`nv' `1' tempvar b`nv' s`nv' gen `b`nv''=. gen `s`nv''=. global B_`nv' = _b[`1'] global B_`nv' = "`1' ${B_`nv'}" mac shift } local j=`nv'+1 global B_`j' count if `dead'==1 & `touse'==1 local ndead = _result(1) /* Set up "windowed" survival times and censoring indicator */ gsort -`touse' `t' gen int `wevents' = sum(`dead'==1) if `touse' /* Calc Cox model for moving windows each containing `window' events. Default window size is 10*(#covariates). Default increment in #events aims for at least 50 estimated b(t)'s. */ tempname oldest est hold `oldest' if `window'<=0 { local window = max(20, int(10*`nv')) } if `increme'<=0 { local increme = max(1, int((`ndead'-`window')/50)) } local midevnt = int(1+`window'/2) local i 0 tempvar touze gen byte `touze'=`touse' while `i'<=`ndead'-`window' { noi di "." _cont local index = `i'+`midevnt' sum `t' if `touze' & (`wevents'>`i') & (`wevents'<=`i'+`window') replace `touze' = `touze' & (`t'>_result(5)) replace `wdead' = cond(`t'>_result(6),0,`dead') cox `t' `vars' `w' if `touze', /* */ `robust' `cluster' `t0' `d' `options' nocoef local j 0 while `j'<`nv' { local j=`j'+1 cap replace `b`j''=_b[`var`j''] if `wevents'==`index' cap replace `s`j''=_se[`var`j''] if `wevents'==`index' } local i=`i'+`increme' } /* Create new variables from coefficients and confidence limits. */ cap drop `coef'* cap drop `se'* local j 0 while `j'<`nv' { local j=`j'+1 mkname1 `var`j'', prefix(`coef') local B $S_1 mkname1 `var`j'', prefix(`se') running `b`j'' `t' if `touse', gen(`B') rep(2) nograph running `s`j'' `t' if `touse', gen(`S') rep(2) nograph locae smooth "smoothed " } else { rename `b`j'' `B' rename `s`j'' `S' } lab var `B' "`var`j'' `smooth'coefficient" lab var `S' "`var`j'' standard error" } est unhold `oldest' noi cox global S_1 `window' global S_2 `increme' global S_E_cmd "stbtcalc" } end ---------------------------------cut here---end-of-stbtcalc.ado-------- ********************************************************************** stbtplot.ado Instructions for use: Cut out the following program and put it in a file called "stbtplot.ado". Be sure that the file is saved as a plain text (ASCII) file and that it has a hard return at the end of the last line in the file. Put stbtplot.ado in your C:\ADO directory (you may have to create this directory) or your Stata working directory (i.e., your current directory). Note: Do not put it in the C:\STATA\ADO directory. Only official Stata ado and hlp files should be placed in this directory. For more information, see [2] ado in Vol. 1 of the Reference manuals. ---------------------------------cut here---start-of-stbtplot.ado------- *! version 1.0.0 27Mar98. program define stbtplot version 5.0 st_is if "$S_E_cmd"!="stbtcalc" { error 301 } local varlist "req min(1)" local if "opt" #delimit ; local options "CI LEvel(int $S_level) Connect(str) Symbol(str) SAving(str) HR TItle(string) L1title(str) B2title(str) noRescale MEan SINgle COEF(str) SE(str) *" ; #delimit cr parse "`*'" local nvar : word count `varlist' if `level'<10 | `level'>99 { error 198 } if "`ci'"!="" & "`single'"!="" { noi di in red "ci and single cannot be used together" exit 198 } local l1 "l1title(`l1title')" if "`rescale'"=="" { local rescale "Rescale"} if "`hr'"!="" { local hr1 "exp(" local hr2 ")" local base 1 tempvar expy qui gen `expy'=. local l1p "Hazard Ratio" } else { local base 0 local l1p "Parameter" } if "`l1title'"=="" { local l1 "`l1p' Estimate" local l1 "l1title(`l1')" } if "`t1title'"!="" {local t1 "t1title(`t1title')"} if "`b2title'"=="" { local b2 : var lab $S_E_depv if "`b2'"=="" { local b2 "$S_E_depv" } local b2 "b2title(`b2')" } else local b2 "b2title(`b2title')" if "`title'"=="" { local title "Beta of t Plot" } if "`saving'"!="" { local save "sav(`saving')" } tempname z scalar `z' = -invnorm((100-`level')/200) capture confirm var _touse if _rc==0 { if "`if'"=="" { local if "if _touse==1" } else local if "`if' & _touse==1" } local tvar $S_E_depv if "`coef'"=="" {local coef BT} if "`se'" =="" {local se BS} parse "`varlist'",parse(" ") tempvar lci uci qui gen `lci'=. qui gen `uci'=. local i 0 if "`single'"=="" { while "`1'"!="" { local i=`i'+1 vlook4 `1' `coef'* local y $S_1 mlook4 `1' local b0: word 2 of $$S_1 if "`ci'"!="" { vlook4 `1' `se'* local s $S_1 } local vlab : var lab `y' local vlab : word 1 of `vlab' local vlab2 : var lab `vlab' if "`valb2'"!="" {local vlab "`vlab2'"} if "`t1title'"=="" {local t1 "t1title(`vlab')"} if "`ci'"!="" { qui replace `lci'=`hr1' `y'-`z'*`s' `hr2' qui replace `uci'=`hr1' `y'+`z'*`s' `hr2' } if "`hr'"!="" { local b0 = exp(`b0') qui replace `expy'=exp(`y') local y "`expy'" } if `nvar'>1 { tempfile gr`i' set graphics off graph `y' `lci' `uci' `tvar', s(iii) c(llll) /* */ ti("`title'") yline(`base',`b0') `l1' `t1' `options' sav(`gr`i'') local graphs "`graphs' `gr`i''" } else { graph `y' `lci' `uci' `tvar', s(iii) c(llll) /* */ ti("`title'") yline(`base',`b0') `l1' `t1' `options' `save' } mac shift } if `nvar'>1 { set graphics on gr using `graphs', `save' } } else { local Bs "`base'" while "`1'"!="" { local i=`i'+1 vlook4 `1' `coef'* local y`i' $S_1 mlook4 `1' local b0: word 2 of $$S_1 local vlab : var lab `y`i'' local vlab : word 1 of `vlab' local vlab2 : var lab `vlab' if "`valb2'"!="" {local vlab "`vlab2'"} if "`hr'"!="" { local b0 = exp(`b0') tempvar expy`i' qui gen `expy`i''=exp(`y`i'') lab var `expy`i'' `vlab' local ys "`ys' `expy`i''" } else { local ys "`ys' `y`i''" } local Bs "`Bs',`b0'" local con "`con'l" local sym "`sym'i" mac shift } graph `ys' `tvar', c(`con') s(`sym') ti("`title'") yline(`Bs') /* */ `l1' `t1' `options' sort `save' } end program define vlook4 version 3.0 local varlist "req ex" parse "`*'" parse "`varlist'", parse(" ") local find "`1'" mac shift while "`1'"!="" { local lbl : variable label `1' local var : word 1 of `lbl' if substr("`var'",1,.)==substr("`find'",1,.) { global S_1 "`1'" exit } mac shift } end program define mlook4 version 3.0 local find "`1'" local i 1 while "${B_`i'}"!="" { local var : word 1 of ${B_`i'} if substr("`var'",1,.)==substr("`find'",1,.) { global S_1 "B_`i'" exit } local i = `i' + 1 } end ---------------------------------cut here---end-of-stbtplot.ado-------- ********************************************************************** stgtcalc.ado Instructions for use: Cut out the following program and put it in a file called "stgtcalc.ado". Be sure that the file is saved as a plain text (ASCII) file and that it has a hard return at the end of the last line in the file. Put stgtcalc.ado in your C:\ADO directory (you may have to create this directory) or your Stata working directory (i.e., your current directory). Note: Do not put it in the C:\STATA\ADO directory. Only official Stata ado and hlp files should be placed in this directory. For more information, see [2] ado in Vol. 1 of the Reference manuals. ---------------------------------cut here---start-of-stgtcalc.ado------- *! version 1.0.1 PR 17Dec97. rev PS 23Dec97 5.25pm * This version looks for a variable called _touse created by stcox * identifying those observations remaining after mark and markout * program define stgtcalc version 5.0 st_is if "$S_E_cmd2"!="stcox" { error 301 } tempname V coef tempvar c matrix `V' = get(VCE) matrix `coef' = get(_b) local vars : colnames (`coef') local t : char _dta[st_t] local dead : char _dta[st_d] local t0 : char _dta[st_t0] local id : char _dta[st_id] local w : char _dta[st_w] local wt : char _dta[st_wt] if "`id'"!="" & "`t0'"==""{ di in red "id not currently supported without t0" exit 198 } if "`w'`wt'"!="" { di in red "weights not currently supported" exit 198 } quietly { capture confirm var _touse if _rc==0 { local if "if _touse==1" } if "`dead'"!="" { local ifdead "if (`dead'!=0 & `dead'!=.)" if "`if'"!="" { local ifdead "`ifdead' & _touse==1" } gen byte `c' = !`dead' `if' } else { gen byte `c' = 0 `if' /* all dead */ local ifdead "`if'" } count `ifdead' local ndead = _result(1) matrix `V' = `V'*`ndead' tempvar hr order mark IND s0 s1 predict `hr' local nobs = _N if "`t0'" != "" { unabbrev `t0', max(1) local t0 "$S_1" local t0ttl "`t0'" capture assert `t0'<`t' `if' if _rc { di in red "t0(`t0') >= `t' in some obs." exit 498 } preserve tempvar index gen long `index' = _n compress `index' sort `index' local fn "$S_FN" tempfile old save "`old'" capture keep `if' keep `vars' `hr' `t' `dead' `t0' `id' $S_E_svn `index' `c' expand 2 gen byte `mark'=-cond(_n<`nobs',2,cond(`c',3,1)) /* -2=enter -3=censor -1=die */ replace `t'=`t0' in 1/`nobs' } else { gen byte `mark'=-cond(`c',3,1) } gen `order' = -`t' sort $S_E_svn `order' `mark' if "${S_E_svn}"!="" { local S_by "by ${S_E_svn} :" } `S_by' gen double `s0'=sum(cond(`mark'==-2,-`hr',`hr')) by $S_E_svn `order' : replace `s0'=`s0'[_N] if `mark'==-1 gen double `s1' = . parse "`vars'", parse(" ") local nv 0 while "`1'"!="" { local nv = `nv'+1 local var`nv' `1' tempvar res`nv' GT`nv' `S_by' replace `s1'=sum(`1'*cond(`mark'==-2,-`hr',`hr')) by $S_E_svn `order' : replace `s1'=`s1'[_N] if `mark'==-1 gen `res`nv'' = cond(`mark'==-1,(`1'-`s1'/`s0'),.) local rvars "`rvars' `res`nv''" local GTvars "`GTvars' `GT`nv''" mac shift } cap drop GT* mxv `rvars' `if', mat(`V') gen(`GTvars') local i 1 while `i'<=`nv' { quietly replace `GT`i'' = `GT`i''+`coef'[1,`i'] mkname1 `var`i'', prefix(GT) rename `GT`i'' $S_1 lab var $S_1 "`var`i'' stand. score res." local i = `i'+1 } if "`t0'" != "" { keep if `t'>`t0' sort `index' merge `index' using `old' drop _merge global S_FN `fn' restore, not } } end *! Based on PDS _mxv version 1.0.1 25May1995 program def mxv version 4.0 local varlist "req ex" local in "opt" local if "opt" local options "Matrix(string) Generate(string)" parse "`*'" tempname nrow scalar `nrow'=rowsof(`matrix') tempname M vec matrix `M' = `matrix' matrix colnames `M' = `varlist' parse "`generate'", parse(" ") quietly { local i 0 while "`1'" !="" { local i = `i'+1 matrix `vec'=`M'[`i',.] matrix score `1' = `vec' `if' `in' mac shift } } end ---------------------------------cut here---end-of-stgtcalc.ado-------- ********************************************************************** stgtplot.ado Instructions for use: Cut out the following program and put it in a file called "stgtplot.ado". Be sure that the file is saved as a plain text (ASCII) file and that it has a hard return at the end of the last line in the file. Put stgtplot.ado in your C:\ADO directory (you may have to create this directory) or your Stata working directory (i.e., your current directory). Note: Do not put it in the C:\STATA\ADO directory. Only official Stata ado and hlp files should be placed in this directory. For more information, see [2] ado in Vol. 1 of the Reference manuals. ---------------------------------cut here---start-of-stgtplot.ado------- *! version 1.0.1 PR 08May98. program define stgtplot version 5.0 st_is *if "$S_E_cmd2"!="stcox" { * error 301 *} local varlist "req min(1)" local if "opt" #delimit ; local options "CI LEvel(int $S_level) Connect(str) Symbol(str) YLIne(str) TItle(string) L1title(str) B2title(str) noRescale MEan *" ; #delimit cr parse "`*'" local nvar : word count `varlist' if `level'<10 | `level'>99 { error 198 } local l1 "l1title(`l1title')" if "`rescale'"=="" { local rescale "Rescale"} if "`l1title'"=="" { local l1 "Parameter Estimate" local l1 "l1title(`l1')" } if "`t1title'"!="" {local t1 "t1title(`t1title')"} if "`b2title'"=="" { local b2 : var lab $S_E_depv if "`b2'"=="" { local b2 "$S_E_depv" } local b2 "b2title(`b2')" } else local b2 "b2title(`b2title')" if "`title'"=="" { local title "Grambsch and Therneau Plot" } if "`yline'"=="" { local yline 0 } local yl "yline(`yline')" tempname z scalar `z' = -invnorm((100-`level')/200) capture confirm var _touse if _rc==0 { if "`if'"=="" { local if "if _touse==1" } else local if "`if' & _touse==1" } local tvar $S_E_depv local prefix GT parse "`varlist'",parse(" ") local i 0 while "`1'"!="" { local i=`i'+1 cap unabbrev "`prefix'`1'" if _rc { di in red "`1' not found" exit 111 } local gvar $S_1 local vlist "`vlist' `gvar'" local stlist "`stlist' `tvar' `gvar'" local vlab`i' : var lab `gvar' local vlab`i' : word 1 of `vlab`i'' local y `gvar' mac shift } quietly { if `nvar'>1 { preserve if "`if'"!="" { keep `if' local if } keep $S_E_depv `vlist' describe, short if _result(1)==0 { error 2000 } local maxobs =_result(4) local reqobs=`nvar'*_result(1) if `maxobs'<`reqobs' { tempfile user save `user' drop _all local reqobs=int(`reqobs'*1.02)+100 set maxobs `reqobs' use `user' erase `user' } tempvar tvar y stack `stlist', into(`tvar' `y') clear lab var `tvar' "Time" tempvar param rename _stack `param' local _stack "`param'" local byst "by `param':" } tempvar Y sm if "`ci'"!="" { tempvar lci uci se } if `nvar'>1 { gen `Y'=. if "`ci'"!="" { gen `lci'=. gen `uci'=. } local j 0 while `j'<`nvar' { local j=`j'+1 local slab "`slab' `j' `vlab`j''" if "`ci'"!="" { running `y' `tvar' if `_stack'==`j', `ci' /* */ gen(`sm') gense(`se') nograph /* */ `mean' replace `lci'=`sm'-`z'*`se' if `_stack'==`j' replace `uci'=`sm'+`z'*`se' if `_stack'==`j' drop `se' } else { running `y' `tvar' if `_stack'==`j', /* */ gen(`sm') /* */ nograph `mean' } replace `Y'=`sm' if `_stack'==`j' drop `sm' } tempname stlab lab def `stlab' `slab' lab val `_stack' `stlab' local by "by(`param') `rescale'" } else { if "`t1title'"=="" { local t1 "t1title(`vlab1')" } if "`ci'"!="" { running `y' `tvar' `if', gen(`Y') `ci' /* */ gense(`se') nograph `mean' gen `lci'=`Y'-`z'*`se' gen `uci'=`Y'+`z'*`se' } else { running `y' `tvar' `if', gen(`Y') /* */ nograph `mean' } } } sort `_stack' `tvar' graph `Y' `lci' `uci' `tvar', s(iii) c(llll) ti("`title'") /* */ xsca(0,.) yline(`yline') `by' `l1' `t1' `options' if `nvar'>1 { if `maxobs'<`reqobs' { drop _all quietly set maxobs `maxobs' } } end ---------------------------------cut here---end-of-stgtplot.ado-------- ********************************************************************** mkname1.ado Instructions for use: Cut out the following program and put it in a file called "mkname1.ado". Be sure that the file is saved as a plain text (ASCII) file and that it has a hard return at the end of the last line in the file. Put mkname1.ado in your C:\ADO directory (you may have to create this directory) or your Stata working directory (i.e., your current directory). Note: Do not put it in the C:\STATA\ADO directory. Only official Stata ado and hlp files should be placed in this directory. For more information, see [2] ado in Vol. 1 of the Reference manuals. ---------------------------------cut here---start-of-mkname1.ado------- *! version 1.0.2 PS/PR 25Apr96. * mkname1 var , prefix() checkin() banned() noclear program define mkname1 version 4.0 local varlist "req ex min(1) max(1)" local options "Prefix(string) Checkin(string) BANned(str) noCLear" parse "`*'" if "`prefix'"=="" | length("`prefix'")>7 { di in red "invalid prefix()" exit 198 } if "`checkin'"!=""{ if index("`checkin'","`varlist'") == 0 { di in red "`varlist' not in check list" exit 111 } } local l=8-length("`prefix'") local i 0 local new ="`prefix'"+substr("`varlist'",1,`l') local rc = index("`banned'","`new'") if "`clear'"!=""{ capture confirm new var `new' local rc2=_rc } else local rc2 0 while (`rc2' | `rc') & length("`new'")<9 { local i=`i'+1 local new "`prefix'`i'" local rc = index("`banned'","`new'") if "`clear'"!=""{ capture confirm new var `new' local rc2=_rc } } if length("`new'")>8 { di in red "Unable to create unique name. " in blue "Try shorter prefix." exit } global S_1 `new' end ---------------------------------cut here---end-of-mkname1.ado--------- ********************************************************************** running.ado Instructions for use: Cut out the following program and put it in a file called "running.ado". Be sure that the file is saved as a plain text (ASCII) file and that it has a hard return at the end of the last line in the file. Put running.ado in your C:\ADO directory (you may have to create this directory) or your Stata working directory (i.e., your current directory). Note: Do not put it in the C:\STATA\ADO directory. Only official Stata ado and hlp files should be placed in this directory. For more information, see [2] ado in Vol. 1 of the Reference manuals. ---------------------------------cut here---start-of-running.ado------- *! version 2.0.0 PS/PR 10Jun97. STB-41 sed9.1 *!! Based on Peter Sasieni's "running4.ado" emailed 09Jun97 program define running version 4.0 local varlist "req ex min(1) max(2)" local if "opt" local in "opt" local weight "aweight" #delimit ; local options "noGraph GEN(string) Knn(str) Double Repeat(int 1) Mean LOGit SPan(real 0) TItle(string) Symbol(str) Connect(str) PEn(str) CI GENSe(str) GENB(str) TWice *"; #delimit cr parse "`*'" local nv = (2-("`mean'"!=""))*(8+16)+9+12*("`weight'"!="")+8*(1-("`twice'"!="")) memchk byte `nv' if _rc==900 {exit} if "`gen'"!="" { confirm new var `gen' } if "`genb'"!="" { if "`twice'`mean'`logit'"!="" { di in red "genb not available with mean, twice or logit" exit 198 } confirm new var `genb' } if "`gense'"!="" { confirm new var `gense' if `repeat'>1 | "`double'`logit'`twice'"!="" { di in red "gense not available with repeat>1, twice or logit" exit 198 } } if "`ci'"!="" { if `repeat'>1 | "`double'`logit'`twice'"!="" { di in red "ci not available with repeat>1, twice or logit" exit 198 } } if `span'!=0 & "`knn'"!="" { di in red "cannot specify both span and knn" exit 198 } if `span'<0 | `span'>2-("`mean'"!="") { di in red "span must be between 0 and " 2-("`mean'"!="") exit 198 } if "`knn'"=="" { local knn 0 /* previous default */ } else { cap confirm var `knn' if _rc { confirm integer num `knn' } else { unabbrev `knn' local kvec $S_1 } } if "`double'"!="" { local repeat=2*`repeat' } local nrep `repeat' if `repeat'>7 { local repeat 7 noisily di in bl "[repeat set to 7]" } parse "`varlist'", parse(" ") local y `1' local x `2' tempvar Y X smooth sy rsy touse kr kl quietly { if "`weight'" != "" { tempvar sw w wt gen `w' `exp' replace `w'=. if `w'<=0 local exp "[`weight'=`w']" } mark `touse' `if' `in' markout `touse' `y' `x' `w' count if `touse' local cnt = _result(1) local IN "in 1/`cnt'" if `cnt'<5 { noisily error 2001 } if `span'!=0 { local knn=(`cnt'*`span'-1)/2 } else { if `knn'==0 {local knn = .5*`cnt'^0.8 } local span=(2*`knn'+1)/`cnt' } if "`kvec'" !="" { sum `kvec' if `touse' if _result(5)<1 | _result(1)!=`cnt' { nois di "knn-variable must be positive & non-missing" exit } local kk int(`kvec'/sqrt(`repeat') +.5) } else { local kk = int(`knn' /sqrt(`repeat') +.5) if `kk'<=0 { noi di in red "Span too small. Increase span or knn." exit 2002 } } if `knn'<=1 & "`ci'"!="" {nois di "[ci not produced when knn=1]"} sum `y' `exp' if `touse' local ycen=_result(3) gen `Y'=`y'-`ycen' if `touse' _crcslbl `Y' `y' if "`x'"!=""{ gen `X'=`x' if `touse' _crcslbl `X' `x' } else { gen `X'=_n if `touse' lab var `X' "n" } set graph off graph `y' `X' `exp' in 1/3, `options' set graph on sort `X' if "`weight'" != "" { gen double `sw'=. replace `w'=. if !`touse' sum `w' noi di in gr "(sum of wgt is " _result(1)*_result(3) ")" replace `w' = `w'/_result(3) gen `wt' = `w' local ks "(`sw'[_n+`kr']-cond(_n>`kl',`sw'[_n-`kl'],0))" local d3 "`sw'[_N]" local wX "`w'*" local wtX "`wt'*" local wn "`wt'`nk'*" local g_swby "by `X':replace `sw'=sum(`w') " local g_wt "by `X':replace `wt'=`sw'[_N]/_N " local gen_swt "replace `sw'=sum(`wt') `IN'" } else { local ks "(`kr'+`kl')" local d3 "_N" } gen `smooth'=`Y' lab var `smooth' "Smooth fit" gen double `sy'=sum(cond(_n<`cnt',(`X'==`X'[_n-1] |`X'==`X'[_n+1]), /* */ `X'==`X'[_n-1])) `IN' gen `rsy' = . local ties=`sy'[`cnt'] local xcen=(`X'[1]+`X'[`cnt'])/2 replace `X'=`X'-`xcen' `IN' sort `X' if `ties'>0 { tempvar Yt `g_swby' /* sw must be defined before d3 is used */ `g_wt' by `X':replace `sy'=sum(`wX'`Y') if `touse' by `X':replace `smooth'=`sy'[_N]/`d3' if `touse' if "`twice'"!=""{ gen `Yt'=`smooth' } } else { local Yt "`Y'" } `gen_swt' gen int `kl' = min(`kk'+1,_n) gen int `kr' = min(`kk',`cnt'-_n) if "`mean'"=="" { local mean "line" tempvar sxy rsxy beta rsx rsxx gen `beta'= . gen double `sxy'=sum(`wtX'`X') `IN' diff `rsx' `sxy' `kr' `kl' `ks' `cnt' gen replace `sxy'=sum(`wtX'((`X')^2)) `IN' diff `rsxx' `sxy' `kr' `kl' `ks' `cnt' gen gen `rsxy'=. } while `repeat'>0 { replace `sy'=sum(`wtX'`smooth') `IN' diff `rsy' `sy' `kr' `kl' `ks' `cnt' replace if "`mean'"=="line" { replace `sxy'=sum(`wtX'`X'*`smooth') `IN' diff `rsxy' `sxy' `kr' `kl' `ks' `cnt' replace replace `beta'=(`rsxy'-`rsx'*`rsy')/(`rsxx'-(`rsx')^2) `IN' replace `smooth'=cond((`rsxx'-`rsx'*`rsx') > max(0, /* */ `rsxx'/1e8), `rsy'+`beta'*(`X'-`rsx'), `rsy') `IN' } else { replace `smooth'=`rsy' `IN' } local repeat=`repeat'-1 } if "`twice'"!="" { tempvar res gen `res' = `Yt'-`smooth' local repeat `nrep' while `repeat'>0 { replace `sy'=sum(`wtX'`res') `IN' diff `rsy' `sy' `kr' `kl' `ks' `cnt' replace if "`mean'"=="line" { replace `sxy'=sum(`wtX'`X'*`res') `IN' diff `rsxy' `sxy' `kr' `kl' `ks' `cnt' replace replace `beta'=(`rsxy'-`rsx'*`rsy')/(`rsxx'-(`rsx')^2) `IN' replace `res'=cond((`rsxx'-`rsx'*`rsx') > max(0, /* */ `rsxx'/1e8), `rsy'+`beta'*(`X'-`rsx'), `rsy') `IN' } else { replace `res'=`rsy' `IN' } local repeat=`repeat'-1 } replace `smooth'=`smooth'+`res' `IN' } if `ties'>0 { by `X':replace `sy'=sum(`smooth') if `touse' by `X':replace `smooth'=`sy'[_N]/_N if `touse' } if "`ci'"!="" | "`gense'"!="" { drop `sy' `rsy' tempvar se sigma2 gen double `se'=(`Y'-`smooth')^2 if `ties'>0 { `g_swby' by `X':replace `se'=sum(`wX'`se') if `touse' by `X':replace `se'=`se'[_N]/`d3' if `touse' } `gen_swt' replace `se'=sum(`wtX'`se') `IN' * Adjustment for estimating parameters made in se not sigma2 * diff `sigma2' `se' `kr' `kl' `ks' `cnt' gen if "`mean'"=="line" { replace `se'=sqrt(`sigma2'* /* */ cond((`rsxx'-`rsx'*`rsx') > max(0,`rsxx'/1e8), /* */ ((`kl'+`kr')/(`kl'+`kr'-2)) /* */ *(1+(`X'-`rsx')^2/(`rsxx'-(`rsx')^2))/`ks' , /* */ 1/(`kl'+`kr'-1) )) `IN' } else { replace `se'=sqrt(`sigma2'/(`kl'+`kr'-1)) `IN' } drop `sigma2' if `ties'>0 { `g_swby' by `X':replace `se'=sum((`se')^2) if `touse' by `X':replace `se'=sqrt(`se'[_N]/_N) if `touse' } } replace `X'=`X'+`xcen' `IN' replace `Y'=`Y'+`ycen' `IN' replace `smooth'=`smooth'+`ycen' `IN' if "`mean'"=="mean" &" `kvec'"=="" { local k2=`cnt'+1-`kk' replace `smooth'=. in 1/`kk' replace `smooth'=. in `k2'/l } if "`logit'"!="" { sum `y' if `touse' if _result(5)<-5e-9 | _result(6)>1+5e-9 { noi di "[Yvar out of range [0,1]. logit not available]" } else { sum `y' if `y'>0 & `y'<1 &`touse' local aL = min(1/(min(`span',1)*`cnt'+1),_result(5)) local aU = max(1-1/(min(`span',1)*`cnt'+1),_result(6)) replace `smooth'=log(cond(`smooth'< `aL', `aL'/(1-`aL'), /* */ cond(`smooth'>`aU',`aU'/(1-`aU'),`smooth'/(1-`smooth')) /* */ )) if `smooth' !=. `IN' local laL = log(`aL'/(1-`aL')) local laU = log(`aU'/(1-`aU')) local adj=(`laU'-`laL')/25 replace `Y'=cond(`y'<`aL', `laL'-(1.3+uniform())*`adj', /* */ cond(`y'>`aU', `laU'+(1.3+uniform())*`adj', /* */ log(`y'/(1-`y')) )) `IN' } } * GRAPH (If required) if "`graph'"!="nograph" { if "`ci'"!="" { local df=`kl'+`kr'-("`mean'"=="line") local t=invt(`df',${S_level}/100) tempvar l u gen `l'=`smooth'-`t'*`se' gen `u'=`smooth'+`t'*`se' local pp "44" } if "`title'" ==""{ local title "Running `mean' smoother" } if "`symbol'" ==""{ if `cnt'<300 { local symbol o } else { local symbol . } } local symbol "`symbol'iii" if "`connect'"==""{ local connect ".lll" } if "`pen'" ==""{ local pen "23`pp'" } drop `kl' `kr' `touse' noisily graph `Y' `smooth' `l' `u' `X' `exp', `options' /* */ ti(`title') s(`symbol') c(`connect') pen(`pen') } if "`gen'`genk'`gense'" !="" {drop `Y' `X'} if "`gen'" !="" {rename `smooth' `gen'} if "`gense'"!="" {rename `se' `gense'} if "`genb'" !="" {rename `beta' `genb'} } global S_1 `knn' global S_2 `span' end prog def diff local new `1' local old `2' local kr `3' local kl `4' local ks `5' local cnt `6' local com `7' `com' `new' =(`old'[_n+`kr'] - cond(_n>`kl',`old'[_n-`kl'],0))/`ks' /* */ in 1/`cnt' end *! version 1.0 17 Jun 1997 program define memchk version 5.0 quietly describe, detail short local width = _result(3) local ws = _result(6) while "`1'"!="" { if "`2'"=="" { error 198 } confirm integer number `2' if "`1'"=="int" { local width = `width' + 2*`2' } else if "`1'"=="byte" { local width = `width' + `2' } else if "`1'"=="long" | "`1'"=="float" { local width = `width' + 4*`2' } else if "`1'"=="double" { local width = `width' + 8*`2' } else if substr("`1'",1,3)=="str" { local len = substr("`1'",4,.) confirm integer number `len' local width = `width' + `len'*`2' } else { error 198 } mac shift 2 } if `width' > `ws' { di _col(8) in red "insufficient memory" exit 900 } end ---------------------------------cut here---end-of-running.ado--------- ********************************************************************** stbtcalc.hlp Instructions for use: Cut out the following help file and put it in a file called "stbtcalc.hlp". Be sure that the file is saved as a plain text (ASCII) file and that it has a hard return at the end of the last line in the file. Put stbtcalc.hlp in your C:\ADO directory (you may have to create this directory) or your Stata working directory (i.e., your current directory). Note: Do not put it in the C:\STATA\ADO directory. Only official Stata ado and hlp files should be placed in this directory. ---------------------------------cut here---start-of-stbtcalc.hlp------ .- help for ^stbtcalc^, ^stbtplot^ (Patrick Royston, Peter Sasieni) .- Time-varying regression coefficients in Cox PH models ----------------------------------------------------- ^stbtcalc^ varlist [^,^ ^w^indow^(^#^)^ ^i^ncrement^(^#^)^ ^sm^ooth stcox_options ] ^stbtplot^ [varlist] [^, ci le^vel^(^#^)^ graph_options ] Description ----------- ^stbtcalc^ calculates local estimates of the regression coefficients in the Cox model defined by ^stset^ with predictors in varlist. The estimates are calculated within moving windows which bracket a fixed number of "deaths" (events). ^stbtplot^ plots the estimates and confidence intervals against the failure times. varlist is (a subset of) the predictors in ^stbtcalc^'s varlist. Options (^stbtcalc^) ------- ^window(^#^)^ sets the window size. The larger the value of #, the smoother the estimated b(t) curve is, but the less far it reaches into the tails of the distribution of failure times. The default value of # is 10 * (number of xvars) or 20, whichever is the larger. ^increment(^#^)^ sets the incremental number of events. The event window is advanced by # failures. If # is set greater than 1, fewer Cox models are fitted and the execution time diminishes, which may be desirable in datasets with a large number of events. The default value of # is ((number of deaths) - (window size))/50 or 1, whichever is the larger. Provided there are enough deaths, this choice guarantees at least 50 local estimates of b(t). ^smooth^ smooths the estimated b(t) curve and its standard error using a running line smoother (see help @running@). stcox_options are any of ^cluster()^, ^hr^, ^robust^. Options (^stbtplot^) ------- ^ci^ causes a ^level()^% confidence band to be plotted around the windowed estimates of b(t). The pointwise confidence intervals are based on the standard error of each b(t). ^level(^#^)^ sets the width of a pointwise confidence band for the estimated b(t). Default: current value of system macro ^$S_level^ (usually 95%). graph_options are any of Stata's ^graph, twoway^ options Remarks ------- ^stbtcalc^ provides a simple method of finding local estimates of the regression coefficient for each covariate. All failures which occur outside the moving time window circumscribed by the minimum and maximum failure times of a fixed number of events become censored. Because of the mode of calculation, no estimates are available in the initial and final phases of follow-up. ^stbtcalc^ saves the local estimates of b(t) to new variables called ^BT^xvar1, ^BT^xvar2 etc, corresponding to the xvars. Standard errors are saved in other new variables called ^BS^xvar1, ^BS^xvar2, ... . If more than one xvar is to be plotted by ^stbtplot^, a panel of plots is produced. The smoothing for b(t) and its standard error (see the ^smooth^ option) uses @running@ with resmoothing option ^repeat(2)^. Example ------- . ^stset survtime cens^ . ^stbtcalc x1 x2 x3 x4 x5^ . ^stbtplot x2^ . ^stbtplot x1 x2 x3 x4, ci xlabel ylabel^ Authors ------- Patrick Royston Imperial College School of Medicine, London FAX: (011)-44-181-383-8573 Peter Sasieni Imperial Cancer Research Fund, London FAX: (011)-44-171-269-3429 Also see -------- Manual: [R] st stcox On-line: help for @stcox@, @stgtcalc@ ---------------------------------cut here---end-of-stbtcalc.hlp-------- ********************************************************************** stbtplot.hlp Instructions for use: Cut out the following help file and put it in a file called "stbtplot.hlp". Be sure that the file is saved as a plain text (ASCII) file and that it has a hard return at the end of the last line in the file. Put stbtplot.hlp in your C:\ADO directory (you may have to create this directory) or your Stata working directory (i.e., your current directory). Note: Do not put it in the C:\STATA\ADO directory. Only official Stata ado and hlp files should be placed in this directory. ---------------------------------cut here---start-of-stbtplot.hlp------ .h stbtcalc ---------------------------------cut here---end-of-stbtplot.hlp-------- ********************************************************************** stgtcalc.hlp Instructions for use: Cut out the following help file and put it in a file called "stgtcalc.hlp". Be sure that the file is saved as a plain text (ASCII) file and that it has a hard return at the end of the last line in the file. Put stgtcalc.hlp in your C:\ADO directory (you may have to create this directory) or your Stata working directory (i.e., your current directory). Note: Do not put it in the C:\STATA\ADO directory. Only official Stata ado and hlp files should be placed in this directory. ---------------------------------cut here---start-of-stgtcalc.hlp------ .- help for ^stgtcalc^, ^stgtplot^ (Patrick Royston, Peter Sasieni) .- Time-varying regression coefficients in Cox PH models (via Grambsch and Therneau residuals) ----------------------------------------------------- ^stgtcalc^ ^stgtplot^ varlist [^, ci le^vel^(^#^)^ ^me^an graph_options ] Description ----------- ^stgtcalc^ calculates local estimates of the regression coefficients in the Cox model most recently fitted for using @stcox@. ^stgtplot^ plots the estimates and confidence intervals against the failure times. varlist is (a subset of) the independent variables in model fitted by ^stcox^. Options ------- ^ci^ causes a ^level()^% confidence band to be plotted around the estimates of the regression coefficients. The pointwise confidence intervals are derived from the running line smoother @running@. ^level(^#^)^ sets the width of a pointwise confidence band for the estimated b(t). Default: current value of system macro ^$S_level^ (usually 95%). ^mean^ uses the ^mean^ option of @running@ to produce a running mean smoother. graph_options are any of Stata's ^graph, twoway^ options Remarks ------- ^stgtcalc^ implements a suggestion of Grambsch and Therneau (1994). For an xvar with estimated regression coefficient b, ^stgtcalc^ calculates g(T) = b + #deaths * V^^-1 * r(T) at each failure time T, where V is the estimated dispersion matrix of the regression coefficients, #deaths is the number of uncensored observations and r(T) is the vector of Schoenfeld residuals at T. The g(T)'s are smoothed and the smoothed values are saved. g(T) is a local estimate of beta, so the shape of the smoothed profile may indicate whether beta depends on time. If so, the dataset has non-proportional hazards and the model will need to be refined. Note that the smoothed g(T) may be unreliable in the initial and final phases of follow-up. This may be partly due to lack of robustness in the smoother. The values of g(T) for each xvar are saved with names ^GT^xvar1, ^GT^xvar2, etc. If more than one xvar is to be plotted by ^stgtplot^, a panel of plots is produced. Example ------- . ^stset survtime cens^ . ^stcox x1 x2^ . ^stgtcalc^ . ^stgtplot x1^ . ^stgtplot x1 x2, ci mean^ Reference --------- Grambsch PM, Therneau TM. Proportional hazards tests and diagnostics based on weighted residuals. Biometrika 81: 515-526 (1994). Authors ------- Patrick Royston Imperial College School of Medicine, London FAX: (011)-44-181-383-8573 Peter Sasieni Imperial Cancer Research Fund, London FAX: (011)-44-171-269-3429 Also see -------- Manual: [R] st stcox On-line: help for @stcox@, @stbtcalc@ ---------------------------------cut here---end-of-stgtcalc.hlp-------- ********************************************************************** stgtplot.hlp Instructions for use: Cut out the following help file and put it in a file called "stgtplot.hlp". Be sure that the file is saved as a plain text (ASCII) file and that it has a hard return at the end of the last line in the file. Put stgtplot.hlp in your C:\ADO directory (you may have to create this directory) or your Stata working directory (i.e., your current directory). Note: Do not put it in the C:\STATA\ADO directory. Only official Stata ado and hlp files should be placed in this directory. ---------------------------------cut here---start-of-stgtplot.hlp------ .h stgtcalc ---------------------------------cut here---end-of-stgtplot.hlp-------- ********************************************************************** running.hlp Instructions for use: Cut out the following help file and put it in a file called "running.hlp". Be sure that the file is saved as a plain text (ASCII) file and that it has a hard return at the end of the last line in the file. Put running.hlp in your C:\ADO directory (you may have to create this directory) or your Stata working directory (i.e., your current directory). Note: Do not put it in the C:\STATA\ADO directory. Only official Stata ado and hlp files should be placed in this directory. ---------------------------------cut here---start-of-running.hlp------- .- help for ^running^ (STB-41: sed9.1) .- Symmetric nearest neighbor smoothing ------------------------------------- ^running^ yvar [xvar] [^if^ exp] [^in^ range] [weight] [^,^ ^ci^ ^d^ouble ^gen(^newvar^)^ ^gense(^sevar^)^ ^genb(^bvar^)^ ^k^nn^(^# | knnvar^)^ ^log^it ^m^ean ^nog^raph ^r^epeat^(^#^)^ ^sp^an^(^#^)^ ^tw^ice graph_options ] ^running^ smooths yvar on xvar. By default the smoothed version is a running line: a running mean is also available. A graph of yvar together with its smooth is plotted against xvar, unless suppressed. If xvar is not provided, yvar is smoothed against the ordered observations. Only analytic weights (^aweight^) are allowed; see help @weights@. Options ------- ^ci^ produces a pointwise confidence interval for the smoothed values of yvar. The width is determined by the current value of the macro $^S_level^. Not available with ^twice^, ^repeat^ or ^logit^. ^double^ doubles the value of ^repeat^. If ^repeat^ is not specified, ^double^ is equivalent to ^repeat(2)^. ^gen(^newvar^)^ creates newvar containing the smoothed values of yvar. Note that newvar will be on a logit scale if ^logit^ is used. ^gense(^sevar^)^ creates sevar containing the pointwise standard error of smoothed values of yvar. Not available with ^twice^, ^repeat^ or ^logit^. ^genb(^bvar^)^ creates bvar containing the local slope estimates. They constitute a local estimate of the derivative of the smoothed values of yvar with respect to xvar. Not available with ^mean^, ^twice^ or ^logit^. ^nograph^ suppresses the graph. ^knn(^# | knnvar^)^ controls the number (K) of Nearest Neighbors used on each side of the smoothed point. You may specify a constant (#) or a variable (knnvar). The value # is stored in $^S_1^. The greater the value, the greater the smoothing. You are not allowed to specify both ^span()^ and ^knn()^. ^logit^ transforms the smooth and plots the y-axis on a logit scale. Zero-one observations are automatically jittered in the vertical scale and are plotted just above and outside the range of the smoothed curve. ^mean^ specifies running-mean least-squares smoothing; default is running-line. ^repeat(^#^)^ specifies the number of times the data are to be smoothed. The default is 1. Increasing ^repeat^ increases the time it takes to calculate the smooth but improves the smooth. ^repeat(2)^ corresponds to "smoothing the smooth". ^span(^#^)^ controls the span or proportion of the data to be used in the symmetric nearest neighbors. If ^span^ is specified, ^knn^ is defined to be (N*^span^-1)/2, where N is the number of observations. If both ^span^ and ^knn^ are specified, an error results. ^span^ must be be in the range (0,2]. (It must be less than 1 when using ^mean^.) Span 2 corresponds to fitting a straight line. Stored in $^S_2^. ^twice^ carries out Tukey's `twicing' procedure whereby residuals from the original fit are smoothed and added back to the fit to obtain the final smooth ("smoothing the rough" or "reroughing" in Tukey's terminology). The result is somewhat rougher than would have been obtained without the application of twicing, but may be a better fit to the data. graph_options are any of the options allowed with ^graph, twoway^; see ^help^ ^graph^. Graph plots yvar followed by its smooth against xvar, the default options are s(oiii) (or s(.iii) with over 299 observations), c(.lll) and p(23) (or p(2344) when ^ci^ is specified). Remarks ------- Subsets of (2*k+1) observations are used for calculating smoothed values for each point in the data except for end points, for which smaller uncentered subsets are used. The subsets consist of the closest k points with xvar values less than or equal to that of the given point, the point itself, and the closest k points with xvar values greater than or equal to the given point. Since the neighborhoods are asymmetric in the tails, the running-mean is subject to bias in the tails. For this reason the smooth is set to missing in the tails if ^mean^ is specified. Other than in the tails, using ^mean^ will produce the same result as using the default smooth whenever the xvar values are evenly spaced. ^repeat(3)^, for instance, first smooths yvar creating yhat1, say; next yhat1 is smoothed creating yhat2, and finally yhat2 is smoothed creating yhat3. Author ------ Peter Sasieni Imperial Cancer Research Fund, London FAX 44-171-269-3429 Patrick Royston Royal Postgraduate Medical School, London FAX 44-181-383-8573 Also see -------- STB: STB-41 sed9.1, STB-24 sed9 Manual: [R] ksm, [R] smooth On-line: help for @ksm@ and @smooth@ ---------------------------------cut here---end-of-running.hlp---------