Statalist


[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

Re: st: repeat: predict after -streg-


From   Steven Samuels <ssamuels@albany.edu>
To   statalist@hsphsun2.harvard.edu
Subject   Re: st: repeat: predict after -streg-
Date   Sat, 14 Jul 2007 07:13:52 -0400

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/




© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index