Second Answer: 2014-03-23
Update: 2014-03-26 I fixed the negative probabilities: I'd made an error in transcribing Emily's code. I also show that nlcom
as suggested on Statalist by Austin Nichols (http://www.stata.com/statalist/archive/2014-03/msg00002.html). I made one correction to Austin's code.
Bootstrapping is still the key to the solution. The target quantities are probabilities calculated by a formula that is based on a nonlinear combination of estimated parameters from streg
. As the estimates are not contained in the matrix This is an ideal situation for bootstrapping. The standard approach is adopted: create a program e(b)
returned after streg
, nlcom
will not estimate the standard errors.myprog
to estimate the parameters; then bootstrap
that program.
In the example, transition probabilities pt for a range of t values are estimated. The user must set the minimum and maximum of the t
range as well as a scalar u
. Of interest, perhaps, is that , since the number of estimated parameters is variable, a foreach
statement is required inside myprog
. Also, bootstrap
requires an argument that consists of a list of estimates returned by myprog
. This list is also constructed in a foreach
loop.
/* set u and minimum and maximum times here */
scalar u = 1
local tmin = 1
local tmax = 3
set linesize 80
capture program drop _all
program define myprog , rclass
syntax anything
streg , nohr dist(weibull)
scalar lambda = exp(_b[ln_p:_cons])
scalar gamma =exp(_b[_t:_cons])
forvalues t = `1'/`2'{
scalar p`t'= 1 - ///
(exp((lambda*((`t'-u)^(gamma)))-(lambda*(`t'^(gamma)))))
return scalar p`t' = p`t'
}
end
webuse cancer, clear
stset studytime, fail(died)
set seed 450811
/* set up list of returned probabilities for bootstrap */
forvalues t = `tmin'/`tmax' {
local p`t' = "p" + string(`t')
local rp`t'= "`p`t''" + "=" + "("+ "r(" + "`p`t''" +"))"
local rlist = `"`rlist' `rp`t''"'
}
bootstrap `rlist', nodots: myprog `tmin' `tmax'
forvalues t = `tmin'/`tmax' {
qui streg, nohr dist(weibull)
nlcom 1 - ///
(exp((exp(_b[ln_p:_cons])*((`t'-u)^(exp(_b[_t:_cons]))))- ///
(exp(_b[ln_p:_cons])*(`t'^(exp(_b[_t:_cons]))))))
}
Results:
Bootstrap results Number of obs = 48
Replications = 50
command: myprog 1 3
p1: r(p1)
p2: r(p2)
p3: r(p3)
------------------------------------------------------------------------------
| Observed Bootstrap Normal-based
| Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
p1 | .7009447 .0503893 13.91 0.000 .6021834 .7997059
p2 | .0187127 .007727 2.42 0.015 .0035681 .0338573
p3 | .0111243 .0047095 2.36 0.018 .0018939 .0203548
------------------------------------------------------------------------------
/* results of nlcom */
------------------------------------------------------------------------------
_t | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_nl_1 | .7009447 .0543671 12.89 0.000 .594387 .8075023
-------------+----------------------------------------------------------------
_nl_1 | .0187127 .0082077 2.28 0.023 .0026259 .0347995
-------------+----------------------------------------------------------------
_nl_1 | .0111243 .0049765 2.24 0.025 .0013706 .0208781
------------------------------------------------------------------------------