version 9.0
/* Name: separately by postoperative risk bcr learning curve.do
Description: multivariable model to test for association betweeen surgeon experience and biochemical recurrence following
radical prostatectomy, separately by pathologic stage.
Plots adjusted 5-year predicted probability of freedom from BCR against surgeon experience, separately by pathologic stage */
* change directory
cd "o:\Outcomes\Andrew\_UROLOGY Biostatistics\_Submitted\Accepted for publication\Bianco surgeons and recurrence\deidentify data for publication\"
* load deidentified data set for analysis
use "master learning curve data set deidentified.dta", clear
* drop patients who received neoadjuvant therapy
quietly drop if nht==1
* create dummy variables for pathologic gleason, which was categorized as <=5 / 6 / 7 / 8 / >=9
g pgrade5 = pathgleason<=5
g pgrade6 = pathgleason==6
g pgrade7 = pathgleason==7
g pgrade8 = pathgleason==8
g pgrade9 = pathgleason>=9
* create variables used in the model (will be manipulated later to plug in the mean values for all covariates)
foreach v of varlist psa ece svi lni pgrade5 pgrade6 pgrade7 pgrade8 pgrade9 {
quietly g model`v' = `v'
}
* create knots for surgeon experience (will be put at the quartiles)
quietly centile experience, c(25 50 75)
local k1 = r(c_1)
local k2 = r(c_2)
local k3 = r(c_3)
* create restricted cubic splines
spbase experience, gen(spexperience) knots(`k1' `k2' `k3')
* stset data
stset ttlastfollowbcr, f(bcr)
* create a program that will be run twice,
* first for patients with organ-confined disease and then for patients with non-organ-confined disease
capture program drop FIT_BY_NOCD
program FIT_BY_NOCD
* the input passed to the program is 0 (for organ-confined) and 1 (for non-organ confined)
local nocd = `1'
if `nocd'==0 {
local nocdtext="ocd"
}
else {
local nocdtext="nocd"
}
* multivariable model, with survival times assumed to follow a log-logistic distribution;
* covariates are PSA, ECE, SVI, LNI, and path Gleason
* specify cluster option using anonsurgeonid, since patients treated by the same surgeon are not independent
streg modelpsa modelece modelsvi modellni modelpgrade6 modelpgrade7 modelpgrade8 modelpgrade9 ///
experience spexperience1 spexperience2 spexperience3 if nocd==`nocd', dist(llog) cluster(anonsurgeonid) nolog
* global p-value for experience
test experience spexperience1 spexperience2 spexperience3
if r(p)<0.0001 {
disp "`nocdtext': Adjusted p-value for the association between experience and recurrence p<0.0001"
}
else {
disp "`nocdtext': Adjusted p-value for the association between experience and recurrence p=" string(r(p), "%9.4f")
}
* calculate adjusted 5-year predicted recurrence-free probability
* save out parameter from model (estimated by data), which is a constant for all patients
scalar gamma = e(gamma)
* set covariates to mean values
foreach v of varlist psa ece svi lni pgrade5 pgrade6 pgrade7 pgrade8 pgrade9 {
* obtain mean values of the covariates in the subgroup of organ-confined or non-organ-confined
quietly sum model`v' if nocd==`nocd'
quietly replace model`v' = r(mean)
}
* calculate the linear prediction and the standard error of the linear prediction
predict xb, xb
predict se, stdp
* calculate predicted survival probability and 95% confidence intervals
* obtain prediction for 5 years, i.e. 60 months
local time = 60
g s`nocdtext' = ( 1 + ( exp(-xb) * `time' ) ^ ( 1 / gamma ) ) ^ -1
g slb`nocdtext' = ( 1 + ( exp(-xb + invnorm(0.975)*se) * `time' ) ^ ( 1 / gamma ) ) ^ -1
g sub`nocdtext' = ( 1 + ( exp(-xb - invnorm(0.975)*se) * `time' ) ^ ( 1 / gamma ) ) ^ -1
* multiply by 100 to obtain estimates in percentages rather than proportions
local vars s slb sub
foreach v of local vars {
quietly replace `v'`nocdtext' = `v'`nocdtext' * 100
}
* set covariates back to their true values
foreach v of varlist psa ece svi lni pgrade5 pgrade6 pgrade7 pgrade8 pgrade9 {
quietly replace model`v' = `v'
}
drop xb se
end
FIT_BY_NOCD 0
FIT_BY_NOCD 1
* plot of predicted 5-year probability of freedom from biochemical recurrence against surgeon experience
twoway (line socd slbocd subocd snocd slbnocd subnocd experience, sort clpat(solid dash dash solid dash dash) ///
clwidth(medthick vthin vthin medthick vthin vthin) clcolor(gs8 gs8 gs8 gs0 gs0 gs0)), ///
ytitle(5-year probability of freedom from BCR (%), margin(medium)) ylabel(50(10)100, angle(horizontal) nogrid) ///
xtitle(Surgeon experience (number of prior surgeries), margin(medium)) xlabel(0(500)2000) ///
legend(off) scheme(s2mono) graphregion(fcolor(white))
* print out some estimates to use in the paper
local vars "ocd nocd"
foreach v of local vars {
disp " "
disp "`v':"
* predicted probability of recurrence at 10 prior cases
quietly sum s`v' if experience==10
local s10 = 100 - r(mean)
disp "Predicted probability of recurrence at 10 prior cases " string(`s10', "%9.1f") "%"
* predicted probability of recurrence at 250 prior cases
quietly sum s`v' if experience==250
local s250 = 100 - r(mean)
disp "Predicted probability of recurrence at 250 prior cases " string(`s250', "%9.1f") "%"
* difference in predicted probability for 250 vs 10 prior cases; confidence interval for difference will be obtained by bootstrap resampling
disp "Difference in predicted probability for 250 vs 10 prior cases " string(`s10' - `s250', "%9.1f") "%"
}