Re: st: repeat: predict after -streg-

[Steven Samuels <ssamuels@albany.edu>]

LINE 0

Here is code that does what Jeph asks.
On Jul 13, 2007, at 2:45 PM, Jeph Herrin wrote:

> I can rank on the hospital effects, but how do I convert the model
> parameters to estimates of the probability of 30-day readmission?
Use the formula  P{time<30} = P{ Z< log(30)-mu)/sigma)  where Z has a  
normal distirbution and mu and sigma are the mean and variance for log 
(time)

> > (log survival time scale or log hazard scale). Check on the  
> > agreement of the rankings from the different specifications. You  
> > may be able to use BIC or some other criterion to choose a "best"  
> > model.
> I did this initially, comparing LLs, AICs, and plots of Cox-Snell
> residuals. The log-normal model carried the day.
/************************CODE  
BEGINS**********************************************************/
/* MAIL MAY ADD GREMLINS TO THIS TEXT.  ZAP THEM IN YOUR EDITOR  
BEFORE RUNNING  */
set more off
capture log close
log using pcox02, replace text

use http://www.stata-press.com/data/r9/catheter ,clear

stset time, fail(infect)

sum age, det
tab female, missing

//Generate fake hospital id
gen hosp_id=113
replace hosp_id=217 if _n>25
replace hosp_id=66 if _n>50
tab hosp_id

/* Generate a unique code numbered 1,2,3,... .  Before this, delete  
any hospitals which will not appear in the regression */
egen hcode= group(hosp_id)  //Give hospitals a temporary sequential code
sum hcode
local n_hosps= `r(max)'  //Count Hospitals
di `n_hosps'
tab hcode, gen(xh)  // Dummy variables for regression

sort hcode
tempfile t1
save `t1'
gen age45=age -45 //center age so that zero is for a real number

local covariates= "age45 female"  //Covariate list
local ncovs: word count `covariates'
local n_predictors =`ncovs'+`n_hosps' //Needed to pick out the sigma
streg xh* `covariates', distribution(lnormal) nocons
matrix beta =e(b)
matrix list beta

di el(beta,1,`n_predictors'+1)

gen sigma=exp(el(beta,1,`n_predictors'+1))    //sigma for the normal  
distribution
di sigma
forvalues i=1/`n_hosps'{
di `i'
gen mu_`i'= el(beta,1,`i')       // mean for hosp_idital i when other  
covariates =0
gen prob30`i'=normal((log(30)-mu_`i')/sigma) //Prob{T<=30} for  
hosp_idital i
di `i' " " round(prob30`i',.01)
}

// Now merge the probabilities back with the original hosp_id numbers
gen dum=1
keep in 1
keep dum prob30*
reshape long prob30, i(dum) j(hcode)
label var prob30 "Est Prob Readmit in <=30 days}"
drop dum
sort hcode
merge hcode using `t1'
bysort hosp_id: keep if _n==1
list hosp_id prob30

/**********************************CODE  
ENDS****************************************/


*
*   For searches and help try:
*   http://www.stata.com/support/faqs/res/findit.html
*   http://www.stata.com/support/statalist/faq
*   http://www.ats.ucla.edu/stat/stata/