version 9.0
/* Name: sensitivity analysis surgeons with at least 100 total cases bcr learning curve.do
Description: sensitivity analysis including only the patients treated by surgeons with at least 100 total cases
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)
}
* KEEP SURGEONS WITH AT LEAST 100 TOTAL CASES (i.e. surgeon experience for the last patient is at least 99)
sort anonsurgeonid
by anonsurgeonid: egen totalcases = max(experience)
keep if totalcases >=99
* 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 %
local vars s slb sub
foreach v of local vars {
quietly replace `v' = `v' * 100
}
* 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") "%"