version 9.0
/* Name: heterogeneity in bcr by surgeon.do
Description: multivariable random-effects model to test for heterogeneity between surgeons in
biochemical recurrence outcomes following radical prostatectomy. */
* 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
* drop surgeons with fewer than 9 cases on the data set - otherwise, the model does not converge,
* and we hypothesize that this is because there are few patients per group.
sort anonsurgeonid
by anonsurgeonid: drop if _N < 9
* 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;
* and shared frailty for patients treated by the same surgeon;
* frailties are assumed to follow an inverse normal distribution
* covariates are year of surgery, PSA, ECE, SVI, LNI, and path Gleason
streg modelsurgyear modelpsa modelece modelsvi modellni modelpgrade6 modelpgrade7 modelpgrade8 modelpgrade9 ///
experience spexperience1 spexperience2 spexperience3, dist(llog) shared(anonsurgeonid) frailty(invgaussian) nolog
* obtain random effects variance and 95% confidence interval
local theta = exp(_b[/ln_the])
local thetalb = exp(_b[/ln_the] - invnorm(0.975) * _se[/ln_the])
local thetaub = exp(_b[/ln_the] + invnorm(0.975) * _se[/ln_the])
* obtain the p-value
if e(p_c)<0.0001 {
local pvalue = "p<0.0001"
}
else {
local pvalue = string(e(p_c), "%9.4f")
}
* display the random effects variance, 95% confidence interval, and p-value
disp "Random effects variance " string(`theta', "%9.3f") "; 95% CI " string(`thetalb', "%9.3f") ///
", " string(`thetaub', "%9.3f") "; `pvalue'"