version 9.0
/* Name: primary analysis bcr learning curve.do
Description: multivariable model to test for association betweeen surgeon experience and biochemical recurrence following
radical prostatectomy. Plots adjusted 5-year predicted probability of freedom from BCR against surgeon experience */
* 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 surgyear psa ece svi lni pgrade5 pgrade6 pgrade7 pgrade8 pgrade9 {
quietly g model`v' = `v'
* obtain mean values of the covariates
quietly sum model`v'
scalar mean`v' = r(mean)
}
* 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)
* multivariable model, with survival times assumed to follow a log-logistic distribution;
* covariates are year of surgery, PSA, ECE, SVI, LNI, and path Gleason
* specify cluster option using anonsurgeonid, since patients treated by the same surgeon are not independent
streg modelsurgyear modelpsa modelece modelsvi modellni modelpgrade6 modelpgrade7 modelpgrade8 modelpgrade9 ///
experience spexperience1 spexperience2 spexperience3, dist(llog) cluster(anonsurgeonid) nolog
* global p-value for experience
test experience spexperience1 spexperience2 spexperience3
if r(p)<0.0001 {
disp "Adjusted p-value for the association between experience and recurrence p<0.0001"
}
else {
disp "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 surgyear psa ece svi lni pgrade5 pgrade6 pgrade7 pgrade8 pgrade9 {
quietly replace model`v' = scalar(mean`v')
}
* 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 = ( 1 + ( exp(-xb) * `time' ) ^ ( 1 / gamma ) ) ^ -1
g slb = ( 1 + ( exp(-xb + invnorm(0.975)*se) * `time' ) ^ ( 1 / gamma ) ) ^ -1
g sub = ( 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' = `v' * 100
}
* plot of predicted 5-year probability of freedom from biochemical recurrence against surgeon experience
twoway (line s slb sub experience, sort clpat(solid solid solid) clwidth(medthick thin thin) clcolor(gs0 gs0 gs0)), ///
ytitle(5-year probability of freedom from BCR (%), margin(medium)) ylabel(75(5)100, angle(horizontal) nogrid) ///
xtitle(Surgeon experience (number of prior surgeries), margin(medium)) xlabel(0(250)1750) ///
legend(off) scheme(s2mono) graphregion(fcolor(white))
* print out some estimates to use in the paper
* predicted probability of recurrence at 10 prior cases
quietly sum s 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 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") "%"