version 9.0
/* Name: surgical margins learning curve.do
Description: multivariable model to test for association betweeen surgeon experience and surgical margins
Plots adjusted predicted probability of positive surgical margin 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 <= 6 / 7 / >=8
g pgrade6 = pathgleason<=6
g pgrade7 = pathgleason==7
g pgrade8 = pathgleason>=8
* 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 pgrade6 pgrade7 pgrade8 {
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 logistic regression model
* 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
logistic sm modelsurgyear modelpsa modelece modelsvi modellni modelpgrade7 modelpgrade8 ///
experience spexperience1 spexperience2 spexperience3, 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 surgical margins p<0.0001"
}
else {
disp "Adjusted p-value for the association between experience and surgical margins p=" string(r(p), "%9.4f")
}
* calculate adjusted predicted 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 pgrade6 pgrade7 pgrade8 {
quietly replace model`v' = scalar(mean`v')
}
* set year of surgery to the rounded year
quietly sum modelsurgyear
quietly replace modelsurgyear = round(r(mean))
* calculate the linear prediction and the standard error of the linear prediction
predict xb, xb
predict se, stdp
* calculate predicted probability and 95% confidence intervals
g phat = exp(xb) / (1 + exp(xb))
g phatlb = exp(xb - invnorm(0.975)*se) / (1 + exp(xb - invnorm(0.975)*se))
g phatub = exp(xb + invnorm(0.975)*se) / (1 + exp(xb + invnorm(0.975)*se))
* multiply by 100 to obtain estimates in percentages rather than proportions
local vars phat phatlb phatub
foreach v of local vars {
quietly replace `v' = `v' * 100
}
* plot of predicted 5-year probability of positive surgical margins against surgeon experience
twoway (line phat phatlb phatub experience, sort clpat(solid solid solid) clwidth(medthick thin thin) clcolor(gs0 gs0 gs0)), ///
ytitle(Probability of a positive surgical margin (%), margin(medium)) ylabel(10(5)50, 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 positive surgical margin at 10 prior cases
quietly sum phat if experience==10
local p10 = r(mean)
disp "Predicted probability of positive surgical margin at 10 prior cases " string(`p10', "%9.1f") "%"
* predicted probability of positive surgical margin at 250 prior cases
quietly sum phat if experience==250
local p250 = r(mean)
disp "Predicted probability of positive surgical margin at 250 prior cases " string(`p250', "%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(`p10' - `p250', "%9.1f") "%"