Re: st: repeat: predict after -streg-

[Jeph Herrin <junk@spandrel.net>]

Steven,

Thanks much for your careful reply. I was trying to
get probabilities for each patient, which was going
nowhere.

cheers,
Jeph


Steven Samuels wrote:
> 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/
*
*   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/