// stgenreg JSS article examples adopath ++ ".\stgenreg_version_0_6_2" use ew_breast_ch7 stset survtime, failure(dead==1) exit(time 5) sts graph, by(dep5) //========================================// // Example 1 /* Weibull model */ streg dep5, dist(weib) estimates store streg stgenreg, loghazard([ln_lambda] :+ [ln_gamma] :+ (exp([ln_gamma]):-1):*log(#t)) nodes(15) ln_lambda(dep5) estimates store stgenreg15 stgenreg, loghazard([ln_lambda] :+ [ln_gamma] :+ (exp([ln_gamma]):-1):*log(#t)) nodes(30) ln_lambda(dep5) estimates store stgenreg30 stgenreg, loghazard([ln_lambda] :+ [ln_gamma] :+ (exp([ln_gamma]):-1):*log(#t)) nodes(50) ln_lambda(dep5) estimates store stgenreg50 stgenreg, loghazard([ln_lambda] :+ [ln_gamma] :+ (exp([ln_gamma]):-1):*log(#t)) nodes(100) ln_lambda(dep5) estimates store stgenreg100 estimates table streg stgenreg15 stgenreg30 stgenreg50 stgenreg100, eq(1,2) stats(ll) se //========================================// // Example 2 /* RCS for log baseline hazard */ cap drop haz* stgenreg, loghazard([xb]) xb(dep5 | #rcs(df(5)) ) nodes(30) predict haz1, hazard ci zeros twoway (rarea haz1_uci haz1_lci _t, sort)(line haz1 _t, sort) /// , ylabel(,format(%3.2f) angle(h)) ytitle("Hazard rate") xtitle("Follow-up time (years)") /// legend(order(1 "95% confidence interval" 2 "Baseline hazard rate" )) yscale(log) /* Now with time-dep effects */ stgenreg, loghazard([xb]) xb(dep5 | #rcs(df(5)) ) nodes(30) predict surv1, survival stgenreg, loghazard([xb]) xb(dep5 | #rcs(df(5)) | dep5:*#rcs(df(3))) nodes(30) predict surv, survival sts gen skm = s, by(dep5) //predict time-dependent hazard ratio ssc install partpred partpred hr, for(dep5 _eq1_cp3*) ci(hr_uci hr_lci) eform twoway (rarea hr_lci hr_uci _t, pstyle(ci) sort) (line hr _t, sort lpattern(solid) lwidth(thick)) if dep5==1 /// ,ytitle("Hazard Ratio") xtitle("Follow-up time (years)") /// yscale(log) twoway (line skm _t if dep5==0,sort connect(stairstep) )( line skm _t if dep5==1,sort connect(stairstep))(line surv1 _t if dep5==0,sort )( line surv1 _t if dep5==1,sort) /// , ylabel(, format(%2.1f) angle(h)) xtitle("Follow-up time (years)") legend(cols(2) ring(0) pos(7) order(1 "Affluent group, KM curve" 2 "Deprived group, KM curve" 3 /// "Affluent group, stgenreg" 4 "Deprived group, stgenreg")) ytitle("Survival") name(g1,replace) /// title("Proportional hazards") twoway (line skm _t if dep5==0,sort connect(stairstep))( line skm _t if dep5==1,sort connect(stairstep))(line surv _t if dep5==0,sort )( line surv _t if dep5==1,sort) /// , ylabel(, format(%2.1f) angle(h)) xtitle("Follow-up time (years)") legend(cols(2) ring(0) pos(7) order(1 "KM curve, dep5 = 0" 2 "KM curve, dep5=1" 3 /// "Predicted survival, dep5=0" 4 "Predicted survival, dep5=1")) ytitle("Survival") name(g2,replace) /// title("Non-proportional hazards") net install grc1leg.pkg grc1leg g1 g2, altshrink nocopies //========================================// // Example 3 /* Generalized Gamma (with proportional hazards) */ local mu [mu] local sigma exp([ln_sigma]) local kappa [kappa] local gamma (abs(`kappa'):^(-2)) local z (sign(`kappa'):*(log(#t):-`mu'):/(`sigma')) local u ((`gamma'):*exp(abs(`kappa'):*(`z'))) local surv1 (1:-gammap(`gamma',`u')):*(`kappa':>0) local surv2 (1:-normal(`z')):*(`kappa':==0) local surv3 gammap(`gamma',`u'):*(`kappa':<0) local pdf1 ((`gamma':^`gamma'):*exp(`z':*sqrt(`gamma'):-`u'):/(`sigma':*#t:*sqrt(`gamma'):*gamma(`gamma'))):*(`kappa':!=0) local pdf2 (exp(-(`z':^2):/2):/(`sigma':*#t:*sqrt(2:*pi()))):*(`kappa':==0) local haz (`pdf1' :+ `pdf2'):/(`surv1' :+ `surv2' :+ `surv3') stgenreg, hazard(exp([xb]):*(`haz')) nodes(30) xb(dep5,nocons) //========================================// // Example 4 //time-varying covariate use prothro,clear stset stop, enter(start) id(id) failure(event=1) list id pro trt _t0 _t _d if id==1 | id==111, noobs sepby(id) stgenreg, loghazard([xb]) xb(pro trt | #rcs(df(3))) nolog